# Gravito-turbulence and dynamo in poorly ionized protostellar discs – I. Zero-net-flux case

A. Riols, W. Xu, G. Lesur, M.W. Kunz, H. Latter

► To cite this version:

A. Riols, W. Xu, G. Lesur, M.W. Kunz, H. Latter. Gravito-turbulence and dynamo in poorly ionized protostellar discs – I. Zero-net-flux case. Monthly Notices of the Royal Astronomical Society, 2021, 506 (1), pp.1407-1426. 10.1093/mnras/stab1637. hal-03267321

HAL Id: hal-03267321

https://hal.science/hal-03267321

Submitted on 12 Aug 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Gravito-turbulence and dynamo in poorly ionized protostellar discs
– I. Zero-net-flux case

A. Riols,<fnref n="1" />* W. Xu,<fnref n="2" /> G. Lesur,<fnref n="1" /> M. W. Kunz<fnref n="2" />,<fnref n="3" /> and H. Latter<fnref n="4" />

Accepted 2021 May 28. Received 2021 May 2; in original form 2021 February 15

## ABSTRACT

In their early stages, protoplanetary discs are sufficiently massive to undergo gravitational instability (GI). This instability is thought to be involved in mass accretion, planet formation via gas fragmentation, the generation of spiral density waves, and outbursts. A key and very recent area of research is the interaction between the GI and magnetic fields in young protoplanetary discs, in particular whether this instability is able to sustain a magnetic field via a dynamo. We conduct 3D, stratified shearing-box simulations using two independent codes, PLUTO and Athena++, to characterize the GI dynamo in poorly ionized protostellar discs subject to ambipolar diffusion. We find that the dynamo operates across a large range of ambipolar Elssaser number Am (which characterizes the strength of ambipolar diffusion) and is particularly strong in the regime Am = 10–100, with typical magnetic to thermal energy ratios of order unity. The dynamo is only weakly dependent on resolution (at least for Am ≲ 100), box size, and cooling law. The magnetic field is produced by the combination of differential rotation and large-scale vertical roll motions associated with spiral density waves. Our results have direct implications for the dynamo process in young protoplanetary discs and possibly some regions of active galactic nucleus discs.

**Key words:** dynamo – instabilities – turbulence – protoplanetary discs.

## 1 INTRODUCTION

Magnetic fields should be ubiquitous in accretion discs and are likely to influence their dynamics and long-term evolution significantly. Disc magnetism is believed to facilitate the accretion process by enabling turbulent angular-momentum transport and by launching large-scale winds (Balbus 2003; Wardle 2007; Armitage 2019). It is also possible that magnetic fields are responsible for the formation of large-scale structures, such as concentric rings (Kunz & Lesur 2013; Bai 2015; Béthune, Lesur & Ferreira 2017; Riols & Lesur 2019; Suriano et al. 2019), which are now widely observed in protoplanetary discs by ALMA and SPHERE (e.g. ALMA Partnership et al. 2015; Avenhaus et al. 2018; Andrews 2020). Despite their perceived importance, direct measurements of magnetic-field strengths in protoplanetary discs are rare. Donati et al. (2005) detected azimuthal fields at a strength of 1 kG in the inner part of the FU-Ori disc using the high-resolution spectropolarimeter ESPaDOnS. In the outer parts of discs, researchers have attempted to measure fields directly from the Zeeman splitting of CN lines (Vlemmings et al. 2019), and have set an upper limit of 20 mG on the toroidal field and 0.8 mG on the vertical field. Indirect techniques using dust polarisation have detected magnetic fields with various intensity and morphologies (Carrasco-González et al. 2010; Stephens et al. 2014; Goddi et al. 2017), but these observations remain disputed (Kataoka et al. 2017;

Tazaki, Lazarian & Nomura 2017). Magnetic fields have been also recorded in meteorite samples, suggesting fields of a few tens of mG at $R > 1$ au in the primitive Solar system (Fu et al. 2014; Cournede et al. 2015; Fu et al. 2020). On the other hand, in active galactic nuclei (AGN) and binary discs, there are some direct detections of a strong or moderate poloidal field (Martí-Vidal et al. 2015), but the existence of magnetic fields is mainly inferred from the jets associated with these objects (see review by Spruit 1996).

A central but difficult question pertains to the origin of these magnetic fields and their sustenance against dissipative effects resulting from the microphysics of the plasma such as non-ideal magnetohydrodynamic (MHD) effects and/or turbulent diffusion. In discs around young stars, one possibility is that the residual field from the primordial nebulae is advected and amplified by the combined effects of gravitational contraction, rotation, and accretion. Indeed, simulations of magnetic star formation incorporating non-ideal MHD effects suggest a field strength ~0.1 G in the central flux tubes of the 'first protostellar core' (e.g. Desch & Mouschovias 2001; Kunz & Mouschovias 2010; Masson et al. 2016; Xu & Kunz 2021), consistent with estimates for the protosolar magnetic field as deduced from meteoritic data (Levy & Sonett 1978). To what extent this magnetic flux is then retained by and transported through the protostellar accretion disc during the Class 0 and I phases depends upon the (unknown) level of turbulent transport and diffusivity in the disc (e.g. Guilet & Ogilvie 2013; Bai & Stone 2017; Leung & Ogilvie 2019).

Another possibility is that any weak seed magnetic field is amplified in the disc itself by a turbulent dynamo (e.g. Pudritz 1981;

<footnote n="1">¹Institut de Planétologie et d'Astrophysique de Grenoble (IPAG), Univ. Grenoble Alpes, CNRS, F-38000 Grenoble, France</footnote>

<footnote n="2">²Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane, Peyton Hall, Princeton, NJ 08544, USA</footnote>

<footnote n="3">³Princeton Plasma Physics Laboratory, PO Box 451, Princeton, NJ 08543, USA</footnote>

<footnote n="4">⁴DAMTP, University of Cambridge, Centre for Mathematical Sciences, Wilberforce Road, Cambridge CB3 0WA, UK</footnote>

<footnote>*E-mail: antoine.riols@univ-grenoble-alpes.fr</footnote>
Balbus & Hawley 1998). For nearly 30 yr, the magneto-rotational instability (MRI; Balbus & Hawley 1991; Hawley, Gammie & Balbus 1995) has been considered as a primary source of turbulence and angular-momentum transport in astrophysical discs. Several studies have pointed out the existence of a sub-critical dynamo involving the MRI that exhibits regular reversals of the large-scale magnetic field, somewhat reminiscent of the 'butterfly' diagram of the solar dynamo (Brandenburg et al. 1995; Hawley, Gammie & Balbus 1996; Lesur & Ogilvie 2008). However, there is today a broad consensus that the MRI dynamo reliably survives only in the very inner parts of protoplanetary discs, within 0.1–1 au where the temperature is high enough to thermally ionize the plasma (Gammie 1996; Sano et al. 2000). Further away, the ionization fraction is too low and the dynamics is regulated by non-ideal MHD effects, such as Ohmic dissipation, the Hall effect, and ambipolar diffusion, which all tend to suppress any form of dynamo process (Fleming, Stone & Hawley 2000; Sano & Stone 2002; Wardle & Salmeron 2012; Bai & Stone 2013; Lesur, Kunz & Fromang 2014; Bai 2015). Even in the inner regions of protoplanetary discs and in more ionized objects like AGN discs, the viability of the MRI dynamo process is not guaranteed, because of the low magnetic Prandtl number that characterizes the plasma in these environments (see Fromang et al. 2007; Balbus & Henri 2008; Käpylä & Korpi 2011; Meheut et al. 2015; Riols et al. 2017b, for more discussions about this issue).

Another potential source of turbulence in accretion discs is the gravitational instability (GI). This instability is thought to manifest during the early phases of protoplanetary discs when their masses are large enough for self-gravity to be dynamically important. Depending on the efficiency of the cooling process, GI either triggers fragmentation of the disc or saturates in a 'gravito-turbulent' state characterized by spiral density waves that are continuously generated and dissipated (Gammie 2001; Rice et al. 2003, 2006; Durisen et al. 2007). Based on measurements of disc masses, recent surveys found that 50 per cent of Class 0 and 10–20 per cent of Class I sources could be unstable to GI (Tobin et al. 2013; Mann et al. 2015). Recent images of large-scale spirals arms in several protoplanetary discs (e.g. Elias 2-27, WaOph 6, and MWC758), and more disordered 'streamer' structures in FU Ori systems, have been attributed to GI activity (Dong et al. 2016; Liu et al. 2016; Pérez et al. 2016; Huang et al. 2018). GI is also believed to occur in AGN discs at distances of ~0.01–0.1 pc, though fragmentation (instead of gravito-turbulence) is expected because the cooling time-scale is much shorter than the dynamical time-scale (Lodato 2007). In particular, this leads to the AGN disc truncation problem (Goodman 2003).

Recently, the ability of GI motions to amplify and sustain a magnetic field has been investigated through local (Riols & Latter 2019, hereafter RL19) and global numerical simulations (Deng, Mayer & Latter 2020). These studies have shown that gravito-turbulence works as a dynamo itself, producing strong azimuthal fields that can reach nearly thermal strengths even in poorly ionized gas characterized by magnetic Reynolds number of order ~10. The magnetic fields generated back-react on the flow significantly and tend to reduce the strength of the spiral density waves. The mechanism behind the amplification of the field, identified by RL19, involves the action of counter-rotating vertical rolls associated with the spiral density waves. These rolls stretch and fold azimuthal field lines (α effect), leading to the generation of a mean radial field, which in turn is sheared to produce toroidal fields (Ω effect). However, the GI-dynamo process has been examined in discs where magnetic dissipation results exclusively from Ohmic resistivity. In Class 0 discs, it is expected that other non-ideal effects, such as ambipolar diffusion, prevail in regions where GI may be active (see Dapp, Basu & Kunz 2012, who have used a detailed chemical equilibrium

model to compute diffusivities in nascent protostellar discs). A natural question that arises is whether or not ambipolar diffusion changes the growth and saturation of the GI dynamo. In particular, is the amplification mechanism similar to that at work in discs regulated by Ohmic dissipation? The answer is non-trivial since there are major differences between the two diffusion processes: unlike Ohmic dissipation, ambipolar diffusion is non-linear, anisotropic, and does not allow reconnection of the field.

Another important question that has been only partially addressed (Riols & Latter 2018a) is the effect of a net vertical field threading the disc on the dynamo. Such a field, if sufficiently strong, could severely impede gravito-turbulence (Riols & Latter 2018a) and have direct consequences for the GI-dynamo process. In addition, such a vertical field is expected to facilitate the launching of a wind (Wardle & Koenigl 1993; Ferreira & Pelletier 1995; Fromang et al. 2013; Lesur, Ferreira & Ogilvie 2013).

In this paper, we focus on the zero-net-flux configuration (without a mean vertical field) and examine the effect of ambipolar diffusion on the gravito-turbulent dynamo. We performed 3D MHD stratified shearing-box simulations with self-gravity and vertical stratification using two different codes, PLUTO and Athena++. (The Hall effect and Ohmic dissipation are omitted.) Broad agreement between the code results is achieved, demonstrating their weak dependence on the implementation of self-gravity and on other numerical details (time-stepping methods, orbital advection algorithm, etc.). We explore different cooling prescriptions (namely, linear cooling and $T^4$ cooling) and consider a wide range of ambipolar Elsasser numbers (Am) from 1 to 200. Given the size of the box necessary to capture GI and thus our limited resolution per *H*, it is not possible to probe the regime of yet larger Am. In Paper II, we examine the effect of a net vertical field on the dynamo behaviour by varying the midplane $\beta$ from $10^6$ to $10^3$.

We find that the gravito-turbulent dynamo in the presence of ambipolar diffusion is very similar to that described by RL19: there is no fundamental difference between the magnetic fields generated in the presence of ambipolar diffusion and those produced subject to Ohmic dissipation. Dynamo amplification is most efficient between $Am = 30$ and $Am = 100$, giving super-equipartition fields (with $\beta \sim 1$) that are predominantly toroidal and whose energy resides primarily at large scales. The dynamo weakens in the regime of large $Am \gtrsim 100$ and for sufficiently small Am. As in RL19, the field is amplified by poloidal velocity rolls associated with large-scale density waves. By using a new technique involving phase-folded simulation snapshots, we provide a clearer view of how gravitationally induced spiral arms amplify magnetic fields.

The structure of the paper is as follows. In Section 2, we describe the model used and the basic equations solved, as well as the numerical setup of our simulations. In Section 3, we study the saturated state of the dynamo as a function of Am and compare the results amongst two different cooling laws. In Section 4, we investigate in detail how spirals generate a mean-field dynamo in the presence of ambipolar diffusion. We study in Section 5 the dependence of this dynamo on numerical details, in particular the resolution and box size. Finally, in Section 6, we discuss the implications of our work for the evolution of protoplanetary and AGN discs.

# 2 MODEL AND NUMERICAL SETUP

## 2.1 Governing equations

We perform numerical simulations of gravito-turbulence and dynamo in a poorly ionized, differentially rotating disc using the local
shearing-box approximation. In this model, a small patch of the disc, co-orbiting with angular velocity $\Omega = \Omega_0 \hat{e}_z$ at a fixed radial location $R_0$, is represented by Cartesian coordinates $(x, y, z)$, with $x = R - R_0$ and $y = R_0\phi$ corresponding to the radial ($R$) and azimuthal ($\phi$) directions of a cylindrical coordinate system. Differential rotation is included through the Coriolis force and a background linear shear $v_0 = -Sx\hat{e}_y$, where $S = -(d\Omega/d\ln R)_{R_0}$ is the local shear frequency; we use $S = (3/2)\Omega_0$, corresponding to Keplerian rotation. We adopt an ideal equation of state for the fluid in the disc, with the pressure $P$ and the mass density $\rho$ related by $P = \rho c_s^2/\gamma$, where $c_s \equiv (\gamma P/\rho)^{1/2}$ is the sound speed and $\gamma$ is the ratio of specific heats (=5/3). The internal energy $U \equiv P/(\gamma - 1)$ of the disc is lost through radiative cooling at the rate $\Lambda^-$. An initially weak magnetic field $B$ with zero net vertical flux is frozen into the minority charged species, which drift relative to the bulk neutral fluid at a rate controlled by collisional drag and parametrized through the field-strength-dependent ambipolar diffusivity $\eta_A$.

