MRI turbulence in vertically stratified accretion discs at large magnetic Prandtl numbers

The discovery of the first binary neutron star merger, GW170817, has spawned a plethora of global numerical relativity simulations. These simulations are often ideal (with dissipation determined by the grid) and/or axisymmetric (invoking ad hoc mean-field dynamos). However, binary neutron star mergers (similar to X-ray binaries and active galactic nuclei inner discs) are characterised by large magnetic Prandtl numbers, $\rm Pm$, (the ratio of viscosity to resistivity). $\rm Pm$ is a key parameter determining dynamo action and dissipation but it is ill-defined (and likely of order unity) in ideal simulations. To bridge this gap, we investigate the magnetorotational instability (MRI) and associated dynamo at large magnetic Prandtl numbers using fully compressible, three-dimensional, vertically stratified, isothermal simulations of a local patch of a disc. We find that, within the bulk of the disc ($z\lesssim2H$, where $H$ is the scale-height), the turbulent intensity (parameterized by the stress-to-thermal-pressure ratio $\alpha$), and the saturated magnetic field energy density, $E_\text{mag}$, produced by the MRI dynamo, both scale as a power with Pm at moderate Pm ($4\lesssim \text{Pm} \lesssim 32$): $E_\text{mag} \sim \text{Pm}^{0.74}$ and $\alpha \sim \text{Pm}^{0.71}$, respectively. At larger Pm ($\gtrsim 32$) we find deviations from power-law scaling and the onset of a plateau. Compared to our recent unstratified study, this scaling with Pm becomes weaker further away from the disc mid-plane, where the Parker instability dominates. We perform a thorough spectral analysis to understand the underlying dynamics of small-scale MRI-driven turbulence in the mid-plane and of large-scale Parker-unstable structures in the atmosphere.


INTRODUCTION
The transport of angular momentum is the key mechanism driving accretion in discs.Turbulence -thought to be driven by the magnetorotational instability (MRI), and potentially also magnetically-driven outflows, facilitate this angular momentum transport.Complicating matters, a characteristic of magnetohydrodynamic (MHD) turbulence is the existence of not one, but two dissipation mechanisms (viscosity and resistivity) and associated length-scales that can drastically alter the nature of MHD turbulence compared to its more familiar hydrodynamic (i.e.Kolmogorov-type) cousin.The ratio of these two diffusivities, known as the magnetic Prandtl number (Pm), is a key parameter in MHD turbulence and in dynamo theory (Rincon 2019), influencing not only the saturated state of the turbulence (with magnetic field growth and associated angular momentum transport), but also the thermal stability of the disc (Balbus & Henri 2008;Potter & Balbus 2017;Kawanaka & Masada 2019).
Different discs are believed to lie in different Pm regimes: the hot, dense, partially neutrino-cooled discs from binary neutron star mergers, and the inner regions of X-ray binaries and active galactic nuclei are thought to lie in the Pm ≫ 1 regime (Balbus & Henri 2008; ★ E-mail: leh50@cam.ac.uk (LEH) Rossi et al. 2008), while cooler, protoplanetary, discs are generally thought to lie in the Pm ≪ 1 regime (Lesur 2021).Global simulations however sophisticated, are almost always either ideal (Kiuchi et al. 2022), and thus characterized by a purely numerical magnetic Prandtl number which is likely of order unity (Minoshima et al. 2015), or otherwise employ mean-field dynamo or viscosity prescriptions (Fujibayashi et al. 2020;Shibata et al. 2021), prescriptions which, while very informative in their own right, are nevertheless ad hoc.Sub-grid modeling attempts to bridge the gap between ideal simulations and those with fully explicit dissipation coefficients, but comes with its own uncertainties (Meheut et al. 2015;Miravet-Tenés et al. 2022).Though it has been known for some time that MRI turbulence is particularly sensitive to the ratio of dissipation scales (Lesur & Longaretti 2007;Pessah & Chan 2008;Simon et al. 2011;Nauman & Pessah 2016;Nauman & Pessah 2018), the MRI in the regime Pm ≫ 1 has only been sparsely investigated in the literature.Recent work (Held & Mamatsashvili 2022;Guilet et al. 2022), however, has begun to forge inroads into this physically interesting, but numerically challenging (Reboul-Salze et al. 2022) regime, though it has so far been limited to unstratified and isothermal disc models.

The need for explicit dissipation coefficients
The magnetorotational instability requires a weak seed field in order to operate, though this need not be a mean (external) field.In fact, the MRI can sustain itself without the presence of a mean field by means of a dynamo, as demonstrated in both unstratified local models (e.g., Hawley et al. 1996;Fromang & Papaloizou 2007;Riols et al. 2017b;Walker et al. 2016;Mamatsashvili et al. 2020), and vertically stratified ones (e.g., Brandenburg et al. 1995;Davis et al. 2010;Shi et al. 2010;Gressel 2010;Salvesen et al. 2016).This setup (referred to as 'zero-net-flux' or ZNF) is particularly attractive because it is agnostic to the strength and geometry of the magnetic field of the central accretor or external medium, both of which can vary considerably from system to system.
Complicating matters, in the absence of explicit dissipation coefficients (viscosity and resistivity) both unstratified and stratified models exhibit the 'convergence problem', a numerical artifact that very likely afflicts ideal global simulations, too.In essence, in ideal MHD simulations, turbulent transport decreases with increasing numerical resolution (Pessah et al. 2007;Fromang & Papaloizou 2007).When vertical stratification is included, which gives the domain a 'global' character in the vertical direction, one encounters the same problem (Bodo et al. 2014;Ryan et al. 2017).Convergence with resolution is restored when dissipation coefficients are explicitly taken into account, however, and also leads to new and interesting physics: in the regime Pm ≲ 1 the zero-net-magnetic flux MRI struggles even to sustain itself (Walker et al. 2016;Riols et al. 2015;Nauman & Pessah 2016;Mamatsashvili et al. 2020), while in the opposite regime turbulent intensity is remarkably sensitive to Pm, exhibiting powerlaw scaling with Pm, at least up moderately large values of Pm ≳ 1 (Potter & Balbus 2017;Held & Mamatsashvili 2022;Guilet et al. 2022).Although a few authors have carried out stratified ZNF simulations at Pm ≳ 1 (Gressel 2010;Davis et al. 2010), the only work we are aware of that investigated the Pm ≫ 1 regime with vertical stratification is Simon et al. (2011), although even they only went up to Pm = 8.In this work, we aim to explore the behavior of the MRI dynamo in stratified discs at much larger magnetic Prandtl numbers (up to Pm=90) than has been previously done.

Magnetic buoyancy instabilities in accretion discs
Real accretion discs are stratified in the direction perpendicular to the disc, and interesting new physics arises when stratification is taken into account.In this case magnetic fields can become buoyant, rising through the disc (Miller & Stone 2000;Blackman & Pessah 2009;Shi et al. 2010;Uzdensky 2013;Dudorov et al. 2019), which is most easily revealed in the well-known 'spacetime diagrams' (Rincon 2019;Lesur 2021).These buoyant fields do not only play a passive role: the zero-net-flux MRI dynamo in stratified discs has been shown to behave much like an  − Ω dynamo (Gressel 2010;Gressel & Pessah 2015;Rincon 2019;Dhang et al. 2023).Finally, the addition of thermodynamics in stratified models further complicates things, leading to scenarios such as the interplay between the MRI and convection (Bodo et al. 2013;Hirose et al. 2014;Scepi et al. 2018;Held & Latter 2021), but comes with its own numerical challenges (Gressel 2013;Held & Latter 2021).
Yet another interesting effect that can occur in stratified discs is that of Parker instability (and associated dynamo).The Parker mode (also known as the undulatory mode) is characterized by perturbations whose wavenumbers are parallel to the magnetic field (k ∥ B) (Pringle & King 2007;Stone & Gardiner 2007;Hillier 2016).The instability involves bending of field lines and plasma slides down these field lines forming overdense regions in the troughs between the crests.Assuming the gravitational acceleration points in the -direction, the linear instability criterion is given by / < 0 (Pringle & King 2007).
The idea that magnetic fields in accretion discs might undergo Parker instability, leading to the generation of vertical field from toroidal and radial fields has been around for a long time (Shu 1974;Tout & Pringle 1992).In early 3D, vertically stratified, zero-net-flux, isothermal simulations the disc was found to be only marginally unstable to Parker instability (Stone et al. 1996), however this was likely due to the restricted vertical size of the domain which encompassed only ±2 scale-heights either side of the mid-plane.Later, simulations in taller boxes (spanning 6-9 scale-heights on either side of the mid-plane) initialized with a relatively strong initial toroidal magnetic field (i.e. with a gas-to-magnetic pressure ratio  ∼ 1-25) did observe the development of Parker instability (in addition to the MRI) both in isothermal discs (Johansen & Levin 2008;Kadowaki et al. 2018), and also in non-isothermal simulations with radiative transfer (Blaes et al. 2007;Shi et al. 2010;Blaes et al. 2011).As we report in this work, we find the emergence of Parker instability even when the disc is initialized with relatively weak ( ∼ 10 3 ) zero-net-magnetic-flux.

