STRATOSPEC is a pseudo-spectral code that solves the incompressible Navier–Stokes equations, either under the Boussinesq approximation (SBO) with the MHD framework, or under the Variable-Density approximation (SVD), where there is no magnetic field. Pseudo-spectral methods have the advantage to reach a spatial accuracy of infinite order. The SBO version has been used recently to investigate turbulent mixing in the Faraday instability [51] and in the magnetic Rayleigh–Taylor instability [52,53]; whereas the SVD version was used to investigate the varying properties of weakly coupled plasma under spherical compression [59].

A classical spectral Fourier collocation method is used with two-third rule dealiasing. The P3DFFT algorithm is used to perform massively parallel Fast Fourier Transforms [60]. The time increment is determined using a third-order, low-storage, strong-stability-preserving Runge–Kutta scheme, with an implicit treatment of diffusive terms. These numerical methods are common for both SBO and SVD.

STRATOSPEC-BOUSSINESQ ( SBO ): Within the incompressible Boussinesq and MHD approximations, the equations solved by STRATOSPEC are the following ones ∂ V ∂ t + V ⋅ ∇ V = − ∇ Π + 1 2 B 2 − g Θ n 3 + B ⋅ ∇ B + ν ∇ 2 V , (1a) ∂ B ∂ t + V ⋅ ∇ B = B ⋅ ∇ V + η ∇ 2 B , (1b) ∂ Θ ∂ t + V ⋅ ∇ Θ = κ ∇ 2 Θ , (1c) ∇ ⋅ V = ∇ ⋅ B = 0 , (1d) where V is the total velocity field and B is the magnetic field scaled as a velocity, defined as B = B ̃ / ρ 0 μ 0 , where B ̃ is the true magnetic field, μ ₀ is the magnetic permeability, ρ ₀ = (ρ ₁ + ρ ₂)/2 is the reference density, with ρ ₁ and ρ ₂ the minimum and maximum values of the density in the simulation domain. Π = P/ρ ₀ is the reduced pressure, ν, η and κ are the kinematic viscosity, magnetic diffusivity and molecular diffusivity respectively, and Θ is the dimensionless scalar field related to the density fluctuations around the reference state ρ ₀ through Θ = log ρ ρ 0 ≃ ρ ̃ ρ 0 , with  ρ = ρ 0 + ρ ̃ and  ρ ̃ ≪ ρ 0 . (2)

Hence, within the Boussinesq approximation, one has simply ρ = ρ ₀(1 + Θ). Moreover, note that within the Boussinesq approximation, the volume and mass fractions α are identical, and related to Θ and ρ through ρ ρ 0 = 1 + 2 A α − 1 2 = 1 + Θ , (3) where A = ( ρ 2 − ρ 1 ) / ( ρ 1 + ρ 2 ) is the Atwood number.

STRATOSPEC-VARIABLE-DENSITY ( SVD ): The Variable-Density approximation is a low Mach number limit for which density fluctuations can be large, in contrast to the Boussinesq approximation. Hence, the scalar field Θ defined in (2) is modified into Θ = log ρ ρ 0 and  Π = P ρ . (4)

Note that the generalized pressure Π is now normalized by the total density for convenience. The hydrodynamic evolution equations are then ∂ V ∂ t + V ⋅ ∇ V = − ∇ Π − Π ∇ Θ − g n 3 + ν ∇ 2 V + 1 3 ∇ D + F ν , (5a) F i ν = ν ∂ Θ ∂ x j ∂ V i ∂ x j + ∂ V j ∂ x i − 2 3 D δ i j , ∂ Θ ∂ t + V ⋅ ∇ Θ = κ ∇ 2 Θ , (5b) D = ∇ ⋅ V = − κ ∇ 2 Θ . (5c)

Note that the advection equation for Θ is formally the same between SBO and SVD, but there is an additional nonlinear term in (5a) involving the pressure, which translates the more complex effects of strong density gradients. In addition, the flow is not incompressible anymore with (5c): this equation shows that the scalar field is driven by the velocity divergence, and that the mixture between two incompressible fluids is compressible. Finally, the Poisson equation to obtain Π is solved using the GMRES algorithm [61].

