Issue 
A&A
Volume 663, July 2022



Article Number  A138  
Number of page(s)  17  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202243219  
Published online  21 July 2022 
Gravitoturbulent dynamo in global simulations of gaseous disks
^{1}
Institut für Astronomie und Astrophysik, Universität Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany
email: william.bethune@unituebingen.de
^{2}
DAMTP, University of Cambridge, CMS, Wilberforce Road, Cambridge CB3 0WA, UK
Received:
28
January
2022
Accepted:
10
June
2022
Context. The turbulence driven by gravitational instabilities (GIs) can amplify magnetic fields in massive gaseous disks. This GI dynamo may appear in young circumstellar disks, whose weak ionization challenges other amplification routes, as well as in active galactic nuclei. Although regarded as a largescale dynamo, only local simulations have so far described its kinematic regime.
Aims. We study the GI dynamo in global magnetohydrodynamic (MHD) models of accretion disks, focusing on its kinematic phase.
Methods. We perform resistive MHD simulations with the PLUTO code for different radiative cooling times and electrical resistivities. A weak magnetic field seeds the dynamo, and we adopt meanfield and heuristic models to capture its essence.
Results. We recover the same induction process leading to magnetic field amplification as previously identified in local simulations. The dynamo is, however, global in nature, connecting distant annuli of the disk via a largescale dynamo mode of a fixed growth rate. This largescale amplification can be described by a meanfield model that does not rely on conventional αΩ effects. When varying the disk parameters we find an optimal resistivity that facilitates magnetic amplification, whose magnetic Reynolds number, ℛ_{m} ≲ 10, is substantially smaller than in local simulations. Unlike local simulations, we find an optimal cooling rate and the existence of global oscillating dynamo modes. The nonlinear saturation of the dynamo puts the disk in a strongly magnetized turbulent state on the margins of the effective range of GI. In our simulations, the accretion power eventually exceeds the threshold required by local thermal balance against cooling, leaving the longterm nonlinear outcome of the GI dynamo uncertain.
Key words: accretion, accretion disks / dynamo / gravitation / magnetohydrodynamics (MHD) / turbulence
© W. Béthune and H. Latter 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1. Introduction
Magnetic fields play a central role in the evolution of accretion disks. They are often considered as essential for mass accretion due to their ability to power jets and disk winds (Blandford & Payne 1982; Pelletier & Pudritz 1992) and to trigger turbulence via the magnetorotational instability (MRI; Balbus & Hawley 1991). Dynamically consequent magnetic fields may be supplied by infalling plasma during the formation of the disk (protostellar ones, for example; see Zhao et al. 2020) or by the overflow of a binary companion (Ju et al. 2016, 2017; Pjanka & Stone 2020), and their accumulation can even lead to magnetically dominated flows (e.g., Tchekhovskoy et al. 2011). Alternatively, magnetic fields may be amplified and sustained by a dynamo process operating inside the disk.
Dynamos associated with the MRI have received the most attention due to their ubiquity in astrophysical disks and their relative strength (Hawley & Balbus 1992; Hawley et al. 1996; Rincon et al. 2007; Lesur & Ogilvie 2008; Gressel & Pessah 2015; Walker et al. 2016; Riols et al. 2017b; Mamatsashvili et al. 2020). Alternative fluid dynamos exist^{1} and become especially important when the ionization fraction is too low to sustain the MRI (e.g., Kawasaki et al. 2021). Recently, turbulence caused by the gravitational instability (GI; Toomre 1964) of massive disks has offered a novel and compelling way to amplify magnetic fields via the socalled GI dynamo (Riols & Latter 2019), not least because it can outcompete the MRI in certain regimes (Riols & Latter 2018a).
The conditions for GI are expected to be met in young circumstellar disks (Durisen et al. 2007; Kratter & Lodato 2016) and active galactic nuclei (AGN; Shlosman & Begelman 1987; Goodman 2003). The instability leads to a selfsustained turbulent state if radiative cooling is sufficiently slow (Gammie 2001; Rice et al. 2003; Lin & Kratter 2016; Booth & Clarke 2019), at which point it facilitates angular momentum transport via spiral density and gravitational perturbations (Paczynski 1978; Lin & Pringle 1987; Balbus & Papaloizou 1999). Because the GI develops in the plane of the disk, much hydrodynamical work has excluded the outofplane dynamics outright by employing twodimensional models, or otherwise neglected its analysis (exceptions include Boley & Durisen 2006; Shi & Chiang 2014; Riols et al. 2017a; Béthune et al. 2021). However, understanding the threedimensional (3D) structure of the flow is absolutely crucial when assessing the dynamo question.
A number of global magnetohydrodynamic (MHD) simulations have described the susceptibility of magnetized disks to fragmentation (Fromang 2005; Forgan et al. 2017; Zhao et al. 2018; Wurster & Bate 2019) or the interplay of the MRI with the GI in disks with a preexisting strong magnetization (Fromang et al. 2004a,b; Deng et al. 2020). But to date, only the local (shearing box) simulations of Riols & Latter (2018a, 2019) have probed the dynamo properties of GI turbulence at low magnetization and in regimes excluding the MRI, thus witnessing the kinematic phase of the GI dynamo and thereby determining its underlying mechanism. However, a problem they faced was the emergence of largescale magnetic fields – relative to the largest turbulent eddies, and indeed the computational domain – which called into question the local approximation itself. It is likely that global (geometrydependent) effects are important and need to be accounted for.
In this paper, we combine the global and the kinematic through MHD simulations of gravitoturbulent disks, thereby presenting clean results on the global GI dynamo process. Unlike previous global studies, we start with magnetic fields that are sufficiently weak such that the MRI is initially ineffective, especially when Ohmic diffusion is added. We use weighted averages to separate the GI dynamo into radially local and global parts. The local part helps us make contact with shearing box results and meanfield dynamo models, whereas the global part is entirely novel and comprises one of the main achievements of the paper. Our focus is on the kinematic dynamo regime, but we present some results on its saturated phase that reinforce the idea that the saturation route ends in a highly magnetized state (plasma β of order unity).
The rest of the paper is subdivided into seven sections. We present our framework in Sect. 2, including the chosen disk model, numerical scheme, and data reduction techniques. Section 3 provides an overview of the disks in which the dynamo operates, including descriptions of the mean flow and magnetic field. In Sect. 4 we decompose the kinematic induction loop into its components and then reconstruct a local geometric interpretation of the dynamo mechanism. We then show how global behaviors stem from the radial diffusion of magnetic energy. In Sect. 5 we establish a reduced dynamo model based on the local mean fields and whose parameters are directly measured from simulations. Section 6 undertakes a survey of various disk parameters to assess the robustness of our previous results. We examine in Sect. 7 the saturation of the GI dynamo as observed in our simulations. Finally, we summarize our results and some of their possible extensions in Sect. 8.
2. Method
Our numerical setup is nearly identical to that of Béthune et al. (2021), the main novelty being the inclusion of a magnetic field. Below we recall the properties of the adopted disk model and of its implementation within the gridbased and shockcapturing code PLUTO. We also define conventions and notations for various quantities used throughout the paper.
2.1. Global disk model and units
We simulated a disk of selfgravitating fluid orbiting a more massive central object. We used both spherical coordinates (r,θ,φ) and cylindrical coordinates (R,φ,z) with the central object at the origin and the rotation axis of the disk as the polar (resp. vertical) direction. We considered only a finite radial extent [r_{in},r_{out}], excluding the central region. For simplicity, we assumed that this frame is inertial by neglecting the reaction of the central object to the gravity of the disk, which is a reasonable approximation for the relatively low disk masses considered (Heemskerk et al. 1992; Michael & Durisen 2010). We modeled radiative energy losses via a simple cooling function that brings the gas temperature toward zero (in the absence of heating) over a prescribed number of orbits. Finally, we assumed that the disk is electrically conducting such that fluid motions can induce and react to magnetic fields according to a resistive MHD description.
Throughout this paper, we take the mass m_{⋆} of the central object as the constant mass unit and denote M ≡ m_{disk}/m_{⋆} the initial mass of the disk. We use the inner radius r_{in} of the simulation domain as the constant length unit. After absorbing the gravitational constant G into the masses, denotes the Keplerian frequency at radius R, and we take Ω_{in} ≡ Ω_{⋆}(r_{in}) as the inverse time unit. The inner Keplerian period is denoted t_{in} ≡ 2π/Ω_{in}. Magnetic fields are measured by their corresponding Alfvén velocities when the gas density ρ = 1 in these units.
2.2. Governing equations
The equations governing the evolution of the gas density ρ, velocity V, pressure P, and magnetic field B are:
We used a constant adiabatic exponent γ = 5/3 to facilitate comparison with Riols & Latter (2019), Deng et al. (2020), and Béthune et al. (2021). The gravitational potential Φ is the sum of a constant central potential −m_{⋆}/r and a timedependent contribution from the disk satisfying Poisson’s equation ΔΦ_{disk} = 4πρ. The cooling function in Eq. (3) takes the form Λ ≃ −(Ω_{⋆}/τ)P, parametrized by a constant dimensionless cooling time^{2}τ. The magnetic field influences the gas’s momentum via the Lorentz force in Eq. (2) and its internal energy via Joule heating in Eq. (3), where J ≡ ∇ × B is the electric current density. For later reference, we denote by ℰ ≡ −V × B the ideal electromotive force (EMF) in Eq. (4). We omitted explicit viscosity in Eqs. (2) and (3) and relied mainly on shocks to thermalize kinetic energy.
We modeled the finite electric conductivity of the gas with an Ohmic resistivity η in Eqs. (3) and (4), which we parametrized by a magnetic Reynolds number ℛ_{m} = Ωh^{2}/η. Since the characteristic thickness h of the disk can vary with its thermal stratification and selfgravity, there are multiple possible definitions for ℛ_{m}. For simplicity, we used the initial profile of the pressure scale height, h_{init}/R ≡ (c_{s}/V_{φ})_{init}, given the isothermal sound speed . The resistivity is therefore constant in time according to the prescribed value of . The actual ratio Ωh^{2}/η may vary with radius and time, but we checked that the characteristic scale height h_{ρ} (see Sect. 2.5 below) remains within ten per cent of the initial hydrostatic scale height h_{init} on average, and they differ by at most 30 per cent locally.
2.3. Numerical integration scheme
We used the finite volume code PLUTO version 4.3 (Mignone et al. 2007) to integrate Eqs. (1)–(4) in conservative form on a static spherical grid. The computational domain (r/r_{in},θ,φ) ∈ [1,32] × [π/2±0.35] × [0,2π] was meshed with 518 × 96 × 512 grid cells, with a logarithmic sampling in radius and a finer meridional resolution near the disk midplane (θ ∈ π/2 ± 0.1).
As Béthune et al. (2021), we used a secondorder RungeKutta time stepping with CFL coefficient 0.3, and a piecewise linear reconstruction (Mignone 2014) of the primitive variables at cell interfaces with the symmetric slope limiter of van Leer (1974). Interface fluxes were computed with the HLLd approximate Riemann solver of Miyoshi & Kusano (2005). At strong discontinuities, we switched to the more dissipative MINMOD slope limiter (Roe 1986) and to a twostate HLL Riemann solver (Harten et al. 1983). The magnetic field obeys the constrained transport formalism (Evans & Hawley 1988) so that ∇ ⋅ B = 0 is guaranteed from the start down to machine precision. We used the UCT_HLL method to reconstruct the EMF on cell edges (Del Zanna et al. 2003; Londrillo & del Zanna 2004). Unacceptable magnetic energy injection occurred at the domain boundaries when using other reconstruction schemes (see Appendix A). Ohmic resistivity was integrated explicitly via supertimestepping (Alexiades et al. 1996), although only the most resistive case (ℛ_{m} = 1) benefited from the associated acceleration.
2.4. Initial and boundary conditions
The initial conditions for the hydrodynamic variables are the same as those of Béthune et al. (2021). The gas initially rotates at , accounting for the enclosed disk mass . The initial surface density Σ ∝ R^{−2} was scaled by the chosen disk mass. The initial sound speed c_{s} ∝ R^{−1/2} gives a roughly constant aspect ratio h/R, and it was adjusted so that the Toomre number Q ≃ Ω_{⋆}c_{s}/πΣ ≈ 1.
We initialized the magnetic field with a net toroidal component confined inside 1.1r_{in} < r < 0.9r_{out} and θ−π/2 < 0.1. This magnetic field is present from the start of our simulations, before the onset of GI turbulence. Compared to the expected dynamo growth times, the initial transition to a quasisteady turbulent state is short and may be neglected. Given how weak the initial magnetic field is, the MRI would operate appreciably for wavenumbers kh_{init} ∼ 10^{6} in ideal MHD (Balbus & Hawley 1992; Ogilvie & Pringle 1996). Such scales are much smaller than our finest grid increment, and numerical diffusion would easily overcome the MRI at our grid scale. The addition of a finite Ohmic resistivity guarantees that the MRI cannot operate before reaching an average B^{2}/P ≳ 10^{−2}.
The boundary conditions for the hydrodynamic variables ρ, V, and P are the same as those of Béthune et al. (2021). We set the tangential components of the magnetic field to zero in the ghost cells of both the radial (r) and meridional (θ) boundaries. When a grid cell adjacent to a boundary supports a nonzero magnetic flux, the component normal to the boundary is automatically determined by ∇ ⋅ B = 0 in the ghost cells.
2.5. Data reduction
Throughout this paper, the word “mean” designates the result of some spatial integration, and fluctuations are defined relative to this average. It is conceptually different from separating the flow into a turbulent component over a quiescent background, as is customary in classic meanfield theory. This distinction should be kept in mind later when we discuss mean flows that originate in turbulent motions, although the average of most quantities does reflect their location and timeindependent properties.
Our simulation diagnostics rely on several averaging procedures to extract clear signals out of turbulence. We use angled brackets ⟨•⟩_{xi} to denote spatial averages over the dimensions x_{i}, omitting the subscripts when unnecessary. We use the subscript ρ to denote a densityweighted meridional and azimuthal average. Since we focus on thin disks, the meridional (θ) and vertical (z) denominations may be used interchangeably.
We used radial averages to maximize the statistical quality of our diagnostics and to extract vertical structures common at all radii. They covered the interval r/r_{in} ∈ [2,8] and were evaluated with a d log r measure corresponding to our grid spacing. To cancel the global radial gradients in the disk, we rescaled some quantities prior to their averaging. We used the densityweighted angular velocity Ω_{ρ} ≡ ⟨V_{φ}/R⟩_{ρ} and isothermal sound speed to define radial profiles used as averaging weights: a typical thickness h_{ρ} ≡ c_{ρ}/Ω_{ρ}, a typical velocity , and a typical magnetic field strength . We then obtained the mean meridional profiles of the reduced velocity v, mass flux q, magnetic field b, electric current density j, and ideal EMF ϵ as:
Each dimensionless ratio on the righthand side effectively fluctuates about the average on the lefthand side by a factor of a few at all radii and azimuths. To produce converged meridional profiles, we further averaged them over time. The results are smooth functions of latitude (θ) representing global and persistent background structures embedded in turbulent fluctuations.
3. Overview of the kinematic dynamo phase
In this section, we describe the broad properties of the turbulent disk and dynamo in a simulation chosen as reference; other disk parameters will be considered in Sect. 6. With a disk mass M = 1/3, a dimensionless cooling time τ = 10, and a magnetic Reynolds number ℛ_{m} = 10, we label this run M3T10R10. It has the same disk mass and cooling time as the reference simulation presented by Béthune et al. (2021), and its moderate magnetic Reynolds number places it in the kinematic dynamo regime identified by Riols & Latter (2019, see their Fig. 2).
3.1. Global, secular growth of magnetic energy
Turbulence takes about 200t_{in} to settle over the whole domain. From this point on, velocity fluctuations evolve on local orbital times and on lengthscales comparable to the local disk thickness, as illustrated in Fig. 1. Meanwhile, Fig. 2 shows the dimensionless magnetization ratio B^{2}/P exhibiting a largescale envelope that grows by 12 orders of magnitude over 900t_{in}.
Fig. 1. Flattened distribution of the gas surface density R^{2}Σ at time 400t_{in} over the whole domain of run M3T10R10; colors range linearly from zero to three times the median of R^{2}Σ. 
Noting that and that the disk supports sonic turbulence, the Alfvén velocity remains negligible compared to both the velocity fluctuations and the sound speed as long as B^{2}/P ≪ 1. In this regime, both the Lorentz force and Joule heating are unimportant, and hence the magnetic induction Eq. (4) reduces to a linear initialvalue problem for B in a given velocity field V. We focus here on this kinematic dynamo phase and postpone the dynamo saturation (B^{2}/P ≳ 1) to Sect. 7.
As long as B is dynamically unimportant, its amplitude should vary exponentially in time. We measured the evolution of during one outer orbit ≈181t_{in} before reaching B^{2}/P = 10^{−2} at any given radius and observed exponential growth. We plot in Fig. 3 the radial profile of the growth rate ω estimated at each radius by fitting the slope of as a function of time.
The solid blue curve in Fig. 3 indicates that the growth rate ω ≈ 1.6 × 10^{−3}Ω_{in} is roughly constant up to r/r_{in} ≲ 12. Beyond r/r_{in} ≳ 12 our measurements of ω become inaccurate because the chosen time interval covers only a few orbital periods, and the dynamo has not settled into a steady growth yet at these radii. If instead we measure the growth rate relative to the local orbital frequency Ω_{ρ}, we obtain a steady increase with radius up to ω/Ω_{ρ} ∼ 1, as shown by the dashed red curve.
Figure 3 reveals that magnetic amplification proceeds with a single growth rate irrespective of the local dynamical time. While significant growth takes hundreds of orbits in the inner parts, it happens on nearly orbital times in the outer parts. Thus, the dynamo mode causally connects distant radii, and its growth is likely limited by its outer extent (see Sect. 4.2). Figure 3 demonstrates a global property of the GI dynamo that shearing box simulations cannot capture. Because the dynamo mode in Fig. 2 exhibits lengthscales much longer than the turbulent driving scale, we also conclude that it is a largescale meanfield dynamo – whose type we determine in Sects. 4 and 5.
Fig. 2. Evolution in time (horizontal axis) of the radial profile (vertical axis) of the disk magnetization B^{2}/P in run M3T10R10. 
Fig. 3. Growth rate of the magnetic field amplitude as a function of radius in run M3T10R10, normalized either by the (fixed) inner Keplerian frequency Ω_{in} (solid blue) or by the (radially varying) densityweighted orbital frequency Ω_{ρ} (dashed red). Only the inner parts r/r_{in} ≲ 12 have reached a steady exponential growth phase. 
3.2. Hydrodynamic (GI) turbulence
In this subsection, we describe the smallscale features of the hydrodynamical turbulence that are important for the morphology and growth of the dynamo field. Figure 4 shows the relative surface density fluctuations
Fig. 4. Equatorial slice in run M3T10R10 at time 400t_{in}, zoomed in on the radial interval r/r_{in} ∈ [2,16]. The color map shows the relative density fluctuations as defined by Eq. (6) while the arrows show the orientation of the magnetic field sampled in the midplane, colored in blue when inward (B_{r} < 0) or green when outward (B_{r} > 0). 
at time 400t_{in} in the equatorial plane of run M3T10R10. This is similar to Fig. 1 after unfolding the azimuthal direction. The bright stripes where the density is maximal correspond to spiral wakes. Because of our choice of initial conditions, these wakes share a similar pitch angle (d log R/dφ) at all radii, thus generating logarithmic spiral arms globally. Due to the relatively short cooling time τ = 10, the GI sustains vigorous turbulence and sharp density contrasts (e.g., Cossins et al. 2009).
To reveal the 3D nature of the flow, we sampled meridional slices of the disk at different angles φ. The left panel of Fig. 5 shows the density and poloidal velocity components in the φ = 3π/8 plane and r/r_{in} ∈ [4,8] interval. The velocity field roughly displays a mirror symmetry about the midplane: the gas simultaneously rises or falls on both sides at any given radius. This symmetry hints at a connection between the midplane disk dynamics and the vertical circulations developing in its corona. The density maxima located at r/r_{in} ≈ 4.0, 5.6, and 6.9 correspond to spiral wakes that are surrounded by varied flow patterns. At r/r_{in} ≈ 4 the velocity converges toward the wake in the midplane and splashes vertically away from it. At r/r_{in} ≈ 5.6 we see the opposite trend: the gas falls vertically toward the wake and evacuates it in the midplane. At r/r_{in} ≈ 6.9 the gas is lifted vertically as it moves radially past the wake.
Fig. 5. Meridional slice at time 400t_{in} and angle φ = 3π/8 in run M3T10R10. Left panel: density distribution flattened by a factor of (R/r_{in})^{3} (color scale) and orientation of the poloidal velocity field, green (resp., blue) being oriented upward (resp., downward). Right panel: toroidal magnetic field relative to the local mean (color scale) and integral curves of the poloidal magnetic field. The density is maximal at r/r_{in} ≈ 4.0, 5.6, and 6.9. 
Distinct poloidal rolls appear at (r/r_{in}, π/2−θ) ≈ (4.5,0.05), (6.4,−0.1), and (6.9,0.05), yet we had difficulty identifying a consistent quadrupolar pattern encompassing wakes, as witnessed in the local simulations of Riols & Latter (2019) and Riols et al. (2021). Instead, such circulations correlate with the gas density in only an average sense, as shown by Béthune et al. (2021, Fig. 14). The discrepancy with local simulations may issue from our different vertical stratifications (see Sect. 3.4.1) or from the symmetry constraints (periodicity and zero net momentum) enforced in shearing boxes. The variety of poloidal flow structures may also reflect different evolutionary stages of the wakes. Béthune et al. (2021, Fig. 12) showed that these wakes are transient and survive only for a fraction of local orbital time. One possible way to envision their destruction is to simply reverse the orientation of the surrounding flow.
3.3. Smallscale magnetic field structure
The velocity fluctuations responsible for producing spiral density maxima also distort the magnetic field on scales much smaller than the global envelope described in Sect. 3.1. The arrows in Fig. 4 represent the orientation of the magnetic field sampled in the equatorial plane. Most arrows point to the right (increasing φ), meaning that the toroidal component B_{φ} > 0 has a prefered sign in the disk midplane. The radial component B_{r} is typically negative inside spiral wakes and positive in the lowdensity region between wakes, in agreement with local simulations (Riols & Latter 2019). While inplane motions cannot lead to dynamo amplification on their own, they do organize the magnetic field along the spiral density arms.
The right panel of Fig. 5 shows the structure of the magnetic field in the same meridional plane as the left panel. The toroidal component B_{φ} is maximal inside the wakes. This follows from the compressive motions that simultaneously accumulate mass and magnetic flux. Omitting resistivity, one might expect the same enhancement in gas density and in toroidal magnetic field for tightly wound spirals. As both increase by a factor of ≈5, we deduce that the accumulation of magnetic flux into spiral wakes proceeds on short timescales compared to the typical resistive diffusion time ∼h^{2}/η over the lengthscale of the spirals.
The poloidal magnetic field lines drawn in Fig. 5 are mostly horizontal in the disk midplane and vertical in its corona. The vertical component is predominantly oriented toward the midplane on both sides, such that the net magnetic flux passing vertically through the disk remains approximately zero at all radii. There is therefore no net vertical magnetic flux generated by turbulence (e.g., by splitting B_{z} into positive and negative patches), nor injection from the radial boundaries of the computational domain. As a consequence, there are no magnetic field lines connecting both vertical boundaries that could carry angular momentum vertically away or drive a magnetized wind.
3.4. Mean vertical structure
Because verticality – with respect to the disk plane – is crucial to the GI dynamo as pictured by Riols & Latter (2019), we take a moment to characterize the mean vertical stratification of the disk in the quasisteady turbulent state of our reference run.
3.4.1. Thermal stratification
Anticipating the role of vertical buoyancy in the GI dynamo, and to facilitate comparison with other work, we first examine the thermal stratification of the disk. Figure 6 shows that the gas temperature T ≡ P/ρ increases by a factor of ≈40 with height in run M3T10R10. Unlike the shearing box simulations of Shi & Chiang (2014) and Riols & Latter (2019), our disks are embedded in a warm corona. Béthune et al. (2021) attributed this feature to their choice of boundary conditions that hinder rotational support and indirectly promote pressure (thermal) support against the gravity of the central object.
Fig. 6. Thermal stratification of the disk. The solid blue line (left axis) shows the gas temperature relative to its midplane value. The dashed red line (right axis) shows the squared buoyancy frequency (Eq. (7)) relative to the local squared orbital frequency. Both profiles were averaged azimuthally, over r/r_{in} ∈ [2,8] and over time from 200t_{in} to 400t_{in}. 
To quantify the buoyant response of the gas to small vertical displacements, we define a squared buoyancy frequency
where s ∝ P/ρ^{γ} measures the gas specific entropy, and where the gravitational potential Φ includes the contribution of the disk. This quantity is also plotted in Fig. 6 after normalizing it by the squared orbital frequency and averaging. We find that everywhere, as expected from the combination of a decreasing density and increasing temperature with height relative to the midplane. The disk is therefore convectively stable according to Schwarzschild’s criterion. In fact, fitting a polytropic relation P ∝ ρ^{Γ} onto the mean vertical profiles gives Γ ≈ 0.7 in run M3T10R10. This stratification is strongly subadiabatic, and even outside the range probed by Riols & Latter (2018b). Based on their results, we would expect the thermal stratification to facilitate the formation of rolls near spiral wakes.
3.4.2. Mean flows
The disk also supports mean gaseous flows varying with height relative to the midplane, as we show in Fig. 7. We distinguish the mean velocity v of the gas (upper panel), which governs the transport of magnetic flux according to Eq. (4), from the mean mass flux q (lower panel) as defined by Eq. (5). One may vanish while the other does not, such that mass and magnetic flux may be transported in different proportions – even in ideal MHD.
Fig. 7. Meridional profiles of the mean flow as defined by Eq. (5), averaged from 200t_{in} to 400t_{in} in run M3T10R10, and focused on the midplane region θ ∈ π/2 ± 0.1. Upper panel: spherical components of the reduced gas velocity. Lower panel: components of the reduced mass flux. The toroidal components are measured on the right vertical axes. 
Looking first at the radial components (solid red), the mean velocity v_{r} (upper panel) is roughly constant and overall negative, suggesting a slow inward advection of (toroidal) magnetic flux. The mean radial mass flux q_{r} (lower panel), on the other hand, reveals mass accretion near z/h_{ρ}≈0.7 and decretion in the midplane. We caution that the disk is not steady on viscous timescales because of our choice of initial conditions, so there is no conflict between mass decretion and outward transport of angular momentum by turbulent stresses.
Considering the toroidal components (dashed blue, right axes), both the gas velocity and mass flux are peaked in the midplane. The orbital velocity v_{φ} (upper panel) decreases with height, and its vertical shear rate dV_{φ}/ dz ≈ −3Ω_{ρ} even exceeds the radial (Keplerian) one. This vertical shear is necessary to satisfy momentum balance in conjunction with the steep thermal stratification (see Fig. 6). We do not expect it to drive a vigorous vertical shear instability because of the slow cooling timescale (τ > 1) and the strongly stabilizing entropy stratification (Urpin & Brandenburg 1998; Nelson et al. 2013).
Finally, the mean meridional mass flux q_{θ} (dotted green, lower panel) indicates that gas is accumulating in the midplane: the disk compresses vertically, possibly in order to keep Q ∼ 1 as mass is transported radially. The meridional mass flux does eventually vanish at large z, assuring us that there are no appreciable mass and thermal energy exchanges between the computational domain and its boundaries. In contrast, the amplitude of the mean meridional velocity v_{θ} (upper panel) increases with height and reaches a plateau beyond the range shown in Fig. 7. This is caused by correlations of the density and velocity fluctuations. On average, the highdensity wakes launch slow vertical updrafts from a small area of the disk; the lofted mass then falls to the midplane through dilute but fast vertical downdrafts. These flows manifest in a nonzero mean vertical velocity upon azimuthal averaging, but with no associated net mass flux.
3.4.3. Largescale magnetic field
Figure 2 shows the radially global structure of the largescale dynamo mode. In this section, we use the weighted radial average (5) to extract its characteristic vertical structure. The resulting profiles are plotted in Fig. 8.
Fig. 8. Meridional profiles of the three components of the reduced magnetic field as defined by Eq. (5) and averaged from 200t_{in} to 400t_{in} in run M3T10R10; the toroidal component b_{φ} is measured on the right axis. 
The radial component b_{r} (solid red) is negative inside z/h_{ρ}≳1.7 and positive further away from the midplane. It reaches at most seven per cent in amplitude, with two minima at z/h_{ρ}≈0.7. The meridional component b_{θ} (dotted green) is positive in the upper half of the disk (θ − π/2 < 0) and negative in the lower half, as observed in the right panel of Fig. 5. It is clearly the weakest component, with an amplitude at most ten times smaller than b_{r}. Finally, the toroidal component b_{φ} (dashed blue, right axis) is positive within 4h_{ρ} from the midplane and only slightly negative beyond, with a maximum amplitude of about 2.3 at the midplane. We conclude that the mean magnetic field is mostly toroidal and confined inside the disk.
According to Eq. (4), the net toroidal flux ∫B_{φ} dS crossing any meridional plane (constant φ) can only change via boundary effects; any dynamo operating inside the domain would generate positive and negative flux in equal proportions. We do observe the growth of a net toroidal flux at all radii, indicating a nonzero circulation of the EMF along the boundary contour. Because we enforce (B_{φ},V_{φ},V_{θ}) = 0 in the ghost cells of the meridional boundaries, the ideal EMF ℰ_{r} should approximately vanish there. The growth of a net positive toroidal flux must then be caused by the Ohmic diffusion of negative B_{φ} < 0 contributions through the boundaries, and by any residual ideal ℰ_{r} since it is not exactly controlled on the boundaries.
4. Analysis of the largescale kinematic dynamo
In this section, we examine the growth of a largescale magnetic field during the B^{2}/P ≪ 1 phase of our reference run M3T10R10. As mentioned in Sect. 3.1, the induction Eq. (4) then reduces to a linear initialvalue problem for the magnetic field B. In Sect. 4.1 we use this radially averaged equation to retrace the feedback loop leading to magnetic amplification. In Sect. 4.2 we insert the (radially) local dynamo process into an idealized global disk model to understand the observed largescale growth behaviors.
4.1. The GI dynamo as a radially local process
We first focus on local aspects of the GI dynamo to make contact with the shearing box simulations of Riols & Latter (2019). We decompose the induction Eq. (4) into elementary terms and then track how each component of the magnetic field acts on the other components to produce an unstable induction cycle.
4.1.1. Mean magnetic field induction
The curl of the ideal EMF appearing in Eq. (4) is equivalent to the divergence of a tensor V ⊗ B − B ⊗ V, leading to the following conservative form in cartesian coordinates:
The component A_{ji} represents the transport of B_{i} along the velocity V_{j}, including advection and compression. The terms S_{ji} represent the stretching of components B_{j} into B_{i} by velocity gradients. The diagonal terms A_{ii} and S_{ii} cancel each other as B_{i} is only sensitive to transverse velocity components V_{j ≠ i}. Using spherical coordinates, we compute the curl of the ideal EMF with the above notation but including the appropriate geometrical factors inside the derivatives. For example, A_{rφ} ≡ −(1/r)∂_{r}[rV_{r}B_{φ}] denotes the radial transport of B_{φ}.
Next, we decompose any flow variable X as an axisymmetric part ⟨X⟩_{φ} plus a fluctuating part X′. After azimuthal averaging, the ideal EMF becomes
We denote by and the transport and stretching terms appearing in the azimuthally averaged Eq. (8) – affecting the mean field ⟨B⟩_{φ} – and arising from the rightmost term of Eq. (9), that is, the correlation of the fluctuations V′ and B′.
To isolate the largescale component of the magnetic field, we also average Eq. (8) in the radial dimension after normalizing it by at every radius. As for the other reduced variables, defined in Eq. (5), this procedure effectively flattens the global gradients in the disk and allows different radii to contribute at the same level to the average.
4.1.2. Mean toroidal magnetic flux
Because the toroidal component B_{φ} is dominant, we take it as the starting point of a feedback loop. The different contributions to the induction of a mean B_{φ} are represented in Fig. 9. The toroidal component is controlled by the largest number of terms, to which we can nevertheless associate simple interpretations.
Fig. 9. Contributions to the induction of a mean toroidal field after normalizing Eq. (8) by and averaging azimuthally, radially over r/r_{in} ∈ [2,8], and over time from 200t_{in} to 400t_{in} in run M3T10R10. 
The A_{θφ} term (dashed, dark green) represents the accumulation of B_{φ} in the midplane due to the meridional velocity. It is dominated by the action of the converging mean flow discussed in Sect. 3.4.2 (v_{θ} in Fig. 7). The Ohmic contribution (dotted black) opposes A_{θφ}, working to reduce B_{φ} in the midplane by spreading it vertically. The symmetry of these two terms translates into an approximate advectiondiffusion balance in the vertical direction.
The other term acting to induce a positive B_{φ} near the midplane is S_{rφ} ≡ (1/r)∂_{r}[rV_{φ}B_{r}] (solid red), corresponding to the stretching of B_{r} into B_{φ} by the differential rotation of the disk. Its fluctuating part in fact reduces to zero, so the growth of a net B_{φ} > 0 essentially comes from the radial shear of a net B_{r} < 0. This is usually refered to as the Ω effect in classic meanfield dynamo theory, and it is clearly operating in our simulations, in agreement with Riols & Latter (2018a, 2019).
Finally, the S_{θφ} term (dotdashed, light green) is approximately mirror symmetric with S_{rφ}, albeit with a smaller amplitude. It represents the stretching of a vertical field (B_{z}) into a toroidal one (B_{φ}) by the gradient ∂_{z}[B_{z}V_{φ}]. Its fluctuation part is comparatively small, so S_{θφ} also acts on the mean fields as an Ω effect. However, it works against the induction of a positive B_{φ} and therefore against the GI dynamo. This term was absent from local simulations because they do not capture the global gradients responsible for a vertical shear in the orbital motion.
4.1.3. Mean radial magnetic flux
We have shown that the mean B_{φ} > 0 appears to be generated primarily by the shear of a mean B_{r} < 0; we now look for the origin of this mean radial field. To do so, we plot the various terms of the averaged induction equation for ⟨B_{r}⟩_{φ} in Fig. 10.
The A_{θr} term (solid, dark green) represents the vertical transport of B_{r} by V_{θ}. Its peak in the midplane is essentially due to the fluctuating part (thin gray). Being opposed to the Ohmic term (dotted black), the turbulent velocity fluctuations appear to play an antidiffusive role on B_{r} in the vertical direction, with an approximate balance similar to that found for B_{φ} in Fig. 9.
The only term inducing ⟨B_{r}⟩_{φ} < 0 on both sides of the midplane, and therefore enhancing the structure shown in Fig. 8, is S_{θr} ≡ (1/rsin(θ))∂_{θ}[sin(θ)V_{r}B_{θ}] (dotdashed, light green). This term corresponds to the stretching of a vertical field (B_{z}) into a radial one (B_{R}) by a vertical gradient ∂_{z}[B_{z}V_{R}]. It is dominated by its fluctuating part , and hence the induction of a mean ⟨B_{r}⟩_{φ} cannot be directly traced back to another meanfield component. In Sect. 5 we nevertheless formulate a closure for the dynamo loop that relies solely on the mean fields.
4.1.4. Fluctuations of the meridional magnetic field
Although the mean meridional component ⟨B_{θ}⟩_{φ} plays no significant role in the induction of ⟨B_{r}⟩_{φ}, the key term does rely on nonaxisymmetric fluctuations . Because ⟨B_{θ}⟩ is so small, essentially measures the energy in these fluctuations. Our final step is to connect this magnetic energy back to the mean fields. To this purpose, we multiply each component A_{jθ} and S_{jθ} of the induction Eq. (8) by 2B_{θ} to obtain (twice) the rate of change of the magnetic energy, that is, . We normalize each term by the local , average them as previously, and plot them in Fig. 11.
Fig. 11. Contributions to the energy stored in after normalizing by at every radius and averaging azimuthally, radially over r/r_{in} ∈ [2,8], and over time from 200t_{in} to 400t_{in} in run M3T10R10. 
Only the S_{φθ} term (dashed blue) is positive throughout and therefore injects energy into the fluctuations of B_{θ}. It corresponds to the generation of a vertical field (B_{z}) from a toroidal one (B_{φ}) by the azimuthal variations of the vertical velocity (∂_{φ}V_{z}). This term closes the loop by generating fluctuations from the mean B_{φ}. These fluctuations can then induce a mean radial field, which then shears back into a mean azimuthal field.
Considering the remaining terms, A_{rθ} (solid orange) corresponds to the radial transport (advection and compression) of and should eventually vanish upon averaging. The S_{rθ} term (dotdashed red) is negative inside the disk, meaning that the stretching of a radial field (B_{R}) into a vertical one (B_{z}) by radial gradients ∂_{R}[B_{R}V_{z}] happens at the expense of magnetic energy.
The energy exchanges of the other components and are also confined inside z/h_{ρ}≲2, confirming that the growth of the disk magnetization follows from a turbulent dynamo inside the disk, with no significant energy input from the domain boundaries. Magnetic energy is mainly injected into but also slightly into , while Ohmic resistivity dissipates most of it into heat and brings the disk close to marginal dynamo stability.
4.1.5. Illustration of the local dynamo process
We have recovered the ingredients of the GI dynamo identified by Riols & Latter (2019) after azimuthal and radial averaging of our global simulation results. We now provide a more geometric interpretation of its successive steps in Fig. 12. As noted in Sect. 3.2, we could not easily identify such regular flow patterns in our simulations, so Fig. 12 should be taken as idealized. In particular, we omit the mean vertical shear ∂_{z}V_{φ} and its ability to stretch B_{z} into B_{φ} (the S_{θφ} term in Fig. 9) in this already cramped representation. Finally, the role of the mean vertical advection of magnetic field is also neglected, even though later in Sect. 5 we suggest that it is potentially important.
Fig. 12. Dynamo mechanism as seen in the frame of a spiral wake; see the text in Sect. 4.1.5 for a stepbystep description. 
In the figure, the thick gray line represents a spiral wake, inclined with respect to the azimuthal (φ) direction. The green arrows at the bottom represent the mean shear flow, while the blue arrows and helices describe nonaxisymmetric velocity fluctuations. The segmented lines represent a magnetic flux tube at four consecutive times. Initially toroidal (yellow), it is first pinched vertically away from the midplane (orange), then twisted into a radial component by the helices (light red), stretched azimuthally by the mean shear flow (dark red), and finally folded to form a loop with two flux tubes bunched in the midplane, reinforcing the azimuthal field we started with.
4.2. Global aspects of the dynamo
By using radially averaged quantities, we lost information concerning the radially global properties of the dynamo, as manifested in Fig. 2. We now examine how magnetic amplification proceeds across the radial span of the disk. We propose a heuristic model to describe global dynamo modes in Sect. 4.2.1 and apply it to estimate turbulent parameters in Sect. 4.2.2.
4.2.1. Heuristic model of a global dynamo
Whatever the equation(s) governing the evolution of , the fact that with a constant ω in Fig. 3 means that we are observing a global eigenmode. It is instructive to build a physically motivated model for the evolution of that could connect the local and global physics of the problem.
Béthune et al. (2021) argued that the properties of GI turbulence were essentially local within this simulation setup. We also know from the simulations of Riols & Latter (2019) and from our analysis of Sect. 4.1 that the GI dynamo can be explained using only radially local quantities. It therefore seems reasonable to assume that magnetic amplification is driven by a local source term operating on local dynamical times, where the dimensionless factor ϖ encapsulates the details of the local dynamo efficiency (i.e., its growth rate).
The second major actor is radial magnetic diffusion, as provided by Ohmic resistivity and turbulence, which provides the “connective tissue” between different radii and thus permits the development of coherent global modes. Our setup is designed such that dimensionless numbers constructed from local quantities – the aspect ratio h/R, cooling timescale τ, Ohmic Reynolds number ℛ_{m} – are independent of radius. We may as well expect the turbulent Reynolds number to be radially constant. The total magnetic diffusivity would then take the form ζΩh^{2}, where ζ denotes an effective magnetic Reynolds number. We therefore propose a simple reactiondiffusion equation to describe the evolution of the magnetic field amplitude:
consisting of local growth plus radial diffusion.
In the next paragraphs, we solve Eq. (10) numerically. But assuming that , and after a change of variable, we can manipulate it into a Schrödinger form (not shown) from which we extract a mode’s turning point R_{tp} – delimiting the region R < R_{tp} where a dynamo mode is localized. An approximate equation for R_{tp} is ω ≃ ϖΩ(R_{tp}), which reveals that the global growth rate of the dynamo mode ω is set by the local forcing rate at the mode’s outermost (and hence slowest) part^{3}.
4.2.2. Fitting global eigenmodes
We validated Eq. (10) a posteriori by comparing its eigenmodes to simulation data, namely the radial profile of at the time that the average disk magnetization reaches B^{2}/P = 10^{−2} anywhere. We computed eigenmodes of Eq. (10) for given pairs of (ϖ,ζ) by successively integrating it in time and rescaling the solution until convergence, imposing at both boundaries^{4} and using the initial simulation profiles of , Ω and h. Noting that the eigenmodes of Eq. (10) are invariant when changing ϖ and ζ in the same proportions, we included the measured growth rate ω/Ω_{in} as an additional optimization constraint. We then used the uncertainty on the measured global growth rate to evaluate uncertainties on the best fit parameters ϖ and ζ. We consistently obtained good fits when the disk supports a steady exponential magnetic growth, as illustrated in Fig. 13 for a subset of simulations.
Fig. 13. Normalized profiles of the magnetic field amplitude measured before dynamo saturation in four simulations featuring a steady exponential growth (solid lines; see legend), and the bestfit eigenmodes from Eq. (10) having the same global growth rate (dashed lines). 
In the reference simulation M3T10R10, the optimal local growth rate ϖΩ ≈ 7.3 × 10^{−2}Ω equals the measured global growth rate ω ≈ 1.6 × 10^{−3}Ω_{in} at a radius r/r_{in} ≈ 14, which we take to be the turning point of the mode. Regarding the effective magnetic Reynolds number ζ ≈ 2.7, it is nearly four times smaller than the prescribed Ohmic Reynolds number ℛ_{m} = 10. This difference is too large to be explained by variations of the disk scale height h (see Sect. 2.2), implying a contribution from turbulent diffusion. We come back to this idea when probing the ideal MHD limit in Sect. 6.1.3.
Although the best fit parameters ϖ and ζ are well constrained by our simulation data, they should be interpreted with caution since Eq. (10) was not derived from first principles. For example, our local growth rates ϖ ∼ 10^{−2} are smaller than those reported by Riols & Latter (2019) in the same parameter regime. If a direct comparison can be made with local simulations, then this difference could result from the new global effects, such as radial transport (ζ in Eq. (10)) and vertical shear (S_{θφ} in Fig. 9).
5. Meanfield model of the GI dynamo
The dynamo process summarized in Sect. 4.1.5 can be understood in terms of stretching, twisting, and folding magnetic field lines. This terminology provides a satisfying geometric interpretation within classic dynamo theory (Moffatt & Dormy 2019; Rincon 2019). It is idealized, however, as the generation of a net poloidal field out of a toroidal one depends on the statistical properties of turbulence (Parker 1955). As a compromise, one may posit a relationship between the turbulent EMF and the local mean fields (Krause & Rädler 2016; Brandenburg 2018; and example applications by Vishniac & Brandenburg 1997; von Rekowski et al. 2003; Fendt & Gaßmann 2018). This approach is justified by the relatively low ℛ_{m} in our simulations.
In this section, we establish such a meanfield closure for the turbulent EMF. The practical goal is to generate ⟨B_{r}⟩_{φ} directly from ⟨B_{φ}⟩. The fundamental goal is to better understand the action of GI turbulence on magnetic fields. We formulate a closure ansatz in Sect. 5.1, constrain its parameters from simulation data in Sect. 5.2, exhibit the resulting meanfield equations in Sect. 5.3, and identify their destabilizing constituents in Sect. 5.4.
5.1. Closure ansatz for the turbulent EMF
We seek a relation between the azimuthally averaged turbulent EMF, ⟨ℰ′⟩_{φ} ≡ ⟨−V′×B′⟩_{φ}, and the mean (axisymmetric) field ⟨B⟩_{φ}. The simplest closure is a proportionality ℰ′∝B, allowing magnetic amplification when coupled to a mean shear – the socalled αΩ dynamo. It has been suggested that the GI dynamo may belong to this family (Riols & Latter 2019; Riols et al. 2021), and we tested this hypothesis against our simulations.
We normalized ℰ′ by its characteristic scale and averaged it over (r,φ) to obtain mean meridional profiles of ϵ′ as defined by Eq. (5). This is similar to Sect. 4.1, where we removed any radial dependence of the induction equation and could thus extract clear profiles. The next step is to connect the reduced EMF ϵ′ to dimensionless vectors derived from the mean magnetic field. This relation must be linear during the kinematic induction stage, leading to the following ansatz:
restricted for simplicity to the reduced magnetic field b and electric current density j. The α′ matrix plays its conventional role of proportionality, whereas β′ can be interpreted as a turbulent resistivity matrix. For simplicity, we assumed these coefficients to be constant, independent of height (unlike, e.g., Gressel & Pessah 2015). Because ∂_{t}⟨B_{R}⟩_{φ} = ∂_{z}⟨ℰ_{φ}⟩, the most straightforward way to induce B_{R} out of B_{φ} is via a nonzero .
5.2. Closure coefficients from simulation data
We computed the profiles of ϵ′, b, and j in every simulation featuring a steady magnetic amplification (see Sect. 6). For each three component of ϵ′, we adjusted the six optimal parameters of Eq. (11) on the interval θ−π/2 ≤ 0.1 close to the midplane. We estimated these parameters in 20 independent simulation profiles evenly spaced from 200t_{in} to 400t_{in}, and then on the corresponding timeaveraged profiles. This procedure helped us identify which parameters are statistically significant: we only retained those incompatible with zero to better than three standard deviations. As a validation, we replaced the ideal EMF by the Ohmic EMF and recovered α′≈0 and , where I is the identity matrix.
Figure 14 shows the timeaveraged turbulent EMF in the reference run M3T10R10, along with the best fit profiles of Eq. (11). We could always obtain good fits for and when including the β′ matrix, but not with the α′ matrix alone. These matrices typically take the following form:
Fig. 14. Meridional profiles of the reduced EMF as defined by Eq. (5). The solid lines are simulation data averaged from 200t_{in} to 400t_{in} in run M3T10R10. The dashed lines are the best fits from Eq. (11) over this interval. The toroidal components are measured on the right axis. 
with , , and . The diagonal coefficient is always compatible with zero, precluding the simplest link from B_{φ} to B_{R}. In fact, the nonzero components of α′ all involve the meridional (θ) dimension. Both and connect to the mean meridional field ⟨B_{θ}⟩_{φ}, which does not appear explicitly in the dynamo process summarized in Sect. 4.1.5. The component is the least well constrained because the quality of the fit for varies across simulations – Fig. 14 shows a rather bad case.
In the turbulent resistivity matrix, we find , which indicates a balance between Ohmic diffusion and turbulent antidiffusion of B_{r} in the vertical direction, in agreement with Fig. 10. The diagonal term is always smaller than , implying an inefficient vertical transport of B_{φ} by turbulence. The offdiagonal terms satisfy and ; they rotate B_{R} into B_{φ} and vice versa, not unlike a (linearized) Hall effect. This process can destabilize axisymmetric perturbations depending on the orientation of the rotation (Balbus & Terquem 2001; Kunz 2008; Squire & Bhattacharjee 2016), which is stabilizing in our case.
5.3. Equations governing the mean fields
To lighten the notation, we use cylindrical coordinates and drop the angled brackets until Sect. 6. Let denote the scaledup turbulent resistivity matrix. Omitting and from Eq. (12) and injecting the remaining coefficients into the azimuthally averaged induction Eq. (4) leads to the following closed system:
In both equations, the first term on the righthand side is due to the mean vertical inflow (v_{θ} in Fig. 7), the second term is vertical diffusion, and the third term is the rotation of magnetic field caused by the offdiagonal turbulent resistivity. The second line of Eq. (14) consists of radial shear (Ω effect) and diffusion.
Approximating and , the radial shear term vanishes in Eq. (14). But we know from Fig. 9 that this term (S_{rφ}) is in fact compatible with a Keplerian shear rate (). This conflict comes from the radial averaging used to define j_{θ}, leading to a degeneracy with b_{φ}. The price of our averaging procedure is that the entire column , and hence we lose radial turbulent diffusion in Eq. (14). To make up for it, the best fit of is biased by an unphysical . We can therefore discard this only remaining α′ coefficient.
Neglecting the effective resistivity in Eq. (13) and the radial gradients of RB_{R} and RB_{φ} in Eq. (14), we finally obtain a radially local system of equations:
whose originality relative to classic meanfield dynamo models stems from the absence of an α effect, the apparently stabilizing pair of offdiagonal resistivity coefficients, and the presence of a compressive mean flow (Rincon 2019).
5.4. Meanfield dynamo route(s)
Multiplying Eq. (15) by B_{R} and Eq. (16) by B_{φ}, we find two possible energy source terms that can drive instability. The first corresponds to the Maxwell stress ∝B_{R}B_{φ} feeding on the mean orbital shear in Eq. (16). This is the classic Ω effect identified in Sect. 4.1.2. The second corresponds to the vertical inflow via ∂_{t}B^{2} = −(∂_{z}V_{z})B^{2}. This apparent mean flow – resulting from correlated velocity and density fluctuations (see Sect. 3.4.2) – formally permits instability even without a mean shear, and is a novelty of our reduced model.
When V_{z} = 0, the offdiagonal coefficients and bring about a pair of vertically propagating waves with a retrograde rotation of the magnetic field (B_{R},B_{φ}). This rotation is apparent in Fig. 12 as a toroidal field line is pinched upward and twisted from B_{φ} > 0 into B_{R} > 0. These terms operate similarly to (but are distinct from) the more common α effect that generates a poloidal field out of a toroidal one.
To explore the instability route associated with V_{z} ≠ 0, we focus on the first two terms on the righthand sides of Eqs. (15) and (16). We computed their eigenmodes as a function of z/h ∈ [0,2] using Chebyshev polynomials over 64 GaussLobatto collocation points. We sought symmetric solutions with respect to the midplane (∂_{z}B_{R} = ∂_{z}B_{φ} = 0 at z/h = 0), with zero net toroidal flux (∫B_{φ} dz = 0), no radial flux at the upper boundary (B_{R} = 0 at z/h = 2), and a constant compression rate ∂_{z}V_{z} ≤ 0. These conditions are meant to represent our simulations; they do allow energy exchanges across boundary at z/h = 2, where (V_{z},B_{φ},∂_{z}B_{R},∂_{z}B_{φ}) ≠ 0. We obtain unstable modes when V_{z} ≠ 0 and the compression rate −∂_{z}V_{z} exceeds the waves’ original frequency (see Appendix B), indicating that the present instability need not be sheardriven. The addition of a radial shear d log Ω/ d log R < 0 in Eq. (16) further destabilizes these modes if the shear rate is small enough, and hence the Ω effect as identified in Sect. 4.1.2 can simultaneously operate.
Finally, we injected the closure coefficients measured in run M3T10R10 into Eqs. (15) and (16) and found the unstable eigenmode drawn in Fig. 15. Its structure qualitatively matches the corresponding simulation profiles drawn in Fig. 8. Magnetic energy is mainly injected into B_{φ}, 80 per cent of which comes from the mean vertical inflow, 17 per cent from the Maxwell stress, and 99.7 per cent is dissipated by Ohmic resistivity. Its growth rate of 5.3 × 10^{−4}Ω makes it only marginally unstable, in line with our simulation diagnostics (see Sect. 4.1.4). In conclusion, the ansatz (Eq. (11)) captures the mean (axisymmetric) field dynamo observed in our simulations without relying on α nor Ω effects. The analysis hence complicates the model constructed in Sect. 4.1.5, which, though relatively transparent, is an idealization that downplays the importance of vertical inflows.
Fig. 15. Eigenmode of (15)–(16) under a prescribed converging flow , satisfying (∂_{z}B_{R},∂_{z}B_{φ}) = 0 at z/h = 0, B_{R} = 0 at z/h = 2, and ∫B_{φ} dz = 0. This solution includes a Keplerian radial shear rate d log Ω/ d log R = −3/2, an Ohmic resistivity η = 0.1Ωh^{2}, and the dynamo coefficients estimated in the reference simulation M3T10R10: , , , and . It is marginally unstable with a growth rate of 5.3 × 10^{−4}Ω. 
6. Varying the disk parameters
We ran simulations with different values of the disk mass M relative to the central object, different dimensionless cooling times τ, and different magnetic Reynolds numbers ℛ_{m}, as listed in Table 1. In this section, we examine how these parameters influence our findings in the kinematic regime of the GI dynamo.
One quantity that was useful in distinguishing different dynamo behaviors across our survey was the ratio of magnetic energy stored in the axisymmetric part of B relative to the total magnetic energy:
Smaller values of χ ∈ [0,1] mean that a larger fraction of magnetic energy is stored in nonaxisymmetric fluctuations, and hence this ratio measures how ordered the magnetic field is.
6.1. Ohmic resistivity
Riols & Latter (2019) reported that a largescale kinematic dynamo only occurs for moderate values of the magnetic Reynolds number (see their Fig. 2). At low Reynolds numbers ℛ_{m} ≲ 2 the dynamo was quenched by resistivity, whereas at large numbers ℛ_{m} ≳ 100 it turned into a more stochastic and smallscale process (Riols & Latter 2018a). We ran simulations with ℛ_{m} ranging from 1 to the “ideal” MHD limit (see Table 1). We recall that the Ohmic resistivity η is constant in time and varies with radius R according to the prescribed value of , compatible with the definition used by Riols & Latter (2019).
6.1.1. High resistivity and optimal Reynolds number
We ran two simulations with an increased Ohmic resistivity compared to the reference case M3T10R10, corresponding to magnetic Reynolds numbers in run M3T10R3 and ℛ_{m} = 1 in run M3T10R1. Both disks supported an exponentially growing magnetic field with global growth rates larger than in the reference run M3T10R10 (see ω/Ω_{in} in Table 1). The global growth rate is maximal in the ℛ_{m} ≈ 3.2 case, confirming the nonmonotonic dependence of the largescale dynamo on resistivity. In comparison, Riols & Latter (2019) found maximal growth rates for ℛ_{m} ≈ 20 regardless of the cooling time τ. This difference of optimal ℛ_{m} may point to different characteristic scales associated with the spiral wakes forming in global versus local simulations. Alternatively, it could result from the radial transport properties in the different flow geometries.
6.1.2. Low resistivity and oscillatory dynamo
We performed one simulation with a larger magnetic Reynolds number ℛ_{m} = 10^{3/2} ≈ 32, labeled M3T10R32. Its dynamo growth rate is smaller than in the reference case M3T10R10, following the trend found at high resistivity above. Unlike other simulations presented so far, Fig. 16 reveals that, during a time period of 1000t_{in}, the orientation of the mean magnetic field flips sign twice while its amplitude grows. This simulation thus demonstrates the existence of an oscillatory dynamo mode that can apparently dominate on large scales. Why they were not reported by Riols & Latter (2019) may be due in part to their initial and boundary conditions, while the oscillations reported by Deng et al. (2020, see their Fig. 5) were attributed to the MRI – outside the kinematic GI dynamo regime. Whether a purely growing mode would take over on longer timescales in run M3T10R32 is unfortunately out of our computational reach.
Fig. 16. Evolution in time (horizontal axis) of the radial profile (vertical axis) of the disk magnetization in the lowresistivity run M3T10R32. The mean magnetic field flips sign twice while growing in amplitude. 
6.1.3. Ideal MHD limit
We ran one simulation with zero explicit Ohmic resistivity (ℛ_{m} → ∞), labeled M3T10Ri in Table 1. While the GI dynamo feeds on spiral motions at the largest scales of the flow, Riols & Latter (2019) also found a smallscale dynamo operating in the ideal MHD limit that remained sensitive to numerical resolution. Being mindful of convergence issues related to the finite numerical dissipation that dominates at the grid scale, we restrict our expectations to a qualitative transition in dynamo behavior.
The meanfield to total magnetic energy ratio, χ ≈ 0.17, takes its smallest value in the ideal MHD limit, translating into a mostly disordered magnetic field. Regardless, we measured a constant (global) magnetic growth rate ω/Ω_{in} ≈ 3.7 × 10^{−3} even in this case. Surprisingly, the growth rate of the weakly resistive run M3T10R32 (ℛ_{m} ≈ 32) is smaller than both the ideal MHD case M3T10Ri and the more resistive reference case M3T10R10. It seems unlikely that our simulations can properly resolve any smallscale dynamo in the ideal MHD limit (Riols & Latter 2018a, 2019), so this ordering may be a fortuitous consequence of the oscillating mode prevailing in run M3T10R32.
We could also fit a global eigenmode of Eq. (10) onto the profile of in run M3T10Ri, as shown in Fig. 13. We thus estimate an effective magnetic Reynolds number ζ ≈ 11, which measures the spread of magnetic energy under the action of turbulence alone. The fact that it is only three times larger than in other runs suggests that the GI dynamo may operate in a broad range of diffusivity conditions, without relying on a large Ohmic resistivity (e.g., Riols et al. 2021). Remarkably, the effective Reynolds number associated with the turbulent transport of angular momentum is also ℛ_{e} = 1/α_{SS} ≈ 10 here, where α_{SS} is the viscosity coefficient of Shakura & Sunyaev (1973) determined by local thermal balance (Gammie 2001). Their ratio gives us a crude but first estimate of the turbulent magnetic Prandtl number 𝒫_{m} ≡ ζ/ℛ_{e} ∼ 1 in ideal gravitoturbulent disks.
6.2. Larger initial disk mass
We varied the disk mass by rescaling the initial gas surface density while keeping Σ ∝ R^{−2}. More massive disks tend to be thicker, such that their Toomre number remains near unity after GI saturation. Run M2T10R10 has an increased disk mass M = 1/2, corresponding to an increased aspect ratio h_{ρ}/R ≈ 0.06, and it also features a steady exponential growth of magnetic energy. Its global magnetic growth rate ω/Ω_{in} ≈ 2.8 × 10^{−3} is about 1.8 times larger than in the reference case M3T10R10. Unfortunately, the different contributions to the induction equation nearly balance each other, such that magnetic growth results from a small residual that is delicate to measure.
The meridional profiles of reduced magnetic field b, as functions of z/h_{ρ}, closely resemble those of run M3T10R10 (see Fig. 8), as does the partition of magnetic energy into , , and . The fractions of meanfield to total magnetic energy χ ≈ 0.6 are compatible in both runs. The different contributions A_{ji} and S_{ji} are also within 20 per cent of the reference case shown in Figs. 9 and 10. When comparing the bestfit eigenmodes of Eq. (10) we find that runs M3T10R10 and M2T10R10 have similar effective Reynolds numbers ζ ≈ 2.6 and that the local growth rate ϖ entirely accounts for the factor of 1.8 of increase in global growth rate. The local dynamo process therefore seems to operate more efficiently in thicker (more massive) disks.
6.3. Longer cooling time
In steady state, the disk must generate enough heat to balance cooling, and hence the prescribed cooling time τ controls the strength of GI turbulence (Gammie 2001). The range of cooling times accessible to global simulations is limited, however. For small values τ ≲ 3, nonmagnetized disks are prone to fragmentation (e.g., Lin & Kratter 2016; Booth & Clarke 2019). For long cooling times, the disk takes longer to reach a quasisteady state, and turbulence may become so weak that numerical diffusion could interfere with GI saturation. For these reasons, we only considered one larger value of τ = 10^{3/2} ≈ 32 in run M3T32R10.
The global magnetic growth rate ω/Ω_{in} ≈ 5 × 10^{−3} is approximately three times larger in run M3T32R10 compared to the reference run M3T10R10. This increased growth rate is at odds with the results of Riols & Latter (2019, see their Fig. 3), who found a monotonic decrease of the dynamo growth rates with increasing τ ∈ [5,100]. Unfortunately, as for the highmass case M2T10R10 above, we could not find convincing evidence of an increased growth rate in the average induction equation. On the one hand, the turbulent velocity dispersion (hence ϵ′) and the mean vertical inflow speed are smaller by a factor of ≈0.6. On the other hand, the local growth rate ϖ ≈ 1.5 × 10^{−2} fitting an eigenmode of Eq. (10) is twice larger than in the reference run.
The fact that magnetic energy grows faster at a larger τ ≈ 32 despite weaker velocity fluctuations is puzzling, but must stem from an increased regularity of the flow, leading to a higher yield of the dynamo process from the available kinetic energy. Although there is no hint of this in the thermal stratification of the disk – which is nearly identical to that shown in Fig. 6 – the meanfield to total magnetic energy ratio χ ≈ 0.75 is indeed substantially larger in run M3T32R10. Because turbulence weakens with increasing τ, we do expect an eventual decay of the magnetic growth rates for τ > 32. As for the case of resistivity above, the GI dynamo may therefore admit an optimal cooling time.
7. Saturation of the GI dynamo
Except for runs in which the GI dynamo was the slowest (M3T10R32 and M3T10R1) the disk always reached a magnetized state such that B^{2}/P ≳ 1. We believe that the two other simulations would reach this state if integrated for longer times because a dynamo operates in both of them and only the Lorentz force ∝B^{2} can bring it out of the linear regime. It is expected that the dynamo saturates with a plasma beta near unity because the turbulent flow is transsonic, and thus the Alfvénic Mach number and plasma beta are of the same order. After saturation, the strong field may support additional instabilities competing with the GI in driving turbulence, or cause the disk to puff up and possibly kill the GI, at least for some period of time.
In this section, we take a brief look at the saturated phase of the GI dynamo as it occurs in our simulations. We decompose the sources of stress in Sect. 7.1 and find a magnetically dominated turbulence whose accretion heating exceeds the imposed cooling. We then look at the Toomre number in Sect. 7.2 and confirm that the GI is marginal at best, and quite likely quenched in this phase of our simulations.
7.1. Accretion power and thermal imbalance
One key quantity to study mass accretion in thin disks is the radial flux of angular momentum. We decompose it as the sum of a hydrodynamic (Reynolds) stress R_{Rφ}, a magnetic (Maxwell) stress M_{Rφ}, and a gravitational stress G_{Rφ} defined by:
(see, for example, Balbus & Papaloizou 1999). To connect with standard viscous disk theory, we normalize these momentum fluxes by the vertically and azimuthally averaged gas pressure, denoting by α_{R, M, G} the respective dimensionless coefficients, and by simply α the sum of all three contributions. If the conversion from orbital to thermal energy happens locally, similarly to a viscous process, then there is a unique value of α such that viscous heating balances the imposed cooling (Gammie 2001). Béthune et al. (2021) verified that this equality holds within their framework of nonmagnetized disks, where shocks provide a heating mechanism without explicit viscosity.
We show in Fig. 17 the different contributions to the radial flux of angular momentum after dynamo saturation in the reference run M3T10R10. The gravitational stress α_{G} ≈ 3 × 10^{−2} is subdominant inside r/r_{in} ≲ 12, and we checked that it keeps decaying outward over time in all the dynamosaturating runs. The magnetic stress corresponds to an effective viscosity α_{M} ≈ 10^{−1} roughly constant in the inner parts of the disk. This is the value taken by α_{G} in purely hydrodynamical simulations. The Reynolds stress is maximal near the inner radial boundary, where its effective viscosity reaches α_{R} ≈ 3 × 10^{−1}, but it decays rapidly with radius. Deng et al. (2020) reported the same trends in their ideal MHD simulations (see their Fig. 12).
Fig. 17. Radial flux of angular momentum as a function of radius after averaging in (θ,φ) and over time from 1136t_{in} to 1200t_{in} in run M3T10R10: total stress (solid black), gravitational stress α_{G} (dashed blue), magnetic stress α_{M} (dotdashed red), and hydrodynamic stress α_{R} (dotted green). The thin horizontal line at α = 10^{−1} corresponds to local thermal balance of viscous heating versus cooling. 
The effective viscosity required for local thermal balance in run M3T10R10 is α = 10^{−1}. Since the total stress is substantially larger than this value, we deduce that the rate of turbulent energy injection is larger than the rate of thermal energy extraction by cooling. Whether the dissipation of turbulent energy remains a local process or not, we expect this energy excess to result in a net heating of the gas inside the computational domain. Indeed, we find that the disk thickness h_{ρ}/R ≲ 7 × 10^{−2} overall increases as the dynamo saturates, especially in the outer parts r/r_{in} ≳ 4, which implies an increase in temperature. We conclude that the disk is out of thermal equilibrium and heats up for at least 200t_{in} in run M3T10R10 (300t_{in} in runs M2T10R10 and M3T10R3).
The disk also appears to thicken in the simulations of Deng et al. (2020, comparing the central and rightmost columns of the bottom row of their Fig. 3). They argued that vertical outflows provide the additional cooling mechanism required to keep the disk steady. This is possible if there is a vertical heat flux making the disk corona hotter than its midplane – as is our case – such that the leaving gas has a high specific internal energy. However, our meridional boundary conditions forbid outflows and thus trap heat inside the domain. Steady thermally driven winds also require sound speeds of Keplerian amplitude at the base of the wind. This is not realized in our simulations (c_{s}/v_{K} ≲ 0.3 according to Fig. 6), and Deng et al. (2020) acknowledged that the “clipped” gas particles considered as outflow might actually fall back to the disk. Whether the disk can achieve a steady state after GI dynamo saturation thus remains uncertain.
7.2. Alternative sources of turbulence
Béthune et al. (2021) measured an average Toomre number
in nonmagnetized disks, when using the densityweighted epicyclic frequency and isothermal sound speed c_{ρ}. Figure 18 shows that the Toomre number Q_{ρ} ≳ 2 after dynamo saturation in run M3T10R10. This is partly due to the increased c_{ρ} and partly to the decreased Σ following radial mass transport. Deng et al. (2020) similarly observed Q increase in their MHD disks relative to purely hydrodynamic ones (see their Fig. 7). Since the magnetic field has reached a thermal strength, it may also transform the linear stability of the disk with respect to GI (e.g., Lin 2014). In the axisymmetric limit, magnetic pressure plays a stabilizing role and enhances the Toomre number by a factor of (Gammie 1996; Kim & Ostriker 2001). The dynamo saturation appears to not only alter the gravitoturbulent flow so as to stop magnetic growth: it also appears to push the (hydrodynamic) GI toward – and possibly beyond – marginal stability.
Fig. 18. Evolution in time (horizontal axis) of the radial profile (vertical axis) of the densityweighted Toomre number during the dynamo saturation in run M3T10R10. 
As a consequence, other mechanisms may be involved in sustaining turbulence after GI dynamo saturation. The mean magnetic field remains predominantly toroidal, but it jumps from before saturation to after. This mean field may support MRI growth on spatially resolved wavelengths (Balbus & Hawley 1992) or magnetic buoyancy may locally overcome the entropy stratification and drive Parker instabilities (Shu 1974; Johansen & Levin 2008). The magnetic field may also be strong enough to allow various global instabilities (Curry & Pudritz 1995, 1996; Das et al. 2018). Evidence for an associated transition in the turbulent (and dynamo) state includes the surge of magnetic energy occurring at t = 1000t_{in} (see Fig. 2). Even if GIs become inefficient according to the Toomre number, this turbulence retains enough compressible and nonaxisymmetric character to support a significant gravitational stress α_{G} in Fig. 17. Seeing how it is sustained for hundreds of orbits, we can exclude the onandoff turbulent cycles suggested by Riols & Latter (2019) over such timescales.
Finally, we note the formation of a series of bands near 1000t_{in} in Fig. 18. They survive for at least 200t_{in} in run M3T10R10 and also appear in the Toomre number of other simulations. They correspond to modulations of the gas surface density, with a weak and anticorrelated modulation of the gas sound speed. However, we verified that there are no genuine gas rings forming in the disk: the density fluctuations are mainly nonaxisymmetric. Such modulations did not appear in purely hydrodynamic disks, nor in the global MHD simulations of Deng et al. (2020). It is unclear whether they result from the particular choices made in our numerical setup, and certainly calls for further longterm simulations outside the scope of this paper.
8. Summary and perspectives
We studied the turbulent amplification of magnetic fields in gravitationally unstable disks by means of global 3D MHD simulations. Our setup was the same as that of Béthune et al. (2021), but seeded with a weak toroidal magnetic field and endowed with Ohmic resistivity, the latter sufficient to suppress the MRI initially. We witnessed the exponential growth of magnetic energy and subsequent saturation of the dynamo in simulations with various values of the imposed cooling time and magnetic Reynolds number. Our main results can be summarized as follows.

