From the butterfly effect to intrinsic randomness: the spontaneous growth of singular shear flows

The butterfly effect is today commonly identified with the sensitive dependence of deterministic chaotic systems upon initial conditions. However, this is only one facet of the notion of unpredictability pioneered by Lorenz, who actually predicted that multiscale fluid flows could spontaneously lose their deterministic nature and become intrinsically random. This effect, which is radically different from chaos, have remained out of reach for detailed physical observations. Here, we substantiate this scenario by showing that it is inherent to the elementary Kelvin--Helmholtz hydrodynamical instability of an initially singular shear layer. We moreover provide evidence that the resulting macroscopic flow displays universal statistical properties that are triggered by, but independent of specific physical properties at micro-scales. This spontaneous stochasticity is interpreted as an Eulerian counterpart to Richardson's relative dispersion of Lagrangian particles, giving substance to the intrinsic nature of randomness in turbulence.

In his seminal works on unpredictability of atmospheric motion, Edward N. Lorenz 1,2 formulated a visionary conjecture on the sensitive dependence of deterministic systems upon initial perturbations or errors. He distinguished between two kinds of unstable behaviors: Either the mismatch between two replica of the same system can be made arbitrarily small by sufficiently reducing their initial discrepancy or, alternatively, the two systems reach diverging states, no matter how small they initially differ from each other. The first type of behavior arises in chaotic systems, wherein small perturbations are amplified exponentially in time. This phenomenon is now celebrated as the butterfly effect 3 and has been widely used to link unpredictability and complexity 4 .
The second behavior corresponds to a much more drastic intrinsic materialization of unpredictability. In Lorenz's own words: certain formally deterministic fluid systems which possess many scales of motion are observationally indistinguishable from indeterministic systems; specifically, that two states of the system differing initially by a small observational error will evolve into two states differing as greatly as randomly chosen states of the system within a finite time interval, which cannot be lengthened by reducing the amplitude of the initial error 2 . According to Lorenz, this second type of unpredictability pertains to fluid flows with a sufficient amount of kinetic energy at small scales, meaning that the velocity field should have a singular spatial structure. As pointed out by Palmer et al. 5 , the dynamical formulation of such fluid systems must be ill-posed: solutions do not depend continuously on initial conditions, therefore, allowing a finite-time separation between initially undistinguishable systems. In spite of being supported by phenomenological arguments [6][7][8][9][10][11][12][13] , a clear physical evidence is still lacking whether such a behavior indeed arises in genuine fluid flows, hereby making them infinitely less predictable than chaos.
In this work, we demonstrate that this scenario is in fact relevant for the fundamental instability in classical fluid mechanics: the Kelvin-Helmholtz (KH) instability which describes the growth of a shear layer from an initially discontinuous velocity profile 14 . This instability is an important constituent in dynamical phenomena 15 ranging from the microworld of quantum fluids 16,17 to macroscopic motions in Earths atmosphere and oceans 18 , extending further to astrophysical problems of supernovae 19 and interstellar clouds 20 . While formally deterministic, the singular shear layer problem is however ill-posed in ideal fluid mechanics 21 . The dominant viewpoint generally attributed to the French mathematician J. Hadamard is that such problems are physically meaningless, unless suitably regularized, so that solutions become uniquely determined and exhibit continuous dependence upon a class of initial data 22,23 . This viewpoint has motivated numerous studies employing analytic perturbation of the vortex sheet initial datum, that carefully analyze finite-time solutions of such deterministic and well-posed problems [24][25][26][27] ; it also motivates the use of KH instability as a test-bed numerical problem 28 , and suggests the deterministic approach to unveil the physics of shear layer flows 29,30 The Hadamard viewpoint surely appears quite at odds with Lorenz's intuition. It also fails to account for the essential footprint of KH flows: their visually striking large-scale features, that allow for identifying easily occurrences of KH instabilities across the physical communities, regardless of the detailed triggering mechanisms. This is the fundamental and long-standing observation that certain macroscopic statistical features of freely evolving shear flows are apparently only mildly dependent upon initial conditions 31,32 . Our objective here is to point out that the classical Kelvin-Helmholtz singular shear layer becomes in fact well-posed and physically relevant when formulated in a probabilistic sense. This implies a change from a deterministic to a probabilistic viewpoint, in order to recouncile Hadamard's to Lorenz's.
To unveil the intrinsically random nature of the shear layer dynamics, we employ a probabilistic approach inspired from earlier studies connected to the theory of spontaneous stochasticity [33][34][35][36][37][38] describing unpredictability of trajectories of Lagrangian particles advected by rough (non-smooth) deter-ministic velocity fields, a key notion beneath the modern view on turbulent mixing [39][40][41][42][43] . This approach makes use of micro-scale regularizing mechanisms, such as viscous dissipation and thermal fluctuations, to select relevant physical solutions at the macro-scale. Inferring from this connection, our main result is to reveal the intrinsic space-time randomness of the shear layer growth, characterized by the existence of a universal non-trivial probability distribution of the velocity field. Therefore, despite unpredictability, the resulting ideal flows possess well-defined statistical properties at finite times, which are triggered by, but not sensitive to the nature of infinitesimal micro-scale details.

