Turbulence in outer protoplanetary disks: MRI or VSI?

The outer protoplanetary disks (PPDs) can be subject to the magnetorotational instability (MRI) and the vertical shear instability (VSI). While both processes can drive turbulence in the disk, existing numerical simulations have studied them separately. In this paper, we conduct global 3D non-ideal magnetohydrodynamic (MHD) simulations for outer PPDs with ambipolar diffusion and instantaneous cooling, and hence conductive to both instabilities. Given the range of ambipolar Els\"{a}sser numbers ($Am$) explored, it is found that the VSI turbulence dominates over the MRI when ambipolar diffusion is strong ($Am=0.1$); the VSI and MRI can co-exist for $Am=1$; and the VSI is overwhelmed by the MRI when ambipolar diffusion is weak ($Am=10$). Angular momentum transport process is primarily driven by MHD winds, while viscous accretion due to MRI and/or VSI turbulence makes a moderate contribution in most cases. Spontaneous magnetic flux concentration and formation of annular substructures remain robust in strong ambipolar diffusion dominated disks ($Am\leq1$) with the presence of the VSI. Ambipolar diffusion is the major contributor to the magnetic flux concentration phenomenon rather than advection.


INTRODUCTION
Turbulence is essential to many fundamental processes in protoplanetary disks (PPDs), including angular momentum transport, dust dynamics, planet migration, and chemical mixing. Turbulence ensues as a consequence of hydrodynamic or magnetohydrodynamic (MHD) instabilities. Seeking observational evidence of these instabilities and distinguishing one from another have none the less been challenging. The unprecedented sensitivity and resolution of Atacama Large Millimeter/submillimeter Array (ALMA) brings a step closer to reveal the origins and nature of the turbulence in PPDs, being capable of probing turbulence strengths (Teague et al. 2016;Flaherty et al. 2017Flaherty et al. , 2018Flaherty et al. , 2020 and gas velocity structures (Teague et al. 2019 in the disk. The magnetorotational instability (MRI; Balbus & Hawley 1991) was considered to be the primary source that drives turbulence and angular momentum transport in PPDs. However, the extremely weak level of ionization inevitably induces nonideal MHD effects, which can weaken or suppress the MRI. Instead, magnetized disk winds are found to dominate the disk accretion (e.g., Bai & Stone 2013;Bai 2017;Béthune et al. 2017;Wang et al. 2019;Gressel et al. 2020). In outer regions of the disk ( 30 AU), numerical simulations demonstrated that the MRI is dampened due to ambipolar diffusion, yield-⋆ E-mail: cc795@cam.ac.uk † E-mail: xbai@mail.tsinghua.edu.cn ing weak levels of turbulence (e.g., Simon et al. 2013a,b;Lesur et al. 2014;Bai 2015;Simon et al. 2018). Latest global 3D non-ideal MHD simulations confirmed this result and further illustrated the coexistence of the weak MRI turbulence and MHD winds (Cui & Bai 2021, hereafter CB21).
As alternative sources of turbulence, purely hydrodynamic mechanisms have caught attention (Lyra & Umurhan 2019;Lesur et al. 2022). Among them, the vertical shear instability (VSI) is the disk analogue of the Goldreich-Schubert-Fricke instability (Goldreich & Schubert 1967). Because the stratification and thermodynamic conditions are very conductive to its onset, the VSI may span a large extent, likely covering the entire outer disks (e.g., Pfeil & Klahr 2019). Its robustness in the non-ideal MHD dominated regions (Latter & Papaloizou 2018;Cui & Bai 2020, hereafter CB20;Cui & Lin 2021;Latter & Kunz 2022) and the potential observational significance (Lin 2019;Flock et al. 2020;Schäfer et al. 2020) render the VSI a very promising hydrodynamic mechanism in PPDs.
The MRI and the VSI could both emerge and sustain in outer disks. However, extant published numerical simulations on them have been performed separately. In this work, we explore the interplay between them. We extend the work of CB21, where they focused on the properties of the MRI by using a cooling timescale comparable to the disk dynamical timescale (τ = 1), and hence inhibits the VSI. We also extend the work of CB20, where they performed simulations on the VSI by 2D non-ideal MHD simulations, and hence precludes the MRI. A three-dimensional, ambipolar diffusion dominated outer disk with instantaneous cooling (τ = 0), conductive to the VSI and the MRI, is now explored.
In addition, CB21 illustrated the presence of magnetic flux concentration, and the subsequent gas annular substructure formation (also known as "zonal flows", Johansen et al. 2009), in ambipolar diffusion dominated outer disks. The magnetic flux concentration phenomenon is proved to be robust; it occurs in situations as diverse as local and global, ideal and non-ideal MHD simulations (e.g., Bai & Stone 2014;Suriano et al. 2018Suriano et al. , 2019Riols & Lesur 2018, 2019Riols et al. 2020;Jacquemin-Ide et al. 2021;Cui & Bai 2021). Local and global simulations illustrate that the magnetic flux concentration could result in dusty rings and gaps formation (Riols et al. 2020;Xu & Bai 2022;Hu et al. 2022). In this work, we further seek to answer that whether the spontaneous magnetic flux concentration phenomenon could exist in VSI turbulent disks.
This paper is organized as follows. In §2, we briefly describe the numerical methods and simulation setup. In §3, we show the simulation results and discuss turbulence properties, angular momentum transport, and mangeitc flux concentraton and annular substructure formation. We finish in §4 with a summary and discussion.

METHODS
The numerical setup of this work closely resembles that of CB21 (see their section 2). The instantaneous cooling is now applied to trigger the VSI. Key elements are outlined here.

Disk Model
The disk model is established such that radial temperature and density profiles are set to power laws by ρ ∝ r −q D and T ∝ r −q T , where the power-law indices are qD = −2, qT = −1, respectively (Bai & Stone 2017, CB21). The gas scale height is H = cs/ΩK, where c 2 s = P/ρ is the squared isothermal sound speed. The disk aspect ratio H/r is a constant given qT = −1. It is set to H/r = 0.1 within ±3.5H across the midplane and transitions to H/r = 0.5 in the atmosphere. The temperature is relaxed to initial equilibrium values over a timescale τ cool . After each simulation time step ∆t, the gas temperature is adjusted by the amount A dimensionless thermal relaxation time is defined as τ ≡ τ cool /P orb , where local orbital time P orb (R) = 2π/ΩK, with Keplerian frequency ΩK = R −3/2 . To trigger the VSI, we adopt a locally isothermal disk with instantaneous cooling τ = 0. The initial magnetic fields are purely poloidal by prescribing an azimuthal vector potential. This yields a purely vertical field at the midplane, and its strength is set according to the midplane plasma beta parameter β0. Ambipolar diffusion is parameterized by ambipolar Elsässer number where vA is the Alfvén speed and ηA is the ambipolar diffusivity. It is taken to be constant within ±3.5H across the midplane, and gradually transitions to Am = 100 from the disk towards the atmosphere.

Numerical Setup
We use the grid-based high-order Godunov MHD code Athena++ (Stone et al. 2020) to conduct global 3D nonideal MHD simulations. We solve the equations of non-ideal MHD using the van Leer time integrator, the HLLD Riemann solver, and the piecewise linear reconstruction. Super timestepping is employed to speed up calculations for diffusive non-ideal MHD physics. The simulations are performed in spherical polar coordinates (r, θ, φ). To better present simulation results, cylindrical coordinate system is used with R = r sin θ and z = r cos θ. In code units, GM = R0 = 1, with M the central star mass and R0 a reference radius at the inner radial boundary. The simulation domain spans over r ∈ [1, 100], θ ∈ [0, π], and φ ∈ [0, π/4]. Three levels of mesh refinement are used as in CB21, where the most refined disk zone achieves 32 cells/H in r, 30 cells/H in θ at the midplane, and 12 cells/H in φ.
At the inner radial boundary, we set temperature to their equilibrium values. Density is extrapolated from the last active zone, assuming ρ ∝ r −q D . The angular velocity is set to the minimum value between equilibrium v φ and the angular velocity of solid body rotation ΩK(R0)R. The other two velocity components are vr = v θ = 0.
Other boundary conditions are standard. At outer radial boundary, the hydrodynamic variables are extrapolated from the last active zone, assuming ρ ∝ r −q D , T ∝ r −q T and v φ ∝ r −1/2 . Radial and meridional velocities vr and v θ are copied directly from the last active zone, except setting vr = 0 when vr < 0. Magnetic field variables of the inner and outer ghost cells are set by Br ∝ r −2 , B θ ∝ const., B φ ∝ r −1 . Our θ domain reaches to the poles, and polar wedge boundary condition is adopted there (Zhu & Stone 2018). Lastly, our φ boundary condition is periodic. Table 1 lists the simulation models. Model name with a 'c' denotes instantaneously cooling in the disk, which are newly conducted. The rest of the simulation models are taken from CB21. Model name with an 'E' indicates those which E φ = 0 is enforced at the inner boundary, in order to stabilize the disk innermost region (CB21). The midplane disk magnetization of the initial poloidal field is set to β0 = 10 4 , where β0 is the ratio of gas pressure to magnetic pressure. For typical disk models (Minimum-mass solar nebula Weidenschilling (1977) or similar models), this level of magnetization yields accretion rates in agreement with observational constraints (Simon et al. 2013b). We employ three levels of ambipolar diffusion strengths, Am ∈ {0.1, 1, 10}. In addition, all simulations are run up to 3000P0, where P0 = P orb (R0) = 2π is the orbital period at the inner boundary.

RESULTS
In this section, the analyses will mostly focus on the turbulence properties and the magnetic flux concentration of simulations with instantaneous cooling τ = 0. The comparison between τ = 0 simulations (this work) and τ = 1 simulations (CB21) will also be covered. In the interest of brevity, we point the reader to CB21 for more detailed technical descriptions, analyses, and discussions on ambipolar diffusion dominated 3D global simulations.

Overview of gas and magnetic field evolution
The evolution of gas and magnetic fields in the disk and wind zone is in a similar manner as of CB21, except that now in the stronger ambipolar diffusion models (Am0.1c and Am1c), the VSI turbulence sets in. We summarize the gas and magnetic fields evolution as follows.
In the disk zone, where ambipolar diffusion dominates, the gas evolution is determined by the MRI or VSI, or a combination of both, depending on the ambipolar diffusion strength (sections 3.2 and 3.3). Both types of turbulence contribute to the angular momentum transport process, while MHD winds is the primary source to drive disk accretion (section 3.4). In the azimuthal plane, vortices are not observed in all of our models (section 3.6).
The initial condition of our simulations sets large-scale purely poloidal magnetic fields threading the disk. Poloidal fields are then wound up by Keplerian rotation to generate toroidal fields, which then dominate the magnetic field strength ( Figure 6). In the disk atmosphere, where ideal MHD dominates, the large-scale poloidal fields result in the launching of magnetized winds. In the case of strong ambipolar diffusion models, poloidal magnetic flux accumulates into flux sheets at different radii, yielding annular substructures in radial surface density profiles (section 3.5).

The dominance between VSI and MRI
Commonly reported in global hydrodynamic and non-ideal MHD simulations, the non-linear saturated state of the VSI is characterized by large-scale coherent vertical oscillations (e.g. Nelson et al. 2013;Cui & Bai 2020), which is inherited from its linear corrugation modes (Nelson et al. 2013;Barker & Latter 2015). Figure 1 shows snapshots of azimuthally averaged vertical velocities vz at t = 2600P0 for different models. The top panels show simulation models with instantaneous cooling τ = 0. For comparison purposes, the bottom panels show the corresponding models from CB21 with the same Am but different cooling timescale τ = 1. Comparing each column in Figure 1, it is clearly seen that the VSI has developed for model Am0.1c and Am1c, while for model Am10Ec the coherent vertical motions do not appear, and the MRI turbulence dominates.
To quantify the effect of the VSI in τ = 0 simulations, we present spatial auto-correlation function (ACF) for vertical velocities in Figure 2: As our simulations are scale free, we adopt a normalized vertical velocity Vz ≡ vz/cs at each radius, and present ACF in natural logarithm of R. The angle brackets · denote the spatial averaging for a box spans over [5,15] 2000,2600]P0. Note that by definition, the ACFs obtained at ∆ ln R = 0 should be strictly equal to unity.
In Figure 2, the left panel presents models with Am = 0.1 ( Table 1). The ACF of model Am0.1c is clearly different from model Am0.1. Model Am0.1c has an alternating positive and negative pattern with strong correlation/anti-correlation amplitudes of |ACF| ∼ 0.1 − 0.2. The characteristic length-scale of this variation is ∆ ln R ∼ 0.25, equivalent to ∆R ∼ 2.5H. This lengthscale should correspond to the VSI corrugation modes wavelength. Indeed, it is roughly consistent with the wavelength seen in top-left panel of Figure 1. For model Am0.1, the oscillatory pattern is not seen, corresponding to the absence of the VSI.
The right panel of Figure 2 shows models with Am = 10. The ACFs have similar patterns for both models, where the amplitudes of ACFs first drop significantly and overshoot becoming negative, then almost vanishing for ∆ ln R 0.2−0.3. The similarity shared by the two curves in Figure 1, together with the lack of oscillatory pattern in ACF, suggest the prevalence of the MRI and the absence of the VSI in Am = 10 models.
The middle panel of Figure 2 shows ACFs for Am = 1 models. As in model Am0.1c, the ACF of model Am1c possesses similar alternating positive and negative patterns, though with lower amplitudes of |ACF| ∼ 0.05. The characteristic length-scale is shorter with ∆ ln R ∼ 0.1, or equivalently ∆R ∼ H, consistent with the wavelength seen in the topmiddle panel of Figure 1. We interpret such features as the competition between the VSI and the MRI turbulence, where larger-scale vertical oscillations of the VSI being partially disrupted by the MRI, leaving lower-amplitude, smaller radial lengthscale oscillations. In other words, our results imply that the VSI and the MRI coexist in model Am1c.

Turbulence levels and stresses
We quantify disk turbulence strengths by computing the root mean square of velocity fluctuations, where v denotes azimuthal averaging. Figure 3 shows the radial profiles of three components of density-weighted turbulent velocities δv/cs, normalized by local sound speed. Spatial averaging is performed over ∆R ∈ [0.95R, 1.05R] and ∆z ∈ [−3H, 3H], and temporal averaging over the last 400P0. For models Am0.1c and Am1c, the presence of the VSI enhances the turbulent velocities to δv/cs ∼ 0.02 − 0.04, comparing to their τ = 1 counterparts with the absence of the VSI. In model Am0.1c, the large amplitudes in meridional turbulent velocities δv θ is associated with VSI coherent vertical oscillations. In model Am10Ec, the turbulent velocities reach δv/cs ∼ 0.1, comparable to model Am10E, as both models are dominated by the MRI turbulence. Figure 4 shows radial profiles of the Shakura-Sunyaev α parameter for Reynolds and Maxwell stress. The Reynolds stress is found to be α ∼ 1 × 10 −4 , α ∼ 5 × 10 −4 , and α ∼ 1 × 10 −2 in the domain R ∈ [5, 15] for models Am0.1c, ). Blue curves denote ACF for simulations models with instantaneous cooling τ = 0. Red curves denote models with τ = 1. From left to right, the three panels present models with Am = 0.1, Am = 1, and Am = 10 in Table 1, respectively.