The GAMERA code [50] solves the MHD equations (ideal, resistive or Hall) in 3D using a finite volume method for curvilinear non-orthogonal geometries. It is mainly used in magnetospheric and heliospheric simulations. It shares the numerical methods and the philosophy of the MHD Lyon–Fedder–Mobarry (LFM) code [48]. In addition, GAMERA was built from scratch for modern numerical architectures. It provides some improvements over LFM numerical schemes (e.g., seventh order upwind or eighth order centered numerical scheme) and the addition of computational advances such as massive hybrid parallelization, loop vectorization or data organization in blocks to optimize execution speed.

In this paper, the ideal compressible MHD equations solved by GAMERA are presented as follows ∂ ρ ∂ t = − ∇ ⋅ ρ V , (6a) ∂ ρ V ∂ t = − ∇ ⋅ ρ V ⊗ V + I ̄ ̄ P − ∇ ⋅ I ̄ ̄ B ̃ 2 2 μ 0 − 1 μ 0 B ̃ ⊗ B ̃ , (6b) ∂ E P ∂ t = − ∇ ⋅ V E P + P − V ⋅ ∇ ⋅ I ̄ ̄ B ̃ 2 2 μ 0 − 1 μ 0 B ̃ ⊗ B ̃ , (6c) ∂ B ̃ ∂ t = − ∇ × E , (6d) ∇ ⋅ B ̃ = 0 , (6e) where E _P = ρu ²/2 + P/(γ − 1) is the plasma energy, E the electric field and ρ the plasma density. The energy equation is formulated using the plasma energy E _P rather than the total energy which makes the system not totally conservative. However this choice was motivated by the fact that this formulation simplifies the numerical calculations with strong background magnetic fields and cold ambient plasmas which can be representative of some regimes present in magnetospheric simulations. In the ideal MHD case, Ohm’s law is simplified to E = − V × B ̃ . (7)

One numerical complexity of the compressible MHD equation system is the potential violation of the magnetic field solenoidal behavior. GAMERA uses the constrained transport algorithm [62] and guarantees that at each iteration the zero magnetic field divergence is maintained to machine precision. The sound speed is given by C s = γ P 0 ρ 0 , (8) with γ = 7/5 (rather than 5/3 in general, e.g., [11]). The Mach number is defined as M = Δ U C s , (9) where ΔU is the imposed initial velocity jump across the shear layer between the two fluids.

The KH configuration is greatly inspired by McNally et al. [54]. It consists of a 2D counter flow of two fluids of densities ρ ₂ and ρ ₁ < ρ ₂, at a fixed pressure P ₀ = 2.5, with zero gravity (g = 0), in a square box of width L = 1 with periodic boundary conditions. The light fluid surrounds the heavy fluid, and streams in the positive x direction with imposed velocity U ₁, while the heavier fluid moves in the negative x direction at velocity U ₂. The shear velocity, ΔU = U ₁ − U ₂ = 1, is used to determine the global Mach number M defined in (9).

The initial profiles of density and velocity are given by: ρ t = 0 = ρ 1 − ρ m e z − 1 / 4 / σ , z ∈ 0 ; 1 / 4 ρ 2 + ρ m e − z + 1 / 4 / σ , z ∈ 1 / 4 ; 1 / 2 ρ 2 + ρ m e z − 3 / 4 / σ , z ∈ 1 / 2 ; 3 / 4 ρ 1 − ρ m e − z + 3 / 4 / σ , z ∈ 3 / 4 ; 1 , (10a) V x t = 0 = U 1 − U m e z − 1 / 4 / σ , z ∈ 0 ; 1 / 4 U 2 + U m e − z + 1 / 4 / σ , z ∈ 1 / 4 ; 1 / 2 U 2 + U m e z − 3 / 4 / σ , z ∈ 1 / 2 ; 3 / 4 U 1 − U m e − z + 3 / 4 / σ , z ∈ 3 / 4 ; 1 , (10b) with ρ _m = (ρ ₁ − ρ ₂)/2 and U _m = (U ₁ − U ₂)/2. In McNally et al. [54], the Atwood number is A = 1 / 3 , which sets the Mach number M = 0.535. The parameter σ = 0.025 is a smoothing parameter that defines the thickness of the initial interface. Smoother profiles have the advantage to ease the convergence in terms of spatial resolution. The initial profiles are illustrated in Figure 1.

