Introduction

In his seminal works on unpredictability of atmospheric motion, Edward N. Lorenz1,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 complexity4.

The second behavior corresponds to a much more 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 arguments6,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.

Physically, this scenario could be the relevant one to address fundamental instabilities in classical fluid mechanics. An example is the Kelvin–Helmholtz (KH) instability, which describes the growth of a shear layer from an initially discontinuous velocity profile14. This instability is an important constituent in dynamical phenomena15 ranging from the micro-world of quantum fluids16,17 to macroscopic motions in Earth’s atmosphere and oceans18, extending further to astrophysical problems of supernovae19 and interstellar clouds20. While formally deterministic, the singular shear layer problem is however ill-posed in ideal fluid mechanics21. 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 data22,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 problems24,25,26,27; it also motivates the use of KH instability as a test-bed numerical problem28, and suggests the deterministic approach to unveil the physics of shear layer flows29,30. The Hadamard viewpoint, however, 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 conditions31,32.

Here we show that the classical Kelvin–Helmholtz singular shear layer becomes in fact well-posed and physically relevant when formulated in a probabilistic sense. Such a change from a deterministic to a probabilistic viewpoint reconciles 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 stochasticity33,34,35,36,37,38 describing unpredictability of trajectories of Lagrangian particles advected by rough (non-smooth) deterministic velocity fields, a key notion beneath the modern view on turbulent mixing39,40,41,42,43. This approach makes use of microscale 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. These properties are triggered by, but not sensitive to the nature of infinitesimal microscale 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 ocean44,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 \(\omega = \partial_x v - \partial_y u\) has the Dirac-delta singularity, \(\omega=U\delta(y)\). This is a steady solution to the inviscid incompressible Euler equation. For the ideal flow, the well-known linear theory14 predicts that a small perturbation with wavelength l has the exponential temporal growth \(\sim \exp (\pi Ut/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 equations46, 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 fluctuations7. Specifically, we consider initial conditions for the vortex sheet with infinitesimal small-scale perturbations of the form

$$t=0:\quad {\omega }_{\varepsilon }(x,y)=\left[1+\varepsilon \eta (x)\right]U\delta (y),$$
(1)

where \(\varepsilon\) is a small parameter and \(\eta(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 regularized dynamics, we focus on two fundamentally distinct formulations: (a) the viscous or hyperviscous (Navier–Stokes) regularization, which is controlled by the small viscosity parameter \(\nu \, > \, 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 \(\nu \propto {N}_{b}^{-1}\). Here, the first case employs a natural dissipative mechanism, 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 uniqueness47,48. Differently, the second case follows the Lagrangian approach in fluid mechanics, also relevant for superfluids49,50, in which one follows point vortices advected by the induced velocity field.

Recent numerical studies32 of the chaotic point-vortex system with finite \(N_b\) indicate a type of statistical universality of the nonlinear 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, \(\nu \to 0\) and \(\varepsilon \to 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 unambiguously defined as a non-trivial stochastic process. Mathematically51,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, \(\omega = U\delta(y)\). Physically, this implies that infinitesimal effects of microscale 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 \(\nu\) here also specifying the noise amplitude as \(\varepsilon \propto \nu\). We estimate the random component of the flow by measuring how two solutions, \((u,v)\) and \((u^{\prime} ,v^{\prime} )\), which differ initially by a very small random perturbation of size \(\varepsilon\), split in time. This is done by introducing the separation energy

$${\mathcal{E}}(t)=\frac{2}{L{U}^{2}}\int_{-L/2}^{L/2}dy{\langle {(u-u^{\prime} )}^{2}+{(v-v^{\prime} )}^{2}\rangle }_{x,\eta },$$

where \(L\) is a prescribed (macroscopic) observation scale and the average \(\langle \cdot \rangle_{x,\eta}\) is over both the x-direction and the statistical realizations of the initial noise. With this choice of normalization, \({\mathcal{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 \(\nu\) and \(\varepsilon\).

Fig. 1: Explosive separation between two realizations of the Navier–Stokes dynamics.
figure 1

Explosive separation is observed by measuring the separation energy \({\mathcal{E}}\) as a function of time t between two solutions of the Navier–Stokes equations, which are ε-close at the initial time. The observational scale is half the size of the computational domain \(L=\pi\) Each line corresponds to a different viscosity \(\nu\) and a different initial separation \(\varepsilon \propto \nu\) In the main figure, quantities are reported on a log–log scale and the dashed black line corresponds to the asymptotic linear law \({\mathcal{E}}\approx 0.14\ t\) emerging in the limit of vanishing viscosity and initial separation. The inset highlights the nonuniversal initial exponential stage. It shows in lin–log scales the data delimited within the black rectangle, along with exponential fitting. To highlight the exponential growth, the separation energy is rescaled by the separation at time \(\tau = 0.0044\) and data are shown as a function of \(t-\tau\).

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 effect3. However, it is apparent that the growth rate depends strongly on the viscous parameter \(\nu\). 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, \(\nu \to 0\) and \(\varepsilon \to 0\). At larger times, the separation energy reaches the asymptotics \({\mathcal{E}}\approx 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 flow53, where the distance between two Lagrangian fluid elements grows algebraically as \(\sim \!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. Fig. 2 shows a typical vorticity distribution for numerical simulations with the two regularization methods; see the Methods section for details of numerical implementations. Panels (a) and (c) correspond to the viscous regularization based on direct numerical simulations of the Navier–Stokes equations in a square domain with periodic boundary conditions. Panels (b) and (d) 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.

Fig. 2: Dynamics of vorticity.
figure 2

Snapshots of the vorticity fields at successive times for a the hyperviscous Navier–Stokes dynamics at \(\nu = 3 \times 10^{-3}\) and b Birkhoff–Rott dynamics using \(2^{17}=131,072\) point vortices. The arrows represent the velocity of the baseflow, characterized by a discontinuity of magnitude \(U\) at \(y=0\), and common to both dynamics. Panels c and d zoom into the structures indicated by the rectangular boxes, revealing the continuous or discrete nature of the vorticity field. The colorscale indicated in Panel c codes the magnitude of the vorticity in the Navier–Stokes case, and is common to panel a and panel c. In all figures, the space variables \(x,y\) and the time variable \(t\) are presented in units of \(L\) and \(L/U\), respectively, where \(L\) is the streamwise extension of the computational domain.

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 streaks28,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 theory56, 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 \(\ell(t)\). The latter describes the width of the vortex layer, and can be conveniently defined as

$$\ell (t)={\left(\frac{1}{U}\int{\left\langle {y}^{2}\omega ({\bf{R}},t)\right\rangle }_{x,\eta }dy\right)}^{1/2},$$
(2)

where \({\bf{R}} = (x,y)\). The average \(\langle \cdot\rangle_{x,\eta}\) is over both the x-direction 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: \(\ell = \alpha U t\), where \(\alpha\) 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 \(\alpha\simeq 0.029\) for both Navier–Stokes and Birkhoff–Rott regularizations; see the inset in Fig. 3.

Fig. 3: Growth of the mixing layer.
figure 3

The growth ot the mixing layer is observed by measuring the evolution of the mixing length \(\ell\) as a function of time. The main figure shows this evolution in log–log scale for both the Navier–Stokes regularization (blue) characterized by the viscosity \(\nu\) and the Birkhoff–Rott regularization (red) characterized by the number of vortices \(N_b\). The black dashed line represents the linear scaling \(\ell \sim 0.029 U t\), where \(U\) is the initial velocity jump. Inset shows the convergence for the universal dimensionless pre-factor \(\alpha=\ell/Ut\) towards its asymptotic value \(\alpha_\infty \simeq 0.029\) (black dashed line).

Statistical universality

We now study multiscale statistical properties of the flow. We start by analyzing the evolution of the normalized vorticity profile,

$${p}_{1}(y,t)=\frac{\ell (t)}{U}{\left\langle \omega ({\bf{R}},t)\right\rangle }_{x,\eta },$$
(3)

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\(\int p_1(y,t) \, dy=\ell(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)\approx {P}_{1}\left(y/\ell (t)\right)\). Panel (b) confirms that the asymptotic \(P_1\), corresponding to vanishing amplitudes of the initial perturbations, is identical for both the Navier–Stokes and Birkhoff–Rott regularizations.

Fig. 4: Statistical universality of the vorticity profile.
figure 4

The vorticity profile within the mixing layer is measured in terms of the one-point statistical quantity \(p_1(y,t)\), which is a function of the distance \(y\) to the initial discontinuity and of the time \(t\). a The main figure shows the vorticity profile \(p_1(y,t)\) as a function of the rescaled variable \(y/\ell(t)\) at various times logarithmically spaced between \(t = 0.02\) and \(t=0.51\) for the Birkhoff–Rott dynamics with the highest number of point-vortices \(N_b\). The inset shows the same quantity but for the Navier–Stokes dynamics with the lowest value of viscosity \(\nu\). In both cases, the collapse of the datasets (up to numerical noise) reveals the asymptotic convergence of the vorticity profile towards an asymptotic profile \(P_1(y/\ell)\). b The same graph, now at fixed time \(t=0.5\) for both the Navier–Stokes and point-vortex simulations, revealing universality of the self-similar profile \(P_1\) with respect to both regularizations.

As a second-order statistical observable, we introduce the isotropic two-point correlator

$${p}_{2}(r,t)=\frac{1}{{U}^{2}}\int{\left\langle \omega ({\bf{R}},t)\omega ({\bf{R}}+{\boldsymbol{\rho }},t)\right\rangle }_{x,\eta }\delta \left(\left\Vert {\boldsymbol{\rho }}\right\Vert -r\right){d}^{2}{\boldsymbol{\rho }}\ dy$$
(4)

written in dimensionless form. It is constructed from the covariance of the vorticity field between the points \({\bf{R}}=(x,y)\) and \({\bf{R}} + {\boldsymbol{\rho}}=(x+\rho_x,y+\rho_y)\) at a given distance \(\left\Vert {\boldsymbol{\rho }}\right\Vert =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 Eq. (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)\approx {P}_{2}\left(r/\ell (t)\right)\). 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\ell(t)\) characterizing a typical distance between nearest vortex blobs in Fig. 2. For very large distances \(r \gg \ell(t)\), it converges to the limiting value \(P_2 \to 2\). This value can be computed by substituting the vortex sheet expression \(\omega=U\delta(y)\) into Eq. (4), which approximates the vortex layer at large scales. It is remarkable that a constant state develops also at small distances \(r \ll \ell(t)\), with the universal value estimated as \(P_2 \to 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.

Fig. 5: Statistical universality of the vorticity correlation.
figure 5

The vorticity correlation is measured in terms of the quantity \(p_2(r,t)\), which depends upon the space separation \(r\) and the time \(t\). a The function \(p_2(r,t)\) is represented in log–lin scales as a function of the rescaled variable \(r/\ell\) (where \(\ell\) is the mixing length) at various times logarithmically spaced between \(t=0.02\) and \(t=0.35\). The main figure corresponds to the point-vortex simulations, while the inset shows the corresponding data for the hyperviscous run. In both cases, the obseved data collapse for \(r/\ell\ge 0.1\) indicates convergence of the large-scale statistics towards a function \(P_2(r/\ell)\). The horizontal dotted lines indicate common asymptotic values, \(P_2\to2\) for large separations and \(P_2 \to 2.4\) for small separations. b The same represention of the function \(p_2(r,t)\) now at the fixed time \(t=0.25\), with both Navier–Stokes and Birkhoff–Rott simulations superimposed, showing universality of the asymptotic profile \(P_2\). The horizontal black dotted line indicates the exact asymptotic value \(P_2\to 2\) and the vertical line indicates the local maximum, pointing to a typical separation between clusters of vortices.

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 multiscale fluid systems irrespective of how small the initial uncertainty2. 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. point-vortex Birkhoff–Rott formulations, we showed that this limiting stochastic process is universal.

The recently developed theory for particles advected by non-smooth deterministic velocity fields33,34,35,36,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 flows57,58,59,60,61, where it is referred to as the Richardson super-diffusion53. 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 particles11. 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 equations62,63, however, is still to be developed51.

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 systems64, as well as for systematic statistical mechanics approaches32. In addition to potential applications in fluid mechanics, which of course include the open problem of developed turbulence65 or the physics of boundary layers43,66, the relevance of the spontaneous stochasticity phenomenon may extend to other nonlinear field theories featuring multiscale dynamics and singularities, e.g., wave turbulence67, nonlinear optics68, and astrophysics69,70.

Methods

Navier–Stokes regularization

Our simulations were performed using the open-source GHOST parallel solver71,72, that employs a pseudo-spectral scheme with standard 2/3 dealiasing on spatial grids with \({N}_{a}^{2}\) points, and second-order Runge–Kutta scheme in time. The timestep is fixed and characterized by Courant number \(\sim 0.2\). Specifically, the Navier–Stokes equations in vorticity form

$${\partial }_{t}\omega +{\bf{u}}\cdot \nabla \omega ={D}_{\nu }\ \ {\rm{with}}\ \ {\bf{u}}=-{\nabla }^{\perp }\left({\Delta }^{-1}\omega \right)$$
(5)

are integrated within a two-dimensional periodic domain \((x,y) \in [0,2\pi]^2\) Results presented in the paper are obtained using a hyperviscous linear dissipation defined through its Fourier-space representation as \({\widehat{D}}_{\nu }=-\nu {k}^{2}{(k/{k}_{\max })}^{6}\widehat{\omega }\), where the regularization wavenumber \({k}_{\max }={N}_{a}/3\) is used; this is a common tool to localize viscous effects at small scales. Simulations with the standard viscous dissipation \(D_\nu = \nu \Delta \omega\) were also carried out employing up to 16, 3842 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 Eq. (1) are implemented along two parallel lines: \(y=0\) and \(y=\pi\) so that the total vorticity is zero. The noise \(\eta(x)\) is generated by a centered uniform stochastic process with spatial correlations \(\left\langle \eta (x)\eta ({x}^{\prime})\right\rangle =\delta (x-{x}^{\prime})/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 Table 1. These parameters follow the rule \({k}_{\max }\propto {N}_{a}\), \(\nu \propto {N}_{a}^{-3/2}\) and \(\varepsilon \propto {N}_{a}^{-3/2}\), so that \(\nu \to 0\) and \(\varepsilon\to 0\) simultaneously with the increasing resolution \(N_a\). To avoid spurious Gibbs oscillations of rough initial conditions, the initial state is diffused by hyperviscosity until the largest-wavenumber \({k}_{\max }\) has negligeable Fourier contribution, comparable to the machine precision. This takes a short initialization time \({t}_{d}\propto {N}_{a}^{-1/2}\), vanishing with increasing resolution. After this initial smoothing, the energy and enstrophy (squared norm of the vorticity) of the initial perturbation in expression Eq. (1) scale as \({\mathcal{E}}({0}^{+}) \sim {N}_{a}^{-3}{\mathrm{log}}\,{N}_{a}\) and \({\mathcal{Z}}({0}^{+}) \sim {N}_{a}^{-1}\): 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  \(\propto N_a\) upon increasing resolution.

Table 1 Parameters for the hyperviscous Navier–Stokes simulations with \(U=1/\sqrt{2}\).

Finally, please note that choosing simultaneous scaling rules for \(\varepsilon \to 0\) and \(\nu \to 0\) is dictated by the fact that these limits in principle do not commute. In particular, letting \(\varepsilon \to 0\) before \(\nu \to 0\), one should expect the process to converge towards the deterministic trivial baseflow.

Birkhoff–Rott regularization

In this case, we use a domain that is 1-periodic in the x-direction and unlimited in the y-direction. The initial condition Eq. (1) is approximated by a periodic discrete row of point-vortices located on the line \(y=0\) at positions \(x_n(0)=n/N_b, n=1,\cdots , {N_{b}}\), and carrying the initial vorticity \({\omega }_{n}=\frac{U}{{N}_{b}}\delta (x-{x}_{n})\delta (y)(1+\varepsilon \ {N}_{b}^{1/2}\ {\eta }_{n})\), where the \(\eta_n\)’s are independent Gaussian variables with unit variance. Using the complex position variables \(z_n=x_n+iy_n\), Biot–Savart law prescribes the advection of each point-vortex as54,73

$${\dot{z}}_{n}^{\star }=\frac{1}{2i} {\sum _{\begin{array}{c}1\le m\le {N}_{b}\\ m\ne n\end{array}}}{\omega }_{m}\ \cot \left[\pi \ ({z}_{n}-{z}_{m})\right],$$
(6)

where the star subscript denotes complex conjugation. We performed simulations of the dynamics Eq. (6) with \(U=1,\ \varepsilon=10^{-5}\), and \(N_b=2^{15}, 2^{16}\) and \(2^{17}\) vortices using the fourth-order Runge–Kutta scheme. The timestep \(\Delta t\) is chosen adaptively to satisfy \(\Delta t\,<\,0.1\ \mathop{\min }_{n\ne m}\left\{| {z}_{n}-{z}_{m}| /| {\dot{z}}_{n}-{\dot{z}}_{m}| \right\}\). 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

$$E=-\frac{1}{4\pi } {\sum _{\begin{array}{c}1\le n,m\le {N}_{b}\\ n\ne m\end{array}}}{\omega }_{n}\ {\omega }_{m}\ {\mathrm{log}}\,\left|\sin \left[\pi \ ({z}_{n}-{z}_{m})\right]\right|,$$
(7)

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 setups32,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 distribution75.

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\pi\) 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.