Angular momentum transport
To illustrate angular momentum transport in our simulations, Figure 5 shows the mass accretion rate driven by MHD winds and by turbulent stresses in units of mass of the sun per year. The mass accretion rate in code units is converted to physical   units by equation (25) The above expression is interpreted as follows. If we consider at the radius R = 1 in code units (R code = 1R0) to correspond to a physical radius of 30 AU, where the surface density is 10 g cm −2 , then the physical accretion rate can be directly read off from equation (4) withṀ code obtained from simulation, and is plotted in Figure 5.
In Figure 5, all three models possess accretion rates on the order of 10 −8 − 10 −7 M⊙ yr −1 , consistent with observed mass accretion rates (Hartmann et al. 1998;Herczeg & Hillenbrand 2008). Mass accretion is predominantly mediated by MHD winds for all three models, as being also witnessed in CB21. For models Am0.1c and Am1c, viscous driven accretion by MRI and VSI turbulence makes minor contributions. For model Am10Ec, Maxwell stress contributes nearly equally to MHD winds due to the vigorous MRI turbulence. Accretion driven by the Reynolds stress is small compared to that of Maxwell stress and MHD winds in all models.

Magnetic flux concentration and ring formation
CB21 demonstrates the presence of spontaneous magnetic flux concentration in strong ambipolar diffusion dominated (Am 3) , weakly MRI turbulent disks. On the other hand, weak ambipolar diffusion (Am = 10) permits rigorous MRI turbulence that disrupts the concentrated flux sheets. Here, we show that the flux concentration can sustain in strong ambipolar diffusion dominated, VSI turbulent disks. For model Am0.1c and Am1c, Figure 6 (top panels) and Figure 7 illustrate that the poloidal magnetic fields assemble into magnetic flux sheets at different radial locations 1 . The flux sheets are nearly stationary for model Am0.1c and slowly propagating radially outwards for model Am1c.
The inhomogeneous distribution of poloidal magnetic fields leads to variations in surface densities. Figure 8 shows the radial profiles of surface density and vertical magnetic field strength for models Am0.1c and Am1c. Overall, the radial variations in surface density profile is comparable to those in CB21. Gas density bumps are developed at radii R = 4, 8 in model Am0.1c and radii R = 4, 7, 10 in model Am1c. They are separated by several gas scale heights ( 5H). The surface density variations, measured by the peak-to-trough ratio minus one, in the bumps are typically modest (of order unity), except for a deep gap at R = 5 in model Am1c. It is not a unique feature for this model but also appeared in B3 and Am3E in CB21 (see their Figure 5), and likely caused by the stochastic nature of the magnetic flux concentration phenomenon.