Under these conditions, the equation governing the disc evolution are

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0, \qquad (1) $$

$$ \frac{\partial(\rho\mathbf{v})}{\partial t} + \nabla \cdot (\rho\mathbf{v}\mathbf{v}) = -2\Omega \times \rho\mathbf{v} - \rho\nabla\Phi - \nabla P + \mathbf{J} \times \mathbf{B}, \qquad (2) $$

$$ \frac{\partial \mathbf{B}}{\partial t} - \nabla \times (\mathbf{v} \times \mathbf{B}) = \nabla \times [\eta_A(\mathbf{J} \times \hat{\mathbf{e}}_b) \times \hat{\mathbf{e}}_b], \qquad (3) $$

$$ \frac{\partial U}{\partial t} + \nabla \cdot (U\mathbf{v}) = -P\nabla \cdot \mathbf{v} - \Lambda^{-}, \qquad (4) $$

where $\mathbf{J} = \nabla \times \mathbf{B}$ is the current density and $\hat{\mathbf{e}}_b \equiv \mathbf{B}/B$ is the magnetic-field unit vector. We decompose the fluid velocity $\mathbf{v}$ into its background and fluctuating parts:

$$ \mathbf{v} = \mathbf{v}_0 + \mathbf{u} = -\frac{3}{2}\Omega_0 x \hat{\mathbf{e}}_y + \mathbf{u}. \qquad (5) $$

The potential $\Phi$ in equation (2) is the sum of the tidal potential in the local frame,

$$ \Phi_c = \frac{1}{2} \Omega_0^2 z^2 - \frac{3}{2} \Omega_0^2 x^2, \qquad (6) $$

and the gravitational potential of the disc itself, $\Phi_s$. The latter obeys the Poisson equation

$$ \nabla^2 \Phi_s = 4\pi G\rho; \qquad (7) $$

its dynamical importance is quantified by the Toomre parameter

$$ Q \equiv \frac{c_s \Omega_0}{\pi G \Sigma}, \qquad (8) $$

where $\Sigma = \int dz \rho$ the surface mass density of the disc.

For most of our simulations we assume a dependence of the radiative cooling on the temperature $T$ that is proportional to $T^4$:

$$ \Lambda^{-} = \frac{U}{\tau_{0}} \left( \frac{T}{T_{0}} \right)^{3}, \qquad (9) $$

where $T_0 \propto c_{s0}^2$ is the initial and uniform background temperature of the disc and $\tau_0$ is a uniform and constant coefficient that represents the effective cooling time-scale when $T = T_0$. This prescription corresponds to an extreme case since cooling is not expected to be steeper than $T^4$. The steepness of this cooling law ensures that a thermal equilibrium between explicit heating and cooling always exists, independently of the location of the vertical boundary. By contrast, most numerical studies of gravito-turbulence in the existing

literature have instead adopted the linear cooling function

$$ \Lambda^{-} = \frac{U}{\tau_{c}}, \qquad (10) $$

with the characteristic time-scale $\tau_c$ referred to as the ‘cooling time’. To make comparison with those studies, and to show that our qualitative results are insensitive to the choice of cooling law, we also performed a few simulations using equation (10). Explicit thermal conduction is neglected.

Finally, we quantify the ambipolar diffusivity $\eta_A$ in the induction equation (3) using the dimensionless ambipolar Elsässer number

$$ \text{Am} \equiv \frac{v_A^2}{\Omega_0 \eta_A}, \qquad (11) $$

where $v_A \equiv B/\sqrt{\rho}$ is the Alfvén speed. When the only charged species in the disc are the ions and the electrons, $\eta_A \approx v_A^2 \tau_{ni}$, where $\tau_{ni}$ is the neutral-ion collision time-scale (proportional to the inverse of the ion mass density). In this case, $\text{Am} \approx (\Omega_0 \tau_{ni})^{-1}$ quantifies the statistical number of times that a neutral particle collides with an ion per orbital time-scale; large values of Am indicate that the neutrals are well coupled to the ions (and thus to the magnetic field). Under more general conditions, $\eta_A$ is obtained by solving a multispecies chemical network in which various ionization sources (X-rays, cosmic rays, radioactive decay, and far-ultraviolet [FUV] from the protostar) compete with charge recombination occurring both in the gas phase and on grain surfaces. For this first study of gravito-turbulence and dynamo in poorly ionized discs, we have opted for simplicity by crudely setting Am = const, effectively assuming a constant density for the charged species. In Appendix C, we adopt an alternative prescription with Am = const in the disc mid-plane and Am = $\infty$ above a specified height, to mimic the presence of a well-ionized surface layer. The resulting simulation demonstrates that our results are similar whether or not such an ionized layer is present.

## 2.2 Numerical methods

We solve equations (1–4) in a 3D shearing-periodic domain of size ($L_x, L_y, L_z$), discretized on a mesh of ($N_x, N_y, N_z$) grid points. To increase the trustworthiness of our results, we have used two different and independent codes, PLUTO (Mignone et al. 2007) and Athena++ (Stone et al. 2020); we perform several comparison runs which are detailed in Section 3. Both codes use conservative finite-volume methods that solve the approximate Riemann problem at each inter-cell boundary and are well adapted to highly compressive flows (including shocks). The Riemann problem is handled by the HLLD solver, which is suitable for MHD. An orbital advection algorithm (similar to that detailed in Mignone et al. (2012)) is used to increase the computational speed and reduce numerical dissipation. The divergence of $\mathbf{B}$ is forced to 0 by a constrained-transport algorithm. In most of our simulations, the electromotive force (EMF) at the edge of each cell is reconstructed by the UCT_CONTACT averaging procedure (in both PLUTO and Athena++). Details about this method can be found in the online PLUTO documentation and in Gardiner & Stone (2005). Both codes conserve the total energy and so the heat equation is not solved directly as in equation (4). Heating by compression, magnetic diffusion, and shocks are all captured implicitly through the conservation of total energy, while cooling is treated as an additional sink term. Finally, the temporal integration is handled by a second-order-accurate Runge-Kutta method in PLUTO, while it is handled by a second-order-accurate van Leer predictor-corrector method in Athena++.
### 2.2.1 Implementation of self-gravity in PLUTO

To compute the self-gravitational potential $\Phi_s$ in PLUTO, we take advantage of the shear-periodic boundary conditions. Following Riols & Latter (2016) and Riols, Latter & Paardekooper (2017a), we first shift back the density in y to the time at which it was last periodic. Then, for each plane at altitude $z_k$, we compute the direct 2D Fourier transform of the density $\hat{\rho}_{k_x, k_y}(z)$, where $k_x$ and $k_y$ are the z-dependent radial and azimuthal wavenumbers. Using equation (7), it is straightforward to show that the Fourier coefficients of the gravitational potential satisfy the Helmholtz equation

$$ \left( \frac{d^2}{dz^2} - k^2 \right) \hat{\Phi}_{k_x, k_y}(z) = 4\pi G \hat{\rho}_{k_x, k_y}(z) \quad (12) $$

with $k^2 = k_x^2 + k_y^2$. This differential equation is solved in the complex plane by means of a fourth-order finite-difference scheme and a direct inversion method. The discretized system takes the form of a linear problem $AX = B$, where X is a vector representing the z-profile of the discretized potential, A is a penta-diagonal matrix, and B is a column vector containing the right hand side of equation (12) and extra coefficients setting the boundary conditions. We use a fast algorithm involving $O(N_Z)$ flops to invert the matrix and obtain the discretized coefficients $\hat{\Phi}_{k_x, k_y}(z_k)$. For each altitude $z_k$, we finally compute the inverse Fourier transform of the potential and shift it back to the initial sheared frame.

The gravitational acceleration $-\rho\nabla\Phi_s$ is then obtained by computing the derivative of the potential in each direction using a fourth-order, cell-centred, finite-difference method. The result is implemented as a source term in the discretized version of the momentum equation (2).

This implementation was tested by Riols et al. (2017a, see their Appendices A and B) using a stratified disc equilibrium and its linear stability.

### 2.2.2 Implementation of self-gravity in Athena++

In Athena++, we solve for $\Phi_s$ using a discrete form of the Poisson equation in which the derivatives in the Laplace operator are computed with second-order finite differences. For example, the (continuous) second derivative $\partial_x^2\Phi_s$ is replaced with the discrete operator $\partial_x^2\Phi_{s,i} = (\Phi_{s,i+1} - 2\Phi_{s,i} + \Phi_{s,i-1})/(\Delta x)^2$, where i is the cell index in the x direction and $\Delta x$ is the grid spacing. $\partial_x^2$ and $\partial_x^2$ share the same eigenfunctions $\exp(ik_x x)$, but the eigenvalue for $\partial_x^2$ is $-\hat{k}_x^2 \equiv 2[\cos(k_x \Delta x) - 1]/(\Delta x)^2$ instead of $-k_x^2$. (Their difference vanishes at long wavelengths, $|k_x \Delta x| \to 0$.) Similar to the PLUTO implementation, to compute $\Phi_s$ we do a shifted 2D Fourier transform in x, y and solve the resulting Helmholtz equation in z, where $k^2$ in equation (12) is replaced by $\hat{k}_x^2 + \hat{k}_y^2$ and $d^2/dz^2$ is computed with a second-order finite difference, to be consistent with the treatment in x, y plane.

Instead of calculating $-\rho\nabla\Phi_s$ at cell centres and directly applying it as a source term, Athena++ calculates the momentum flux due to self gravity using the (second-order accurate) gravitational stress evaluated at cell interfaces. This treatment better conserves momentum, and is more accurate than applying the gravitational acceleration as a source term whenever the density is discontinuous or under-resolved.

We have tested this implementation of self gravity in Athena++ using various stratified disc equilibria and their linear stability (see appendix A).

## 2.3 Numerical setup and free parameters

In our codes, the unit of time is $\Omega^{-1}$ the orbital frequency and the unit of length is $H_0 \equiv c_{s0}/\Omega_0$ the standard disc scale height (corresponding to the actual scale height of an isothermal disc without self-gravity; the true scale height with self-gravity is smaller than this value). For most of our simulations, we use a box of size $L_x = L_y = 20H_0$ and $L_z = 9H_0$. The horizontal size is chosen so that three or four spiral structures may be captured in the box at the same time, with the caveat that the shearing-box approximation, which assumes $L_x, L_y \ll$ disc radius, is only marginally justified at such large box sizes. The vertical size is chosen so that the energy that escapes through the vertical boundary is not too large in comparison to the energy lost through radiative cooling. Our fiducial numerical resolution is 256 $\times$ 256 $\times$ 144 cells, which corresponds to 13 cells per $H_0$ in the horizontal directions and 16 cells per $H_0$ in the vertical direction. (This horizontal box size and resolution are identical to those of RL19.) In Section 5, we study the dependence of our results on resolution and box size.

The boundary conditions are periodic in y and shear-periodic in x. In the vertical direction, for hydrodynamic variables (velocity, density, and pressure), we use a simple outflow boundary condition and suppress inflow by setting $v_z = 0$ in the ghost zones if the local velocity points towards the disc. For the magnetic field, we use a vertical boundary condition in which the field is 'mirrored' into the ghost zones with $B_x$ and $B_y$ changing sign (such that $B_x = B_y = 0$ on the boundary). We checked that using the mirror symmetry or forcing $B_x = B_y = 0$ in all ghost zones gives the same results. For the self-gravitating potential, PLUTO assumes that

$$ \frac{d}{dz} \hat{\Phi}_{k_x, k_y} (\pm L_z/2) = \mp k \hat{\Phi}_{k_x, k_y} (\pm L_z/2). \quad (13) $$

This condition is an approximation of the Poisson equation in the limit of low density and ensures that the potential does not diverge at infinity. Athena++ simply sets $\hat{\Phi}_{k_x, k_y} = 0$ in the ghost zones for $k_x \neq 0$ and $k_y \neq 0$. This boundary condition for the potential matches that of PLUTO at $z \to \pm\infty$ (if the box had infinite size in z).

To compensate for the loss of mass via outflow, the gas is replenished near the mid-plane so that the total mass in the box (and therefore $\Sigma$) is maintained constant during the time of each simulation. When adding mass, we change the density of all cells by a constant factor, while fixing momenta and total energy density.

To avoid excessively small time steps, especially when a magnetic field is present, a finite density floor $\rho_{floor} = 10^{-4}$ is used. Such a density floor can be problematic in a relatively tall box, because the gravitational acceleration (accounting for both self gravity and the vertical component of stellar gravity) acting on the floor density creates unphysical acceleration that continually pulls material towards the mid-plane and generates artificial heating. We alleviate this problem by calculating the gravitational force as $-(\rho - \rho_{floor})\nabla\Phi$ instead of $-\rho\nabla\Phi$. Large values of magnetic diffusivity can also cause prohibitive time steps, and so we have capped $\eta_A$ so that

$$ \eta_A = \min \left\{ \eta_{\text{cap}}, \frac{1}{\text{Am}} \frac{B^2}{\rho \Omega} \right\}, \quad (14) $$