The initial perturbation is inspired from Nykyri et al. [63], with the initial vertical velocity profile being a sum of m _max contributions as follows V z x = V 0 2 m max ∑ m min m max sin m 2 π L x + ϕ m , (11) where V ₀ = 0.01 and each ϕ _m is a random phase that is fixed for all the GAMERA and STRATOSPEC simulations. Two initial conditions are considered in the next sections. The single-mode configuration (SM), with m _min = m _max = 2 (thus only the m = 2 mode is present) and ϕ ₂ = 0, so that (11) reduces to the initial condition used in McNally et al. [54]. In the multi-mode configuration (MM), we set m _min = 1 and m _max = 25.

For MHD simulations, an initial uniform magnetic field is imposed parallel to the flow, written as B _x (t = 0) = B ₀.

As mentioned in the introduction, the objective of the present study is to compare GAMERA and STRATOSPEC in the inviscid incompressible limit. This amounts to decrease the Mach number in GAMERA, which solves the compressible Euler equations, and to decrease the diffusion coefficients in STRATOSPEC, which solves the Boussinesq Navier–Stokes equations.

To do so, we first consider the hydrodynamic case in section 3, and both the single-mode and multi-mode initial conditions. We choose to work with small density contrasts between the counter flowing fluids to approach the Boussinesq limit, where the MHD module is available in STRATOSPEC. Hence, starting from the McNally et al. [54] configuration, the Atwood number is decreased from A = 0.33 to A = 0.05 , and this value is kept throughout the study. Still, the large density contrast case for A = 0.33 is addressed briefly in the Supplementary Figure S3A to show that the criteria chosen for A = 0.05 are still relevant. In addition, even smaller density contrasts are also investigated in the Supplementary Figure S2A, namely, A = 0.01 , which require more constrained parameters.

In the following analysis, the results of the simulations are compared at three levels: i) qualitatively, ii) 1D profiles in z averaged along the flow direction x, and iii) spectral analysis of the perturbed quantities. Qualitatively, plots of the full density fields are compared at three different dimensionless times, t = 1, 2 and 3 (with t = t ^⋆ΔU/L, with t ^⋆ the dimensioned time), to assess whether the structures are well reproduced within the two codes GAMERA and STRATOSPEC. Then, averages are performed, at t = 3, along the horizontal, homogeneous, periodic direction x to produce 1D profiles that depend only on the vertical coordinate z. The average of some quantity A along the flow direction is defined as ⟨ A ⟩ = 1 L ∫ 0 L A x , z d x , (12) and the average in the inhomogeneous direction reads a β ′ a β ′ ̄ = 1 Δ L ∫ L bottom L top ⟨ a β ′ a β ′ ⟩ d z , with  β = x  or  β = z , (13) with a′ = A − ⟨A⟩ and ΔL = L _top − L _bottom. For the SM initial condition, due to the top/bottom symmetry of the flow, the average is performed from L _bottom = L/2 to L _top = L, whereas for the MM initial condition, L _bottom = 0. Finally, spectral scale-by-scale comparisons are also provided at t = 3. We write a ̂ as the Fourier transform of the fluctuating field a′ and k the modulus of the wavevector k . The various spectra are defined as follows E a a k , t = ∫ S k a ̂ k , t a ̂ − k , t d 2 k , (14) where a is either v _x , v _z or ρ and S _k is the spherical shell of radius k = | k |. The one-point global variance can then be obtained either from the spherically-averaged spectra or the 1D profiles through ∫ E a a k , t d k = 1 Δ L ∫ ⟨ a ′ z , t a ′ z , t ⟩ d z . (15)

This methodology is used for both the hydrodynamic case in section 3 and the MHD case in section 4.

The density contrast is set to A = 0.05 , with the averaged density being ρ ₀ = (ρ ₁ + ρ ₂)/2 = 1. To evaluate the effects of compressibility for a given Atwood number, the Mach number M is decreased in GAMERA by gradually increasing the reference pressure from P ₀ = 2.5 to P ₀ = 80. Results are shown in Figure 2A.

We note that for the largest value of the Mach number, M = 0.535, there are strong discrepancies. The densities of the unmixed fluids depart from their rest values. This can be seen on the mean density ⟨ρ⟩ in Figure 2A, in particular far away from the shear layers, e.g., at z ≃ 0.5 and z ≃ 1.0, and consequently, also on the density variance ⟨ρ′²⟩. Regarding the velocity correlations, the effect of the reference pressure is less pronounced such that the horizontal mean velocity ⟨V _x ⟩ is hardly affected. For all quantities, a satisfactory convergence is reached for M ≤ 0.189. Results cannot be distinguished for M = 0.134 and M = 0.095.