Motivation and outline
Our aim in this paper is to build on our earlier work (Held & Mamatsashvili 2022) on the saturation and energetics of zero-net-flux MRI turbulence (alternatively referred to as the 'MRI dynamo') in the regime of large magnetic Prandtl number, mainly by taking into account more realistic disc physics, in particular the effect of vertical stratification.To facilitate comparison to our unstratified, isothermal, simulations, here we also adopt an isothermal equation of state.
The structure of the paper is as follows.We outline our governing equations, numerical algorithms, and key parameters and diagnostics in Section 2. In Section 3 we present various results in physical space such as the scaling of turbulent intensity and magnetic field strength with magnetic Prandtl number, and also the vertical structure of the disc, including evidence for the Parker instability in the disc atmosphere.In Section 4 we perform a detailed spectral analysis of the turbulence, focusing on how energy transfers in spectral space contribute to the dominant dynamics in different vertical parts of the disc.Finally, we present our conclusion in Section 5.As some of our key results related to magnetic buoyancy (Parker instability) occur predominantly in the atmosphere of the disc, we also investigated the effects of changing the vertical boundary conditions and of increasing the vertical box size in Appendices A1 and A2, respectively.
Readers primarily interested in how turbulent transport and magnetic field strength scale with Pm should jump straight to Section 3, in particular Figure 1.Readers primarily interested in the key results of our spectral analysis and how these relate to the self-sustenance mechanism of the dynamo at different heights in the disc should jump to Section 4.3.

Governing equations
We work in the shearing box approximation (Hawley et al. 1995;Latter & Papaloizou 2017), which treats a local region of a disc as a Cartesian box located at some fiducial radius  =  0 and orbiting with the angular frequency of the disc at that radius Ω 0 ≡ Ω( 0 ).
A point in the box has Cartesian coordinates (, , ) along the radial, azimuthal/toroidal, and vertical directions, respectively.In this rotating frame, the equations of non-ideal MHD are (1) with the symbols taking their usual meanings.We close the system with the equation of state for an isothermal gas  =  2   where  2  is the constant sound speed.
All our simulations are vertically stratified and the effective gravitational potential is embodied in the tidal acceleration g eff = 2Ω 2 0 e  − Ω 2 0 e  (third term on the right-hand side of Equation 2), where  is the dimensionless shear parameter  ≡ −  ln Ω/ ln  | = 0 .For Keplerian discs  = 3/2, a value we adopt throughout this paper.
To control the magnetic Prandtl number (see Section 2.2) we employ explicit diffusion coefficients.The viscous stress tensor is given by T ≡ 2S, where  is the kinematic viscosity, and S ≡ (1/2) [∇u + (∇u) T ] − (1/3)(∇ • u)I is the traceless shear tensor (Landau & Lifshitz 1987).The explicit magnetic diffusivity is denoted by : it is related to the resistivity  via  ≡ / 0 , where  0 is the permeability of free space (note that from now on we will use the terms resistivity and magnetic diffusivity interchangeably).Note that in a real disc the microscopic viscosity and resistivity (and thus the magnetic Prandtl number) will depend on temperature and density (see Rossi et al. (2008) and Kawanaka & Masada (2019) for the explicit dependence of Pm on temperature and density in a neutrino-cooled disc).However, as the simulations discussed in this paper are isothermal, we keep the viscosity and resistivity fixed in space and time in any given simulation.

Important parameters
The magnetic Reynolds number compares inductive to resistive effects and is given by where   is the isothermal speed of sound,  =   /Ω 0 is the scaleheight (see Section 2.3.3 for definitions), and  is the magnetic diffusivity.
The Reynolds number compares inertial to viscous forces and is given by where  is the kinematic viscosity.Finally, the ratio of Rm to Re defines the magnetic Prandtl number, which serves as the key control parameter in our simulations.

Code
For our simulations we use the conservative, finite-volume code PLUTO (Mignone et al. 2007).We employ the HLLD Riemann solver, 2nd-order-in-space linear interpolation, and the 2nd-orderin-time Runge-Kutta algorithm.In addition, in order to enforce the condition that ∇ • B = 0, we employ Constrained Transport (CT), and use the UCT-Contact algorithm to calculate the EMF at cell edges (Gardiner & Stone 2005).To allow for longer time-steps, we take advantage of the FARGO scheme (Mignone et al. 2012).When explicit resistivity  and viscosity  are included, we further reduce the computational time via the Super-Time-Stepping (STS) scheme (Alexiades et al. 1996).Ghost zones are used to implement the boundary conditions.
We use the built-in shearing box module in PLUTO (Mignone et al. 2012).Rather than solving Equations ( 1)-(3) (primitive form), PLUTO solves the governing equations in conservative form.

Initial conditions
All our simulations are initialized from an equilibrium exhibiting a Gaussian density profile: where  0 is the mid-plane density at initialization, and  0 is the scaleheight at the mid-plane at initialization (formally defined below).
The background velocity is given by u 0 = −Ω 0  e  .At initialization we usually perturb all the velocity components with random noise exhibiting a flat power spectrum.The perturbations u have maximum amplitude of about 5 × 10 −2  0 , unless stated otherwise.Here  0 is the sound speed at initialization.All simulations are initialized with commonly used zero-net-flux (ZNF) magnetic field configuration B 0 =  0 sin (2/  )e  , where   is the radial box size.We define the field strength at initialization  0 through the plasma beta parameter  0 ≡ 2 0  0 / 2 0 , where  0 =  0  2 0 is the pressure at the mid-plane.We set  0 ≡ 1000 in all our simulations.

Units
Note that from this point onwards, all quantities are given in terms of dimensionless (code) units.Time units are selected so that Ω 0 = 1.The length unit is chosen so that the initial sound speed  0 = 1, which in turn defines a reference scale-height  0 ≡  0 /Ω 0 = 1.Finally the mass unit is set by the initial mid-plane density, which is  0 = 1.Magnetic field is expressed in units of  0 √  0  0 .Pressure, stresses, and energy densities are expressed in units of  2 0  0 .Note that we will occasionally drop the subscript on Ω 0 ,  0 , and  0 in the text from this point onwards.

Box size and resolution
The majority of our simulations are run at a resolution of   ×   ×   = 512×512×1024 in a box of size [  ,   ,   ] = [4, 4, 8] (i.e.128 cells per scale-height).In our unstratified paper we found that a resolution of 128 cells/ was sufficient to resolve the resistive scale at Rm = 18750 (see Appendix A of Held & Mamatsashvili (2022)).We have also run select simulations at lower resolutions of 32 cells/ and 64 cells/.To investigate the dependence on the vertical box size, we have repeated our fiducial run (Pm = 4, box size [4, 4, 8]) in a taller box of size [4, 4, 10] (keeping the resolution per scale-height fixed at 128 cells per ).We discuss this box size study in Appendix A2.All the simulations described in this paper are listed in Tables B1-B3 in Appendix B.