where $\eta_{\text{cap}} = 3$ (in units of $c_{s0}^2/\Omega$) for Am = 1 and $\eta_{\text{cap}} = 1$ for larger Am. For the Athena++ simulations that employ double resolution (512 $\times$ 512 $\times$ 288), ambipolar diffusion is calculated with a super-time-stepping algorithm (Alexiades, Amiez & Gremaud 1996) to reduce cost. We have tested at our fiducial resolution that super-time-stepping does not affect the results.
<table><caption>**Table 1.** List of hydrodynamical and zero-net-flux MHD simulations, combining Athena++ (AT-) and PLUTO (PL-) runs. The column 'transient' indicates the time before the simulation reaches a quasi-steady turbulent state. The parameters $τ_0$, $τ_c$, and Am are defined respectively in (9), (10), and (11).</caption><thead><tr><th>Run</th><th>Resolution</th><th>L<sub>z</sub><br>(in H<sub>0</sub>)</th><th>Duration<br>(in Ω<sup>-1</sup>)</th><th>Transient<br>(in Ω<sup>-1</sup>)</th><th>Cooling</th><th>τ<sub>0</sub> or τ<sub>c</sub><br>(in Ω<sup>-1</sup>)</th><th>Am</th></tr></thead><tbody><tr><td>AT-GHydro</td><td>256 × 256 × 144</td><td>9</td><td>500</td><td>100</td><td>T<sup>4</sup></td><td>10</td><td>—</td></tr><tr><td>AT-GHydro-lin</td><td>256 × 256 × 144</td><td>9</td><td>500</td><td>100</td><td>Linear</td><td>11.5</td><td>—</td></tr><tr><td>AT-ZNF-ideal</td><td>256 × 256 × 144</td><td>9</td><td>2000</td><td>500</td><td>T<sup>4</sup></td><td>10</td><td>—</td></tr><tr><td>AT-ZNF-ideal-lin</td><td>256 × 256 × 144</td><td>9</td><td>1000</td><td>450</td><td>Linear</td><td>9.5</td><td>—</td></tr><tr><td>AT-ZNF-ideal-HR</td><td>512 × 512 × 288</td><td>9</td><td>400</td><td>100</td><td>T<sup>4</sup></td><td>10</td><td>—</td></tr><tr><td>AT-ZNF-Am100</td><td>256 × 256 × 144</td><td>9</td><td>500</td><td>200</td><td>T<sup>4</sup></td><td>10</td><td>100</td></tr><tr><td>AT-ZNF-Am100-HR</td><td>512 × 512 × 288</td><td>9</td><td>700</td><td>0</td><td>T<sup>4</sup></td><td>10</td><td>100</td></tr><tr><td>AT-ZNF-Am10</td><td>256 × 256 × 144</td><td>9</td><td>450</td><td>200</td><td>T<sup>4</sup></td><td>10</td><td>10</td></tr><tr><td>AT-ZNF-Am10-lin</td><td>256 × 256 × 144</td><td>9</td><td>100</td><td>20</td><td>Linear</td><td>7</td><td>10</td></tr><tr><td>AT-ZNF-Am10-HR</td><td>512 × 512 × 288</td><td>9</td><td>100</td><td>0</td><td>T<sup>4</sup></td><td>10</td><td>10</td></tr><tr><td>AT-ZNF-Am1</td><td>256 × 256 × 144</td><td>9</td><td>500</td><td>200</td><td>T<sup>4</sup></td><td>10</td><td>1</td></tr><tr><td>PL-GHydro</td><td>256 × 256 × 144</td><td>9</td><td>250</td><td>50</td><td>T<sup>4</sup></td><td>10</td><td>—</td></tr><tr><td>PL-GHydro-Lz12</td><td>256 × 256 × 144</td><td>12</td><td>470</td><td>50</td><td>T<sup>4</sup></td><td>10</td><td>—</td></tr><tr><td>PL-GHydro-lin</td><td>256 × 256 × 144</td><td>9</td><td>400</td><td>0</td><td>Linear</td><td>11.5</td><td>—</td></tr><tr><td>PL-ZNF-ideal</td><td>256 × 256 × 144</td><td>9</td><td>800</td><td>250</td><td>T<sup>4</sup></td><td>10</td><td>—</td></tr><tr><td>PL-ZNF-ideal-lin</td><td>256 × 256 × 144</td><td>9</td><td>660</td><td>200</td><td>Linear</td><td>11.5</td><td>—</td></tr><tr><td>PL-ZNF-Am200</td><td>256 × 256 × 144</td><td>9</td><td>300</td><td>50</td><td>T<sup>4</sup></td><td>10</td><td>200</td></tr><tr><td>PL-ZNF-Am100</td><td>256 × 256 × 144</td><td>9</td><td>480</td><td>260</td><td>T<sup>4</sup></td><td>10</td><td>100</td></tr><tr><td>PL-ZNF-Am100-Lz12</td><td>256 × 256 × 144</td><td>12</td><td>400</td><td>175</td><td>T<sup>4</sup></td><td>10</td><td>100</td></tr><tr><td>PL-ZNF-Am100-lin</td><td>256 × 256 × 144</td><td>9</td><td>350</td><td>0</td><td>Linear</td><td>4</td><td>100</td></tr><tr><td>PL-ZNF-Am100-HR</td><td>512 × 512 × 288</td><td>9</td><td>160</td><td>100</td><td>T<sup>4</sup></td><td>10</td><td>100</td></tr><tr><td>PL-ZNF-Am30</td><td>256 × 256 × 144</td><td>9</td><td>240</td><td>0</td><td>T<sup>4</sup></td><td>10</td><td>30</td></tr><tr><td>PL-ZNF-Am10</td><td>256 × 256 × 144</td><td>9</td><td>460</td><td>200</td><td>T<sup>4</sup></td><td>10</td><td>10</td></tr><tr><td>PL-ZNF-Am10-Lz12</td><td>256 × 256 × 144</td><td>12</td><td>200</td><td>20</td><td>T<sup>4</sup></td><td>10</td><td>10</td></tr><tr><td>PL-ZNF-Am10-lin</td><td>256 × 256 × 144</td><td>9</td><td>100</td><td>0</td><td>Linear</td><td>7</td><td>10</td></tr><tr><td>PL-ZNF-Am3</td><td>256 × 256 × 144</td><td>9</td><td>200</td><td>20</td><td>T<sup>4</sup></td><td>10</td><td>3</td></tr></tbody></table>

Our initial conditions depend on whether or not a magnetic field is present. For our purely hydrodynamical simulations, we start from a self-gravitating disc equilibrium in which the strength of self gravity is specified by $Q_0 = 1$. The equilibrium is initially isothermal with temperature $T_0 = c_{s0}^2/γ$ (in code units). The density in the mid-plane is fixed to 1 so that the surface density is $Σ = 0.917$. Random, non-axisymmetric, finite-amplitude density, and velocity perturbations are injected to trigger a gravito-turbulent state. For our MHD runs, we use a gravito-turbulent steady state taken from a hydrodynamical simulation as the initial condition and add in a weak seed toroidal field. In this case, the seed field is of the form

$$B_y = B_{y0} \sin(2\pi z / L_z) \qquad (15)$$

and its strength is characterized by the initial mid-plane

$$\beta_{y0} = \frac{2P_0}{B_{y0}^2} = \frac{2\rho_0 c_{s0}^2}{\gamma B_{y0}^2}. \qquad (16)$$

The list of runs is shown in Table 1 with the different input parameters specified in the columns.

## 2.4 Diagnostics

### 2.4.1 Averaged quantities

To analyse the statistical behaviour of the turbulent flow, we define the standard and weighted box averages

$$\langle X \rangle = \frac{1}{V} \int_V dV X \quad \text{and} \quad \langle X \rangle_W = \frac{1}{V \langle \rho \rangle} \int_V dV \rho X, \qquad (17)$$

respectively, where $V = L_x L_y L_z$ is the volume of the box. For example, $E_k = \frac{1}{2} \langle \rho u^2 \rangle$ and $E_m = \frac{1}{2} \langle B^2 \rangle$ are the box-averaged kinetic and

magnetic energies, respectively. We also define the horizontally averaged vertical profile of a variable:

$$\bar{X}(z) = \frac{1}{L_x L_y} \iint dxdy X. \qquad (18)$$

This averaging will be useful when discussing mean-field equations in the context of dynamo (see Section 2.4.4).

There are three important box-averaged quantities with which we characterize the gravito-turbulent steady state. First, we define the 2D Toomre parameter

$$Q_W = \frac{\langle c_s \rangle_W \Omega_0}{\pi G \Sigma}, \qquad (19)$$

where $Σ = L_z \langle \rho \rangle$ is the average surface density of the disc. This quantity is directly relevant for characterizing the stability of GI spiral modes. Secondly, we follow the custom of quantifying the efficiency of the angular-momentum transport in terms of an $α$ parameter, although one which combines the contributions from the Reynolds ($H_{xy}$), gravitational ($G_{xy}$), and Maxwell ($M_{xy}$) stresses:

$$\alpha = \frac{\langle H_{xy} + G_{xy} + M_{xy} \rangle}{\langle P \rangle}, \qquad (20)$$

where

$$H_{xy} = \rho u_x u_y, \quad G_{xy} = \frac{1}{4\pi G} \frac{\partial \Phi}{\partial x} \frac{\partial \Phi}{\partial y}, \quad \text{and} \quad M_{xy} = -B_x B_y.$$

Finally, the average disc thermodynamics is described by the average effective cooling time-scale

$$\tau_{\text{eff}} = \frac{\langle U \rangle}{\langle \Lambda^{-1} \rangle}. \qquad (21)$$
For simulations with a linear cooling law, $\tau_{\text{eff}} = \tau_c$, while for simulations with $T^4$ cooling (equation 9),

$$ \tau_{\text{eff}} = \tau_0 \langle U \rangle \left\langle U \left( \frac{T}{T_0} \right)^3 \right\rangle^{-1} . \qquad (22) $$

This last relation implies that the average effective cooling time can be different from $\tau_0$ depending on how the temperature deviates from the isothermal condition $T = T_0$. This is important in the context of making direct comparisons between runs with $T^4$ cooling and runs with linear cooling. To facilitate such comparisons, we first compute the turbulent steady state using $T^4$ cooling, measure $\tau_{\text{eff}}$ via equation (22), and then set $\tau_c = \tau_{\text{eff}}$ in a complementary run using the linear cooling law (10).

### 2.4.2 Fourier decomposition and spectra

To analyse the structure of the flow and, in particular, the role of spiral waves and small-scale dynamics in the dynamo process, we employ a 2D Fourier decomposition of the fields. Namely, we denote by $\hat{u}^e(k_x, k_y, z)$ and $\hat{B}^e(k_x, k_y, z)$ the horizontal 2D decomposition (in Eulerian wavenumbers $k_x$ and $k_y$) of the turbulent velocity and magnetic field, for a given altitude $z$, and averaged over a given period of time. In this mathematical representation, the fields are a sum of a mean ($k_x = k_y = 0$) mode, axisymmetric modes with $k_y = 0$ ($k_x \neq 0$), and non-axisymmetric modes with $k_y \neq 0$ (commonly referred as ‘shearing waves’). These Fourier decompositions are calculated via a method described by Riols et al. (2017a, section 2.5.2) using the FFT algorithm. The 2D kinetic and magnetic energy power spectra are then defined as

$$ E_K(k_x, k_y, z) = \frac{1}{2} \sum H_0^{-1} |\hat{u}^e(k_x, k_y, z)|^2, \qquad (23) $$

$$ E_M(k_x, k_y, z) = \frac{1}{2} |\hat{B}^e(k_x, k_y, z)|^2. \qquad (24) $$

### 2.4.3 Phase-folding technique

As shown by Riols & Latter (2018b), when 3D dynamics are allowed in simulations of gravito-turbulence, the customary spiral density waves are accompanied by characteristic hydrodynamical motions such as vertical rolls that play an important role in amplifying magnetic fields. So far the analysis of these motions, as well as of the geometry of the magnetic around the spiral arms, has focused on single snapshots of the flow and field (e.g. see fig. 15 of RL19). However, the strong gravito-turbulence makes it difficult to interpret these snapshots. To filter out these turbulent fluctuations and see more clearly how spiral waves (and their accompanying vertical rolls) shape the magnetic field, we phase-fold our data as follows. First, we take a series of snapshots from moments when the fields are exactly periodic in $x$ (typically once every $2\Omega^{-1}$). From a Fourier decomposition of the 2D surface density, we then identify the most prominent spiral mode over the course of the simulation, say, $k_y = 2\pi/L_y$ and $k_x = k_{x_s}$. For each snapshot, we calculate the phase of this mode and shift all the fields in $x$ to ensure that this phase is reduced to zero (in that way we align the phase of the spiral waves). We then stack the shifted fields upon one another to obtain a time average. We perform an additional average of the data over the azimuthal ($y$) direction along the spiral structure, i.e. along $(k_{x_s}/k_{x_0})x + y = \text{const}$ (this is done by shifting the data in $x$ and averaging in $y$). The result is a phase-folded poloidal map ($x, z$) of the flow, where systematic perturbations associated with the most prominent spiral mode are

preserved but turbulent fluctuations uncorrelated with the spiral mode are suppressed by the averaging procedure.

### 2.4.4 Mean-field equations and EMFs

One of our principal goals is to characterize the ability of gravito-turbulent motions to supply and sustain a mean field dynamo. Therefore, a key quantity is the mean magnetic field $\vec{B}(z)$, averaged in the $x$ and $y$ directions, corresponding to the $k_x = k_y = 0$ Fourier component of $\vec{B}$. The equations governing its radial ($x$) and toroidal ($y$) projections are:

$$ \frac{\partial \vec{B}_x}{\partial t} = -\frac{\partial}{\partial z} (\vec{\mathcal{E}}_y + \vec{\mathcal{E}}_y^{\text{NI}}), \qquad (25) $$