A smaller Atwood number is briefly discussed in the Supplementary Figure S2B, and we show that a smaller Mach number must be chosen in GAMERA to ensure a satisfactory comparison with STRATOSPEC. Convergence in terms of spatial resolution for GAMERA is also addressed in the Supplementary Figure S1. We show that resolution with 10,24² and 20,48² points yield almost similar results, so that 20,48² points are chosen for the simulations yielding the main comparisons of this study.

In order to approach the Euler limit (i.e., vanishing bulk viscosity) with the Boussinesq version of STRATOSPEC (SBO), the diffusion coefficients ν = κ are decreased. Conjointly, the number of points is increased to ensure a sufficient spatial resolution. Mean fields in Figure 2B are only slightly affected by a decrease of the diffusion coefficients compared to second-order correlations. The density variance is significantly increased in the turbulent mixing zone when the kinematic viscosity ν is lowered. On the contrary, both horizontal and vertical kinetic energies, along with the vertical mass flux, are less impacted.

In conclusion, for the Atwood number discussed in this study, A = 0.05 , the inviscid incompressible limit is reached for a resolution of 20,48² points, a Mach number M = 0.134 for GAMERA, and diffusion coefficients of ν = κ = 1.27 × 10⁻⁶ for STRATOSPEC. These parameters are retained throughout the article.

The density field for the three codes (GAM, SBO and SVD) is shown for the single-mode perturbation (m = 2) at three different times in Figure 3.

There is an excellent qualitative agreement, with the presence of rolling structures as the consequence of the mean shear. There are no small-scale shear instabilities: indeed, smaller vortices are observed for large density contrasts due to baroclinic torque [64]. They can be seen in the test case of McNally et al. [54] at A = 1 / 3 , and are reproduced with GAMERA for the same Atwood number in Supplementary Figure S3A.

A slight difference can be observed at t = 3 when looking at the bottom vortices, which are slightly shifted towards the left in the GAMERA simulation (top line) compared with the STRATOSPEC-BOUSSINESQ simulation (middle), where the vortices are fixed at the same location. We evaluate a vortex speed along the x-direction of ≃ − 0.012 in GAMERA. The same motion is confirmed with the STRATOSPEC-VARIABLE-DENSITY simulation (bottom line), showing that the advection of the vortices is a non-Boussinesq effect, due to non-zero density contrasts, and not a compressibility effect [65]. Note that if the thickness of the shear layer is neglected, the phase velocity of the KHI can be evaluated as v _ph = (ρ ₁ U ₁ + ρ ₂ U ₂)/(ρ ₁ + ρ ₂) in the linear regime [1], that would provide a speed of −0.025 in our case. The observed discrepancy can be explained either by finite-thickness effects, or by the fact that the vortices are clearly in the nonlinear phase.

Horizontally averaged profiles are compared in Figure 4A. Looking at GAM (black) and SBO (red) first, slight differences in intensity can be observed for the mean density ⟨ρ⟩ at the center of the shear layer, as well as for the density variance ⟨ρ′²⟩. The horizontal kinetic energy ⟨ v x ′ 2 ⟩ profiles show some differences around the local maxima, which is the consequence of the Boussinesq simulation SBO yielding perfectly symmetric statistics between the light and heavy fluids. The SVD simulation (blue) in Figure 4A allows one to disentangle the origins of this asymmetry. Indeed, for ⟨ρ⟩, GAM and SVD are superimposed, indicating that the difference with SBO at the center of the mixing region is mainly a non-Boussinesq effect rather than a viscous effect. The same conclusion holds for ⟨ v x ′ 2 ⟩ .

In contrast, SBO and SVD are superimposed for ⟨ρ′²⟩, which shows that this difference with GAM is attributable to non-zero diffusion coefficients in the STRATOSPEC simulations. We investigate further the asymmetry in the GAMERA simulation with a smaller Atwood number in the Supplementary Figure S2B. The asymmetry is less pronounced when lowering the Atwood number, provided the Mach number is jointly decreased, leading to a good agreement between GAM, SBO and SVD. Conversely, at a much larger Atwood number, we show that diffusion smooths the perturbations and regularizes the small-scale shear instabilities in STRATOSPEC (see the Supplementary Figure S3A).