Mechanism: advection or ambipolar diffusion?
Various physical mechanisms have been proposed for magnetic flux concentration in non-ideal MHD simulations of PPDs. We point the reader to a summary in section 6.2 of CB21. As these mechanisms involve distinct physical processes and are still under debate, here we do not seek to verify them. Instead, we aim to distinguish the role between advection and (ambipolar) diffusion of the magnetic field in flux concentration phenomenon in our simulations.
The rate of change of magnetic flux is described by the Faraday's Law. In its integral form, the rate of change of poloidal magnetic flux is related to the toroidal electric field (e.g., Lubow et al. 1994;Bai & Stone 2017), where the poloidal magnetic flux is Br(r, θ) r 2 sin θ dθ (6) and the toroidal electric field is On the right hand side, the first term corresponds to the advection of magnetic field B with the fluid at velocity v, and the second term arises from the diffusion of magnetic filed through the fluid due to ambipolar diffusion; ηA is the ambipolar diffusivity, JAD = −(J × b) × b is the effective current density associated with ambipolar diffusion, where J = c∇ × B/4π is the current density, and b = B/B is the unit vector of magnetic field. The angle brackets · denote azimuthal averaging, and above equations are written in Gaussian units. Taking model Am0.1c as an example, Figure 9 shows the rate of change of poloidal magnetic flux, represented by −RE φ , at a evolutionary time of t = 2000P0 averaged over ±200P0. Over this period, magnetic flux continues to build up at two radial intervals, R ∼ 3 − 5 and R ∼ 8 − 11, as seen in Figure 7. Flux accumulation (depletion) requires a positive (negative) radial gradient in −RE φ . Indeed, the right panel of Figure 9 reveals prominent features in total −RE φ at those two radial intervals, which reflect the local increase or decrease in magnetic flux, for which local increase corresponds to magnetic flux concentration. The similarity between the middle and right panels indicates that ambipolar diffusion is the major contributor to the accumulation of poloidal magnetic flux, whereas contribution from advection is minor as evident from the left panel. This result is also valid for model Am1c and for models in CB21 with Am 3.