$$ \frac{\partial \vec{B}_y}{\partial t} = -S\vec{B}_x + \frac{\partial}{\partial z} (\vec{\mathcal{E}}_x + \vec{\mathcal{E}}_x^{\text{NI}}), \qquad (26) $$

where $\vec{\mathcal{E}}(z) = \overline{u \times B}$ is ideal component of the horizontally averaged EMF and

$$ \vec{\mathcal{E}}^{\text{NI}}(z) = \overline{\eta_{\Lambda}(\mathbf{J} \times \hat{e}_b \times \hat{e}_b)} \qquad (27) $$

its non-ideal component. The ideal EMF can be decomposed into a sum of a ‘spiral’ component $\vec{\mathcal{E}}^* = \hat{u}^* \times \hat{B}^*$, where the ‘★’ superscript refers to the phase-folded profile of the fields (see Section 2.4.3) and a turbulent component defined by $\vec{\mathcal{E}}^I = \overline{u \times B - u^* \times B^*}$.

We remind the reader that the mean radial and toroidal fields in the box, respectively $\langle B_x \rangle$ and $\langle B_y \rangle$, equal to the $z$-average of $\vec{B}_x(z)$ and $\vec{B}_y(z)$, are free to evolve because of the open boundaries at the top and bottom of the computational domain.

## 3 PROPERTIES OF THE GRAVITO-TURBULENT DYNAMO

### 3.1 Initial flows: hydrodynamical simulations

A customary way to study the dynamo process in simulations is to compute the (unmagnetized) hydrodynamic state of fully developed turbulence, introduce a zero-net flux magnetic seed field, and see whether that field grows or decays over time. Following this scheme, we first simulate purely hydrodynamical gravito-turbulence, with the aim of comparing two different cooling laws (linear and $T^4$) as well as both codes, PLUTO and Athena++ (a similar code comparison was done by Riols et al. (2017a) using PLUTO and RODEO). The key parameters that control the strength of the turbulence are the cooling times $\tau_0$ (when $T^4$ cooling is used; see equation (9)) and $\tau_c$ (when linear cooling is used; see equation (10)). For $T^4$ cooling, we use the fiducial $\tau_0 = 10\Omega^{-1}$ (runs AT-GHydro and PL-GHydro). For linear cooling, we choose $\tau_c = 11.5\Omega^{-1}$, which matches the $\tau_{\text{eff}}$ measured in the $T^4$ runs (runs AT-GHydro-lin and PL-GHydro-lin, see Section 2.4.1 and Table 2).

Table 2 summarizes some box-averaged quantities obtained with PLUTO and Athena++. We find excellent agreement between the two codes for the hydrodynamical runs, despite the different methodology used to compute self-gravity. In particular, in both cases, the Toomre parameter settles at $Q_W = 1.12$ and the ratio of the gravitational and Reynolds stresses is around 2.5. The box-averaged quantities do not seem to depend too much on the cooling law, although the kinetic energy and Reynolds stress are slightly larger in the linear cooling runs.

From energy conservation and the fact that outflows are observed to be very weak in the purely hydrodynamical case, one may expect
<table><caption>**Table 2.** Box-averaged and flux quantities measured in the simulations of Table 1. $τ_{eff}$ is the effective cooling time defined in (21). $Q_w$ is the Toonne parameter defined in (19). $U$ is the internal energy, $E_k$ and $E_m$ the kinetic and magnetic energies defined in Section 2.4.1. $H_{xy}$, $G_{xy}$, $M_{xy}$ are defined in Section 2.4.1 and correspond to the Reynolds, gravitational, and magnetic stresses. $\alpha$ is the transport coefficient defined in 20 and $\dot{m}_w$ is the vertical mass loss rate per unit area.</caption><thead><tr><th>Run</th><th>τ<sub>eff</sub></th><th>Q<sub>w</sub></th><th>⟨U⟩</th><th>E<sub>kl</sub>(⟨U⟩)</th><th>E<sub>ml</sub>(⟨U⟩)</th><th>⟨H<sub>xy</sub>⟩/(⟨U⟩)</th><th>⟨G<sub>xy</sub>⟩/(⟨U⟩)</th><th>⟨M<sub>xy</sub>⟩/(⟨U⟩)</th><th>α</th><th>ṁ<sub>w</sub>/(ΣΩ)</th></tr></thead><tbody><tr><td>AT-GHydro</td><td>11.91</td><td>1.116</td><td>0.1173</td><td>0.2067</td><td>0</td><td>0.01646</td><td>0.03848</td><td>0</td><td>0.08241</td><td>-6.67 × 10<sup>-6</sup></td></tr><tr><td>AT-GHydro-lin</td><td>11.5</td><td>1.146</td><td>0.1257</td><td>0.2851</td><td>0</td><td>0.01860</td><td>0.0400</td><td>0</td><td>0.08790</td><td>1.14 × 10<sup>-4</sup></td></tr><tr><td>AT-ZNF-ideal</td><td>9.38</td><td>1.175</td><td>0.1297</td><td>0.1838</td><td>0.07449</td><td>0.01822</td><td>0.03105</td><td>0.02212</td><td>0.10708</td><td>1.15 × 10<sup>-4</sup></td></tr><tr><td>AT-ZNF-ideal-lin</td><td>9.5</td><td>1.442</td><td>0.2034</td><td>0.1092</td><td>0.1535</td><td>0.01806</td><td>0.01161</td><td>0.05048</td><td>0.1202</td><td>1.01 × 10<sup>-3</sup></td></tr><tr><td>AT-ZNF-ideal-HR</td><td>13.14</td><td>1.126</td><td>0.1184</td><td>0.1410</td><td>0.03550</td><td>0.01277</td><td>0.02471</td><td>0.01321</td><td>0.07603</td><td>6.01 × 10<sup>-5</sup></td></tr><tr><td>AT-ZNF-Am100</td><td>4.31</td><td>1.376</td><td>0.1767</td><td>0.1596</td><td>0.6007</td><td>0.04185</td><td>0.007985</td><td>0.1258</td><td>0.2634</td><td>2.11 × 10<sup>-3</sup></td></tr><tr><td>AT-ZNF-Am100-HR</td><td>5.67</td><td>1.297</td><td>0.1574</td><td>0.1472</td><td>0.3358</td><td>0.03043</td><td>0.01328</td><td>0.08580</td><td>0.1943</td><td>9.17 × 10<sup>-4</sup></td></tr><tr><td>AT-ZNF-Am10</td><td>7.01</td><td>1.229</td><td>0.1418</td><td>0.1627</td><td>0.2854</td><td>0.02367</td><td>0.02842</td><td>0.04959</td><td>0.1525</td><td>5.13 × 10<sup>-4</sup></td></tr><tr><td>AT-ZNF-Am10-lin</td><td>7</td><td>1.230</td><td>0.1503</td><td>0.2305</td><td>0.2841</td><td>0.0250</td><td>0.03167</td><td>0.05415</td><td>0.1662</td><td>1.04 × 10<sup>-3</sup></td></tr><tr><td>AT-ZNF-Am10-HR</td><td>6.98</td><td>1.235</td><td>0.1427</td><td>0.1636</td><td>0.2603</td><td>0.02493</td><td>0.02667</td><td>0.04889</td><td>0.1507</td><td>3.97 × 10<sup>-4</sup></td></tr><tr><td>AT-ZNF-Am1</td><td>10.89</td><td>1.130</td><td>0.1204</td><td>0.2537</td><td>0.02023</td><td>0.01696</td><td>0.03887</td><td>0.004507</td><td>0.0905</td><td>1.94 × 10<sup>-5</sup></td></tr><tr><td>PL-GHydro</td><td>11.41</td><td>1.119</td><td>0.1180</td><td>0.2417</td><td>0</td><td>0.01673</td><td>0.04171</td><td>0</td><td>0.0876</td><td>-1.7 × 10<sup>-6</sup></td></tr><tr><td>PL-GHydro-Lz12</td><td>10.62</td><td>1.13</td><td>0.09024</td><td>0.248</td><td>0</td><td>0.01944</td><td>0.04363</td><td>0</td><td>0.0946</td><td>-4.38 × 10<sup>-8</sup></td></tr><tr><td>PL-GHydro-lin</td><td>11.5</td><td>1.144</td><td>0.1255</td><td>0.2996</td><td>0</td><td>0.0199</td><td>0.0409</td><td>0</td><td>0.0911</td><td>1.10 × 10<sup>-4</sup></td></tr><tr><td>PL-ZNF-ideal</td><td>10.56</td><td>1.153</td><td>0.1246</td><td>0.1832</td><td>0.0473</td><td>0.01679</td><td>0.03118</td><td>0.01576</td><td>0.0955</td><td>7.75 × 10<sup>-5</sup></td></tr><tr><td>PL-ZNF-ideal-lin</td><td>11.5</td><td>1.170</td><td>0.1307</td><td>0.213</td><td>0.0372</td><td>0.01707</td><td>0.03053</td><td>0.01342</td><td>0.09153</td><td>1.15 × 10<sup>-4</sup></td></tr><tr><td>PL-ZNF-Am200</td><td>7.41</td><td>1.234</td><td>0.1425</td><td>0.1712</td><td>0.1504</td><td>0.02269</td><td>0.02494</td><td>0.04565</td><td>0.139</td><td>2.81 × 10<sup>-4</sup></td></tr><tr><td>PL-ZNF-Am100</td><td>3.93</td><td>1.398</td><td>0.1821</td><td>0.1700</td><td>0.6149</td><td>0.04505</td><td>0.006020</td><td>0.1397</td><td>0.2861</td><td>1.97 × 10<sup>-3</sup></td></tr><tr><td>PL-ZNF-Am100-Lz12</td><td>3.78</td><td>1.398</td><td>0.1367</td><td>0.1777</td><td>0.6574</td><td>0.04502</td><td>0.0074</td><td>0.1379</td><td>0.2854</td><td>8.32 × 10<sup>-4</sup></td></tr><tr><td>PL-ZNF-Am100-lin</td><td>4</td><td>1.562</td><td>0.2195</td><td>0.1802</td><td>0.6095</td><td>0.04740</td><td>0.0045</td><td>0.1537</td><td>0.3083</td><td>3.57 × 10<sup>-3</sup></td></tr><tr><td>PL-ZNF-Am100-HR</td><td>4.81</td><td>1.354</td><td>0.1707</td><td>0.1536</td><td>0.3877</td><td>0.0331</td><td>0.0118</td><td>0.1030</td><td>0.2218</td><td>-</td></tr><tr><td>PL-ZNF-Am30</td><td>6.20</td><td>1.2849</td><td>0.1537</td><td>0.1188</td><td>0.6269</td><td>0.02848</td><td>0.00680</td><td>0.08688</td><td>0.1832</td><td>1.21 × 10<sup>-3</sup></td></tr><tr><td>PL-ZNF-Am10</td><td>7.21</td><td>1.224</td><td>0.1403</td><td>0.1601</td><td>0.2864</td><td>0.02307</td><td>0.02601</td><td>0.04993</td><td>0.1485</td><td>4.32 × 10<sup>-4</sup></td></tr><tr><td>PL-ZNF-Am10-Lz12</td><td>7.46</td><td>1.220</td><td>0.1045</td><td>0.1582</td><td>0.2845</td><td>0.02152</td><td>0.0249</td><td>0.04725</td><td>0.1405</td><td>3.17 × 10<sup>-4</sup></td></tr><tr><td>PL-ZNF-Am10-lin</td><td>7</td><td>1.222</td><td>0.1476</td><td>0.2075</td><td>0.2899</td><td>0.02456</td><td>0.02456</td><td>0.05657</td><td>0.1585</td><td>-</td></tr><tr><td>PL-ZNF-Am3</td><td>9</td><td>1.180</td><td>0.1235</td><td>0.17652</td><td>0.07058</td><td>0.01727</td><td>0.03389</td><td>0.01472</td><td>0.1</td><td>-</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Top: snapshot of mid-plane density for AT-GHydro ($T^4$ cooling) at an arbitrary time in steady state. Density perturbation is dominated by GI-driven spirals density waves. Bottom: same but for run AT-GHydro-lin (linear cooling). Small-scale fluctuations seem more visible in the latter case.</figcaption></figure>

the time-averaged stress to follow the Gammie (2001) relation,

$$ \alpha \approx \frac{1}{q\Omega(\gamma - 1)\tau_{\text{eff}}} = \frac{1}{\Omega\tau_{\text{eff}}} \quad (\gamma = 5/3, q = 3/2) \quad (28) $$

The values of $\alpha$ are shown in Table 2; they follow the Gammie (2001) relation within a relative error of less than 5 per cent.

For both the linear- and $T^4$-cooling runs, the turbulence is super-sonic, highly compressible, and characterized by large-scale spiral density waves, which are particularly strong for the values of $\tau_0$ and $\tau_c$ considered. A snapshot of the density structure is shown in Fig. 1 for the runs AT-GHydro and AT-GHydro-lin. We clearly identify three or four spiral waves with a dominant $k_y = k_{y0} = 2\pi/L_y$ mode. The shape of these spiral structures and their amplitude (or contrast) do not seem to depend significantly on the cooling law. Fig. 2 shows the kinetic-energy spectra $E_k(k_y)$ averaged in time, $z$, and $k_x$ for the linear- and $T^4$-cooling runs. We find that the two spectra are very similar, except perhaps that small-scale motions are more prominent in the linear-cooling run. On large scales the power spectrum agrees closely with $k^{-2}$ but for smaller scales ($k_y H_0 \gtrsim 4$) the kinetic energy is steeper. These slopes are similar to those found by Booth &

<figure><img src="image_4.png" /><figcaption>**Figure 2.** 1D kinetic energy spectra $E_k$, averaged in time, $z$, and $k_x$ for AT-GHydro and AT-GHydro-lin. The two are very similar, except slightly higher amplitudes at small scale (high $k_y$) for the linear cooling run.</figcaption></figure>