Boundary conditions and mass source term
We use standard shear-periodic boundary conditions (BCs) in the -direction (see Hawley et al. 1995), and periodic boundary conditions in the -direction.In the vertical direction, we keep the ghost zones associated with the thermal variables in isothermal hydrostatic equilibrium, in the manner described in Zingale et al. (2002).For the velocity components we use standard outflow boundary conditions in the vertical direction, whereby the vertical gradients of all velocity components are zero (numerically, we set variables in the ghost zones equal to those in the active cells bordering the ghost zones).For the magnetic field we employ 'vertical field' boundary conditions, also known as 'pseudo-vacuum' boundary conditions.Explicitly these are defined as setting   = 0,   = 0, and   / = 0 at the vertical boundaries.
As we are interested in the dynamics in the disc atmosphere || > 2, we also investigate the effect of using different vertical boundary conditions for the magnetic field.Altogether we explored three different types of boundary condition in the vertical direction: vertical field zBCs, outflow zBCs, and perfect conductor zBCs.These results are described in Appendix A1 and the corresponding simulations are listed in Table B3 in Appendix B.
Finally, to prevent mass-loss through the vertical boundary from depleting the mass in the box, we employ a simple mass source term.(This mimics what occurs in a real disc or in a global disc simulation, where mass lost from any given annulus through outflows is replenished by accretion of material from a neighboring annulus.)At the end of the th step, we subtract the total mass in the box at the end of that step   from the total mass in the box at initialization  0 .This mass difference Δ  ≡  0 −   is added back into the box with the same profile used to initialize the density (cf.Equation 7).Thus the total mass in the box remains constant in any given simulation.

Diagnostics
Below we define various diagnostics in physical space.For diagnostics in Fourier space the reader should refer to Section 4.1.

Averaged quantities
The volume-average of a quantity  is denoted ⟨⟩ and is defined as where  is the volume of the box.Note that occasionally we average only over a part of the box instead of the entire domain, for example || < 2 to capture diagnostics in the 'bulk' of the disc, or || > 2 to capture diagnostics in the 'atmosphere' of the disc.In this case the region we average over is stated in the text.We are also interested in averaging certain quantities (e.g.magnetic energy density or turbulent stresses) over time.The temporal average of a quantity  is denoted ⟨⟩  and is defined as where we integrate from some initial time   to some final time   and Δ ≡   −   .
The horizontal average of a quantity  is denoted ⟨⟩   and is defined as Horizontal averages over different coordinate directions (e.g. over the -and -directions) are defined in a similar manner.

Reynolds and magnetic stresses and transport 𝛼-parameter
In accretion discs, the radial transport of angular momentum is related to the -component of the total stress in which    ≡     is the Reynolds stress, where   ≡   + Ω is the perturbation of the y-component of the total velocity   about the background Keplerian flow  0 = −Ω 0  and    ≡ −    is the magnetic (Maxwell) stress.Note that since we will exclusively refer to fluctuating part of velocity u below, from now on we will drop the  and simply refer to the perturbed as u, which should not be confused with the total velocity.The total stress is related to the classical dimensionless angular momentum transport parameter .This can be defined either by normalizing the total stress by the volume-averaged gas pressure ⟨⟩ or alternatively by normalizing by the mid-plane pressure at initialization where we remind the reader that, in code units, the mid-plane sound speed and mid-plane density at initialization are simply  0 = 1 and  0 = 1.Note that due to the effects of vertical stratification (for which the pressure, and therefore density, decreases monotonically from the mid-plane), averaging alpha over the entire box can result in alpha defined by Equation 12being as much as 1.5-3 times that defined by Equation 13.Different authors use different definitions for  in the literature, and so caution is needed when comparing values from different sources (e.g., Pessah et al. 2008;Held & Latter 2018).

Magnetic buoyancy
To determine whether (and where) the magnetic field becomes buoyantly unstable, we calculate the vertical profile of the square of the magnetic buoyancy frequency  2  , which we define as follows (Shi et al. 2010): where Regions where  2  < 0 are unstable to the Parker instability, which In the bottom panel we distinguish between  (solid triangles, cf.Equation 12).and  0 (empty triangles, cf.Equation 13).Note that in the stratified simulations, spatial averages were taken over the bulk of the disc (|  | < 2), only.Unstratified (stratified) simulations were run in boxes of size 4 × 4 × 4 and 4 × 4 × 8, respectively.All simulations were run at a resolution of 128 cells-per-scale-height  and at at fixed magnetic Reynolds number of Rm = 18750.The abscissa is in log scale to base 4.
causes toroidal magnetic field lines to rise and escape from disc to the atmosphere.

RESULTS IN PHYSICAL SPACE
In this section we present various diagnostics from our vertically stratified simulations of MRI-turbulence in the large magnetic Prandtl number regime (Pm ≥ 4).We restrict our attention here to physical space, and defer a spectral analysis of the results to Section 4. Readers who are primarily interested in how turbulent transport and magnetic field strength scale with the magnetic Prandtl number may wish to focus their attention on 3.1 (in particular Figures 1 and 2), only.
All the simulations described in this section have been run in a box of size [  ,   ,   ] = [4, 4, 8], and at a resolution of   ×   ×   = 512 × 512 × 1024, corresponding to 128 cells per scale-height  in all three directions.The magnetic Reynolds number Rm is fixed at Rm = 18750 across all runs (this value stems from a resolution study we carried out in our unstratified runs: it is high enough to achieve Reynolds numbers in the hundreds at large Pm, but low enough for the resistive scale to be well-resolved.See Appendix A of Held & Mamatsashvili (2022)).We increase the magnetic Prandtl number from simulation to simulation by decreasing the Reynolds number Re. 2 2 We run each simulation for 200 orbits, which gives the turbulence ample time to settle into quasi-steady state.All timeaverages are taken between orbit 100 and orbit 200.For the Pm = 4 case, we have checked that the key results are qualitatively similar regardless of choice of vertical boundary conditions or of the vertical box size (see Appendix A1 and A2, respectively).All the simulations are listed in Table B1 in Appendix B.

Scalings of 𝛼 and 𝐸 𝑚𝑎𝑔 with Pm
The key results are plotted in Figure 1, where we show the dependence of the (time-and volume-averaged) total saturated magnetic energy  mag = (1/2) 2 and turbulent transport  parameter on magnetic Prandtl number.The green stars show the results of our previous unstratified simulations (Held & Mamatsashvili 2022) at fixed Rm = 18750 in boxes of size 4 × 4 × 4 and identical resolution 128 cells per scale-height.All other data points are from our new, vertically stratified, simulations.Note that for the stratified simulations (denoted by triangles) we have restricted volume-averages to the bulk of the disc (|| < 2).In the bottom panel of this figure, for comparison we also show the two variants of  defined in Section 2.4.2:empty triangles represent  (the total stress normalized by the volume-averaged pressure ⟨⟩; cf.Equation 12), while solid triangles represent  0 (the total stress normalized by the mid-plane pressure at initialization  0 ; cf.Equation 13).Due to the stratification (vertical pressure gradient),  is larger than  0 by a factor of around 1.70.Note that in unstratified boxes  is also sensitive to the box size and vertical-to-radial box aspect ratio, which we discuss in more detail in our unstratified work (Held & Mamatsashvili 2022).
When we restrict our averages to the bulk of the disc (|| < 2), which contains 95% of the total mass in the box, the results of the stratified simulations for both  mag -Pm and -Pm scaling are in excellent agreement with those of our earlier unstratified simulations.At intermediate values of Pm we find power law scaling with slopes of  ∼ 0.74 and  ∼ 0.71 for each quantity, respectively, in agreement to two decimal places with the results of the unstratified runs, while for Pm ≳ 32 we observe the onset of a plateau.Thus, we find that in the bulk of the disc stratification does not appreciably change dependence of the turbulent dynamics on the magnetic Prandtl number.As we show below, this conclusion changes as we move away from the bulk of the disc where the MRI dominates the dynamics, and into into the atmosphere of the disc where another process due to magnetic buoyancy -Parker instability -comes into play.