We conclude from the mean fields ⟨ρ⟩ and ⟨V _x ⟩ that the large-scale flows are well reproduced by the three codes. The second-order correlations ⟨ρ′²⟩, ⟨ v x ′ 2 ⟩ , ⟨ v z ′ 2 ⟩ and ⟨ v z ′ ρ ′ ⟩ illustrate that some differences persist at the smallest scales, which are mainly due to viscosity and diffusion in STRATOSPEC.

Spectra of the second-order moments, which include both the information at large and small scales, are shown in Figure 4B and support this point. Excellent agreement is found between the codes for the large scales (small k) while there are more discrepancies at small scales (large k). Most of the small-scale differences can be explained by a stronger dissipation in STRATOSPEC due to non-zero diffusion coefficients. The comparison with the Variable-Density formulation of STRATOSPEC also shows that differences at intermediate scales are only due to non-Boussinesq effects, as already pointed out in Figure 4A.

We continue the analysis of the hydrodynamic KHI by considering the multi-mode case (MM) given in Eq. 11 with m _min = 1 and m _max = 25 modes. The density fields of the GAMERA and STRATOSPEC simulations are first presented in Figure 5. The agreement is excellent, showing that nonlinearities are well captured by both codes. Contrary to the SM case, MM perturbations create a top/bottom asymmetry, already visible at t = 2. Indeed, the merging process, where small vortices are absorbed by bigger ones, results in larger vortices of various shapes and does not proceed at the same rate at the two shear layers.

The asymmetry of the flow requires computation of the 1D horizontally-averaged profiles over the whole domain z ∈ [0, 1] as done in Figure 6A. The overall agreement is again satisfactory. The strong asymmetry between the top and bottom parts in the horizontal kinetic energy ⟨ v x ′ 2 ⟩ is well captured and is attributed to the two bottom vortices which have different heights.

Finally, the density and velocity spectra for the MM case are presented in Figure 6B. Similar to the SM perturbation, large scales are superimposed while more discrepancies persist at small scales. Nevertheless, the overall agreement is better as the flow becomes less sensitive to diffusion at smaller scales due to the interaction of several modes.

The effects of gradually increasing the mean magnetic field magnitude from B ₀ = 0 to B ₀ = 0.2 is shown in Figure 7A. Compared with the hydrodynamic case (B ₀ = 0), the vortices are progressively stretched and flatten when B ₀ increases, with the same orientation compared to the interface. Till B ₀ = 0.1 (corresponding to M a ‖ = 10 ), vortices fold the original interface. For the largest intensity, B ₀ = 0.2, exactly at the threshold M a ‖ = 5 for having folded structures [35], the vortices have almost vanished and only thin filaments of density close to ρ ₀ remain. This case will be analyzed further below.

Conversion, through the mean shear and magnetic field, of vertical kinetic energy v z ′ 2 ̄ into horizontal kinetic energy v x ′ 2 ̄ , both computed according to (13), can be analyzed by considering their time evolution shown in Figure 7B. Beyond the good agreement between GAM and SBO simulations, we observe an increasing horizontal kinetic energy at the expense of the vertical kinetic energy as the magnitude of B ₀ increases. This corresponds to the damping of the vortices and their stretching in the flow direction. The induced magnetic energies in Figure 7B decrease with increasing B ₀ because structures are more aligned with the flow and less distorted. Note that the magnetic energy in the flow direction b x ′ 2 ̄ is always stronger than in the transverse direction.

We pursue the investigation of the single-mode magnetized KHI for the case B ₀ = 0.2. The density fields are shown at three different times in Figure 8, with a nice agreement between the STRATOSPEC and GAMERA simulations. Differences with the hydrodynamic case are significant. The magnetic field acts as a tension that tends to stabilize the interface and prevent the rolling of the structures. At t = 2, unlike the hydrodynamic case (Figure 3), there is almost no rolling. At t = 3, only thin filaments survive. Indeed, being exactly at the threshold for rolling, the magnetic tension seems to be sufficiently strong for hindering the folding at t = 2, and even for causing its regression at t = 3. This is particularly visible in the bottom row of Figure 8, where magnetic field lines have been drawn in white for the GAMERA simulation. The line corresponding to the bluish plasma, clearly separating the yellow and dark-blue regions, is almost folded at t = 2 but goes back to a nearly straight line at t = 3. We also note that there is a very good correspondence between the field lines and the density structures that, in the incompressible limit, are generated by the sole advection. This means that the ideal MHD frozen-in law is well respected during the simulation duration. In fact, there is no sign of magnetic island along the magnetic inversion layers that have been created in correspondence with the thin filaments. Numerical diffusivity is thus sufficiently low for preventing the development of Type II VIR [36].