Clarke (2019). In both cases, spiral waves are accompanied by strong axisymmetric structures with $k_y = 0$ and $k_x = k_{x0} = 2\pi/L_x$, which have been shown to trigger small-scale inertial modes via a parametric instability in GI turbulence (Riols et al. 2017a).

The left and central panels of Fig. 3 shows the vertical profiles of the rms radial and vertical velocities, $\sqrt{\overline{u_x^2}}$ and $\sqrt{\overline{u_z^2}}$ respectively, for run PL-GHydro ($T^4$ cooling; solid line) and run PL-GHydro-lin (linear cooling; dashed line). Around the midplane and below $z = 2.5H_0$, the two runs share the same behaviour. In the mid-plane ($z = 0$), the velocity fluctuations are dominated by radial motions. However, higher up in the disc corona, the vertical velocity component becomes comparable. These profiles issue from the combination of the vertical roll motions that accompany the spiral waves at $z \lesssim H_0$ (see Riols & Latter 2018b), and small-scale inertial modes that attack these spirals at all altitudes (but with some predominance at $z \gtrsim H_0$). The main difference between the two runs occurs near the box boundaries at $z \gtrsim 2.5H_0$. In the case of $T^4$ cooling, the vertical rms velocity fluctuation is $\sim 0.05c_{s0}$, while for the case of linear cooling it is $\sim 10$ times higher. As a result, the linear-cooling run exhibits significantly more outflow (although it remains small overall), with mass-loss rates $\dot{m}/(\Sigma\Omega) \approx 1.1 \times 10^{-4}$ (see Table 2). We checked that this difference is obtained both with PLUTO and with Athena++. Finally, the right-hand panel of Fig. 3 shows that the thermal structure of the disc is different depending on the cooling prescription. The corona is almost two times hotter in the case of a linear cooling. However, these differences are only minor, since the results (in particular the dynamo behaviour) depend only weakly on the cooling law (see Section 3.4).

## 3.2 Gravito-turbulent dynamo in ideal MHD

Having established consistency between the two codes and determined a converged state of hydrodynamical gravito-turbulence, we now inject an initially weak toroidal seed field into the steady state and evolve the system using ideal MHD, i.e. without ambipolar diffusion. In this case, magnetic diffusion is purely numerical and occurs on the grid scale. The geometry of the initial field is set by equation (15) with $\beta_{y0} \ge 10^4$. This large value of $\beta_{y0}$ ensures that the Lorentz force is negligible during the first few hundred $\Omega_0^{-1}$, at least on large scales, and thus any dynamo action is purely kinematic. We focus first on the $T^4$ cooling runs (AT-ZNF-ideal and PL-ZNF-ideal).
<figure><img src="image_5.png" /><figcaption>Figure 3. Vertical profile of rms radial (left), vertical (centre) velocities, and sound speed (right). Quantities are averaged in x, y and time from the PLUTO hydrodynamical runs.</figcaption></figure>

We find that the magnetic field grows rather slowly, on time-scales of 250 $\Omega_0^{-1}$ (PLUTO) to 500 $\Omega_0^{-1}$ (Athena++) before reaching a quasi-steady state. In saturation, the magnetic energy remains smaller than the kinetic energy, with $E_m = 0.047\langle U \rangle$ for PLUTO and $E_m = 0.074\langle U \rangle$ for Athena++. The ideal MHD run obtained with Athena++ also supports a greater Maxwell stress than the PLUTO run. The discrepancy has two possible origins. First, the dynamo saturation and growth rates are likely sensitive to the amount of numerical dissipation (see RL19), and the time-integrator and the orbital advection algorithm are different between the two codes. In particular, the ideal MHD saturated state depends on resolution (see Section 5) and on the numerical scheme used to reconstruct the EMF at the cell edges (each scheme producing different amount of numerical dissipation, see Appendix B). Secondly, we do not integrate the simulations long enough for the quantities to be statistically meaningful. In particular, the magnetic energy $E_m$ undergoes significant and aperiodic fluctuations on long time-scales $\sim 1000 \Omega^{-1}$, suggesting that the uncertainty related to a finite averaging window could be large. However, we emphasize that this bursty behaviour concerns only the ideal runs; when including ambipolar diffusion, the magnetic energy much more regular and smooth in time.

## 3.3 Gravito-turbulent dynamo with ambipolar diffusion

We now explore the behaviour of the dynamo when ambipolar diffusion is included. Unlike Ohmic dissipation, ambipolar diffusion is a non-linear process with an associated electric field that is cubic in B. Therefore, we expect that in the kinematic regime (at least when B is very weak) ambipolar diffusion plays a negligible role in initial growth of the magnetic field, compared to numerical diffusion, and the results will be almost identical to those of ideal-MHD. Thus the calculation of growth rates during the kinematic phase, such as performed by RL19, is not meaningful in the presence of ambipolar diffusion only. Instead, we focus on the level of magnetic energy attained at saturation. To calculate these saturated states, we start the simulation either from the ideal MHD state or from an already existing state having ambipolar diffusion. In this way, we scan a large number of Am (ambipolar Elssaser number) ranging from 1 to 200. Table 2 summarizes the different box-averaged quantities obtained at saturation for various Am.

Before going into the box average analysis, we start by showing a view of the flow and the topology of the magnetic field. Fig. 4 shows

a snapshot of one of the simulations with Am = 10. We see that magnetic structures reside at large scales and are correlated with the spiral density waves (we will return to this point later in Section 4).

We now first focus on the $T^4$ cooling runs for which $\tau_0 = 10 \Omega^{-1}$. Fig. 5 (top) shows the magnetic energy $E_m$ and kinetic energy $E_k$ (normalized to thermal energy $\langle U \rangle$) as a function of Am. Circular markers correspond to Athena++ runs, while diamond markers correspond to PLUTO runs. When explicit diffusion is added, both codes agree well with one another, suggesting that numerical diffusion play a lesser role than in the ideal case (even for Am = 100). As discussed later in Section 5, results are weakly dependent on resolution, although the grid is probably influencing the small-scale dynamics (especially for Am = 100).

We see that the dynamo is strongest when Am = 30 and Am = 100, for which we have $E_m \sim 0.6\langle U \rangle$. In this regime, the magnetic energy is typically 3–5 times larger than the turbulent kinetic energy. This corresponds to a regime of super-equipartition fields, for which we expect important feedback of the Lorentz force on the gas flow. Indeed, the kinetic energy is reduced by a factor of almost 2 compared to the case of unmagnetized gravito-turbulence, while the gravitational stress is reduced by a factor 8 (see bottom panel of Fig. 5). However, Fig. 5 (bottom) shows that the transport coefficient $\alpha$ is maximal when Am = 100, for which it is three times larger than in the hydrodynamic case. Energy is mainly extracted from the shear through the Maxwell stress in this regime.

The drop in gravitational stress and kinetic energy is mainly because GI is weakened by strong magnetic fields. One way this occurs is through magnetic tension; tangled magnetic fields can disrupt the hydrodynamical flow and in particular slow down the motions associated with spiral waves. Another way is by decreasing the GI growth rate: the effective Toomre parameter, including magnetic pressure, is

$$Q_{\text{eff}} \approx Q_{\text{hydro}} \sqrt{1 + \frac{1}{\beta_t}}, \qquad (29)$$

where $\beta_t = \langle P \rangle / E_m$ (Kim & Ostriker 2001). This is 40 per cent greater than the hydrodynamic Q for $\beta_t \sim 1$. Moreover, in the presence of magnetic fields, there is an increase of temperature and, consequently, an increased Q. To illustrate this effect, we plot in Fig. 6 the evolution of the box-average Toomre parameter $Q_W$ as a function of time for Am = 100. As the field is amplified and reaches $\beta_t \sim 1$, $Q_W$ increases by 20 per cent for $T^4$ cooling (and even by a larger factor
<figure><img src="image_6.png" /><figcaption>Figure 4. Snapshots taken at $z = 0$ for AT-ZNF-Am10 at an arbitrary time in steady state. The dominant spiral structures have $n_x \approx 3-4$.</figcaption></figure>

for linear cooling, see next section). By combining the effect of magnetic pressure and the rise of thermal pressure, the effective Toomre parameter in the highly magnetized state (in particular Am = 100) can be 75 per cent higher than in pure hydrodynamics, causing a drastic reduction of GI activity.

Fig. 5 indicates that the dynamo is weakened for Am $\lesssim$ 10. When Am = 1, the magnetic energy has been reduced by a factor 30 compared to Am = 100. Because ambipolar diffusion depends nonlinearly on B, it is possible that the dynamo survives no matter how small Am, but will support an infinitely small amount of magnetic energy at saturation. In the other limit, for Am $\gtrsim$ 100, we see from Fig. 5 that the strength of the dynamo is reduced as Am is increased. This behaviour is reminiscent of simulations with Ohmic dissipation, and is due to the kink instability developing at small scales, which can weaken the dynamo (see discussion in RL19).

To summarize, our results show that there is little difference in terms of box-average quantities between discs regulated by ambipolar diffusion and those controlled by Ohmic dissipation (RL19). The dynamo strength peaks at intermediate Am $\sim$ 30–100. To be more quantitative and allow comparison with Ohmic runs, we calculate the effective magnetic Reynolds number defined as :

$$Rm_{\text{eff}} = \left\langle \frac{\Omega H^2}{\eta_A} \right\rangle = \left\langle \frac{c_{s0}^2}{v_A^2} \right\rangle \text{Am} \quad (30)$$

When Am = 30 and Am = 100, we find $Rm_{\text{eff}} \sim 0.8$Am. For comparison, dynamo activity was observed to be strongest when Rm $\sim$ 20–100 in the case of Ohmic dissipation (RL19). The dynamo

behaves like a slow dynamo, with magnetic energy decreasing with increasing Am above a critical value. However, we emphasize that we may not be able to probe the behaviour of the dynamo accurately at larger Am at our fiducial resolution, because for Am $\gtrsim$ 200, the minimum length scale of magnetic perturbation that is not strongly damped by ambipolar diffusion (i.e. damping rate at scale $\lambda$, $\eta_A/\lambda^2$, is $\lesssim \Omega$) becomes comparable to or smaller than the grid scale. It is likely that, in this case, the grid (rather than explicit diffusion) controls the small-scale dynamics of the field.

### 3.4 Dependence on the cooling law

So far we have analysed the MHD simulations with $T^4$ cooling. In this section, we study how those results depend on the cooling law and make comparison with the linear cooling runs. We start with the case of ideal MHD, without explicit diffusion. The cooling time $\tau_c$ in Athena++ is chosen so that it is equal to the average effective cooling time-scale $\tau_{\text{eff}}$ measured in the ideal $T^4$ cooling run, rounded off to the next 0.5. We found that the magnetic energy is twice larger when a linear cooling is assumed, and it exceeds the kinetic energy (see Table 2). As explained in Section 3.2, this difference may be due to the relatively short averaging window (note also that the linear and $T^4$ cooling runs are not averaged on the same time-scale). In PLUTO the cooling time is chosen to be equal to the effective $\tau_{\text{eff}}$ corresponding to the hydrodynamical $T^4$ run, which appears to be close to the effective $\tau_{\text{eff}} = 10.56\Omega^{-1}$ measured in the ideal MHD run. Table 2 shows that the magnetic energy and turbulent properties
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Top: box averaged magnetic energy (red) and kinetic energy (blue), normalized to the averaged thermal energy, as functions of Am for $\tau_0 = 10 \Omega^{-1}$. Bottom: transport efficiency $\alpha$ as a function of Am. All quantities are measured in the saturated state. The diamond and circle markers correspond, respectively, to the PLUTO and Athena++ runs. For comparison, we plot at the right edge of each panel the same quantities for the ideal case (without explicit ambipolar diffusion).</figcaption></figure>

(stress, $Q_W$, etc.) at saturation depend weakly on the choice of the cooling law.

We then compare $T^4$ and linear cooling when explicit ambipolar diffusion is added. A comparison is undertaken for Am = 10 and Am = 100 (see Table 2). The cooling time $\tau_c$ is again chosen so that it is equal to the average $\tau_{eff}$ measured in the corresponding $T^4$ cooling run. For Am = 10 we then have $\tau_c = 7\Omega^{-1}$ and for Am = 100 we have $\tau_c = 4\Omega^{-1}$. Table 2 shows that there is a good agreement in terms of magnetic energy and stress between the $T^4$ and the linear cooling runs. One notable difference is perhaps for Am = 100 where $Q_W$ is slightly larger (and consequently $G_{xy}$ smaller) when a linear cooling law is used.

In addition to that, Fig. 6 (bottom) shows that the Toomre parameter exhibits strong oscillations in the case of linear cooling, which are not present in the $T^4$ cooling run. These oscillations are correlated with the variations of the magnetisation (or $\beta_t$) and were interpreted in Riols & Latter (2018a) and RL19 as a thermal cycle involving a switch on and off of the GI activity. One possibility is that this cycle is initiated by a thermal instability, which is suppressed by the steep dependence of the $T^4$ cooling law. Further investigations will be necessary to understand the nature of this cycle.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Time-evolution of turbulent beta plasma $\beta_t = \langle P \rangle / E_m$ (top) and Toomre parameter (bottom) for Am = 100. Blue/solid curves are for the $T^4$ cooling run with $\tau_0 = 10 \Omega^{-1}$ (PL-ZNF-Am100), while the red/dashed curves are for the linear cooling run (PL-ZNF-Am100-lin), restarting from a $T^4$ cooling state.</figcaption></figure>

## 3.5 Fourier modes and butterfly diagrams