Vertical dependence of scaling
A new result is that the scaling of turbulent transport and saturated magnetic energy with magnetic Prandtl number is stronger in the bulk of the disc than in the atmosphere, as seen in Figure 2, where we plot the time-and horizontally-averaged magnetic energy density ⟨⟨ mag ⟩   ⟩  .Different colors correspond to horizontal-averages taken at fixed || = {0, 1, 2, 3, 3.5}, respectively.The scaling clearly becomes weaker further away from the mid-plane: for example we measure a slope of  ∼ 0.84 at the mid-plane (|| = 0), but a much more shallow slope of  ∼ 0.36 at || = 3.5.Thus our results foreshadow that different dynamics/dynamo processes (mainly MRI and Parker instability) are at play in different parts of the disc, as discussed in greater detail in Section 3.2 below.Note that these processes have also been observed in earlier stratified simulations (Johansen & Levin 2008;Blaes et al. 2007;Shi et al. 2010;Blaes et al. 2011;Kadowaki et al. 2018), although those simulations were ideal and initialized with relatively strong ( 0 ∼ 1 − 25) net-toroidalmagnetic flux (which aids the development of the Parker instability), in contrast to our non-ideal, zero-net-flux runs with relatively weak initial magnetic field ( 0 = 103 ).We analyze the dynamics of these instabilities in more detail in Fourier space in Section 4.

Dependence of results on Pm
To determine whether the MRI is well-resolved we also measured the MRI quality-factor in the Pm = 4 run.We define this diagnostic as   ≡  MRI,z /Δ, where Δ is the grid size in the vertical direction. 3 At the mid-plane we find   ∼ 26 (  ∼ 0.2), while  = 2 we find   ∼ 68 (  ∼ 0.5).Thus the critical MRI wavelength is very well resolved within the bulk of the disc, and also fits comfortably in the box.

Structure of the flow
Next we consider the flow field in the -plane, which is shown Figure 3 from a simulation characteristic of the power-law scaling region (Pm = 4, top row) and also from a simulation characteristic of the plateau region (Pm = 90, bottom row).The density is shown in the left-hand panels: due to the inclusion of stratification, the bulk of the disc is concentrated between ±2, surrounded by a tenuous atmosphere.As in our unstratified runs, we observe density waves superimposed on the turbulence at lower Pm but not at higher Pm (Held & Mamatsashvili 2022;Shi et al. 2016).
The magnetic field is shown in the middle and right-hand panels of Figure 3 for   and   , respectively.A small-scale, turbulent magnetic field dominates the flow within ± of the mid-plane, while in the atmosphere (|| ≳ 2.2) the field is still turbulent, but characterized by noticeably larger structures.Further into the atmosphere still, the character of the magnetic field changes and we observe only large-scale structures (of size ∼ 0.8).In particular,   is characterized by streaky, vertically-oriented structures.This indicates that in the disc atmosphere the magnetic field structure is not really characteristic of MRI turbulence, but rather an ordered, vertical field has formed.This is due to magnetic buoyancy, specifically the Parker instability, which dominates MRI in the atmosphere, as we show in Section 3.2.4.Finally, we observe larger structure in both the velocity (not shown) and magnetic fields as Pm is increased, as expected given the low Reynolds number (Re = 208) at Pm = 90. 4Reassuringly, the flow field near the mid-plane at both values of Pm resembles that observed in our unstratified simulations (see Figure 3 of (Held & Mamatsashvili 2022))

Vertical disc structure
To characterize the vertical structure of the disk, we look at diagnostics such as the time-and horizontally-averaged vertical profiles of density, plasma beta, and  at different Pm (not shown).As we increase Pm, we observe the largest changes to the disk structure within the bulk of the disk, as foreshadowed in Figure 2. In the power-law scaling region at low Pm, the density remains nearly unchanged from  its initial Gaussian profile, while the plasma beta parameter drops from  ∼ 55 at the mid-plane to  ∼ 5 at || ∼ 4.The entire disc remains gas-pressure-dominated at all values of Pm that we investigated, although magnetic pressure contributes significantly to the total pressure as Pm is increased ( ≲ 10 at the mid-plane at Pm = 90).Meanwhile, the  parameter is largest at the mid-plane and drops monotonically away from the mid-plane.The profiles of , plasma beta, and magnetic field are relatively uniform within the bulk of the disk at Pm = 4, but becomes increasingly peaked around  = 0 as Pm increases.The latter behavior is likely due to the increased importance of buoyancy in the disc region as Pm is increased, which results in the buoyantly unstable region extending closer to the mid-plane (see Figure 6).
In Figure 4 we show vertical profiles of all horizontally-averaged root-mean-square (rms) velocity and magnetic field components at Pm = 4.The profiles were time-averaged between around orbit 100 and orbit 200.At the mid-plane the dominant velocity components are the radial and azimuthal components (top panel).All velocity components increase away from the mid-plane (essentially scaling with the Alfvén speed   ≡ / √ ), and at || ≳ 3 the vertical component   (black curve) increases rapidly, quickly becoming the dominant component.This suggests that outflows dominate the dynamics in the upper reaches of the disc, probably due to the emergence of the large-scale vertical field   in this region, as seen in its -slice in the top-right-hand panel of Figure 3. Turning to the magnetic field vertical profiles (bottom panel of Figure 4), the toroidal component   (red curve) dominates over the radial (blue curve) and vertical (black curve) components at nearly all .The toroidal field decreases monotonically with height away from the mid-plane, while the radial and vertical field are relatively constant within the bulk of the disc.In the atmosphere the toroidal and radial field drop rapidly, however, becoming comparable to the vertical field around ±3.This suggests a magnetically-driven outflow dominates the dynamics away from the bulk of the disc even in the present case of zero-net-vertical magnetic flux, an outcome that otherwise generally requires sufficiently strong net-vertical-flux (Lesur 2021).

Spacetime diagrams
To examine the temporal behaviour of the magnetic field in more detail, we next turn to spacetime diagrams of toroidal field   (see Figure 5), which are shown between orbit 360 and orbit 400 from our low resolution runs (32 cells per scale-height ).The top panel is from a simulation at Pm = 4, while the bottom panel is from a simulation at Pm = 90.We observe the characteristic 'butterfly' pattern in   , which has been observed in previous stratified MRIturbulence simulations (e.g., Davis et al. 2010;Gressel 2010;Simon et al. 2011;Bodo et al. 2014;Salvesen et al. 2016).The toroidal field changes sign every 5-10 orbits, and this period does not appear to be sensitive to the magnetic Prandtl number.The field reversals are most evident outside ± of the mid-plane.Closer inspection reveals that field reversals are already present very close to the mid-plane at low Pm, where the disc is otherwise highly turbulent, but the butterfly pattern inside the disc || < 2 is disrupted at Pm = 90.

Parker instability
To investigate whether the disc is unstable to the Parker instability we plot the square of the magnetic buoyancy frequency  2  (Equation 14).We find that, within the bulk of the disc (|| ≲ 2),  2  ≈ 0 (bottom panel of Figure 6), and thus the bulk of the disc is marginally stable to Parker instability.Outside this region  2  < 0, and increases in absolute value with height.This shows that the atmosphere of the disc is magnetically buoyantly unstable.
Another characteristic of the Parker instability is the bending (undulation) of magnetic field lines, and associated over-and underdensities in the gas density as plasma slides along the field lines away from the crests and collects in the troughs.See, for example, Figure 1 of Johansen & Levin (2008) which is taken from a 2D simulation in an azimuthally extended box with strong net-toroidal-flux.In our 3D zero-net-magnetic-flux simulations, toroidal field is produced mostly by the shear and the flow field is much more complex than in 2D simulations.Thus such clear undulations of the field lines with over-densities in the troughs of the undulations is difficult to detect.Nevertheless, we have observed many examples of bending (and rising) field lines in our simulation.In Figure 6 we show an example from one snapshot at orbit 135.The colorplot shows the density in the -plane between  = 2.65 (upper edge of the disc) and  = 3.3.Note that in this part of the domain the square of the magnetic buoyancy frequency is negative, and so we expect the flow to be Parker unstable.The black lines are the magnetic field lines, while the yellow arrows indicate the velocity field.An upwelling of the magnetic field can clearly be seen between  = −0.5and  = 2.The velocity (yellow arrows) points along the magnetic field streamline on either side of the undulation, and two overdense regions can be seen in the troughs.Thus we have direct visual evidence of Parker instability in our vertically stratified zero-net-flux shearing box simulations.

