© W. Béthune and H. Latter 2022

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¹ 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 so-called 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 self-sustained 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 out-of-plane dynamics outright by employing two-dimensional 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 three-dimensional (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 large-scale 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 (geometry-dependent) 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 mean-field 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 grid-based and shock-capturing 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 self-gravitating 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, ${\mathrm{\Omega }}_{\star }\equiv \sqrt{m_{\star }/R^3}$ 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:

$$\begin{array}{cc}& {\partial }_t\rho =-\mathrm{\nabla }·\left(\rho \mathit{V}\right),\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill & \rho ({\partial }_t\mathit{V}+\mathit{V}·\mathrm{\nabla }\mathit{V})=-\mathrm{\nabla }P-\rho \mathrm{\nabla }\mathrm{\Phi }+\mathit{J}\times \mathit{B},\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill & {\partial }_{t}P+\mathit{V}·\mathrm{\nabla }P=-\gamma P\mathrm{\nabla }·\mathit{V}+(\gamma -1)\eta J^2+\mathrm{\Lambda },\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill & {\partial }_t\mathit{B}=\mathrm{\nabla }\times (\mathit{V}\times \mathit{B}-\eta \mathit{J}).\hfill \end{array}$$(4)

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 time-dependent 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²τ. 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²/η. Since the characteristic thickness h of the disk can vary with its thermal stratification and self-gravity, 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 $c_{\mathrm{s}}\equiv \sqrt{P/\rho }$. The resistivity is therefore constant in time according to the prescribed value of ${\mathcal{R}}_m\equiv {\mathrm{\Omega }}_{\star }{h}_{\text{init}}^2/\eta$. The actual ratio Ωh²/η 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 second-order Runge-Kutta 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 two-state 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 super-time-stepping (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 $V_{\phi }=\sqrt{[m_{\star }+m_{\mathrm{disk}}\left(R\right)]/R}$, accounting for the enclosed disk mass $m_{\mathrm{disk}}\left(R\right)={\int }_{r{}_{\text{in}}}^{R}2\mathrm{\pi }\mathrm{\Sigma }x\text{d}x$. The initial surface density Σ ∝ R⁻² 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 $B_{\phi }^2={10}^{-12}P$ 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 quasi-steady 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⁶ 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²/P ≳ 10⁻².

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 mean-field 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 time-independent 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 density-weighted 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 density-weighted angular velocity Ω_ρ ≡ ⟨V_φ/R⟩_ρ and isothermal sound speed $c_{\rho }\equiv \sqrt{{\langle P\rangle }_{\theta ,\phi }/{\langle \rho \rangle }_{\theta ,\phi }}$ to define radial profiles used as averaging weights: a typical thickness h_ρ ≡ c_ρ/Ω_ρ, a typical velocity $\overline{V}\equiv {\mathrm{\Omega }}_{\rho }R$, and a typical magnetic field strength $\overline{B}\equiv \sqrt{{\langle B^2\rangle }_{\theta ,\phi }}$. 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:

$$\begin{array}{c}\hfill (\mathit{v},\mathit{q},\mathit{b},\mathit{j},\mathit{ϵ})\equiv {\langle (\frac{\mathit{V}}{\overline{V}},\frac{\rho \mathit{V}}{\mathrm{\Sigma }{\mathrm{\Omega }}_{\rho }},\frac{\mathit{B}}{\overline{B}},\frac{\mathit{J}}{\overline{B}/h_{\rho }},\frac{\mathcal{E}}{\overline{V}\overline{B}})\rangle }_{\phi ,r}.\end{array}$$(5)

Each dimensionless ratio on the right-hand side effectively fluctuates about the average on the left-hand 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²/P exhibiting a large-scale envelope that grows by 12 orders of magnitude over 900t_in.

Fig. 1. Flattened distribution of the gas surface density R2Σ at time 400tin over the whole domain of run M3T10R10; colors range linearly from zero to three times the median of R2Σ.

Noting that $P=\rho c_{\text{s}}^2$ and that the disk supports sonic turbulence, the Alfvén velocity $V_{\mathrm{A}}\equiv \sqrt{B^2/\rho }$ remains negligible compared to both the velocity fluctuations and the sound speed as long as B²/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 initial-value problem for B in a given velocity field V. We focus here on this kinematic dynamo phase and postpone the dynamo saturation (B²/P ≳ 1) to Sect. 7.

As long as B is dynamically unimportant, its amplitude $\overline{B}$ should vary exponentially in time. We measured the evolution of $\overline{B}$ during one outer orbit ≈181t_in before reaching B²/P = 10⁻² 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 $\frac{1}{2}log\left({\overline{B}}^2\right)$ as a function of time.

The solid blue curve in Fig. 3 indicates that the growth rate ω ≈ 1.6 × 10⁻³Ω_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 large-scale mean-field 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 B2/P in run M3T10R10.

Fig. 3. Growth rate of the magnetic field amplitude $\overline{B}\equiv \sqrt{{\langle B^2\rangle }_{\theta ,\phi }}$ as a function of radius in run M3T10R10, normalized either by the (fixed) inner Keplerian frequency Ωin (solid blue) or by the (radially varying) density-weighted orbital frequency Ωρ (dashed red). Only the inner parts r/rin ≲ 12 have reached a steady exponential growth phase.

3.2. Hydrodynamic (GI) turbulence

In this subsection, we describe the small-scale 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

$$\begin{array}{c}\hfill \stackrel{\sim }{\rho }\equiv \frac{{\langle \rho -{\langle \rho \rangle }_{\theta ,\phi }\rangle }_{\theta }}{{\langle \rho \rangle }_{\theta ,\phi }}\end{array}$$(6)

Fig. 4. Equatorial slice in run M3T10R10 at time 400tin, zoomed in on the radial interval r/rin ∈ [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 (Br < 0) or green when outward (Br > 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 400tin and angle φ = 3π/8 in run M3T10R10. Left panel: density distribution flattened by a factor of (R/rin)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/rin ≈ 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. Small-scale 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 low-density region between wakes, in agreement with local simulations (Riols & Latter 2019). While in-plane 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²/η 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 quasi-steady 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/rin ∈ [2,8] and over time from 200tin to 400tin.

To quantify the buoyant response of the gas to small vertical displacements, we define a squared buoyancy frequency

$$\begin{array}{c}\hfill N_z^2\equiv -\frac{1}{\gamma }\left({\partial }_z\mathrm{\Phi }\right)({\partial }_{z}log\left[s\right]),\end{array}$$(7)

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 $N_z^2\ge 0$ 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 sub-adiabatic, 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 200tin to 400tin 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 high-density 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. Large-scale magnetic field

Figure 2 shows the radially global structure of the large-scale 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 200tin to 400tin 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 large-scale kinematic dynamo

In this section, we examine the growth of a large-scale magnetic field during the B²/P ≪ 1 phase of our reference run M3T10R10. As mentioned in Sect. 3.1, the induction Eq. (4) then reduces to a linear initial-value 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 large-scale 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:

$${\mathrm{\partial }}_t{B}_i={\displaystyle \sum _j}\underset{A_{ji}}{\underbrace{-{\mathrm{\partial }}_j\left[V_j{B}_i\right]}}+\sum _j\underset{S_{ji}}{\underbrace{{\mathrm{\partial }}_j\left[B_j{V}_i\right]}}+\mathrm{Ohm}.$$(8)

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_rB_φ] 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

$$\begin{array}{c}\hfill {\langle \mathcal{E}\rangle }_{\phi }=-{\langle \mathit{V}\rangle }_{\phi }\times {\langle \mathit{B}\rangle }_{\phi }-{\langle {\mathit{V}}^{\prime }\times {\mathit{B}}^{\prime }\rangle }_{\phi }.\end{array}$$(9)

We denote by $A_{\mathit{ji}}^{\prime }$ and $S_{\mathit{ji}}^{\prime }$ 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 large-scale component of the magnetic field, we also average Eq. (8) in the radial dimension after normalizing it by ${\mathrm{\Omega }}_{\rho }\overline{B}$ 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 ${\mathrm{\Omega }}_{\rho }\overline{B}$ and averaging azimuthally, radially over r/rin ∈ [2,8], and over time from 200tin to 400tin 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 advection-diffusion 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 $S_{r\phi }^{\prime }$ 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 mean-field dynamo theory, and it is clearly operating in our simulations, in agreement with Riols & Latter (2018a, 2019).

Finally, the S_θφ term (dot-dashed, 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_zV_φ]. Its fluctuation part $S_{\theta \phi }^{\prime }$ 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.

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

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 $A_{\theta r}^{\prime }$ (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_rB_θ] (dot-dashed, 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_zV_R]. It is dominated by its fluctuating part $S_{\theta r}^{\prime }$, and hence the induction of a mean ⟨B_r⟩_φ cannot be directly traced back to another mean-field 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 $S_{\theta r}^{\prime }$ does rely on nonaxisymmetric fluctuations $B_{\theta }^{\prime }$. Because ⟨B_θ⟩ is so small, $\langle B_{\theta }^2\rangle$ 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, ${\partial }_t{B}_{\theta }^2$. We normalize each term by the local ${\mathrm{\Omega }}_{\rho }{\overline{B}}^2$, average them as previously, and plot them in Fig. 11.

Fig. 11. Contributions to the energy stored in $B_{\theta }^2$ after normalizing ${\partial }_t{B}_{\theta }^2$ by ${\mathrm{\Omega }}_{\rho }{\overline{B}}^2$ at every radius and averaging azimuthally, radially over r/rin ∈ [2,8], and over time from 200tin to 400tin 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 $B_{\theta }^{\prime }$ 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 $B_{\theta }^2$ and should eventually vanish upon averaging. The S_rθ term (dot-dashed 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_RV_z] happens at the expense of magnetic energy.

The energy exchanges of the other components $B_r^2$ and $B_{\phi }^2$ 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 $B_{\phi }^2$ but also slightly into $B_r^2$, 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 ∂_zV_φ 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 step-by-step 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 $\overline{B}(R,t)$, the fact that ${\partial }_t\overline{B}\simeq \omega \overline{B}$ 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 $\overline{B}$ 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 $\varpi \mathrm{\Omega }\overline{B}$ 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², where ζ denotes an effective magnetic Reynolds number. We therefore propose a simple reaction-diffusion equation to describe the evolution of the magnetic field amplitude:

$$\begin{array}{c}\hfill \frac{\partial \overline{B}}{\partial t}\simeq \varpi \mathrm{\Omega }\overline{B}+\frac{1}{R}\frac{\partial }{\partial R}\left(R\zeta \mathrm{\Omega }{h}^2\frac{\partial \overline{B}}{\partial R}\right),\end{array}$$(10)

consisting of local growth plus radial diffusion.

In the next paragraphs, we solve Eq. (10) numerically. But assuming that $\overline{B}\propto {\text{e}}^{\omega t}$, 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³.

4.2.2. Fitting global eigenmodes

We validated Eq. (10) a posteriori by comparing its eigenmodes to simulation data, namely the radial profile of $\overline{B}$ at the time that the average disk magnetization reaches B²/P = 10⁻² anywhere. We computed eigenmodes of Eq. (10) for given pairs of (ϖ,ζ) by successively integrating it in time and rescaling the solution until convergence, imposing $\overline{B}=0$ at both boundaries⁴ and using the initial simulation profiles of $\overline{B}$, Ω 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 $\overline{B}$ measured before dynamo saturation in four simulations featuring a steady exponential growth (solid lines; see legend), and the best-fit 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⁻²Ω equals the measured global growth rate ω ≈ 1.6 × 10⁻³Ω_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⁻² 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. Mean-field 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 mean-field 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 mean-field 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 so-called α-Ω 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 $\overline{V}\overline{B}$ 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:

$$\begin{array}{c}\hfill {\mathit{ϵ}}^{\prime }={\alpha }^{\prime }·\mathit{b}+{\beta }^{\prime }·\mathit{j},\end{array}$$(11)

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 ${\alpha }_{\phi \phi }^{\prime }$.

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 time-averaged 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 ${\beta }^{\prime }\approx {\mathcal{R}}_m{}^{-1}I$, where I is the identity matrix.

Figure 14 shows the time-averaged turbulent EMF in the reference run M3T10R10, along with the best fit profiles of Eq. (11). We could always obtain good fits for ${ϵ}_r^{\prime }$ and ${ϵ}_{\phi }^{\prime }$ when including the β′ matrix, but not with the α′ matrix alone. These matrices typically take the following form:

$$\begin{array}{c}\hfill {\alpha }^{\prime }\simeq \left(\begin{array}{ccc}0& {\alpha }_{r\theta }^{\prime }& 0\\ {\alpha }_{\theta r}^{\prime }& 0& 0\\ 0& {\alpha }_{\phi \theta }^{\prime }& 0\end{array}\right),{\beta }^{\prime }\simeq \left(\begin{array}{ccc}{\beta }_{\mathit{rr}}^{\prime }& 0& {\beta }_{r\phi }^{\prime }\\ 0& 0& 0\\ {\beta }_{\phi r}^{\prime }& 0& {\beta }_{\phi \phi }^{\prime }\end{array}\right),\end{array}$$(12)

Fig. 14. Meridional profiles of the reduced EMF as defined by Eq. (5). The solid lines are simulation data averaged from 200tin to 400tin 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 ${\alpha }_{r\theta }^{\prime }>0$, ${\alpha }_{\theta r}^{\prime }\approx 1$, and ${\alpha }_{\phi \theta }^{\prime }<0$. The diagonal coefficient ${\alpha }_{\phi \phi }^{\prime }$ 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 ${\alpha }_{r\theta }^{\prime }$ and ${\alpha }_{\phi \theta }^{\prime }$ connect to the mean meridional field ⟨B_θ⟩_φ, which does not appear explicitly in the dynamo process summarized in Sect. 4.1.5. The ${\alpha }_{\theta r}^{\prime }$ component is the least well constrained because the quality of the fit for ${ϵ}_{\theta }^{\prime }$ varies across simulations – Fig. 14 shows a rather bad case.

In the turbulent resistivity matrix, we find ${\beta }_{\phi \phi }^{\prime }\approx -{{\mathcal{R}}_m}^{-1}$, 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 ${\beta }_{\mathit{rr}}^{\prime }\lesssim 0$ is always smaller than ${\beta }_{\phi \phi }^{\prime }$, implying an inefficient vertical transport of B_φ by turbulence. The off-diagonal terms satisfy ${\beta }_{r\phi }^{\prime }\approx -{{\mathcal{R}}_m}^{-1}$ and ${\beta }_{\phi r}^{\prime }\gtrsim 0$; 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 ${\eta }^{\prime }\equiv {\mathrm{\Omega }}_{\rho }{h}_{\rho }^2{\beta }^{\prime }$ denote the scaled-up turbulent resistivity matrix. Omitting ${\alpha }_{r\theta }^{\prime }$ and ${\alpha }_{\phi \theta }^{\prime }$ from Eq. (12) and injecting the remaining coefficients into the azimuthally averaged induction Eq. (4) leads to the following closed system:

$$\begin{array}{cc}& {\partial }_t{B}_R\simeq -{\partial }_z\left[V_z{B}_R\right]+(\eta +{\eta }_{\phi \phi }^{\prime }){\partial }_z^2{B}_R-{\eta }_{\phi R}^{\prime }{\partial }_z^2{B}_{\phi },\hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill & {\partial }_t{B}_{\phi }\simeq -{\partial }_z\left[V_z{B}_{\phi }\right]+(\eta +{\eta }_{\mathit{RR}}^{\prime }){\partial }_z^2{B}_{\phi }-{\eta }_{R\phi }^{\prime }{\partial }_z^2{B}_R\hfill \\ \hfill & +{\partial }_R([V_{\phi }+{\alpha }_{\mathit{zR}}^{\prime }\overline{V}]B_R+(\eta /R){\partial }_R\left[R{B}_{\phi }\right]).\hfill \end{array}$$(14)

In both equations, the first term on the right-hand 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 off-diagonal turbulent resistivity. The second line of Eq. (14) consists of radial shear (Ω effect) and diffusion.

Approximating $\overline{V}\simeq V_{\phi }$ and ${\alpha }_{\mathit{zR}}^{\prime }\simeq -{\alpha }_{\theta r}^{\prime }\approx -1$, 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 (${\alpha }_{\mathit{zR}}^{\prime }=0$). 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 $({\beta }_{r\theta }^{\prime },{\beta }_{\theta \theta }^{\prime },{\beta }_{\phi \theta }^{\prime })\approx 0$, and hence we lose radial turbulent diffusion in Eq. (14). To make up for it, the best fit of ${ϵ}_{\theta }^{\prime }$ is biased by an unphysical ${\alpha }_{\theta r}^{\prime }\ne 0$. We can therefore discard this only remaining α′ coefficient.

Neglecting the effective resistivity $\eta +{\eta }_{\phi \phi }^{\prime }\approx 0$ in Eq. (13) and the radial gradients of RB_R and RB_φ in Eq. (14), we finally obtain a radially local system of equations:

$$\begin{array}{cc}& {\partial }_t{B}_R\simeq -{\partial }_z\left[V_z{B}_R\right]-{\eta }_{\phi R}^{\prime }{\partial }_z^2{B}_{\phi },\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill & {\partial }_t{B}_{\phi }\simeq -{\partial }_z\left[V_z{B}_{\phi }\right]-{\eta }_{R\phi }^{\prime }{\partial }_z^2{B}_R+(\eta +{\eta }_{\mathit{RR}}^{\prime }){\partial }_z^2{B}_{\phi }\hfill \\ \hfill & +\left(\frac{\mathrm{d}log\mathrm{\Omega }}{\mathrm{d}logR}\right)\mathrm{\Omega }{B}_R,\hfill \end{array}$$(16)

whose originality relative to classic mean-field dynamo models stems from the absence of an α effect, the apparently stabilizing pair of off-diagonal resistivity coefficients, and the presence of a compressive mean flow (Rincon 2019).

5.4. Mean-field 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_RB_φ 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 ∂_tB² = −(∂_zV_z)B². 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 off-diagonal coefficients ${\eta }_{R\phi }^{\prime }$ and ${\eta }_{\phi R}^{\prime }$ 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 right-hand sides of Eqs. (15) and (16). We computed their eigenmodes as a function of z/h ∈ [0,2] using Chebyshev polynomials over 64 Gauss-Lobatto collocation points. We sought symmetric solutions with respect to the midplane (∂_zB_R = ∂_zB_φ = 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 ∂_zV_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_φ,∂_zB_R,∂_zB_φ) ≠ 0. We obtain unstable modes when V_z ≠ 0 and the compression rate −∂_zV_z exceeds the waves’ original frequency (see Appendix B), indicating that the present instability need not be shear-driven. 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⁻⁴Ω 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 $V_z/\overline{V}=-0.11z/h$, satisfying (∂zBR,∂zBφ) = 0 at z/h = 0, BR = 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Ωh2, and the dynamo coefficients estimated in the reference simulation M3T10R10: ${\beta }_{\mathit{RR}}^{\prime }=-0.02$, ${\beta }_{R\phi }^{\prime }=-0.11$, ${\beta }_{\phi R}^{\prime }=7.5\times {10}^{-3}$, and ${\beta }_{\phi \phi }^{\prime }=-0.1$. 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:

$$\begin{array}{c}\hfill \chi \left(r\right)\equiv {\langle \frac{{\langle \mathit{B}\rangle }_{\phi }·{\langle \mathit{B}\rangle }_{\phi }}{{\langle \mathit{B}·\mathit{B}\rangle }_{\phi }}\rangle }_{\theta }.\end{array}$$(17)

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 large-scale 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 small-scale 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 ${\mathcal{R}}_m\equiv {\mathrm{\Omega }}_{\star }{h}_{\text{init}}^2/\eta$, 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 ${\mathcal{R}}_m=\sqrt{10}\approx 3.2$ 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 large-scale 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 low-resistivity 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 small-scale 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 mean-field 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⁻³ 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 small-scale 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 $\overline{B}$ 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⁻². 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⁻³ 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 $B_r^2$, $B_{\theta }^2$, and $B_{\phi }^2$. The fractions of mean-field 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 best-fit 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 quasi-steady 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⁻³ 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 high-mass 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⁻² 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 mean-field 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²/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² 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:

$$\begin{array}{cc}& {\mathrm{R}}_{R\phi }\equiv \rho (v_R-{\langle v_R\rangle }_{\rho })(v_{\phi }-{\langle v_{\phi }\rangle }_{\rho }),\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill & {\mathrm{M}}_{R\phi }\equiv -B_R{B}_{\phi },\hfill \end{array}$$(19)

$$\begin{array}{cc}\hfill & {\mathrm{G}}_{R\phi }\equiv \frac{\left({\partial }_R{\mathrm{\Phi }}_{\mathrm{disk}}\right)\left({\partial }_{\phi }{\mathrm{\Phi }}_{\mathrm{disk}}\right)}{4\mathrm{\pi }RG}\hfill \end{array}$$(20)

(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⁻² is subdominant inside r/r_in ≲ 12, and we checked that it keeps decaying outward over time in all the dynamo-saturating runs. The magnetic stress corresponds to an effective viscosity α_M ≈ 10⁻¹ 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⁻¹, 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 1136tin to 1200tin in run M3T10R10: total stress (solid black), gravitational stress αG (dashed blue), magnetic stress αM (dot-dashed 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⁻¹. 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⁻² 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

$$\begin{array}{c}\hfill Q_{\rho }\equiv \frac{{\kappa }_{\rho }{c}_{\rho }}{\mathrm{\pi }G\mathrm{\Sigma }}\approx 1\end{array}$$(21)

in nonmagnetized disks, when using the density-weighted epicyclic frequency ${\kappa }_{\rho }\equiv \sqrt{(2{\mathrm{\Omega }}_{\rho }/R)\text{d}\left[R^2{\mathrm{\Omega }}_{\rho }\right]/\text{d}R}$ 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 $\sqrt{1+B^2/P}$ (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 density-weighted 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 ${\langle B_{\phi }\rangle }_{\theta ,\phi }\approx 0.33\overline{B}$ before saturation to $\approx 0.86\overline{B}$ 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 on-and-off 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 long-term 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.

1. During its kinematic phase, the GI dynamo grows a large-scale 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.

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

3. The kinematic dynamo can be described via simple mean-field models upon azimuthal averaging. While its mechanism obeys the classic stretch-twist-fold terminology, its effect on the mean fields formally differs from traditional α-Ω models.

4. 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 long-term 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 GI-induced 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).

Acknowledgments

W.B. acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG) through Grant KL 650/31-1, 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 Baden-Württemberg through bwHPC and the DFG through grant no INST 37/935-1 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

Appendix A: Choice of EMF reconstruction scheme

Integrating the magnetic induction Eq. (4) over a bounded surface 𝒮 and applying Stokes’ theorem gives

$$\begin{array}{c}\hfill {\partial }_t{\displaystyle {\int }_{\mathcal{S}}\mathit{B}·\mathrm{d}\mathit{S}=-{\oint }_{\partial \mathcal{S}}\mathit{E}·\mathrm{d}\mathit{\ell },}\end{array}$$(A.1)

where E is the total EMF as seen in the chosen reference frame. The method of constrained transport implements this equation by taking the surface-integrated magnetic flux as variable to evolve in time (Evans & Hawley 1988). In a finite-volume framework, the different components of the discretized magnetic field therefore live on their respective interfaces and can be called “face-centered”. On the other hand, the circulation of the EMF involves “edge-centered” variables.

There are different ways to reconstruct the edge-centered EMF given the values of the face-centered and cell-centered 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²/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 ⟨B2⟩/⟨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/tin ≈ 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 ${\langle B_{\phi }\rangle }_{\theta ,\phi }/\overline{B}\approx 0$ 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 $B_{\theta }^2$ and ${\partial }_t{B}_{\theta }^2$ 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²/P ≲ 10⁻² 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 boundary-driven amplification, though at its slower pace.

Appendix B: Minimal mean-field instability

Equations (15) & (16) can be further reduced to the minimal form:

$$\begin{array}{cc}\hfill {\partial }_t{B}_R\simeq & -{\partial }_z\left[V_z{B}_R\right]-{\eta }^{\prime }{\partial }_z^2{B}_{\phi },\hfill \end{array}$$(B.1)

$$\begin{array}{cc}\hfill {\partial }_t{B}_{\phi }\simeq & -{\partial }_z\left[V_z{B}_{\phi }\right]+{\eta }^{\prime }{\partial }_z^2{B}_R,\hfill \end{array}$$(B.2)

where, for simplicity, we took the same off-diagonal 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 ω₀ ≈ 5.59. Unstable (real positive) eigenvalues appear when the amplitude of the compression rate |∂_zV_z| becomes comparable to the fundamental frequency ω₀, 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 −∂zVz relative to the fundamental mode frequency ω0 obtained at Vz = 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 R2Σ at time 400tin over the whole domain of run M3T10R10; colors range linearly from zero to three times the median of R2Σ.
In the text

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

	Fig. 3. Growth rate of the magnetic field amplitude $\overline{B}\equiv \sqrt{{\langle B^2\rangle }_{\theta ,\phi }}$ as a function of radius in run M3T10R10, normalized either by the (fixed) inner Keplerian frequency Ωin (solid blue) or by the (radially varying) density-weighted orbital frequency Ωρ (dashed red). Only the inner parts r/rin ≲ 12 have reached a steady exponential growth phase.
In the text

	Fig. 4. Equatorial slice in run M3T10R10 at time 400tin, zoomed in on the radial interval r/rin ∈ [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 (Br < 0) or green when outward (Br > 0).
In the text

Fig. 5. Meridional slice at time 400tin and angle φ = 3π/8 in run M3T10R10. Left panel: density distribution flattened by a factor of (R/rin)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/rin ≈ 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/rin ∈ [2,8] and over time from 200tin to 400tin.
In the text

	Fig. 7. Meridional profiles of the mean flow as defined by Eq. (5), averaged from 200tin to 400tin 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 200tin to 400tin 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 ${\mathrm{\Omega }}_{\rho }\overline{B}$ and averaging azimuthally, radially over r/rin ∈ [2,8], and over time from 200tin to 400tin 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 $B_{\theta }^2$ after normalizing ${\partial }_t{B}_{\theta }^2$ by ${\mathrm{\Omega }}_{\rho }{\overline{B}}^2$ at every radius and averaging azimuthally, radially over r/rin ∈ [2,8], and over time from 200tin to 400tin 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 step-by-step description.
In the text

	Fig. 13. Normalized profiles of the magnetic field amplitude $\overline{B}$ measured before dynamo saturation in four simulations featuring a steady exponential growth (solid lines; see legend), and the best-fit 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 200tin to 400tin 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 $V_z/\overline{V}=-0.11z/h$, satisfying (∂zBR,∂zBφ) = 0 at z/h = 0, BR = 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Ωh2, and the dynamo coefficients estimated in the reference simulation M3T10R10: ${\beta }_{\mathit{RR}}^{\prime }=-0.02$, ${\beta }_{R\phi }^{\prime }=-0.11$, ${\beta }_{\phi R}^{\prime }=7.5\times {10}^{-3}$, and ${\beta }_{\phi \phi }^{\prime }=-0.1$. 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 low-resistivity 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 1136tin to 1200tin in run M3T10R10: total stress (solid black), gravitational stress αG (dashed blue), magnetic stress αM (dot-dashed 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 density-weighted Toomre number during the dynamo saturation in run M3T10R10.
In the text

	Fig. A.1. Evolution of the average magnetization ratio ⟨B2⟩/⟨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/tin ≈ 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 −∂zVz relative to the fundamental mode frequency ω0 obtained at Vz = 0; blue (resp. red) crosses represent the real (resp. imaginary) part of the eigenvalue.
In the text