During its kinematic phase, the GI dynamo grows a largescale and predominantly toroidal magnetic field. It is a global process facilitated by a substantial amount of resistivity (ℛ_{m} ≲ 10), though it may also operate in the ideal MHD limit using the effective resistivity of GI turbulence.

The feedback loop leading to runaway magnetic amplification is the same as in the local simulations of Riols & Latter (2019), for which vertical (outofplane) motions are crucial. The vertical shear resulting from a global thermal stratification typically opposes the dynamo.

The kinematic dynamo can be described via simple meanfield models upon azimuthal averaging. While its mechanism obeys the classic stretchtwistfold terminology, its effect on the mean fields formally differs from traditional αΩ models.

Magnetic amplification stops near equipartition of the magnetic and turbulent energy densities. A saturated state of strong magnetization (i.e., plasma β ∼ 1) is likely a robust result of the dynamo mechanism, and will be shared by both circumstellar and AGN settings. Here the magnetic stress becomes the primary source of turbulent energy, heating the disk faster than it can cool and bringing it away from GI, in agreement with Deng et al. (2020).
Our focus through most of this paper was on the kinematic (linear) phase of the GI dynamo. Although our simulations reach nonlinear dynamo saturation, we believe that this final phase is largely determined by boundary conditions and limited integration times, and hence predictions regarding actual astrophysical objects might be premature at this point. The heat accumulating inside the computational domain may power outflows in more extended domains (Deng et al. 2020); the subsequent evacuation of magnetic flux may allow a new and hot steady state or overshoot and bring the disk back to the kinematic regime in a cyclic fashion (Riols & Latter 2019). This calls for longterm simulations of the nonlinear dynamo phase in extended domains.
Another route open for investigation concerns the microphysics. We have opted for the simplest modeling of the gas resistivity and radiative energy losses, fixing the magnetic Reynolds number ℛ_{m} and cooling time τ to be constant, though in real disks both vary. On the one hand, a realistic treatment of radiative transfer can help clarify the role of GIinduced magnetic fields on disk fragmentation (Deng et al. 2021). On the other, circumstellar disks are notoriously nonideal electrical conductors, and the existence and properties of a GI dynamo should be further explored in realistic Hall and ambipolar MHD regimes (Riols et al. 2021).
For an example of waveinduced dynamo, see Vishniac et al. (1990).
We denote the dimensionless cooling time by τ, as Gammie (2001) and Riols & Latter (2019), in place of β used by Béthune et al. (2021).
The growth of the global MRI is limited in an analogous way, though in that context, it is magnetic tension rather than diffusion that links disparate radii into a coherent global mode (Latter et al. 2015).
Except in run M3T10Ri for which we imposed at R = r_{in} only to better match simulation profiles; see Fig. 13.
Acknowledgments
W.B. acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG) through Grant KL 650/311, an instructive discussion with Yannick Ponty about numerical dynamo experiments, and a most uplifting visit to the DAMTP. H.N.L. acknowledges funding from STFC grant ST/T00049X/1. The authors acknowledge support by the High Performance and Cloud Computing Group at the Zentrum für Datenverarbeitung of the University of Tübingen, the state of BadenWürttemberg through bwHPC and the DFG through grant no INST 37/9351 FUGG. They also thank the anonymous reviewer, Antoine Riols, and François Rincon for their constructive feedback on the submitted version of the paper. The preparation of this paper was overshadowed by mentor and friend Willy Kley’s passing away; we carried his positive attitude throughout.
References
 Alexiades, V., Amiez, G., & Gremaud, P.A. 1996, Commun. Numer. Methods Eng., 12, 31 [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1992, ApJ, 400, 610 [CrossRef] [Google Scholar]
 Balbus, S. A., & Papaloizou, J. C. B. 1999, ApJ, 521, 650 [Google Scholar]
 Balbus, S. A., & Terquem, C. 2001, ApJ, 552, 235 [NASA ADS] [CrossRef] [Google Scholar]
 Balsara, D. S., & Spicer, D. S. 1999, J. Comput. Phys., 149, 270 [NASA ADS] [CrossRef] [Google Scholar]
 Béthune, W., Latter, H., & Kley, W. 2021, A&A, 650, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 [CrossRef] [Google Scholar]
 Boley, A. C., & Durisen, R. H. 2006, ApJ, 641, 534 [Google Scholar]
 Booth, R. A., & Clarke, C. J. 2019, MNRAS, 483, 3718 [Google Scholar]
 Brandenburg, A. 2018, J. Plasma Phys., 84, 735840404 [Google Scholar]
 Cossins, P., Lodato, G., & Clarke, C. J. 2009, MNRAS, 393, 1157 [Google Scholar]
 Curry, C., & Pudritz, R. E. 1995, ApJ, 453, 697 [NASA ADS] [CrossRef] [Google Scholar]
 Curry, C., & Pudritz, R. E. 1996, MNRAS, 281, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Das, U., Begelman, M. C., & Lesur, G. 2018, MNRAS, 473, 2791 [NASA ADS] [CrossRef] [Google Scholar]
 Del Zanna, L., Bucciantini, N., & Londrillo, P. 2003, A&A, 400, 397 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Deng, H., Mayer, L., & Latter, H. 2020, ApJ, 891, 154 [Google Scholar]
 Deng, H., Mayer, L., & Helled, R. 2021, Nat. Astron., 5, 440 [NASA ADS] [CrossRef] [Google Scholar]
 Durisen, R. H., Boss, A. P., Mayer, L., et al. 2007, in Protostars and Planets V, eds. B. Reipurth, D. Jewitt, & K. Keil, 607 [Google Scholar]
 Evans, C. R., & Hawley, J. F. 1988, ApJ, 332, 659 [Google Scholar]
 Fendt, C., & Gaßmann, D. 2018, ApJ, 855, 130 [NASA ADS] [CrossRef] [Google Scholar]
 Forgan, D., Price, D. J., & Bonnell, I. 2017, MNRAS, 466, 3406 [NASA ADS] [CrossRef] [Google Scholar]
 Fromang, S. 2005, A&A, 441, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fromang, S., Balbus, S. A., & De Villiers, J.P. 2004a, ApJ, 616, 357 [NASA ADS] [CrossRef] [Google Scholar]
 Fromang, S., Balbus, S. A., Terquem, C., & De Villiers, J.P. 2004b, ApJ, 616, 364 [NASA ADS] [CrossRef] [Google Scholar]
 Gammie, C. F. 1996, ApJ, 462, 725 [NASA ADS] [CrossRef] [Google Scholar]
 Gammie, C. F. 2001, ApJ, 553, 174 [Google Scholar]
 Gardiner, T. A., & Stone, J. M. 2005, J. Comput. Phys., 205, 509 [NASA ADS] [CrossRef] [Google Scholar]
 Goodman, J. 2003, MNRAS, 339, 937 [NASA ADS] [CrossRef] [Google Scholar]
 Gressel, O., & Pessah, M. E. 2015, ApJ, 810, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Harten, A., Lax, P. D., & Leer, B. V. 1983, SIAM Rev., 25, 35 [Google Scholar]
 Hawley, J. F., & Balbus, S. A. 1992, ApJ, 400, 595 [NASA ADS] [CrossRef] [Google Scholar]
 Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1996, ApJ, 464, 690 [NASA ADS] [CrossRef] [Google Scholar]
 Heemskerk, M. H. M., Papaloizou, J. C., & Savonije, G. J. 1992, A&A, 260, 161 [NASA ADS] [Google Scholar]
 Johansen, A., & Levin, Y. 2008, A&A, 490, 501 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ju, W., Stone, J. M., & Zhu, Z. 2016, ApJ, 823, 81 [NASA ADS] [CrossRef] [Google Scholar]
 Ju, W., Stone, J. M., & Zhu, Z. 2017, ApJ, 841, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Kawasaki, Y., Koga, S., & Machida, M. N. 2021, MNRAS, 504, 5588 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, W.T., & Ostriker, E. C. 2001, ApJ, 559, 70 [NASA ADS] [CrossRef] [Google Scholar]
 Kratter, K., & Lodato, G. 2016, ARA&A, 54, 271 [NASA ADS] [CrossRef] [Google Scholar]
 Krause, F., & Rädler, K. H. 2016, Meanfield Magnetohydrodynamics and Dynamo Theory (Elsevier) [Google Scholar]
 Kunz, M. W. 2008, MNRAS, 385, 1494 [NASA ADS] [CrossRef] [Google Scholar]
 Latter, H. N., Fromang, S., & Faure, J. 2015, MNRAS, 453, 3257 [Google Scholar]
 Lesur, G., & Ogilvie, G. I. 2008, A&A, 488, 451 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lin, M.K. 2014, ApJ, 790, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, M.K., & Kratter, K. M. 2016, ApJ, 824, 91 [Google Scholar]
 Lin, D. N. C., & Pringle, J. E. 1987, MNRAS, 225, 607 [Google Scholar]
 Londrillo, P., & del Zanna, L. 2004, J. Comput. Phys., 195, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Mamatsashvili, G., Chagelishvili, G., Pessah, M. E., Stefani, F., & Bodo, G. 2020, ApJ, 904, 47 [NASA ADS] [CrossRef] [Google Scholar]
 Michael, S., & Durisen, R. H. 2010, MNRAS, 406, 279 [Google Scholar]
 Mignone, A. 2014, J. Comput. Phys., 270, 784 [Google Scholar]
 Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [Google Scholar]
 Miyoshi, T., & Kusano, K. 2005, J. Comput. Phys., 208, 315 [NASA ADS] [CrossRef] [Google Scholar]
 Moffatt, K., & Dormy, E. 2019, Selfexciting Fluid Dynamos (Cambridge University Press) [CrossRef] [Google Scholar]
 Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610 [Google Scholar]
 Ogilvie, G. I., & Pringle, J. E. 1996, MNRAS, 279, 152 [NASA ADS] [CrossRef] [Google Scholar]
 Paczynski, B. 1978, Acta Astron, 28, 91 [NASA ADS] [Google Scholar]
 Parker, E. N. 1955, ApJ, 122, 293 [Google Scholar]
 Pelletier, G., & Pudritz, R. E. 1992, ApJ, 394, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Pjanka, P., & Stone, J. M. 2020, ApJ, 904, 90 [NASA ADS] [CrossRef] [Google Scholar]
 Rice, W. K. M., Armitage, P. J., Bate, M. R., & Bonnell, I. A. 2003, MNRAS, 339, 1025 [Google Scholar]
 Rincon, F. 2019, J. Plasma Phys., 85, 205850401 [NASA ADS] [CrossRef] [Google Scholar]
 Rincon, F., Ogilvie, G. I., & Proctor, M. R. E. 2007, Phys. Rev. Lett., 98, 254502 [NASA ADS] [CrossRef] [Google Scholar]
 Riols, A., & Latter, H. 2018a, MNRAS, 474, 2212 [NASA ADS] [CrossRef] [Google Scholar]
 Riols, A., & Latter, H. 2018b, MNRAS, 476, 5115 [Google Scholar]
 Riols, A., & Latter, H. 2019, MNRAS, 482, 3989 [Google Scholar]
 Riols, A., Rincon, F., Cossu, C., et al. 2017a, A&A, 598, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Riols, A., Latter, H., & Paardekooper, S. J. 2017b, MNRAS, 471, 317 [Google Scholar]
 Riols, A., Xu, W., Lesur, G., Kunz, M. W., & Latter, H. 2021, MNRAS, 506, 1407 [NASA ADS] [CrossRef] [Google Scholar]
 Roe, P. L. 1986, Ann. Rev. Fluid Mech., 18, 337 [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33 [NASA ADS] [Google Scholar]
 Shi, J.M., & Chiang, E. 2014, ApJ, 789, 34 [Google Scholar]
 Shlosman, I., & Begelman, M. C. 1987, Nature, 329, 810 [Google Scholar]
 Shu, F. H. 1974, A&A, 33, 55 [NASA ADS] [Google Scholar]
 Squire, J., & Bhattacharjee, A. 2016, J. Plasma Phys., 82, 535820201 [NASA ADS] [CrossRef] [Google Scholar]
 Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2011, MNRAS, 418, L79 [NASA ADS] [CrossRef] [Google Scholar]
 Toomre, A. 1964, ApJ, 139, 1217 [Google Scholar]
 Urpin, V., & Brandenburg, A. 1998, MNRAS, 294, 399 [Google Scholar]
 van Leer, B. 1974, J. Comput. Phys., 14, 361 [NASA ADS] [CrossRef] [Google Scholar]
 Vishniac, E. T., & Brandenburg, A. 1997, ApJ, 475, 263 [NASA ADS] [CrossRef] [Google Scholar]
 Vishniac, E. T., Jin, L., & Diamond, P. 1990, ApJ, 365, 648 [NASA ADS] [CrossRef] [Google Scholar]
 von Rekowski, B., Brandenburg, A., Dobler, W., Dobler, W., & Shukurov, A. 2003, A&A, 398, 825 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Walker, J., Lesur, G., & Boldyrev, S. 2016, MNRAS, 457, L39 [Google Scholar]
 Wurster, J., & Bate, M. R. 2019, MNRAS, 486, 2587 [NASA ADS] [Google Scholar]
 Zhao, B., Caselli, P., Li, Z.Y., & Krasnopolsky, R. 2018, MNRAS, 473, 4868 [Google Scholar]
 Zhao, B., Tomida, K., Hennebelle, P., et al. 2020, Space Sci. Rev., 216, 43 [CrossRef] [Google Scholar]
Appendix A: Choice of EMF reconstruction scheme
Integrating the magnetic induction Eq. (4) over a bounded surface 𝒮 and applying Stokes’ theorem gives
where E is the total EMF as seen in the chosen reference frame. The method of constrained transport implements this equation by taking the surfaceintegrated magnetic flux as variable to evolve in time (Evans & Hawley 1988). In a finitevolume framework, the different components of the discretized magnetic field therefore live on their respective interfaces and can be called “facecentered”. On the other hand, the circulation of the EMF involves “edgecentered” variables.
There are different ways to reconstruct the edgecentered EMF given the values of the facecentered and cellcentered MHD variables, with different accuracy and stability properties (Balsara & Spicer 1999; Londrillo & del Zanna 2004; Gardiner & Stone 2005). Throughout this paper we used the UCT_HLL scheme of Del Zanna et al. (2003) implemented in PLUTO 4.3. However, we also tried the less dissipative UCT0 and UCT_CONTACT schemes of Gardiner & Stone (2005) and observed spurious magnetic amplification in both cases.
Figure A.1 shows how the average disk magnetization evolves over time in a simulation equivalent to our reference case (see Fig. 2) but using the UCT0 reconstruction scheme. After about 25t_{in}, the ratio B^{2}/P increases by more than six orders of magnitude in less than ten inner orbits. The corresponding growth rate is too fast to be associated with any fluid dynamo, leaving boundary effects as the only possibility.
Fig. A.1. Evolution of the average magnetization ratio ⟨B^{2}⟩/⟨P⟩ in radius (vertical axis) and time (horizontal axis) in a simulation equivalent to our reference case M3T10R10 (see Fig. 2) but using the UCT0 reconstruction scheme. This ratio increases by more than six orders of magnitude in less than ten inner orbits at time t/t_{in} ≈ 25. 
We set the tangential components of the magnetic field to zero in the ghost cells, and we verified that there is no tangential magnetic flux injection appearing in the meridional profiles of b as defined by (5) and measured at time 100t_{in}. However, it is the ∇ ⋅ B = 0 condition that controls the component normal to the boundaries. We found that the average proportion of toroidal field drops to as soon as magnetic energy goes up. We deduce that magnetic energy is injected in B_{θ} fluctuations from the meridional boundaries, and confirm this after examining and in the equivalent of Fig. 11.
If magnetic energy was injected from the domain boundaries in our simulations with UCT_HLL, it did so with no significant consequence. Interestingly, the magnetization saturates below B^{2}/P ≲ 10^{−2} in Fig. A.1, and the reduced magnetic field b displays the same meridional profiles as in Fig. 8. These two points suggest that the GI dynamo may still be operating underneath the boundarydriven amplification, though at its slower pace.
Appendix B: Minimal meanfield instability
Equations (15) & (16) can be further reduced to the minimal form:
where, for simplicity, we took the same offdiagonal resistivities η′> 0 in both equations. We can then set η′ = 1 and focus on z ∈ [0,1] by rescaling space and time without loss of generality. We look for eigenmodes satisfying the conditions described in Sect. 5.4 under a prescribed velocity V_{z} ∝ z. In the case of V_{z} = 0, the fundamental mode has a purely imaginary eigenvalue with frequency ω_{0} ≈ 5.59. Unstable (real positive) eigenvalues appear when the amplitude of the compression rate ∂_{z}V_{z} becomes comparable to the fundamental frequency ω_{0}, as illustrated in Fig. B.1. For larger values of the compression rate, unstable eigenvalues acquire a nonzero imaginary part and therefore correspond to growing oscillations — similar to Fig. 16.
Fig. B.1. Eigenvalues of Eqs. (B.1) & (B.2) for unstable modes with η′ = 1 over z ∈ [0,1], satisfying the conditions described in Sect. 5.4, and under different compression rates −∂_{z}V_{z} relative to the fundamental mode frequency ω_{0} obtained at V_{z} = 0; blue (resp. red) crosses represent the real (resp. imaginary) part of the eigenvalue. 
All Tables
All Figures
Fig. 1. Flattened distribution of the gas surface density R^{2}Σ at time 400t_{in} over the whole domain of run M3T10R10; colors range linearly from zero to three times the median of R^{2}Σ. 

In the text 
Fig. 2. Evolution in time (horizontal axis) of the radial profile (vertical axis) of the disk magnetization B^{2}/P in run M3T10R10. 

In the text 
Fig. 3. Growth rate of the magnetic field amplitude as a function of radius in run M3T10R10, normalized either by the (fixed) inner Keplerian frequency Ω_{in} (solid blue) or by the (radially varying) densityweighted orbital frequency Ω_{ρ} (dashed red). Only the inner parts r/r_{in} ≲ 12 have reached a steady exponential growth phase. 

In the text 
Fig. 4. Equatorial slice in run M3T10R10 at time 400t_{in}, zoomed in on the radial interval r/r_{in} ∈ [2,16]. The color map shows the relative density fluctuations as defined by Eq. (6) while the arrows show the orientation of the magnetic field sampled in the midplane, colored in blue when inward (B_{r} < 0) or green when outward (B_{r} > 0). 

In the text 
Fig. 5. Meridional slice at time 400t_{in} and angle φ = 3π/8 in run M3T10R10. Left panel: density distribution flattened by a factor of (R/r_{in})^{3} (color scale) and orientation of the poloidal velocity field, green (resp., blue) being oriented upward (resp., downward). Right panel: toroidal magnetic field relative to the local mean (color scale) and integral curves of the poloidal magnetic field. The density is maximal at r/r_{in} ≈ 4.0, 5.6, and 6.9. 

In the text 
Fig. 6. Thermal stratification of the disk. The solid blue line (left axis) shows the gas temperature relative to its midplane value. The dashed red line (right axis) shows the squared buoyancy frequency (Eq. (7)) relative to the local squared orbital frequency. Both profiles were averaged azimuthally, over r/r_{in} ∈ [2,8] and over time from 200t_{in} to 400t_{in}. 

In the text 
Fig. 7. Meridional profiles of the mean flow as defined by Eq. (5), averaged from 200t_{in} to 400t_{in} in run M3T10R10, and focused on the midplane region θ ∈ π/2 ± 0.1. Upper panel: spherical components of the reduced gas velocity. Lower panel: components of the reduced mass flux. The toroidal components are measured on the right vertical axes. 

In the text 
Fig. 8. Meridional profiles of the three components of the reduced magnetic field as defined by Eq. (5) and averaged from 200t_{in} to 400t_{in} in run M3T10R10; the toroidal component b_{φ} is measured on the right axis. 

In the text 
Fig. 9. Contributions to the induction of a mean toroidal field after normalizing Eq. (8) by and averaging azimuthally, radially over r/r_{in} ∈ [2,8], and over time from 200t_{in} to 400t_{in} in run M3T10R10. 

In the text 
Fig. 10. Same as Fig. 9 but for the induction of a mean radial field. 

In the text 
Fig. 11. Contributions to the energy stored in after normalizing by at every radius and averaging azimuthally, radially over r/r_{in} ∈ [2,8], and over time from 200t_{in} to 400t_{in} in run M3T10R10. 

In the text 
Fig. 12. Dynamo mechanism as seen in the frame of a spiral wake; see the text in Sect. 4.1.5 for a stepbystep description. 

In the text 
Fig. 13. Normalized profiles of the magnetic field amplitude measured before dynamo saturation in four simulations featuring a steady exponential growth (solid lines; see legend), and the bestfit eigenmodes from Eq. (10) having the same global growth rate (dashed lines). 

In the text 
Fig. 14. Meridional profiles of the reduced EMF as defined by Eq. (5). The solid lines are simulation data averaged from 200t_{in} to 400t_{in} in run M3T10R10. The dashed lines are the best fits from Eq. (11) over this interval. The toroidal components are measured on the right axis. 

In the text 
Fig. 15. Eigenmode of (15)–(16) under a prescribed converging flow , satisfying (∂_{z}B_{R},∂_{z}B_{φ}) = 0 at z/h = 0, B_{R} = 0 at z/h = 2, and ∫B_{φ} dz = 0. This solution includes a Keplerian radial shear rate d log Ω/ d log R = −3/2, an Ohmic resistivity η = 0.1Ωh^{2}, and the dynamo coefficients estimated in the reference simulation M3T10R10: , , , and . It is marginally unstable with a growth rate of 5.3 × 10^{−4}Ω. 

In the text 
Fig. 16. Evolution in time (horizontal axis) of the radial profile (vertical axis) of the disk magnetization in the lowresistivity run M3T10R32. The mean magnetic field flips sign twice while growing in amplitude. 

In the text 
Fig. 17. Radial flux of angular momentum as a function of radius after averaging in (θ,φ) and over time from 1136t_{in} to 1200t_{in} in run M3T10R10: total stress (solid black), gravitational stress α_{G} (dashed blue), magnetic stress α_{M} (dotdashed red), and hydrodynamic stress α_{R} (dotted green). The thin horizontal line at α = 10^{−1} corresponds to local thermal balance of viscous heating versus cooling. 

In the text 
Fig. 18. Evolution in time (horizontal axis) of the radial profile (vertical axis) of the densityweighted Toomre number during the dynamo saturation in run M3T10R10. 

In the text 
Fig. A.1. Evolution of the average magnetization ratio ⟨B^{2}⟩/⟨P⟩ in radius (vertical axis) and time (horizontal axis) in a simulation equivalent to our reference case M3T10R10 (see Fig. 2) but using the UCT0 reconstruction scheme. This ratio increases by more than six orders of magnitude in less than ten inner orbits at time t/t_{in} ≈ 25. 

In the text 
Fig. B.1. Eigenvalues of Eqs. (B.1) & (B.2) for unstable modes with η′ = 1 over z ∈ [0,1], satisfying the conditions described in Sect. 5.4, and under different compression rates −∂_{z}V_{z} relative to the fundamental mode frequency ω_{0} obtained at V_{z} = 0; blue (resp. red) crosses represent the real (resp. imaginary) part of the eigenvalue. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.