SPECTRAL ANALYSIS
To understand the energetics and dynamics that underpin the turbulence, following Held & Mamatsashvili (2022), we now present a detailed spectral analysis of the results first discussed above in physical space.We begin by transforming the induction equation in Fourier (wavenumber) space in Section 4.1.In Section 4.2 we present a spectral analysis at fixed Pm = 4, including energy spectra and the spectra of the dynamical terms governing magnetic field evolution at different heights in the disc.In Section 4.4 we discuss how the results change with the magnetic Prandtl number.The reader who is only interested in how the main physical processes (MRI and Parker instability) sustain the turbulence/dynamo should jump to the summary in Section 4.3.

Governing equations
We start the spectral analysis by decomposing velocity and magnetic field perturbations into spatial Fourier modes in the radial − and toroidal −coordinates, along which the perturbations are shearperiodic and periodic, respectively.Along the vertical -coordinate, there is no periodicity due to the presence of stratification, so we do not Fourier transform in .Thus, we have (Riols et al. 2017a) where  ≡ (u, B) and f ≡ ( ū, B) denote the corresponding Fourier amplitudes.The grid in Fourier (k-)space is determined by the horizontal sizes of the flow domain   and numerical resolution   , where  ∈ {, }, such that the cell sizes in Fourier space are given by Δ  = 2/  , and hence the wavenumbers run through the values   =   Δ  , where   = 0, ±1, ±2, ..., ±  /2 is an integer.Performing a Fourier decomposition in horizontal slices at each height  allows us to explore different dynamical regimes discussed in Section 3 -MRI-turbulence dominating in the bulk of the disc at || ≲ 2 and a magnetic buoyancy, or Parker instability regime dominating in the upper layers || ≳ 2. 5  Substituting Equation ( 15) into induction Equation (3), we obtain the equations for the spectral magnetic field components, 5 Note, we refer here to which process dominates the dynamics: the MRI is present at all , and is thus still active in the atmosphere.
where in these and other spectral equations below we have chosen to denote Fourier transforms of the products of any two quantities with the subscript k instead of over-bars (which should not cause confusion) and Δ k = − 2  −  2  + 2 / 2 is the 2D Laplace operator in the (  ,   )-plane.The Fourier transform of the nonlinear terms -the products of velocity and magnetic field components (    ) k , where ,  ∈ {, , }, come from the Fourier transform of the electromotive force u × B in Equation ( 3) and are given by convolutions in Fourier space (Mamatsashvili et al. 2020) They describe the net effect of the nonlinear triadic interactions of a mode with k with two other modes k ′ and k − k ′ .Note that stratification does not explicitly appear in the induction equation ( 3).However, it influences the vertical velocity   in the momentum equation ( 2) and, via the latter, the magnetic field dynamics.Thus, to characterize the effect of stratification, we have separated out the contributions of   and   in the nonlinear terms in the magnetic field Equations ( 16)-( 18).Multiplying both sides of these equations by the complex conjugates B *  , B *  , and B *  , respectively, we obtain the equations for the squared moduli of the spectral magnetic field components, which are central in the spectral analysis.Their right-hand sides contain terms of linear and nonlinear origin.The linear terms are: (i) the shear-induced drift, −  /  , which simply advects the field for non-axisymmetric (  ≠ 0) modes along the   -axis without any new energy production, (ii) The Maxwell stress multiplied by the shear parameter , which is responsible for the energy exchange between (toroidal) magnetic field and the background disc flow.In this case, the Maxwell stress mediates the growth of toroidal field via stretching of radial field by Keplerian shear, which is a main part of the (non-modal) MRI process (Herault et al. 2011;Squire & Bhattacharjee 2014;Mamatsashvili et al. 2013;Gogichaishvili et al. 2017;Riols et al. 2017b).It is the main term supplying (injecting) energy into turbulence.
(iii) Ohmic dissipation terms, (iv) The nonlinear transfer terms are:

, and N (𝑣𝑎) 𝑦
, where  ∈ , , , which correspond to the nonlinear terms in the original equations (Equations 16-18).As mentioned above, they are classified into the terms explicitly containing the vertical velocity   and magnetic field   , and the terms independent of these components.The explicit forms and physical meanings of these nonlinear terms are as follows: The nonlinear induction-advection terms depending only on horizontal velocity, (  ,   ), and magnetic field (  ,   ) are: which describe, respectively, the production of B from B and vice versa due to stretching by horizontal velocity variation (shear) along -and -directions (the first terms in the brackets) and their horizontal advection/transport (second terms in the brackets).
The nonlinear terms containing vertical field   are: which describe, respectively, the production of B and B from B due to stretching by the horizontal velocity variation (shear) along vertical -axis, and which describes the advection of vertical field by horizontal velocity.Finally, the nonlinear terms explicitly depending on   are: which describe the vertical advection of the horizontal field components B and B , respectively, by   , hence characterizing the effect of buoyancy on the magnetic field, and which describes production of the vertical field B from the horizontal field ( B , B ) due to stretching of the latter by horizontal shear of   .The spectral equations ( 19)-( 21) are similar to those in the unstratified case in Mamatsashvili et al. (2020) and Held & Mamatsashvili (2022), except that here we have isolated the nonlinear terms dependent on   and   and hence affected by stratification.

Fiducial simulation at fixed Pm = 4
Following our previous studies on the spectral dynamics of MRIturbulence (Gogichaishvili et al. 2017;Mamatsashvili et al. 2020;Held & Mamatsashvili 2022), we first focus on the fiducial simulation Pm4Rm18750Re4687 with Pm = 4 and perform a detailed analysis of the energy spectra and the dynamical terms in the (  ,   )-plane both near the disc mid-plane at  = 0 and in the atmosphere at  = 3.2 in order to understand the sustaining dynamics of turbulence/dynamo at different heights.Then, we explore how these spectra and the dynamical balances change with Pm.Below, we average all the spectral quantities in time from 100 to 200 orbits when the turbulence is already well in a quasi-steady state.We also average over a small height interval Δ = 0.4, around both  = 0 and  = 3.2, the two locations at which we perform a spectral analysis.
Using the spectral analysis we can study, apart from other modes, the dynamics of the special   =   = 0 mode, referred to as the dynamo mode, because in physical space it in fact represents a solution that is uniform in  and  and depends only on .The dynamo mode gives rise to the large-scale horizontally-averaged magnetic field, which exhibits the well-known "butterfly"-like variation in the (, )-plane in vertically stratified discs (see Figure 5 and, e.g.Davis et al. 2010;Shi et al. 2010;Gressel 2010;Simon et al. 2011;Guan & Gammie 2011;Salvesen et al. 2016).This dynamo field has usually been studied in physical space using a conventional mean-field  − Ω approach (e.g., Brandenburg et al. 1995; Johansen & Levin 2008;  Gressel 2010; Gressel & Pessah 2015; Shi et al. 2016; Brandenburg  2018; Gressel & Pessah 2022), 6 In this case, the details of nonlinear mode interactions responsible for the sustenance and dynamics of the dynamo mode and the associated large-scale magnetic field are packed into  and  tensorial parameters.On the other hand, the spectral approach that we use in this section, i.e. analysing in detail the nonlinear energy transfers in Fourier space, provides deeper insights into the self-sustenance mechanism of the dynamo in stratified discs.due to buoyancy (right column) in the (  ,   )-plane.The bottom row shows these terms at  = 0 (the mid-plane), while the top shows the terms at  = 3.2 (the disc atmosphere).The spectra were taken from a simulation run at Pm = 4 and Rm = 18750 with a resolution of 128/.At the mid-plane (bottom row), the dominant process is the production of the radial field from the toroidal field due to the non-linear transverse cascade, which is described by the positive values of  (ℎ)  > 0 (yellow/red areas in the left-hand column).The radial field is also produced from the vertical field, which is described by the positive values of  ()  > 0 (yellow/red areas in the middle column).By contrast, the vertical advection term is negative,  ()  < 0, at small wavenumbers (blue areas in the right-hand column), implying removal of large-scale field by vertical buoyancy.In the atmosphere (top row), on the other hand, the efficiency of the transverse cascade is relatively small and the radial field is mainly produced from the vertical one.