Results
The singular shear layer. The classical mathematical formulation of the KH instability refers to perturbations of the idealized interface between two parallel streams with the velocity difference U . We will focus on two-dimensional incompressible formulation, as it is more accessible for accurate numerical simulations; in applications, this formulation refers to large-scale motions in the atmosphere and ocean 44,45 . Such streams are defined by the constant horizontal velocities: u(x, y) = −U/2 in the upper half-plane y > 0 and u(x, y) = U/2 in the lower half-plane y < 0; the vertical speed v(x, y) is everywhere zero. This flow can be seen as a vortex sheet localized on the horizontal axis y = 0, where the vorticity field ω = ∂v/∂x − ∂u/∂y has the Dirac-delta singularity, ω = U δ(y). This is a steady solution to the inviscid incompressible Euler equation. For the ideal flow, the well-known linear theory 14 predicts that a small perturbation with wavenlength l has the exponential temporal growth ∼ exp(πU t/l). The growth rate becomes infinite when the perturbation scale vanishes, implying the explosive breakdown of linear theory and signaling the ill-posedness of the underlying equations 46 , as required for Lorenz's intrinsic unpredictability. In order to select physically relevant solutions, we employ simultaneously two different small-scale mechanisms.
The first mechanism refers to a vanishingly small perturbation, which can be viewed as the effect of a small noise induced, e.g., by microscopic fluctuations 7 . Specifically, we consider initial conditions for the vortex sheet with infinitesimal small-scale perturbations of the form where ε is a small parameter and η(x) is a perturbation profile generated by white noise. The second mechanism desingularizes the equations of motion at vanishingly small scales, such that the resulting regularized equations can be used to evolve the flow. In order to explore the physical properties of the stochastic regularization, we focus on two fundamentally distinct formulations: (a) the viscous or hyperviscous (Navier-Stokes) regularization, which is controlled by the small viscosity parameter ν > 0 and (b) the point-vortex (Birkhoff-Rott) approximation, where the regularization is determined by the number of vortices N b and controlled by the small parameter ν ∝ N −1 b . Here, the first case employs a natural dissipative mechanism, FIG. 1. Explosive separation between two realizations of the Navier-Stokes dynamics measured as the energy E(t) of the difference between two solutions, which are ε-close at the initial time. The values of viscosity and initial separation corresponding to the three graphs are given in Tab. I, and the observational scale is half the size of the computational domain L = π. The main figure is shown in logarithmic scales and demonstrates the asymptotic algebraic law E ≈ 0.14t emerging in the limit of vanishing viscosity and initial separation. The inset, shown in linear horizontal and logarithmic vertical scales, zooms into the non-universal initial exponential stage. and the flow is considered in the so-called Eulerian formulation. For deterministic initial data, adding a vanishingly small viscosity is known to select weak solutions of the 2D Euler equations but is not enough to ensure uniqueness 47,48 . Differently, the second case follows the Lagrangian approach in fluid mechanics, also relevant for superfluids 49,50 , in which one follows point vortices advected by the induced velocity field. Recent numerical studies 32 of the chaotic point-vortex system with finite N b indicate a type of statistical universality of the non-linear dynamics with respect to the initial data over a wide range of timescales. While finite amplitude of any random initial perturbation trivially entails statistical description, robustness of the macroscopic features with respect to the perturbation size points towards a deeper phenomenon. In this work we unveil the origin of dynamical universality, showing that for vanishingly small values of viscosity and noise, the KH instability develops into a universal and spontaneously stochastic flow.
By spontaneous stochasticity, we mean that such a flow is intrinsically random. Universality signifies that the statistical properties of the flow are independent of the type of regularization, except at very small scales and times. The mathematical formulation of this phenomenon corresponds to the vanishing viscosity and noise, ν → 0 and ε → 0, where the two limits are taken simultaneously in a suitable fashion, such that physically meaningful norms of the initial disturbance (such as the kinetic energy) vanish; see the Methods section for more details. This limiting procedure grants the originally ill-posed inviscid problem with a well-posed statistical interpretation where the asymptotic flow is unambigu- ously defined as a non-trivial stochastic process. Mathematically 51,52 , one expects that each realization of this limiting process solves the incompressible Euler equations in a weak sense for the initial ideal vortex sheet, ω = U δ(y). Physically, this implies that infinitesimal effects of micro-scale regularization and noise do not select a unique deterministic KH flow, but rather a well-defined statistical solution, with universal non-stationary probability distribution.
The spontaneous butterfly effect. In order to reveal how the intrinsic randomness in the vortex layer emerges at small times, we first focus on the Navier-Stokes dynamics, with the viscous parameter ν here also specifying the noise amplitude as ε ∝ ν. We estimate the random component of the flow by measuring how two solutions, (u, v) and (u , v ), which differ initially by a very small random perturbation of size ε, split in time. This is done by introducing the separation energy where L is a prescribed (macroscopic) observation scale and the average · x,η is over both the x-direction and the statistical realizations of the initial noise. With this choice of normalization, E varies between 0 and 1, and those two extremal values are reached when the two solutions, respectively, coincide or fully decorrelate over the observational window. The results are presented in Fig. 1, where the separation energy is plotted for the three different sets of parameters ν and ε.
Initially the flows are very close and diverge exponentially as seen in the inset of Fig. 1. This exponential growth is governed by a positive Lyapunov exponent, which is the distinctive feature of the usual butterfly effect 3 . However, it is apparent that the growth rate depends strongly on the viscous parameter ν. In fact, the main part of Fig. 1 demonstrates that the exponential growth is a transient behavior. The transient region is shifted to smaller and smaller times in the combined limit of infinitesimal viscosity and initial perturbation, ν → 0 and ε → 0. At larger times, the separation energy reaches the asymptotics E ≈ 0.14t, which is independent of the initial separation and has the order of the total energy within the flow. This behavior quantifies spontaneous stochasticity for infinitesimal viscosity and initial disturbance: solutions which are initially undistinguishable and deterministic become distinct and random at finite times. It is analogous to Richardson's relative dispersion in fully developed turbulent flow 53 , where the distance between two Lagrangian fluid elements grows algebraically as ∼ t 3/2 , independently of how close they initially are. The explosive separation of Fig. 1 can be seen as the Eulerian counterpart to this phenomenon as it pertains to the full description of the flow. This type of randomness is spontaneous since it builds up in infinitesimal time and requires only infinitesimal perturbations to be revealed.
Time development of the vortex layer. We will now argue that while the growth of the vortex layer is an ill-posed dynamics from the deterministic viewpoint, it becomes well-posed in a probabilistic sense, meaning that macroscopic statistical features of the flow are independent of the noise and of the type of regularization. Figure 2 shows a typical vorticity distribution for numerical simulations with the two regularization methods; see the Methods section for details of numerical implementations. The panels on the left correspond to the viscous regularization based on direct numerical simulations of the Navier-Stokes equations in a square domain with periodic boundary conditions. The panels on the right correspond to simulations of the Birkhoff-Rott dynamics. Visually, the vorticity distributions obtained for each type of regularization look akin to each other on macroscopic scales. It is only upon zooming into the microscopic structure, that the distinction between the viscous and the point-vortex regularization becomes apparent, as the first is continuous and the second is discrete.
On macroscopic scales, the development of the shear layer consists in a cascade process of collisions and the subsequent mergers of smaller vortex blobs into larger ones; see Fig. 2.
Here, the blobs are small concentrated regions of high vorticity, surrounded by halos of lower-vorticity streaks 28,31,54,55 . In the collision process, a part of the vorticity is scattered away from the blobs, and it may be absorbed by the same or other blobs at later times. Measurements suggests that the vorticity is divided into two approximately equal parts corresponding to the blobs and to the background flow, while the enstrophy is mostly carried by the blobs.
We now verify numerically that the vortex layer can be qualified as being universal and spontaneously stochastic, namely that the increase of numerical resolution along with decreasing viscosity and noise yield a unique stochastic solution at finite times. By the scaling theory 56 , one may expect that statistical properties of the macroscopic dynamics depend only on the two independent quantities, namely, the velocity jump U and the mixing length (t). The latter describes the width of the vortex layer, and can be conveniently defined as where R = (x, y). The average · x,η is over both the xdirection and the statistical realizations of the initial noise. For the Birkhoff-Rott regularization, this quantity represents the standard-deviation of point-vortices across the vortex layer. From dimensional analysis, we expect that the mixing layer grows linearly in time: = αU t, where α is a dimensionless coefficient. Figure 3 presents the numerical results shown in log-log scale, which verify this asymptotic linear scaling law. The vortex layer is formed after a short transient time, and this transient time decreases as we decrease the regularization parameter and noise. The subsequent evolution yields the universal pre-factor α ≈ 0.029 for both Navier-Stokes and Birkhoff-Rott regularizations; see the inset in Fig. 3.
Statistical universality. We now study multi-scale statistical properties of the flow. We start by analyzing the evolution of the normalized vorticity profile, obtained from the vorticity distribution averaged with respect to the x-direction and the statistical realizations. In the Birkhoff-Rott case, p 1 relates to the one-point distribution of point-vortices along the y direction so that p 1 (y, t) dy = (t). Numerical vorticity profiles are shown in Fig. 4. Panel (a) highlights the self-similarity of this one-point statistical quantity, that is p 1 (y, t) ≈ P 1 (y/ (t)). Panel (b) confirms that the asymptotic profile P 1 , corresponding to vanishing amplitudes of the initial perturbations, is identical for both the Navier-Stokes and Birkhoff-Rott regularizations.
As a second-order statistical observable, we introduce the isotropic two-point correlator (4) written in dimensionless form. It is constructed from the covariance of the vorticity field between the points R = (x, y) and R+ρ = (x+ρ x , y+ρ y ) at a given distance ρ = r. The function p 2 (r, t) is strongly dominated by the relative positions of concentrated blobs of vorticity (see Fig. 2), for which the integrated product of vorticities in (4) is maximal. Results of numerical simulations are shown in Fig. 5, where Panel (a) displays the asymptotic self-similar relation p 2 (r, t) ≈ P 2 (r/ (t)). Panel (b) demonstrates that the asymptotic form of p 2 is universal: at a given finite time, it collapses onto a single universal function P 2 in the limit of vanishing regularization and noise for both Navier-Stokes and Birkhoff-Rott formulations. One can see that the profiles for different regularizations are distinct only at small scales, which reflects the different small-scale structures visible in the bottom plots of Fig. 2. The function P 2 has a pronounced maximum around r = 1.4 (t) characterizing a typical distance between nearest vortex blobs in Fig. 2. For very large distances r (t), it converges to the limiting value P 2 → 2. This value can be computed by substituting the vortex sheet expression ω = U δ(y) into (4), which approximates the vortex layer at large scales. It is remarkable that a constant state develops also at small distances r (t), with the universal value estimated as P 2 → 2.4; see Fig. 5. This value characterizes the uniform statistical distribution of the vorticity at distances which are small compared to the mixing length, but still larger than the regularization scale. In our view, this constant asymptotic value reflects the self-similar nature of the reconnection process, where vortex blobs constantly attract each other, before merging into larger blobs.