As in RL19, we checked whether the GI dynamo in the presence of ambipolar diffusion preferentially amplifies large or small-scale fields. For that purpose, we show in Fig. 7 (top left panel) the evolution of magnetic energy associated with different Fourier modes in the runs PL-ZNF-Am10 (Am = 10, $T^4$ cooling). The calculation is restricted to the midplane $z = 0$, as it contains most of the magnetic energy. Clearly the dominant component of the magnetic field is the $k_x = k_y = 0$ mode (mean field $\bar{B}$), which is essentially toroidal. The dynamo in this regime is thus large scale. We see that the next important component is the $k_y = 2\pi/L_y$ mode associated with the large-scale spiral waves. Although significantly smaller than the mean field (its energy being about eight times smaller), this component is necessary for the large-scale dynamo process. As we will show, it contributes to the turbulent EMF and helps to regenerate the mean field against ambipolar diffusion. Finally, Fig. 7 suggests that the axisymmetric modes ($k_y = 0$) with $k_x > 0$ are quite negligible in the magnetic budget.

In the lower left panels of Fig. 7, we show the time-evolution of the vertical profiles of the mean magnetic field $\bar{B}_x(z)$ (centre) and $\bar{B}_y(z)$ (bottom). The large-scale field remains confined within the mid-plane region of the disc and is mainly toroidal with $\bar{B}_y/\bar{B}_x \simeq 20$ at $z = \pm H_0$. Both poloidal and toroidal components have the same polarity in the mid-plane, but $\bar{B}_x$ changes its sign at $z \simeq \pm 0.2 H_0$. Note that the evolution of the radial field show strong variations in the mid-plane which are probably associated with the intermittent nature of GI activity. Finally an important feature is that the large-scale field does not flip (the polarity remains fixed) and is symmetric about the mid-plane (i.e. it has an 'even' geometry). All these properties are similar to those obtained in simulations with explicit Ohmic resistivity (RL19). Note that we obtain a similar field geometry for all zero-net flux simulations with ambipolar diffusion (even for Am =
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Top: time evolution of magnetic energy contained within the mid-plane ($z=0$), projected on to different Fourier harmonics, for the $T^4$ cooling law with $\tau_0 = 10\Omega^{-1}$. Blue/solid lines correspond to the magnetic energy of the mean field ($k_x = 0, k_y = 0$), while red/dashed lines correspond to its toroidal projection. Purple/dashed lines represent the magnetic energy projected on to the axisymmetric modes ($k_x \neq 0, k_y = 0$). Green and cyan solid lines correspond to the energy of non-axisymmetric components, respectively $k_y = 2\pi/L_y$ and $k_x = 4\pi/L_y$. Lower panels: space-time $(t, z)$ diagram of the mean $\bar{B}_x$ and $\bar{B}_y$. The left-hand panels are for $Am = 10$ (run PL-ZNF-Am10) and the right-hand panels are for the ideal case (run PL-ZNF-ideal).</figcaption></figure>

200). Also the profile of $\bar{B}_x(z)$ and $\bar{B}_y(z)$ remains similar when increasing the vertical box size (see Section 5).

As a comparison, we computed the same diagnostics for the ideal run PL-ZNF-ideal. Fig. 7 (bottom right) shows that the mean field flips many times during the course of the simulation but its evolution is irregular (it can be either antisymmetric or symmetric about the mid-plane). The period of the flipping is about $200 \Omega^{-1}$ in the saturated state. Fig. 7 (top right) shows that, unlike the case $Am = 10$, the mean component of the field $\bar{B}$ in the mid-plane is not the dominant one. Even when it is symmetric about the mid-plane (or at least when it is not changing polarity at $z=0$), its energy is comparable to the non-axisymmetric mode $k_y = 2\pi/L_y$.

## 4 GENERATION OF MEAN FIELD BY SPIRALS

In this section, we investigate how the mean field is generated and sustained by GI motions. Our focus will be on the role of spiral density waves; as has been shown by RL19, spirals are the key to explaining the generation of $\bar{B}$ in discs regulated by Ohmic dissipation. One would like to know whether or not these structures are also important when ambipolar diffusion dominates.

### 4.1 EMF budget

We show in Fig. 8 the vertical profiles of the EMFs and the terms involved in the evolution of $\bar{B}_x$ (top) and $\bar{B}_y$ (bottom) for $Am = 10$ and $T^4$ cooling law (run AT-ZNF-Am10). This includes the spiral and turbulent part of the mean ideal EMF ($\bar{\mathcal{E}}^*$ and $\bar{\mathcal{E}}^t$ defined in Section 2.4.4), as well as the shear stretching (the ‘omega effect’) and the ambipolar EMF $\bar{\mathcal{E}}^{NI}$ (see equations (25 and 26)). Recall that the spiral term $\bar{\mathcal{E}}^*$ is computed using the phase-folded technique (see Section 2.4.3), using a spiral pattern with $k_{x_s} = 4k_{x_0}$ with $k_{x_0} = 2\pi/L_x$ (this corresponds to the most prominent spiral mode seen in the Athena++ and PLUTO runs). All the terms are averaged in time, during the saturated phase (~250 $\Omega^{-1}$), and over the horizontal plane.

Let us begin with the toroidal component, $\bar{B}_y$. Above the mid-plane for $|z| > 0.25H_0$, it is clear from the bottom right panel of Fig. 8 that $\bar{B}_y$ is generated mainly via the background shear flow (red line). In the central panels, the radial turbulent EMF (orange) has a negative gradient in $z$, like ambipolar diffusion (green), and therefore acts as a turbulent diffusion on the mean field. In the mid-plane ($|z| < 0.25H_0$), however, the configuration is different since
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Time-averaged vertical profile of mean magnetic field $\bar{B}_x$, $\bar{B}_y$; mean electromotive force $\bar{\mathcal{E}}_x$, $\bar{\mathcal{E}}_y$; and $d\bar{B}_x/dt$, $d\bar{B}_y/dt$ for run AT-ZNF-Am10 in its steady state. The behaviours are characteristic of runs with relatively high Am. We separate the EMF into several components: $\mathbf{u} \times \mathbf{B}$ in the $n_x = 4$, $n_y = 1$ phase folded profile, which mainly accounts for the contribution from $n_x = 4$, $n_y = 1$ spirals waves (blue curves); the remaining $\mathbf{u} \times \mathbf{B}$, which mainly accounts for turbulent EMF (orange curves); and the non-ideal EMF, which in this case accounts for ambipolar diffusion. We use the same colour code for $d\bar{B}/dt$ due to these components of the electric field, and the red curve in the $d\bar{B}_y/dt$ panel corresponds to the generation of $\bar{B}_y$ through the shearing of $\bar{B}_x$. The divergences seen at the boundaries are probably due to the boundary conditions, but they do not affect the dynamo process, which happens mainly near the mid-plane (see Section 5.2).</figcaption></figure>

$\bar{B}_y$ and $\bar{B}_x$ have the same polarity, and therefore the omega effect tends to destroy the large-scale field (this explains the little dip in the $\bar{B}_y$ profile at $z=0$). We note that the EMF budget for $\bar{B}_y$ is actually similar to that obtained with Ohmic resivity (see fig. 8 of RL19). Despite its negative contribution in the mid-plane, the omega effect is essential in the vertically averaged budget, and therefore a poloidal field is required for the large-scale dynamo to work.

We then study how the mean $\bar{B}_x$ is produced. The top panels of Fig. 8 shows the spiral, turbulent, and ambipolar EMFs profiles associated with the radial field $\bar{B}_x$. For $|z| < 0.25H_0$ (mid-plane region), the top right panel shows that $\bar{B}_x > 0$ is generated mainly via the spiral EMF in $y$ (blue) and dissipated through the non-ideal ambipolar EMF (green) and the turbulent EMF (orange). Above the mid-plane, for $|z| > 0.25H_0$ and $|z| \lesssim 1.5H_0$, $\bar{B}_x < 0$ and the gradients of $\bar{\mathcal{E}}_y^t$, $\bar{\mathcal{E}}_y^*$, and $\bar{\mathcal{E}}_y^{N1}$ are just reversed, meaning again that the spiral EMF helps to sustain the mean radial field against ambipolar and turbulent diffusion.

In summary, we have shown that the mean field dynamo relies on the production of $B_y$ from $B_x$ by the omega effect (primarily) and on the generation of $B_x$ by relatively large-scale fluctuations issuing from the spiral density waves.

## 4.2 Dynamo process and phase-folded profiles

As demonstrated in the previous subsection, spiral motions are key to regenerating the mean field $\bar{B}_x$. But how exactly do they operate? In the presence of Ohmic dissipation, RL19 proposed a physical picture that highlights the role of large-scale vertical rolls associated with spiral density waves in the generation of $\bar{B}_x$ (see their sections 5.1 and 5.2). These rolls are a fundamental feature of linear density waves in stratified discs and arise from a baroclinic effect associated with the thermal stratification of the disc (Riols & Latter 2018b).

The mechanism of radial field generation from a toroidal field takes place through the following steps (see in particular fig. 12 of RL19):

(i) First, spirals horizontally compress and tilt the initially toroidal field. This creates a negative (positive) $B_x$ inside (outside) of the spiral wave. On average, no mean radial field ($\bar{B}_x$) is produced via this process.

(ii) At the centre of the spiral, the field is pinched and lifted by the vertical velocity associated with the roll structure. This carries some magnetic flux to higher altitude and creates a $\bar{B}_z$ component with opposite polarity on both sides of the spiral.

(iii) As the rolls 'turn over' above the disc at $z \sim H_0$, the field is then stretched in the radial direction and folded, producing a net $\bar{B}_x$ with opposite sign to $B_y$ in the corona (through the EMF term $\tilde{u}_x \bar{B}_z$).

One may wonder if the mechanism described above and detailed in RL19 still hold in discs dominated by ambipolar diffusion. To this end, we examined the gas motions and the magnetic topology around spiral waves for Am = 10. Instead of looking at fixed snapshots, as in RL19, we compute the 'mean' motion and magnetic field associated with the spirals by phase folding the data (see Section 2.4.3 for details about this technique). This procedure offers a visualization of the flow in a poloidal plane $(x, z)$, in which perturbations associated with the most prominent spiral mode are preserved but turbulent and uncorrelated fluctuations are suppressed. The result is shown in Fig. 9. The most prominent spiral mode in the Athena++ and PLUTO runs is the mode $k_y = 2\pi/L_y$ and $k_x = 4k_{x0} = 8\pi/L_x$ (see Fig. 1, which shows four spiral arms in the box). First, we clearly identify in the top left panel of Fig. 9 the pattern of four counter-rotating, large-scale rolls associated with each density maximum. These rolls are located within $z \pm H_0$, but the radial motions directed outward from the density wave persist at larger $z$ in the corona (see also $v_x$ panel).
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Phase-folded profile in the steady state of run AT-ZNF-Am10. Turbulent fluctuations are suppressed by phase-folding multiple snapshots, and the velocity rolls associated with the spiral density waves and the twisting of field lines due to the velocity rolls are clearly visible, despite their amplitudes being often much smaller than those of the turbulent fluctuations.</figcaption></figure>

Fig. 9 shows that the toroidal field $B_y$ in the mid-plane is compressed within the spiral density maxima (see $B_y$ panel) and twisted along the direction of the spirals, which produces a negative $B_x$ component inside the spirals and a positive component outside (see $b_x$ panel). Fig. 4 (top right panel) shows a snapshot of the magnetic-field lines projected in the mid-plane and confirms this picture. This affirms step (i). The vertical motions associated with the velocity rolls of the spirals then bring the negative $B_x$ upwards into the corona (see $B_x$ panel). This leads to a mean $\bar{B}_x > 0$ in the mid-plane and a mean $\bar{B}_x < 0$ above, consistent with Fig. 7. Another way to check this result is to see that, at the location of the density maxima where $B_x < 0$, the vertical velocity of the rolls is always directed outwards (positive $\tilde{u}_z$ if $z > 0$ and negative if $z < 0$). In the inter-arm regions, we observe the opposite configuration. A clear correlation is then obtained between $u_z$ and $B_x$, which generates a positive gradient of $-\mathcal{E}_y$ and then a positive $\bar{B}_x$ near the mid-plane. The $B_z$ panel of Fig. 9 shows that, for $z > 0$, a positive vertical field is produced on the left side of the spiral wave front and a negative $B_z$ on the other side of the wave front. This configuration can be explained simply by the fact that the field line, which crosses the spiral wave at some angle,

is pinched and lifted up at the centre of the spiral by the rolls. All of these elements confirm that step (ii) is at work in our simulations. Finally, the $B_x$ panel of Fig. 9 shows red ‘mushroom-like’ structure in the contours of $B_x$, reflected on either side of the midplane. These structures are similar to those observed in the snapshots of fig. 15 of RL19 and indicate that the field that reaches $z \simeq H_0$ is stretched horizontally by the counter-rotating motions of each roll. We see that, in the corona, the vertical component of the field is correlated with the rolls’ horizontal motions (in regions where $u_x < 0$, $B_z$ is positive, and vice versa). Such correlations lead to a negative EMF $-\mathcal{E}_y = \tilde{u}_x \bar{B}_z$, resulting in a negative gradient for the EMF and therefore a negative $\bar{B}_x$, as expected from Fig. 7. The physical processes described in step (iii) seem to be at work in our simulations.

While generating $\bar{B}_x$ with an opposite sign to $\bar{B}_y$ in the corona, the dynamo mechanism above tends to generate $\bar{B}_x$ with the same sign as $\bar{B}_y$ at the mid-plane, which is visible in Fig. 8. This effect, together with the shear of $\bar{B}_x$ to generate $\bar{B}_y$, leads to cyclic flipping of the polarity in the absence of diffusion, and may explain the dependence of field geometry (flipping or fixed) on the level of diffusion seen in Fig. 7 and RL19.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Spectra of magnetic energy for runs AT-ZNF-Am10 (fiducial resolution) and AT-ZNF-Am10-HR (double resolution). Increasing resolution leads to more small-scale (high $k_y$) fluctuations. However, the amplitudes of large-scale (low $k_y$) modes, which dominate the magnetic energy, are similar.</figcaption></figure>