Energy spectra
The time-averaged magnetic energy spectrum   in the (  ,   )plane is shown in Figure 7 at  = 0 and 3.2 (the kinetic energy spectrum   has a similar structure and is not shown here).It has a typical anisotropic structure in the (  ,   )-plane due to the background shear, with nearly the same inclination towards the  axis at both heights.Such an anisotropic spectrum is typical of MRI-turbulence and has also been observed in spectral analyses of unstratified MRI-turbulence (Lesur & Longaretti 2011;Murphy & Pessah 2015;Gogichaishvili et al. 2017;Held & Mamatsashvili 2022).With increasing , the spectrum becomes more concentrated at smaller wavenumbers, i.e., larger-scale modes become dominant over smaller-scale ones (see also the 1D spectra in Figure 11).
This behavior of   with height is related to different characteristic (correlation) length-scales of the turbulence near the mid-plane and in the atmosphere.Near the mid-plane this length   ∼   /Ω 0 (Walker et al. 2016), which is much less than the scale-height, since   ≪   at  = 0 and hence the spectrum extends to higher wavenumbers ∼ 1/  ≫ 1/.On the other hand, in the buoyancydominated atmosphere, the dynamics is mainly governed by Parker instability, whose characteristic length-scale is set through the gravitational acceleration  = Ω 2 0 , Alfvén   and sound   speeds, (Parker 1967;Blaes et al. 2007) where it has been assumed that in the atmosphere   ∼   .Therefore, the energy spectra at  = 3.2 extends to wavenumbers ∼ 1/  much smaller than 1/  at the midplane, since the ratio   /  ∼   / ,=0 ≫ 1.This means that turbulent structures are of much smaller scale at the midplane than in the atmosphere consistent with the spatial distributions of the variables in Figure 3 (see also Blackman & Pessah 2009).

Spectra of the dynamical terms
We now turn to the analysis of the spectral dynamics of turbulence and consider each magnetic field component separately: The dynamics of the radial field   is governed by the nonlin- , and N () , which are shown in Figure 8 in the (  ,   )-plane both at the disc mid-plane  = 0 and in the atmosphere  = 3.2.The radial component plays a central role in the self-sustaining process of MRI-turbulence since at all heights, the energy-carrying toroidal field can be produced only from the radial field via stretching by Keplerian shear (see below).At the midplane, like the energy spectra, all the nonlinear terms are strongly anisotropic in Fourier space due to the Keplerian shear, depending on the polar angle  = (  /( 2  +  2  ) 1/2 ) of the horizontal wavevector k.The main consequence of this anisotropy for N (ℎ)  is the redistribution of power over wavevector orientation (-angle) in the (  ,   )-planethe nonlinear transverse cascade (Mamatsashvili et al. 2014;Gogichaishvili et al. 2017;Mamatsashvili et al. 2020), which transfers the radial field energy, | B | 2 /2, from "giver" wavenumbers for which N (ℎ)  < 0 (blue) to "receiver" wavenumbers for which N (ℎ)  > 0 (yellow and red), as seen in the left-hand column of Figure 8.We refer to this region in spectral space encompassing  (see equation 24), so the role of the nonlinear transverse cascade is in fact to continually replenish and amplify the radial field from the toroidal one at the "receiver" wavenumbers where N (ℎ)  > 0. In the atmosphere, the vital area in the (  ,   )-plane shrinks to smaller wavenumbers and the transverse cascade weakens.
The radial field is also generated from the vertical field via N ()  at those wavenumbers where it is positive, N ()  > 0 (yellow/red areas in the middle column of Figure 8).At the mid-plane, this process is less intensive but still comparable to the transverse cascade, whereas it dominates the latter in the atmosphere.However, the situation is different for the k = 0 dynamo mode, where N .From these four terms the dominant ones are M and the vertical transport term N ()  which are shown in Figure 9.They exhibit a similar anisotropy due to the shear as do their counterparts for the radial field.In the vital area, only the Maxwell stress is positive, M > 0 (left-hand column) and hence amplifies (injects) toroidal field energy | B | 2 /2 by stretching the radial field due to shear.At the mid-plane, the amplification of B is due to the non-modal MRI-growth process and hence spans a much broader range of wavenumbers than in the atmosphere, where it is related to Parker instability and instead is concentrated at the smallest wavenumbers peaking around k = 0.In the vital area, N ()  < 0 (right-hand column), draining the largescale toroidal field energy due to buoyancy (see also Blackman & Tan (2004) and Blackman & Pessah (2009)).
The dynamics of the vertical field   is governed by the nonlinear terms N (ℎ)  and N ()  , which are shown in Figure 10, having a similar anisotropic structure in the (  ,   )-plane due to the shear.
In the vital area, N (ℎ)  is positive at intermediate wavenumbers (lefthand column), and generates the vertical field from the horizontal (mainly toroidal) one.This process is, however, due to two different sources: non-modal MRI process near the mid-plane at  = 0 (where buoyancy is weak) and Parker instability in the atmosphere, where magnetic buoyancy, having  2  < 0, plays a major role (Figure 6).Consequently, at  = 3.2, the term N (ℎ)  operates at wavenumbers comparable to 1/  , which corresponds to the length-scale of the Parker instability   ∼  2  / as given by Equation (32) above.We estimated this length to be   ∼ 0.3 from the simulations and showed for reference as a circle of radius 1/  ∼ 3 in the map of N (ℎ)  in Figure 10, which covers the location of the peaks of this term.At k ≈ 0, N (ℎ)  = 0, implying that there is no production of the large-scale vertical field.By contrast, the horizontal advection term  Finally, the drift terms due to shear  in equations ( 19)-( 21), do not produce (inject) energy into the modes, and serve only to balance the joint action of the dynamical terms in quasi-steady state, so we do not show the drift terms here.They cause the radial wavenumber   of non-axisymmetric (  ≠ 0) modes to change in time and cross the vital area.As a result, the growth of these modes acquires a transient, or non-modal character (Balbus & Hawley 1992;Mamatsashvili et al. 2013;Squire & Bhattacharjee 2014).7