Discussion
Opposed to the exponential growth of observational errors in chaotic systems, now known as the butterfly effect, Lorenz anticipated that finite-time randomness can emerge in some formally deterministic multi-scale fluid systems irrespective of how small the initial uncertainty 2 . In our work, we demonstrated that this scenario is relevant for the flows generated by the Kelvin-Helmholtz instability of a singular shear layer. Being ill-posed for the ideal fluid mechanics, this system triggers a spontaneous stochastic process for the vorticity field, when considered in the combined limit of vanishing microscale regularization and vanishing random perturbation. In particular, considering hyperviscous Navier-Stokes vs. pointvortex Birkhoff-Rott formulations, we showed that this limiting stochastic process is universal.
The recently developed theory for particles advected by non-smooth deterministic velocity fields 33-37 demonstrate similar spontaneously stochastic properties when infinitesimal noise is introduced: fluid parcels separate in finite time no matter how close they are initially. Such an approach proved to be relevant for particles in turbulent flows [57][58][59][60][61] , where it is referred to as the Richardson super-diffusion 53 . Within this modern perception, Lorenz's scenario of intrinsic randomness can be seen as the Eulerian counterpart of spontaneous stochasticity applied to the flow fields rather than to individual particles 11 . In our point of view, this extension is rather natural, because the Eulerian description of fluid dynamics involves the transport of active quantities, such as momentum or vorticity. A comprehensive mathematical framework of this effect dealing with statistical solutions to fluid equations 62,63 , however, is still to be developed 51 .
Whether the finite-time unpredictability is caused by a small though finite "butterfly" as in chaos or by an arbitrarily small perturbation as in spontaneous stochasticity, the distinction may be illusive for an unprepared observer. However, these two phenomena represent fundamentally different physical and mathematical mechanisms that bridge the deterministic classical world to genuinely random systems. As such, spontaneous stochasticity could then possibly provide the salient mechanism that generates intrinsic randomness in natural phenomena. It provides a conceptual framework that justifies the relevance of stochastic modeling for apparently deterministic fluid dynamical systems 64 , as well as for systematic statistical mechanics approaches 32 . In addition to potential applications in fluid mechanics, which of course include the open problem of developed turbulence 65 or the physics of boundary layers 43,66 , the relevance of the spontaneous stochasticity phenomenon may extend to other nonlinear field theories featuring multi-scale dynamics and singularities, e.g., wave turbulence 67 , nonlinear optics 68 and astrophysics 69,70 .