# 5 DEPENDENCE ON RESOLUTION AND BOX SIZE

## 5.1 Resolution

To study whether our fiducial resolution is sufficient to capture all of the relevant dynamics, we re-perform the set of zero-net-flux, $T^4$ runs at double resolution (512 × 512 × 288) for ideal MHD, $Am = 100$ and $Am = 10$ with Athena++, and $Am = 100$ with PLUTO. In Athena++, we restart the simulations from an MHD quasi-steady turbulent state at fiducial resolution and interpolate the magnetic field so that $\nabla \cdot \mathbf{B} = 0$ is preserved. (Because the code stores magnetic fluxes through cell interfaces, a simple linear interpolation of magnetic flux density on to new interfaces guarantees that the total magnetic flux through a cell remains zero when doubling resolution.) In PLUTO, we start from a gravito-turbulent hydrodynamical state and add a strong toroidal magnetic perturbation in the mid-plane of the form $B_y = 0.2 \exp(-z^2/2H_0^2)$.

When $Am = 10$, the fiducial-resolution run AT-ZNF-Am10 and the double-resolution run AT-ZNF-Am10-HR show very similar box-averaged quantities (Table 2). The magnetic energy spectra for the two runs are also very similar, although the double-resolution runs show more small-scale fluctuations (Fig. 10). To confirm that these additional small-scale fluctuations do not affect large-scale dynamics, we calculate the rate of energy generation and transport in the steady state for the two simulations (Fig. 11). Most of the transport rates are barely affected by the increase in resolution. In particular, the similarity in the grid-scale dissipation of magnetic energy between the two resolutions, despite the different level of small-scale fluctuations, suggests that the energy cascade rate is set by large-scale processes that are relatively insensitive to the grid resolution.

When $Am = 100$, the box-averaged quantities in the fiducial- and double-resolution runs (with Athena++) appear to be quite different but, as Fig. 12 shows, this is mainly because the double-resolution run contains a flip in the sign of the large-scale magnetic field while the fiducial run does not. If we exclude the flip from the averaging, the averaged quantities of the two runs becomes similar (the magnetic energy in the double-resolution run is ≈70 per cent of that obtained at the fiducial resolution). The increased tendency for

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Average rate of energy generation and transport in steady state in runs AT-ZNF-Am10 (top) and AT-ZNF-Am10-HR (bottom). The energies are normalized by the thermal energy, and energy generation and transport rates are normalized by the thermal energy × Ω. Shear generates kinetic energy through Reynolds and gravitational stresses, and magnetic energy through Maxwell stress. We split the conversion of magnetic to thermal energy into two channels: numerical dissipation at the grid scale and ambipolar heating.</figcaption></figure>

flipping, on the other hand, is probably because of the decreased numerical dissipation. Previously, we showed that non-ideal diffusivity tends to suppress flipping, and it is likely that numerical dissipation can produce a similar effect.

For ideal MHD, the fiducial- and double-resolution runs give significantly different box-averaged quantities, especially for the magnetic energy. The lack of numerical convergence is to be expected, since we believe small-scale kink instabilities, which are sensitive to the level of numerical dissipation in the absence of non-ideal effects, affect the strength of the dynamo (RL19).

## 5.2 Vertical box size

Another important test, especially for dynamo problems, is the dependence of the dynamo properties on the vertical box size. Our fiducial simulations are run with $L_z = 9H_0$, and so we extend the box size to $L_z = 12H_0$, keeping the same resolution per $H_0$. Using PLUTO, we first re-computed the hydrodynamical states with the extended box assuming $T^4$ cooling. We found very similar box-averaged quantities between $L_z = 9H_0$ and $L_z = 12H_0$ (Table 2). We then performed two MHD runs with ambipolar diffusion, having $Am = 10$ and $Am = 100$ (using PLUTO). To avoid the long kinematic phase, we introduced a strong magnetic perturbation at the beginning of the simulation of the form $B_y = 0.1 \exp(-z^2/8H_0^2)$. Fig. 13 (top
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Evolution of magnetic energy and butterfly diagrams for $\bar{B}_x$, $\bar{B}_y$ in runs AT-ANF-Am100 (fiducial resolution) and AT-ANF-Am100-HR (double resolution). The large difference in magnetic energy in Table 2 is mainly because the double-resolution run contains a flip of magnetic polarity, during which the magnetic energy is low due to the lack of a strong mean field. Before and after the flip, the magnetic energy remains similar to that obtained at fiducial resolution.</figcaption></figure>

panels) shows the vertical profiles of $\bar{B}_x$ and $\bar{B}_y$ for $L_z = 9H_0$ and $L_z = 12H_0$, averaged during the saturated state. Clearly, for both Am = 10 and Am = 100, the magnetic profiles are almost independent of vertical box size. The bottom panels of Fig. 13 show the ideal part of the Poynting flux as a function of z for both $L_z = 9H_0$ and $L_z = 12H_0$. There is again little difference between the two, indicating that the amount of magnetic energy that leaves the disc surface is almost identical. We also checked that the EMF profiles are not affected by increasing $L_z$. All of these results suggest that the dynamo properties are independent of vertical box size.

## 6 DISCUSSION AND CONCLUSIONS

In conclusion, we conducted a series of 3D MHD shearing-box simulations of self-gravitating discs, stratified in the vertical direction and with zero net flux, to understand whether GI is able to generate and sustain a large-scale magnetic field. The physical and numerical setup is similar to that used in RL19. The main difference is that we examined the effect of ambipolar diffusion (rather than Ohmic dissipation) on the dynamo, a non-ideal effect that is more likely to impact the GI-unstable regions of poorly ionized protoplanetary discs. We examined the fundamental properties of the dynamo and its saturation for a wide range of ambipolar Elsasser numbers (Am). We also studied the dependence of the dynamo on the cooling law, resolution and box size. We obtained the following results:

(i) The GI dynamo in the presence of ambipolar diffusion is optimal for ambipolar Elsasser numbers between Am = 30 and 100. In that regime, the magnetic field reaches thermal strength, with $\beta \sim 1$. The field is mainly large scale (approximately box size in the horizontal direction) and dominated by the mid-plane toroidal component. At larger Am, the magnetic energy weakens, which suggests that the dynamo is ‘slow’. These results are reminiscent of the case with Ohmic dissipation (RL19). At small Am (Am < 20), the field strength decreases and the magnetic to thermal energy ratio scales as:

$$E_m/E_T \sim 0.02 \text{ Am}^{1.15} \qquad (31)$$

Therefore the field strength is expected to vary according to the following power law:

$$B \sim 0.15 \text{ Am}^{0.57} (\Sigma H)^{1/2} \Omega \qquad (32)$$

(ii) The mechanism of field generation is similar to that described in RL19 and relies on the vertical circulations (poloidal ‘rolls’) associated with spiral density waves. We used a new phase-folding procedure that made it especially clear to see these characteristic motions and the consequent field topology.

(iii) The saturation level of the dynamo is independent of the cooling law, except that we obtained more variability (in particular, thermal cycles) when a linear cooling law is employed. The origin of these cycles is still not well understood.

(iv) Finally, the dynamo is independent of the code used and robust to various numerical details (e.g. vertical box size, EMF reconstruction algorithm). We found that resolution does not affect the results provided that Am is sufficiently small (Am $\lesssim$ 100). However, in the ideal-MHD limit (without explicit diffusion), the results are not converged with resolution, suggesting that diffusion plays a crucial role in the amplification of the field. The behaviour of the dynamo as Am tends to infinity remains difficult to predict.

The dynamo process that we study in this paper has potential direct implications for magnetic-field generation in the outer regions (typical beyond 20 AU) of young and massive protoplanetary discs. It is expected that in these regions the Class 0/I discs are subject to ambipolar diffusion and self-gravity, though the values of Am in the mid-plane of such discs are highly uncertain. Indeed, detailed ionization maps of young Class 0/I disc have not been produced in the literature, on contrast to classical T Tauri discs. Nevertheless, the chemical models of Masson et al. (2016) suggest a range of Am between 0.02 and 20 in the outer regions ($R > 50$ au) of nascent protostellar discs during the protostellar core collapse (note that their definition of Am is different from ours and so we have converted their Am to the one we use in this paper). According to our Fig. 5, this lies in the regime of weak dynamo, but the range of expected values for the magnetic-field strength remains fairly broad. In the best case (Am $\sim$ 20), the magnetic energy is roughly half of the thermal energy, which implies field strengths of a few tens of mG at 30 au (see discussion of RL19, for details of this estimation). This is comparable with values of the field estimated from meteoritic samples, although these values are generally given at $R \lesssim 10$ au (Fu et al. 2020). Overall, given the uncertainty over Am, it is currently a challenge to predict the intensity of the magnetic field generated by GI in Class 0/I discs.

Our work is also relevant to understanding the generation of magnetic fields in AGN discs. First, given their small thickness, AGN discs are susceptible to GI (Menou & Quataert 2001; Goodman 2003; Lodato 2007) beyond ~0.01 pc, and therefore could be subject to the GI dynamo. For the relatively low densities in these regions,
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Top panels: vertical profiles of the mean $\bar{B}_x$ and $\bar{B}_y$ averaged in time. Bottom panels: vertical profiles of the ideal Poynting flux $B^2v_z - B_z(B \cdot v)$. The right-hand panels are for Am = 10 and the left-hand panels are for Am = 100. For each panel, we compare the fiducial box with $L_z = 9H_0$ (runs PL-ZNF-Am10 and PL-ZNF-Am100) with the taller box with $L_z = 12H_0$ (runs PL-ZNF-Am10-Lz12 and PL-ZNF-Am100-Lz12).</figcaption></figure>

ambipolar diffusion can be more important than resistive diffusion, but for reasonable mass accretion rates models predict that Am ≫ 100 (Menou & Quataert 2001). Our simulations cannot probe this regime but suggest that the dynamo would not be optimal. Such values of Am remain poorly constrained and depend on the efficiency of thermal ionization and photo-ionization. It can be also influenced by the geometry of the disc in the case it is warped. In any case, AGNs (as well as protoplanetary discs) are believed to be threaded by a large-scale vertical magnetic field (Prior & Gourgouliatos 2019), which could significantly alter the dynamo properties. Future work is needed to understand the impact of this vertical field. This will be the focus of paper II.

Finally, it is worth mentioning that our work might be applied to galactic dynamos. Indeed radio observations revealed that some spiral galaxies (M81, M51, and M33) are dominated by bisymmetric-spiral magnetic fields. These magnetic spiral patterns are strongly correlated with the large-scale galactic spiral arms. Our result suggests that spiral density waves can directly produce magnetic fields in these objects as well as the observed alignment between magnetic and density spiral patterns (see Fig. 4). This scenario could compete with to the standard mean field theory which is usually invoked to explain the dynamo in galactic discs (Chiba & Tosa 1990; Moss 1997; Elstner et al. 2000; Rüdiger & Hollerbach 2004; Gressel et al. 2008).

## ACKNOWLEDGEMENTS

AR and GL have received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 815559 (MHDiscs)). High-performance computing resources were pro-

vided by the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by Grenoble research communities; and the PICSciE-OIT TIGRESS High Performance Computing Center and Visualization Laboratory at Princeton University.

## DATA AVAILABILITY

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

## REFERENCES

* Alexiades V., Amiez G., Gremaud P.-A., 1996, Commun. Numer. Methods Eng., 12, 31

* ALMA Partnership et al., 2015, ApJ, 808, L3

* Andrews S. M., 2020, ARA&A, 58, 483

* Armitage P. J., 2019, Saas-Fee Advanced Course, 45, 1

* Avenhaus H. et al., 2018, ApJ, 863, 44

* Bai X.-N., 2015, ApJ, 798, 84

* Bai X.-N., Stone J. M., 2013, ApJ, 767, 30

* Bai X.-N., Stone J. M., 2017, ApJ, 836, 46

* Balbus S. A., 2003, ARA&A, 41, 555

* Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214

* Balbus S. A., Hawley J. F., 1998, Rev. Mod. Phys., 70, 1

* Balbus S. A., Henri P., 2008, ApJ, 674, 408

* Béthune W., Lesur G., Ferreira J., 2017, A&A, 600, A75

* Booth R. A., Clarke C. J., 2019, MNRAS, 483, 3718

* Brandenburg A., Nordlund A., Stein R. F., Torkelsson U., 1995, ApJ, 446, 741

* Carrasco-González C., Rodríguez L. F., Anglada G., Martí J., Torrelles J. M., Osorio M., 2010, Science, 330, 1209

* Chiba M., Tosa M., 1990, MNRAS, 244, 714
Cournede C., Gattacceca J., Gounelle M., Rochette P., Weiss B. P., Zanda B., 2015, *Earth Planet. Sci. Lett.*, 410, 62

Dapp W. B., Basu S., Kunz M. W., 2012, A&A, 541, A35

Del Zanna L., Bucciantini N., Londrillo P., 2003, A&A, 400, 397

Deng H., Mayer L., Latter H., 2020, ApJ, 891, 154

Desch S. J., Mouschovias T. C., 2001, ApJ, 550, 314

Donati J.-F., Paletou F., Bouvier J., Ferreira J., 2005, *Nature*, 438, 466

Dong R., Vorobyov E., Pavlyuchenkov Y., Chiang E., Liu H. B., 2016, ApJ, 823, 141

Durisen R. H., Boss A. P., Mayer L., Nelson A. F., Quinn T., Rice W. K. M., 2007, Protostars Planets V, 607

Elstner D., Otmianowska-Mazur K., von Linden S., Urbanik M., 2000, A&A, 357, 129

Ferreira J., Pelletier G., 1995, A&A, 295, 807

Fleming T. P., Stone J. M., Hawley J. F., 2000, ApJ, 530, 464