Summary of dynamo self-sustenance mechanisms at different heights
We now summarize the self-sustenance processes of the turbulence and dynamo at different vertical locations in the disc (at the mid-plane  = 0, and in the atmosphere  = 3.2).Let us briefly remind the reader of our key diagnostics: we have Fourier transformed the induction equation, and so obtained various terms governing the evolution of each spectral magnetic energy component (see . By plotting the 2D spectra of these dynamical terms at the disc mid-plane and in atmosphere we can construct a picture of the key processes sustaining each field component (and thus the dynamo and turbulence) in different parts of the disc, which we describe below.

Dynamical processes in the bulk of the disc.
The dynamical balances sustaining the turbulence near the mid-plane ( = 0) primarily involve the radial and toroidal field components, since these two carry the most energy (see bottom panel in Figure 4).In this case, a kind of anisotropic inverse cascade induced by the shear -more accurately known as a nonlinear transverse cascadeplays a key role, seeding the radial field B from the toroidal field B at small to intermediate wavenumbers (i.e.|  | ≲ 30, |  | ≲ 20, at  = 0).This process is characterized by the nonlinear term in the radial field equation ( 19), which is anisotropic in the (  ,   )-plane due to shear, and which exhibits regions of "giver" wavenumbers where N (ℎ)  < 0 (blue regions in left-hand column of Figure 8) and "receiver" wavenumbers where N (ℎ)  > 0 (yellow and red regions).Physically, this results in a transfer of radial field energy (| B | 2 /2) over angles (i.e., transversely) across the (  ,   )plane, from a select set of larger ("giver") wavenumbers to another set of smaller ("receiver") ones.Thus, the nonlinear transverse cascade replenishes the radial field at those wavenumbers where N (ℎ)  > 0. At small wavenumbers, the radial field is also amplified through induction from the vertical field, since the nonlinear term responsible for this process is positive, N ()  > 0 (see middle column of Figure 8).The "receiver" modes then drift along the   -axis through the injection area, that is, the region in Fourier space where the Maxwell and magnetic   (red) energy spectra at Pm = 4 (solid) and Pm = 90 (dash-dot) at fixed Rm = 18750, both at the mid-plane  = 0 (bottom) and in the atmosphere  = 3.2 (top).As Pm increases, the kinetic energy spectra increases mostly at small  and becomes steeper due to increased viscosity.The spectral magnetic energy, on the other hand, increases at all  with Pm: at the mid-plane it shifts to lower , while in the atmosphere it neither changes shape nor shifts, but increases only in magnitude by about the same factor at all .
stress is positive (M > 0) and hence injects magnetic energy into perturbations (see left-hand column of Figure 9).In this case, the energy injection is mainly due to MRI, resulting in the amplification of the toroidal field energy, | B | 2 /2, during the time that the modes drift across the injection area. 8This amplified toroidal field, which provides the dominant contribution to N (ℎ)  (see Equation 24), in turn produces the radial field via the transverse cascade process, as explained above, thereby closing the main self-sustenance cycle.(All the other non-linear terms in the toroidal field energy equation (Equation 20) are negative at small and intermediate wavenumbers, where most of the energy resides, as seen in Figure 9, and thus act as sinks opposing the self-sustenance.)This self-sustaining scheme at the mid-plane (where buoyancy is unimportant) is similar to that found in our unstratified studies Mamatsashvili et al. (2020); Held & Mamatsashvili (2022). 8The growth of perturbations due to MRI during some finite time interval is also known as the non-modal MRI (Squire & Bhattacharjee 2014; Gogichaishvili et al. 2017).

Dynamical processes in the disc atmosphere.
In the disc atmosphere (|| ≳ 2), the shear-modified Parker instability (Foglizzo & Tagger 1995;Johansen & Levin 2008) dominates over the MRI, as seen in Figure 6 showing the increasingly negative magnetic buoyancy (squared),  2  < 0, and hence more intensive Parker instability with increasing .Unlike at the mid-plane where the toroidal and radial fields dominates, in the atmosphere it is the toroidal and vertical field components that carry the most energy, and hence play the main role in the self-sustenance of the dynamo.The radial field is smaller, though still important (Figure 4).In Figure 8 we show that at  = 3.2 the dominant process producing the radial field is induction from the vertical field described by N ()  > 0, while the transverse cascade is subdominant.This radial field is subsequently stretched out by the Keplerian shear and thus produces and amplifies the toroidal field via the positive Maxwell stress M > 0 (Figure 9).Finally, the toroidal field generates the vertical field for non-axisymmetric modes due to Parker instability (mediated by the positive nonlinear term N (ℎ)  > 0 in the vital area, see Figure 10), thereby closing the self-sustenance cycle in the atmosphere.Note that despite the significance of the Parker instability in the atmosphere, the role of shear is still important, as it maintains the dominant toroidal field (which then becomes Parker-unstable) from the radial one.
Thus, we have identified two distinct but co-existing self-sustaining processes for the dynamo in vertically stratified discs.The first process operates primarily in the bulk of the disc (|| ≲ 2), and is based on the interplay between the non-modal MRI, which amplifies the toroidal field   from the radial field   , and the non-linear transverse cascade, which provides positive feedback, producing radial field from toroidal field.The second process operates in the disc atmosphere (|| ≳ 2), and is based on the interplay between Parker instability (which produces vertical field   from toroidal field   ), non-linear transfers (which generate radial field   from the vertical field), and Keplerian shear (which stretches the radial field into the toroidal field).The self-sustaining process near the mid-plane extends over a much wider range of length-scales than self-sustaining process in the atmosphere: the latter process is localized to large scales (comparable to the size of the system).This picture is consistent with one of the first disc dynamo models proposed by Tout & Pringle (1992) based on the interplay between MRI and Parker instability (see also Blackman & Tan 2004;Gressel 2010).However, a key difference between our model and theirs is that here these two instabilities determine the dynamo action at different heights: the MRI dynamo dominates near the mid-plane, whereas the shear-modified Parker instability dominates in the atmosphere.

Dependence of the results on Pm
Having outlined the turbulence sustenance mechanisms at different heights in the disc, let us now analyse how the dynamical terms underlying these mechanisms vary with magnetic Prandtl number Pm, both at the mid-plane and in the atmosphere.To capture the dynamics in the two regions of interest in the regime of large Pm (Figure 1), we focus on Pm = 4 (start of the power-law scaling region) and Pm = 90 (plateau region).The Pm-dependence of MRI-turbulence and clarification of its physical nature in unstratified disc models were first elucidated in Mamatsashvili et al. (2020) and Held & Mamatsashvili (2022).Here, we generalize these studies to the vertically stratified case.For comparison, we compute 1D ring-averaged spectra of the energies and the dynamical terms, i.e.
∫ 2  0 (..) , as a function of (right-hand column) at Pm = 4 (solid lines) and 90 (dash-dot lines) at the mid-plane  = 0 (bottom row), and in the atmosphere  = 3.2 (top row), for Rm = 18750, at a resolution of 128/.With increasing Pm, all these terms generally increase in magnitude more significantly at the mid-plane than in the atmosphere.At the mid-plane they shift to smaller , whereas in the atmosphere they retain their shape (peaks) as a function of .The main nonlinear terms N sustaining the turbulence at  = 0 and  = 3.2, respectively, exhibit a similar trend: they become increasingly negative at 20 ≲  ≲ 50 as Pm increases from 4 to 90, indicating a decrease in the direct cascade, i.e. a decrease in the transfer of energy to higher  and hence loss by resistive dissipation, and this in turn results in an overall increase in the turbulence level with Pm.

Energy spectra
Figure 11 shows the ring-averaged kinetic and magnetic energy spectra (compensated by ) at Pm = 4 and 90.These spectra do not have a clear power-law form because of the strong anisotropy in the (  ,   )-plane (Figure 7), as also pointed out in Murphy & Pessah (2015).On a promising note, we do find pretty clear peaks in all the spectra (for the kinetic energy spectrum at Pm = 90 the peaks are shallow, but still visible), both at the mid-plane and in the atmosphere, whereas Nauman & Blackman (2014) did not find such peaks in their vertically stratified ZNF simulations (they employed the same box size as we have, but much lower resolution).This led them to conclude that they were not capturing the outer scale of the turbulence in the atmosphere, whereas it follows from the energy spectra in Figure 7 that our box does encompass the outer scale (see below) in the atmosphere.
The ring-averaged kinetic energy spectra (blue curves) increase with increasing Pm at lower  (large scales) both at the mid-plane and in the atmosphere, that is, the velocity of large-scale modes increases.The spectra reach peaks at  , = 11 at the mid-plane and  , = 6.4 in the atmosphere for Pm = 4, which move to smaller  , = 3.1 and  , = 4.7, respectively, for Pm = 90 (here the subscript  stands for 'maximum ', and 'd' and 'a' denote mid-plane and atmosphere, respectively).After the peak, the kinetic energy decreases with  more gently at Pm = 90 than at Pm = 4. Thus, the effect of increasing Pm is to "stretch out" the kinetic energy spectrum by extending the viscous range to smaller : the viscous The ring-averaged magnetic energy spectra (red curves), on the other hand, increase at all  as Pm increases, however, they do so differently at the mid-plane and in the atmosphere.At the mid-plane, smaller  (larger scales) gain more power than higher  (smaller scales), so that the spectral peak shifts to lower  with increasing Pm.On the other hand, in the atmosphere, the spectrum retains its shape (i.e., does not shift in ) as Pm rises from 4 to 90 and mainly increases only in magnitude by about the same factor (of around 3.3) at all .At Rm = 18750 used here, the resistive   ∼ √ Rm = 137 (Held & Mamatsashvili 2022) is larger than the viscous   .

Spectra of the dynamical terms
Next, we examine the 1D ring-averaged dynamical terms and their dependence on Pm, which are plotted in Figure 12 at Pm = 4 and Pm = 90 both at the mid-plane and in the atmosphere.It is seen in this figure that, like the energies, the dynamical terms generally increase in magnitude several times with increasing Pm (compare solid and dash-dot lines for Pm = 4 and Pm = 90, respectively), however, their behaviour with Pm differs at the mid-plane and in the atmosphere.
Dynamics at the mid-plane of the disc.At the mid-plane  = 0, the ring-averaged spectra of the dynamical terms shift to smaller  as Pm is increased (bottom row of Figure 12), exhibiting a similar trend found in our unstratified simulations (Mamatsashvili et al. 2020;Held & Mamatsashvili 2022).In those studies, we interpreted this scales back to large-scales anisotropically in Fourier space.This transverse cascade becomes increasingly efficient at reseeding radial magnetic field as Pm is increased.
In the disc atmosphere, on the other hand, the dynamo primarily involves the interplay of toroidal and vertical magnetic field components, and the dependence on Pm is weaker.The radial field is generated primarily from vertical field due to vertical gradients in the (perturbed) velocity.This field is then sheared out to produce toroidal field, and vertical field is generated primarily from toroidal field by means of magnetic buoyancy due to the Parker instability, thus closing the loop.This process is reminiscent of, though not entirely the same as, the idea put forward by Tout & Pringle (1992) which involved the interplay between MRI and Parker instability, though our analysis shows that, while both instabilities play a role in the dynamo, they dominate in different parts of the disc.
By taking into account the effects of buoyancy and vertical stratification in this work (the second in a series of papers on the MRI in the regime of large magnetic Prandtl number), we have moved one step closer to modeling the conditions under which the MRI occurs in certain types of astrophysical objects, such as discs from binary neutron star mergers and proto-neutron stars.We have shown the turbulent transport is very sensitive to magnetic Prandtl number, at least at intermediate values of Pm, which could render the flow thermally unstable (Balbus & Henri 2008;Potter & Balbus 2017;Kawanaka & Masada 2019).Thus future work should take into account the effects of thermodynamics, neutrino cooling, and temperature and and density-dependent diffusion coefficients, in order to probe the thermal stability of discs in the large-Pm regime.

Figure 1 .
Figure 1.Top: scaling of time-and-spatially averaged magnetic energy density with magnetic Prandtl number Pm in isothermal shearing box simulations of the MRI with zero-net-magnetic flux.Bottom: scaling of alpha viscosity parameter with Pm.The symbols have the following meanings: green stars denote unstratified simulations from (Held & Mamatsashvili 2022), orange triangles denote vertically stratified simulations from this work.In the bottom panel we distinguish between  (solid triangles, cf.Equation12).and  0 (empty triangles, cf.Equation13).Note that in the stratified simulations, spatial averages were taken over the bulk of the disc (|  | < 2), only.Unstratified (stratified) simulations were run in boxes of size 4 × 4 × 4 and 4 × 4 × 8, respectively.All simulations were run at a resolution of 128 cells-per-scale-height  and at at fixed magnetic Reynolds number of Rm = 18750.The abscissa is in log scale to base 4.

Figure 3 .
Figure 3. Structures of density  (left column) and magnetic field components   (middle column) and   (right column) in the -plane from two simulations at Pm = 4 (top row) and Pm = 90 (bottom row).The magnetic Reynolds number is fixed at Rm = 18750 in all simulations.

√
)  (see footnote 6 ofHeld & Mamatsashvili (2022) for a derivation).Even in the the plateau region where Re reaches its minimum value (Re = 208),   ∼ 0.4H.Thus the viscous scale is less than  in all our simulations.