Azimuthal vortices
In 3D hydrodynamic simulations of the VSI, vortices in the R−φ plane are commonly reported (e.g., Richard et al. 2016;Manger et al. 2020). On the other hand, surface density variations induced by the magnetic flux concentration could trigger Rossby wave instability (RWI) and subsequently produce vortices. However, seen in the right column of Figure 10, no vortices are found in all of our simulation models, regardless of the value of Am, despite that the azimuthal resolution is sufficiently fine to resolve VSI and RWI induced vortices (e.g., Manger et al. 2020;Li et al. 2001). We suspect the absence of vortices might be caused by 1) the limited azimuthal domain which requires high surface density variations to trigger RWI (Ono et al. 2016). Taking case (iii) in Ono et al. (2016) as an example, which is closest to our simulation model, their Figure 7 demonstrates that to excite m = 7 or higher-order RWI modes, the surface density variation shall satisfy peakto-trough ratio much larger than unity, exceeding that found for model Am0.1c in section 3.5; 2) the presence of magnetic field may prevent the rolling up of vortex sheets by magnetic tension (Lyra & Klahr 2011). Future studies should address this issue.

CONCLUSIONS AND DISCUSSION
This paper extends the previous work of CB21 to examine the interplay between VSI and MRI in ambipolar diffusion dominated outer PPDs, via global 3D non-ideal MHD simulations. Given the range of ambipolar Elsässer numbers explored, our main findings are: . Top panels: surface densities Σ as a function of radius at t = 0, 600, 1200, 1800, 2400, 3000P 0 from light to dark colors. Bottom panels: vertical magnetic fields normalized by initial values Bz/B z,t=0 as a function of radius at t = 1200P 0 (red) and t = 3000P 0 (black). Similar plot for simulation models with τ = 1 can be found in Figure 5 of CB21.
Angular momentum transport is dominated by magnetized disk winds in all models, where viscous accretion via turbulence plays a moderate role in most cases. Formation of annular substructures by spontaneous magnetic flux concentration remains robust in the VSI turbulent disks (model Am0.1c and Am1c). The magnetic flux concentration is primarily attributed to the ambipolar diffusion.
The dominance between the VSI and the MRI with the presence of magnetic fields has been studied in linear theory. In the ideal MHD limit, the MRI linear modes grow faster than the VSI by a factor of ∼ R/H (Latter & Papaloizou 2018), while non-ideal MHD effects can rescue the VSI and dampen the MRI (Cui & Lin 2021;Latter & Kunz 2022). Our simulation results show consistency with the linear theory for small (Am = 0.1) and large (Am = 10) ambipolar Elsässer numbers. For the intermediate case of Am = 1, our results for the first time demonstrate the co-existance of the VSI and the MRI turbulence in ambipolar diffusion dominated outer disks.
We employ a disk equilibrium temperature gradient to be T ∝ r −q T , with power-law index qT = 1. In a realist disk, this power law is expected to be shallower qT ≈ 1/2 (Chiang & Goldreich 1997), which will result in weaker VSI turbulent velocities δv and smaller VSI-induced α values (Nelson et al. 2013;Manger et al. 2021). Moreover, instantaneous cooling (τ = 0) is employed in all models of this work. The cooling timescale in a realistic disk depends on two processes: the thermal coupling between the gas and dust, and the optical depth (Bae et al. 2021). Taken these processes into account, it is found that greater than ∼ 10 AU, the fast cooling timescales could permit the VSI growth . But these timescales are also finite, which will result in weaker δv and smaller α compared to the instantaneous cooling adopted in this work (Nelson et al. 2013;Manger et al. 2021).
The coherent vertical oscillations driven by corrugation modes is a unique feature of the VSI, which can help differentiate itself from other dynamical processes in PPDs (e.g., Flock et al. 2020;Barraza-Alfaro et al. 2021;Blanco et al. 2021). Although in model Am0.1c the dominant radial lengthscale of the corrugation modes is relatively large, in model Am1c the interplay between VSI and MRI results in less coherent oscillations with narrower radial lengthscales of ∼ H, making it difficult to be identified and spatially resolved. Further complexity can arise from the inertial wave parametric instability that attacks the VSI corrugation modes, producing even smaller-scale turbulence in disks Svanberg et al. 2022).
This work along with CB21 demonstrates the robustness of magnetic flux concentration phenomenon and annular substructure formation in ambipolar diffusion dominated regions. To confront observations with theory, dust particles shall be incorporated in our simulations to study its dynamics under the influence of magnetic flux concentration. In MRIturbulent local shearing box simulations, Xu & Bai (2022) found that magnetic flux concentration leads to dust clumping at pressure maxima. It would be intriguing to test this finding in 3D global simulations, possibly with the presence of VSI turbulence. Moreover, the VSI and MRI may possess non-Kolmogorov energy spectra, and it could modify the grain collisional velocities calculated from classic Kolmogorov turbulence (Ormel & Cuzzi 2007;Gong et al. 2020Gong et al. , 2021.
Future study should focus on relaxing the limited azimuthal domain, in particular to examine the existence and evolution of vortices. It would also be useful to further incorporate dust particles in our global non-ideal MHD simulations and to explore dust trapping in pressure bumps induced by the magnetic flux concentration. Lastly, it is necessary to employ sufficient fine resolution to resolve the parametric instability that attacks the VSI corrugation modes and drives small-scale inertial wave turbulent cascade.

DATA AVAILABILITY
The data underlying this article will be shared on reasonable request to the corresponding author.