Methods
Navier-Stokes regularization. Our simulations were performed using the open-source GHOST parallel solver 71,72 , that employs a pseudo-spectral scheme with standard 2/3 dealiasing on spatial grids with N 2 a points, and are integrated within a two-dimensional periodic domain (x, y) ∈ [0, 2π] 2 . Results presented in the paper are obtained using a hyperviscous linear dissipation defined through its Fourier-space representation as Dν = −νk 2 (k/kmax) 6 ω, where the regularization wavenumber kmax = Na/3 is used; this is a common tool to localize viscous effects at small scales. Simulations with the standard viscous dissipation Dν = ν∆ω were also carried out employing up to 16, 384 2 grid points: they lead to the same conclusions though for a shorter interval of scales due to numerical limitations.
To comply with periodicity, the initial conditions (1) are implemented along two parallel lines: y = 0 and y = π, so that the total vorticity is zero. The noise η(x) is generated by a centered uniform stochastic process with spatial correlations η(x)η(x ) = δ(x − x )/3. Simulations are stopped at sufficiently small times to avoid interaction between the layers and nonlocal effects of a finite domain. Then, the two resulting vortex layers are analyzed separately in the respective domains.
We performed nine simulations for each of the three parameter sets from Tab. I. These parameters follow the rule kmax ∝ Na, ν ∝ N −3/2 a and ε ∝ N −3/2 a , so that ν → 0 and ε → 0 simultaneously with the increasing resolution Na. To avoid spurious Gibbs oscillations of rough initial conditions, the initial state is diffused by hyperviscosity until the largestwavenumber kmax has negligeable Fourier contribution, comparable to the machine precision. This takes a short initialization time t d ∝ N −1/2 a , vanishing with increasing resolution. After this initial smoothing, the energy and enstrophy (squared norm of the vorticity) of the initial perturbation in expression (1) scale as E(0 + ) ∼ N −3 a log Na and Z(0 + ) ∼ N −1 a : Both quantities vanish upon increasing resolution. The initial (regularized) conditions therefore converge to the straight (unperturbed) vortex sheet in both energy and enstrophy metrics. We observe that our choice of scaling guarantees that the most unstable linear mode in our simulation is ∝ Na upon increasing resolution.
Finally, please note that choosing simultaneous scaling rules for ε → 0 and ν → 0 is dictated by the fact that these limits in principle do not commute. In particular, letting ε → 0 before ν → 0, one should expect the process to converge towards the deterministic trivial base flow. Birkhoff-Rott regularization. In this case, we use a domain that is 1periodic in the x-direction and unlimited in the y-direction. The initial condition (1) where the star subscript denotes complex conjugation. We performed simulations of the dynamics (6) with U = 1, ε = 10 −5 , and N b = 2 15 , 2 16 and 2 17 vortices using the fourth-order Runge-Kutta scheme. The timestep ∆t is chosen adaptively to satisfy ∆t < 0.1 min n =m {|zn − zm|/|żn −żm|}. This sort of Courant-Friedrichs-Lewy upper-bound is enough to ensure stability and accuracy of the numerical scheme. Indeed, with such a choice, the kinetic energy, defined as is conserved with a high-enough accuracy. At the end of the simulations, the relative error in energy remains below 10 −4 . Qualitatively similar CFL criteria were employed in previous simulations of point-vortex systems in similar set-ups 32,74 . While this CFL criterion in principle does not preclude divergence of dynamical trajectories due to numerical noise, it is enough to ensure convergence in a statistical sense. Upon increasing resolutions, the space-time distributions of point-vortices indeed converge towards a non-trivial measure.
In mathematical terms, the observed statistical universality reflects convergence in distribution 75 .
Finite-size effects. Formally, the scaling theory for the singular vortex sheet addresses flows evolving within an unbounded domain. In our numerical simulations, this domain has a finite size in the x-direction, taken as L = 2π for the Navier-Stokes and L = 1 for the Birkhoff-Rott models. In all figures, we displayed the results in dimensionless form using t * = L/U as the time unit, for which the mixing layer invades roughly 3% of the computational domain.

Data availability
The data that support the findings of this study are available from the corresponding author on request.

Code availability
The simulation and post-processing codes that have been used to produce the results of this study are available from the corresponding author on request.