Figure 4 .
Figure 4. Vertical disc structure (vertical profiles of time-and-horizontally averaged quantities) from the fiducial simulation at Pm = 4. Top: root-meansquare velocity components.Bottom: root-mean-square magnetic field components.The vertical dashed lines at  = 0.0 and  = 3.2 show locations at which we performed the spectral analysis (see Section 4.2).

Figure 5 .
Figure 5. Spacetime diagrams for the horizontally-averaged toroidal magnetic field   taken from a simulation at the onset of the power-law scaling region at Pm = 4 (top) and from a simulation in the plateau region at Pm = 90 (bottom).The data are taken from low resolution runs (32 cells per scaleheight ).Note, to facilitate easier comparison, we set the colorbar limits of the top panel to match those in the lower panel, where the toroidal field is somewhat stronger.

Figure 6 .
Figure 6.Evidence of Parker instability.Top: flow field from a snapshot in the -plane at magnetic Prandtl number Pm = 4.The colorplot shows the density, black lines show magnetic field lines, and orange arrows indicate the velocity field.The upwelling of magnetic field between  = −0.5and  = 2.0 with plasma flowing down from the crest of the upwelling along the magnetic field lines and into the troughs is characteristic of structures expected from the Parker instability.Bottom: vertical profiles of magnetic buoyancy frequency squared at different Pm in the upper half-plane.(Note that the upper -cut-off of  = 3.3 on the -axis of the bottom panel has been chosen to coincide with that on the -axis of the top panel.)

Figure 7 .
Figure 7. Magnetic energy spectra in the (  ,   )-plane near the mid-plane at  = 0 (bottom panel) and in the atmosphere at  = 3.2 (top panel) for Pm = 4, Rm = 18750 and resolution 128/.The spectra are anisotropicstrongly inclined towards the   -axis -due to Keplerian shear, and become concentrated towards smaller wavenumbers with increasing height.

Figure 8 .
Figure 8. Spectra of the nonlinear terms governing the radial magnetic field   : induction-advection terms  (ℎ)  (left column) and  ()  (middle column) together with the vertical transport term  ()

Figure 9 .
Figure 9. Spectra of the dynamical terms governing the toroidal magnetic field   : Maxwell stress M (left column) and the vertical transport term  ()  (right column), which are the dominant terms for this field component, in the (  ,   )-plane at  = 0 (bottom) and  = 3.2 (top) for the same run as in Figure 8.The positive M (yellow/red) at small wavenumbers ensures amplification of toroidal field energy | B | 2 /2.This process at  = 0 is a part of the non-modal MRI-growth and extends over a broader range of wavenumbers than at  = 3.2, where it is localized about k = 0 and related to Parker instability.The term  () is negative at small wavenumbers near the mid-plane and in the atmosphere, thus acting as a sink for the toroidal field, which indicates removal of this field component due to buoyancy.
small and intermediate wavenumbers (i.e.|  | ≲ 30, |  | ≲ 20, at  = 0) where radial field is produced, as the vital area of the turbulence.Since the toroidal field   is the dominant field component, it gives the largest contribution to N (ℎ) the vital area both at  = 0 and  = 3.2 (blue areas in the right-hand column), removing the large-scale radial field due to buoyancy at all .The dynamics of the toroidal field   is governed by the Maxwell stress M and the non-linear terms N

Figure 10 .
Figure 10.Spectra of the dynamical terms for the vertical magnetic field   : induction term hand column) and the horizontal advection term  ()  (right-hand column) in the (  ,   )-plane at  = 0 (bottom row) and  = 3.2 (top row) for the same simulation as in Figure 8.At small wavenumbers, the vertical field is produced from the horizontal (i.e.toroidal) field due to Parker instability at those wavenumbers where  (ℎ)  > 0 (yellow/red).These wavenumbers are comparable to the characteristic wavenumber of the Parker instability, 1/  ∼ 3 (see the text), as shown by a circle with this radius for reference (top-left hand panel).By contrast, the advection term  ()  < 0 acts as a sink (blue regions).is negative N ()  < 0 in the vital area at all  (right-hand column), transferring the vertical field energy to larger wavenumbers.

Figure 11 .
Figure11.Comparison of the ring-averaged compensated kinetic   (blue) and magnetic   (red) energy spectra at Pm = 4 (solid) and Pm = 90 (dash-dot) at fixed Rm = 18750, both at the mid-plane  = 0 (bottom) and in the atmosphere  = 3.2 (top).As Pm increases, the kinetic energy spectra increases mostly at small  and becomes steeper due to increased viscosity.The spectral magnetic energy, on the other hand, increases at all  with Pm: at the mid-plane it shifts to lower , while in the atmosphere it neither changes shape nor shifts, but increases only in magnitude by about the same factor at all .