Fromang S., Papaloizou J., Lesur G., Heinemann T., 2007, A&A, 476, 1123

Fromang S., Latter H., Lesur G., Ogilvie G. I., 2013, A&A, 552, A71

Fu R. R. et al., 2014, *Science*, 346, 1089

Fu R. R., Kehayias P., Weiss B. P., Schrader D. L., Bai X.-N., Simon J. B., 2020, *J. Geophys. Res. (Planets)*, 125, e06260

Gammie C. F., 1996, ApJ, 457, 355

Gammie C. F., 2001, ApJ, 553, 174

Gardiner T. A., Stone J. M., 2005, *J. Comput. Phys.*, 205, 509

Goddi C., Surcis G., Moscadelli L., Imai H., Vlemmings W. H. T., van Langevelde H. J., Sanna A., 2017, A&A, 597, A43

Goodman J., 2003, MNRAS, 339, 937

Gressel O., Elstner D., Ziegler U., Rüdiger G., 2008, A&A, 486, L35

Guilet J., Ogilvie G. I., 2013, MNRAS, 430, 822

Hawley J. F., Gammie C. F., Balbus S. A., 1995, ApJ, 440, 742

Hawley J. F., Gammie C. F., Balbus S. A., 1996, ApJ, 464, 690

Huang J. et al., 2018, ApJ, 869, L43

Käpylä P. J., Korpi M. J., 2011, MNRAS, 413, 901

Kataoka A., Tsukagoshi T., Pohl A., Muto T., Nagai H., Stephens I. W., Tomisaka K., Momose M., 2017, ApJ, 844, L5

Kim W.-T., Ostriker E. C., 2001, ApJ, 559, 70

Kunz M. W., Lesur G., 2013, MNRAS, 434, 2295

Kunz M. W., Mouschovias T. Ch., 2010, MNRAS, 408, 322

Lesur G., Ogilvie G. I., 2008, A&A, 488, 451

Lesur G., Ferreira J., Ogilvie G. I., 2013, A&A, 550, A61

Lesur G., Kunz M. W., Fromang S., 2014, A&A, 566, A56

Leung P. K. C., Ogilvie G. I., 2019, MNRAS, 487, 5155

Levy E. H., Sonett C. P., 1978, in Gehrels T., Matthews M. S., eds, IAU Colloq. 52: Protostars and Planets. University of Arizona Press, Tucson, USA, p. 516

Liu H. B. et al., 2016, *Sci. Adv.*, 2, e1500875

Lodato G., 2007, *Nuovo Cimento Rivista Serie*, 30, 293

Mann R. K., Andrews S. M., Eisner J. A., Williams J. P., Meyer M. R., Di Francesco J., Carpenter J. M., Johnstone D., 2015, ApJ, 802, 77

Martí-Vidal I., Muller S., Vlemmings W., Horellou C., Aalto S., 2015, *Science*, 348, 311

Masson J., Chabrier G., Hennebelle P., Vaytet N., Commerçon B., 2016, A&A, 587, A32

Meheut H., Fromang S., Lesur G., Joos M., Longaretti P.-Y., 2015, A&A, 579, A117

Menou K., Quataert E., 2001, ApJ, 552, 204

Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, ApJS, 170, 228

Mignone A., Flock M., Stute M., Kolb S. M., Muscianisi G., 2012, A&A, 545, A152

Moss D., 1997, MNRAS, 289, 554

Pérez L. M. et al., 2016, *Science*, 353, 1519

Prior C., Gourgouliatos K. N., 2019, A&A, 622, A122

Pudritz R. E., 1981, MNRAS, 195, 881

Rice W. K. M., Armitage P. J., Bate M. R., Bonnell I. A., 2003, MNRAS, 339, 1025

Rice W. K. M., Lodato G., Pringle J. E., Armitage P. J., Bonnell I. A., 2006, MNRAS, 372, L9

Riols A., Latter H., 2016, MNRAS, 460, 2223

Riols A., Latter H., 2018a, MNRAS, 474, 2212

Riols A., Latter H., 2018b, MNRAS, 476, 5115

Riols A., Latter H., 2019, MNRAS, 482, 3989

Riols A., Lesur G., 2019, A&A, 625, A108

Riols A., Latter H., Paardekooper S.-J., 2017a, MNRAS, 471, 317

Riols A., Rincon F., Cossu C., Lesur G., Ogilvie G. I., Longaretti P.-Y., 2017b, A&A, 598, A87

Rüdiger G., Hollerbach R., 2004, The Magnetic Universe: Geophysical and Astrophysical Dynamo Theory. Wiley-VCH, Weinheim, Germany

Sano T., Stone J. M., 2002, ApJ, 577, 534

Sano T., Miyama S. M., Umebayashi T., Nakano T., 2000, ApJ, 543, 486

Simon J. B., Lesur G., Kunz M. W., Armitage P. J., 2015, MNRAS, 454, 1117

Spruit H. C., 1996, preprint pp astro-ph/9602022

Stephens I. W. et al., 2014, *Nature*, 514, 597

Stone J. M., Tomida K., White C. J., Felker K. G., 2020, ApJS, 249, 4

Suriano S. S., Li Z.-Y., Krasnopolsky R., Suzuki T. K., Shang H., 2019, MNRAS, 484, 107

Tazaki R., Lazarian A., Nomura H., 2017, ApJ, 839, 56

Tobin J. J. et al., 2013, ApJ, 779, 93

Vlemmings W. H. T. et al., 2019, A&A, 624, L7

Wardle M., 2007, Ap&SS, 311, 35

Wardle M., Koenig A., 1993, ApJ, 410, 218

Wardle M., Salmeron R., 2012, MNRAS, 422, 2737

Xu W., Kunz M. W., 2021, MNRAS, 502, 4911

## APPENDIX A: TEST OF ATHENA++
### IMPLEMENTATION

In this appendix, we present two tests of our numerical setup in Athena++ to verify that self gravity is implemented correctly. These tests are similar to those in Appendix A and B of Riols et al. (2017a), which were used to test self gravity in PLUTO.

We first examine the stability of the background hydrostatic equilibrium. We initialize a un-magnetized disc in hydrostatic equilibrium with uniform temperature and $Q = 1$. The setup (including vertical resolution) is identical to the initial condition of our hydrodynamic simulations, except that we do not introduce any initial perturbation and we use a much smaller domain in $x$ and $y$ so that the disc is stable against GI. We evolve this equilibrium state for $100 \Omega^{-1}$ for different box sizes $L_z$. The amplitude of fluctuations in the velocity at $100 \Omega^{-1}$, shown in Fig. A1, increases slightly for a taller box, but remains well below the sound speed in all cases. The $L_z = 6H_0$ result is similar to Fig. A2 of Riols et al. (2017a). Excluding the floor density from the gravitational force is necessary for a tall box, otherwise the velocity perturbation is much larger in the atmosphere.

<figure><img src="image_16.png" /><figcaption>**Figure A1.** Vertical velocity perturbation for different vertical box sizes after evolving the background disc equilibrium for $t = 100 \Omega^{-1}$ using Athena++.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure A2.** Growth rate of the most unstable axisymmetric mode with $k_x = 2\pi/5H_0$ for an isothermal disc with different $Q$. The blue curve is the analytic result from Riols et al. (2017a) and orange crosses are measured from Athena++ simulations.</figcaption></figure>

Our second test measures the linear growth rate of the fastest growing axisymmetric mode and compares it with analytic result from Riols et al. (2017a). For simplicity, we consider a disc with an isothermal equation of state. Fig. A2 shows the growth rate of modes with $k_x = 2\pi/5H_0$ for different $Q$; the growth rates measured from simulations agree very well with analytic predictions.

## APPENDIX B: DEPENDENCE ON THE EMF RECONSTRUCTION SCHEME

We study in this appendix the dependence of the dynamo behaviour on the reconstruction scheme of the EMF. In PLUTO and Athena++, staggered components of the magnetic field are evolved as surface-averaged quantities and are updated using Stokes' theorem. This requires constructing an averaged EMF along the face edges using a reconstruction or averaging technique from the face centres. Different techniques are employed, all of them having different levels of numerical diffusion. The UCT HLL method uses a 2D Riemann solver based on a four-state HLL flux function (see Del Zanna, Bucciantini & Londrillo 2003). It is generally the most diffusive technique. The UCT0 and UCT CONTACT employs the face-to-edge integration procedures proposed by Gardiner & Stone (2005), where EMF derivatives are averaged from neighbouring zones (UCT0) or selected according to the sign of the contact mode (UCT CONTACT). These methods are known to have reduced diffusion. Fig. B1 shows the evolution of magnetic energy computed in PLUTO without explicit diffusion using different EMF reconstruction schemes. In case of linear cooling, all methods seem to converge towards the same state with $E_m \lesssim 0.01$. However, in the case of $T^4$ cooling, we see an important gap between the UCT HLL (blue line) and UCT CONTACT (orange line). In particular, the UCT HLL scheme seems to produce a much stronger field, in agreement with the fact that the GI dynamo is boosted in the presence of a small amount of dissipation. UCT0 is not part of this test because the method appears to be numerically unstable for this problem. This is somewhat expected, given that this reconstruction scheme is the one which is the least dissipative, and which is known to produce spurious oscillations even in the linear regime (Gardiner & Stone 2005). The green line in Fig. B1 corresponds to a simulation using the UCT CONTACT method and starting from the strong dynamo state obtained with UCT HLL. We see that the magnetic energy

<figure><img src="image_18.png" /><figcaption>**Figure B1.** Time-evolution of magnetic energy for different EMF reconstruction schemes and cooling laws (top: linear cooling, bottom: $T^4$ cooling). In the bottom panel, the orange line corresponds to a simulation with UCT CONTACT initialized from a seed field, while the green curve has the same reconstruction scheme but starts using a state from the UCT HLL run.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure B2.** Same as the bottom panel of Fig. B1, but with $A_m = 100$. There is no longer a significant difference between UCT HLL and UCT CONTACT.</figcaption></figure>

decays rapidly to the state corresponding to the orange curve. This rules out the existence of a bi-stable dynamo state and indicates that the non-linear saturation of the ideal simulation is dependent on the reconstruction scheme (and in particular sensitive to the level of numerical diffusion induced by these scheme).

The result is expected to be much less sensitive to the reconstruction scheme when ambipolar diffusion is included, since in this case the dissipation of field at small scale is no longer dominated by numerical dissipation. Indeed, at $A_m = 100$, the magnetic energy
saturates to the same value for both UCT HLL and UCT CONTACT (Fig. B2).

APPENDIX C: DYNAMO WITH IONIZED LAYER

So far, all of our simulations have employed a simple prescription with uniform ion density and Elasasser number Am in the box. However it is expected that, when the star is formed (Class I discs), it radiates FUV and potentially ionizes the surface of the disc. We then test the dynamo behaviour when such an ionized layer is present. This allows us to check also the robustness of the dynamo regarding the physical environment present above the disc surface and near the box vertical boundary. To simplify the problem, we use the same model as in Simon et al. (2015), which captures the essential physics of disc ionization, namely, the effect of FUV in the disc corona. In this model, Am is constant in the mid-plane (we chose Am = 10) and increases abruptly up to a certain height $z_{io}$. This height corresponds to the altitude at which FUV can penetrate, i.e. when the horizontally averaged density integrated from the vertical boundary, $\Sigma_i(z)$, becomes larger than $\Sigma_{ic} = 0.01 \text{ g cm}^{-2}$. The profile of Am is then given by:

$$ Am = Am_0 + 3.3 \times 10^{12} x_i(z) \frac{\rho}{\rho_0} (R_0/1 \text{ au})^{-5/4}, \quad (C1) $$

where $Am_0 = 10$ is the constant mid-plane value, $R_0 = 50$ au is the radius of interest, and

$$ x_i(z) = 10^{-1} \exp \left[ - \left( \frac{\Sigma_i(z)}{\Sigma_{ic}} \right)^4 \right]. \quad (C2) $$

is the ionization fraction. Note that the ionization fraction considered in the corona is much higher than what one might expect in a real protoplanetary disc (typically $x_i \sim 10^{-5}$), but this value ensures that we are almost in the ideal regime with Am $\gg 10$. To convert numerical $\Sigma_i$ into g cm$^{-2}$, we have used a typical radial profile of Class 0–I discs with $\Sigma = 1500 (R/1 \text{ au})^{-1} \text{ g cm}^{-2}$. Using such a prescription, the transition between the disc and the ionized layer takes place at $z = 3H_0$. Above this height, Am reaches values $\sim 10^5$ to $\sim 10^6$.

The simulation with an ionized layer is initialized using a state taken from run PL-ZNF-Am10 and is run for a time $200 \Omega^{-1}$. Fig. C1 (top) shows the vertical profiles of the mean magnetic field components averaged in time over the saturated state, for the case with an ionized layer and the case without (run PL-ZNF-Am10). We see that the profiles are almost identical in both cases, suggesting that the large-scale GI dynamo is not altered by the presence of an FUV layer. We also plot in the bottom panels the vertical profiles of the rms magnetic fluctuations, which provides an idea of the vertical distribution of the turbulent field. Again, there is no major difference between the case with and without FUV layer, except that turbulent magnetic fluctuations are slightly stronger for $z > 3H_0$ in the case with an FUV layer.

<figure><img src="image_20.png" /><figcaption>Figure C1. Top: vertical profile of mean magnetic field components in the radial, azimuthal, and vertical directions. Bottom: root-mean-square magnetic fluctuations in each direction, where fluctuations are defined as $B - \bar{B}$. The blue/solid line is for the simulation with an ionized layer and Am = 10 in the mid-plane. The purple/dashed line is for the simulation with uniform Am = 10.</figcaption></figure>

This paper has been typeset from a TeX/Latex file prepared by the author.