Regarding the horizontally-averaged profiles, there is an excellent agreement between GAM and SBO for the mean fields in Figure 9. For ⟨ρ⟩, there is almost no mixing, with a flat profile between ρ ₁ and ρ ₂, whereas a wider mixing region with a mean density ρ ₀ exists in the hydrodynamic case (see Figure 4A). This is consistent with the magnetic field preventing mixing through the development of small scales. The mean horizontal velocity ⟨V _x ⟩ is distorted due to the thin elongated filaments detaching from the interface. The induced mean magnetic field ⟨B _x ⟩ has significantly evolved compared with its initial constant state. It is vanishing in the inner regions where the thin filaments form, at z ≃ 0.68 and z ≃ 0.82, and much stronger inside the thin mixing layer. Regarding second-order correlations, we now also include the horizontal and vertical magnetic energies, ⟨ b x ′ 2 ⟩ and ⟨ b z ′ 2 ⟩ , in Figure 9. There is an overall good agreement for all statistics between the GAM and SBO simulations, with sharp variations being well captured. We also confirm this good agreement for a mean magnetic field twice as small (B ₀ = 0.1) for which vortices are much more distorted (see 1D profiles in Supplementary Figure S4.

Spectra of the various correlations are presented in Figure 10. Similar to the purely hydrodynamic case, we find excellent agreement at large scales. However, smaller scales (k/(2π) > 30) are found to be quite different for E _ρρ . This is attributed to more mixing in STRATOSPEC due to non-zero diffusion coefficients, which significantly reduces the variance. Since the MHD case is less distorted (the mean magnetic field smooths the small scales), the velocity and magnetic spectra remain very close, with an agreement at both large and small scales. The spectral index in the inertial range is roughly between −2 and −2.5.

We address the MM perturbation in the MHD case with B ₀ = 0.2 in Figure 11. Even though small scales are suppressed quite early by the mean magnetic field, some large vortices survive. The top/bottom asymmetry in the shear layers is quite visible.

The rolling is less pronounced than in the purely hydrodynamic case (see Figure 5) but it is more active than in the SM counterpart (Figure 8) with rolls being able to develop. It is worth noticing that the non-linear evolution of the KHI can be strongly influenced by the number of modes that have been perturbed at t = 0, even if some of them are growing slowy, or not growing at all, during the linear phase. This point has been raised in Matsumoto and Seki [66], where the vortex merging at large scales strongly depensd on the perturbation of stable small-scale modes. Similarly, SM or MM initial conditions lead to a completely different non-linear evolution in Nakamura and Fujimoto [67] and in Faganello et al. [68].

Like in SM case, we note a very good agreement between field lines and density structures (Figure 11, bottom row). Moreover, no sign of Type II VIR is present, even if very thin magnetic inversion layers have been generated by the vortex motion, like the one at x ∼ 0.15, x ∼ 0.3, at t = 3.

The agreement between GAM and SBO for the horizontally-averaged profiles of density, velocity and magnetic fields is quite remarkable, as shown in Figure 12. Small differences in amplitude for the various variances are observed and are attributed to non-zero viscosity in STRATOSPEC simulations. The top/bottom asymmetry is well recovered for all quantities including the variances. Modulations of the induced horizontal mean magnetic field ⟨B _x ⟩ are less intense than in the SM case.

Finally, spectra of the various second-order correlations are shown in Figure 13. Compared with the MM purely hydrodynamic case (Figure 6B), the magnetic field effect and smoother profiles lead to an agreement up to larger k, namely, k/(2π) > 200. Similar to the SM perturbation, all spectra of GAM and SBO are in excellent agreement (up to k/(2π) > 200), except for the density variance departing at k/(2π) ≃ 30. The agreement at large and intermediate scales is quite convincing. The spectral index in the inertial range is once again between −2 and −2.5: we make no further comment on this since 2D turbulence induced by single-mode or multi-mode perturbations is quite singular.