Helical and nonhelical large-scale dynamos in thin accretion discs

The dynamics of accreting and outgoing flows around compact objects depends crucially on the strengths and configurations of the magnetic fields therein, especially of the large-scale fields that remain coherent beyond turbulence scales. Possible origins of these large-scale magnetic fields include flux advection and disc dynamo actions. However, most numerical simulations have to adopt an initially strong large-scale field rather than allow them to be self-consistently advected or amplified, due to limited computational resources. The situation can be partially cured by using sub-grid models where dynamo actions only reachable at high resolutions are mimicked by artificial terms in low-resolution simulations. In this work, we couple thin-disc models with local shearing-box simulation results to facilitate more realistic sub-grid dynamo implementations. For helical dynamos, detailed spatial profiles of dynamo drivers inferred from local simulations are used, and the nonlinear quenching and saturation is constrained by magnetic helicity evolution. In the inner disc region, saturated fields have dipole configurations and the plasma $\beta$ reaches $\simeq 0.1$ to $100$, with correlation lengths $\simeq h$ in the vertical direction and $\simeq 10h$ in the radial direction, where $h$ is the disc scale height. The dynamo cycle period is $\simeq 40$ orbital time scale, compatible with previous global simulations. Additionally, we explore two dynamo mechanisms which do not require a net kinetic helicity and have only been studied in shearing-box setups. We show that such dynamos are possible in thin accretion discs, but produce field configurations that are incompatible with previous results. We discuss implications for future general-relativistic magnetohydrodynamics simulations.


INTRODUCTION
Accretion discs around compact objects power the most luminous objects in the Universe; yet, a thorough understanding of how the gravitational energy of the in-falling matter is converted into radiation is still lacking.In recent years, a common consensus has been built that magnetic fields play central roles in the dynamics of accretion discs, including initiating turbulent flows, enhancing angular momentum transport, and launching jets (Balbus & Hawley 1991;Narayan et al. 2003;Yuan & Narayan 2014).In cases where the central engine is a black hole, jet launching through the Blandford-Znajek mechanism depends on the net magnetic flux accumulated on the horizon (Tchekhovskoy et al. 2011).
General-relativistic magnetohydrodynamics (GRMHD) simulations have become increasingly feasible to study accretion physics in the vicinity of black holes.Given that it is still costly to resolve turbulent flows to the extent where the magnetic field amplification therein, i.e., the dynamo processes, becomes resolution-independent, it has been a common practice to start simulations with gas already threaded by some large-scale poloidal magnetic fluxes.Yet, simulation results ⋆ Email address for correspondence: hongzhe.zhou@sjtu.edu.cnrely on how such initial magnetic fields are arranged (e.g., Beckwith et al. 2008;Barkov & Baushev 2011;Narayan et al. 2012): single-loop poloidal fields lead to strong net fluxes on the horizon and therefore more powerful jets, and fields with alternating polarity lead to episodic and less powerful jets.It remains an active research topic how such large-scale fields should be properly modeled and how the accretion dynamics is affected.
There are two possible origins for the large-scale fields in accretion discs: (i) the large-scale fields are amplified ab initio in the accretion flow, i.e., through large-scale dynamo (LSD) actions, and (ii) such large-scale fields are advected from the outer disc region.For the latter scenario to work, the diffusion time scale of the large-scale fields has to be longer than the advection time scale.Shearing-box simulations suggest that the turbulent magnetic Prandtl number is of order unity (Guan & Gammie 2009), seemingly implying inefficient advection in thin discs, but the corona field (Beckwith et al. 2009), disc vertical structure (Guilet & Ogilvie 2012, 2013), or disc winds (Cao & Spruit 2013) may help alleviate the problem.For thick discs, advection can be more efficient (Cao 2011;Ressler et al. 2020;Dhang et al. 2023b).
In this work, we focus on the alternative explanation of disc magnetization, i.e., through LSDs.Ab initio amplifica-tion of magnetic fields that are coherent above the turbulence scale have been observed in both local shearing-box (Brandenburg et al. 1995;Bai & Stone 2013) and global Newtonian or post-Newtonian (Hogg & Reynolds 2016;Rodman & Reynolds 2023) simulations, but much less often reported in the costly GRMHD simulations (Liska et al. 2020).
Mean-field electrodynamics is commonly used to understand the dynamics of large-scale magnetic fields in plasmas.In such a framework, turbulent fields are averaged out and their mean dynamics contribute to the evolution of the mean magnetic field (for a comprehensive review, see Brandenburg & Subramanian 2005).The coherent effect on the mean magnetic field from random turbulent fields is manifested as the turbulent electromotive force (EMF), E, in the induction equation.Mathematically, E = ⟨u × b⟩, where ⟨• • •⟩ means an ensemble average and u and b are turbulent velocity and magnetic fields, respectively.A common practice is to expand E in terms of the gradients of the mean magnetic field B, where αij and η ijk are tensorial turbulent transport coefficients.A key question is then how these dynamo coefficients can be written in terms of other mean quantities in the system, such as density stratification and differential rotation.A direct way is to measure these turbulent transport coefficients in numerical simulations, and sophisticated tools have been invented to do the job.The test-field method (Brandenburg et al. 2008;Rheinhardt & Brandenburg 2010;Käpylä et al. 2022) has been applied in local MRI simulations (Gressel 2010), and the projection method (Brandenburg & Sokoloff 2002) can be applied to more complicated geometries, such as in thick-disc simulations (Dhang et al. 2020).
One of the practical purposes of a comprehensive understanding of dynamo processes is for better sub-grid dynamo models.Three-dimensional (3D) simulations of accretion discs are capable of dynamos, but are computationally expensive to achieve a resolution-independent regime.In contrast, two-dimensional (2D) simulations are less costly, but intrinsically cannot host a dynamo (Cowling 1933).In order to mimic dynamo processes that are available only at high resolution, sub-grid dynamo models are used in 2D or lowresolution 3D simulations, where a turbulent EMF is artificially added to the induction equation to amplify magnetic fields.
The implementation and application of helical sub-grid dynamo effects have been explored by several groups (von Rekowski et al. 2003;Bucciantini & Del Zanna 2013;Bugli et al. 2014;Stepanovs et al. 2014;Sądowski et al. 2015;Fendt & Gaßmann 2018;Dyda et al. 2018;Tomei et al. 2020;Vourellis & Fendt 2021).However, previous sub-grid models have not exploited the results from analytical theories and local shearing-box simulations: a shortcoming of these works is that both the nonlinear feedback from the magnetic field to the αij coefficient (i.e., dynamo quenching) and the saturation strength employ parameterised forms.In the present work, we improve the previous helical dynamo models in the following ways: (i) the spatial profile of αij is derived based on previous analytical and test-field simulations results, and (ii) the quenching prescription of αij and saturation level of the large-scale field are determined by magnetic helicity evolution.As we shall see, the new model yields dynamo waves and cycle periods in agreement with both local and global direct numerical simulations (DNS) of accretion-disc dynamos.
In addition to the helical dynamo effect, a new class of shear dynamos has been noticed in shearing-box simulations (Brandenburg 2005;Yousef et al. 2008), where nonhelical turbulence and a background shear flow are the only ingredients necessary for the growth of a large-scale magnetic field.Due to such rather simple requirement, it is expected that shear dynamos may widely operate in differentially rotating flows, such as accretion discs.How shear dynamos should be modeled is still under debate (Zhou & Blackman 2021), and two theoretical models have been proposed: (i) the shear-current effect (SCE; Rogachevskii & Kleeorin 2003;Squire & Bhattacharjee 2015), and (ii) the incoherent-α effect (IAE; Vishniac & Brandenburg 1997;Sridhar & Singh 2014;Jingade et al. 2018).Both theories have only been demonstrated in shearing-box setups.In this work, we implement these nonhelical dynamo mechanisms in sub-grid models, and investigate their dynamo capabilities in the disk geometry.
The rest of this work is organized as follows.In section 2 we introduce the axisymmetric model which couples a thin disc with helical or nonhelical large-scale dynamos.In section 3 we introduce the simulation setup that numerically solves the derived models, and the results are presented and discussed in section 4. We conclude in section 5.

AXISYMMETRIC DISC-DYNAMO MODEL
In this section we introduce the mean-field model that applies to a geometrically thin, optically thick disc.We consider a fully turbulent disc which results from the magnetorotational instability (MRI; Velikhov 1959;Chandrasekhar 1961;Balbus & Hawley 1991), whose turbulent time scale is comparable to the local Keplerian time scale and hence much shorter than that of the LSD.Focusing on the LSD processes, we neglect the back-reaction of the LSD-generated magnetic field onto the turbulence and disc structure in this work.In section 5 we discuss the possible outcomes of such feedbacks and the possibility of a more complete disc-dynamo model.
Throughout this work, the physical quantities mentioned (e.g., gas density, pressure, magnetic fields, etc.) all refer to the large-scale ones, i.e., those survive after averaging out the turbulent fluctuations, unless otherwise specified.We use (r, θ, ϕ) for spherical coordinates, and (ϖ, ϕ, z) for cylindrical coordinates.

The thin accretion disc
We consider a Shakura-Sunyaev type disc filled with polytropic gas, where P is the gas pressure, ρ is the density, γ is the ratio of specific heats, and a subscript 0 indicates quantities measured at the disc mid-plane at r = r0.The thermal sound speed is and hence P0 = ρ0c 2 s0 .The vertical balance between the gravitational force and the gas pressure gradient near the disc mid-plane gives where Ω is the Keplerian rotation rate.We denote ϵ = h/ϖ as the dimensionless disc scale height, and define |θ − π/2| ≥ 2ϵ as the corona region.h is also taken to be the typical scale of the growing mode of the LSDs, and the corresponding wave number is k1 = h −1 .The turbulence time scale is roughly τ = Ω −1 .We denote the turbulence length scale by l, which also gives the turbulence velocity scale u = l/τ = lΩ.The typical wave number of turbulent fluctuations is then k2 = l −1 .
The turbulent viscosity is where αSS is the Shakura-Sunyaev parameter.The turbulent viscosity should be of order ν ∼ ul = l 2 Ω.Combining with equations ( 4) and ( 5), we have the estimates (Blackman 1998) and which are useful when deriving dynamo coefficients in terms of cs and Ω.
We ignore the feedback from large-scale fields on the turbulence, and use αSS = 0.3 (King et al. 2007) which remains constant and uniform.

Mean-field dynamo terms
The evolution equation for the large-scale magnetic field can be derived by averaging the induction equation over a suitable scale that is much larger than the turbulence scale.The meanfield equation is (for a derivation see, e.g., Brandenburg & Subramanian 2005) where B is the large-scale magnetic field, U is the mean flow including the accretion inflow and possible large-scale outflows like disc winds, but excluding the Keplerian motion, and E is the turbulent EMF.The last term in equation ( 8) generates toroidal fields from the poloidal ones through differential rotation (the Ω effect), and hence for a successful dynamo, the turbulent EMF must contain terms that generates poloidal fields from the toroidal ones.In this work we consider the turbulent EMFs of the form where αij and ηij are turbulent transport coefficients, and J = ∇ × B/µ0 is the mean current density.
The diagonal components of ηij contribute to the turbulent diffusion of large-scale magnetic fields.We consider isotropic diffusion, and the diagonal components of ηij are Here the first factor, αSScsh, is the same as the turbulent viscosity, and the second factor in the square brackets phenomenologically captures the vertical structure of η, inspired by test-field results in stratified shearing-box simulations (Gressel 2010).The factor of 3/8 ensures that the average of η over one disc scale height is αSScsh, and hence we consider a turbulent magnetic Prandtl number of order unity.
The coupling of large-scale fields and turbulence will lead to modifications of η, which we do not include in this work.Some dynamo models using parametrized forms of η include Zanni et al. (2007) and Stepanovs & Fendt (2014).
The remaining components in αij and ηij vary in different models and are derived in the next few subsections.In the follows, we describe in detail three dynamo models: (i) the αΩ dynamo, (ii) the incoherent α dynamo, and (iii) the shearcurrent dynamo.The αΩ model is the conventional helical dynamo that has been widely applied to explain stellar and galactic magnetization, whereas the latter two need no net kinetic helicity and are less explored, particularly in a thin-disc geometry.For each, we derive explicit forms of the turbulent EMF ( 9) with physics-or simulation-motivated quenching formulae, expressed in terms of cs and Ω that can be readily implemented in an accretion disc context.A summary of dynamo models is in table 1, and numerical solutions of the models are given in the next section.

The αΩ model
We consider a minimalist helical dynamo model where the offdiagonal components of ηij vanish, and hence the turbulent EMF takes the form From stratified shearing-box simulations, all components of αij have roughly the same order of magnitude (Gressel 2010), but the α ϕϕ component is the most relevant as it is responsible for generating poloidal magnetic fields from toroidal ones.Hence we only consider α ϕϕ in what follows, making the model an αΩ-type dynamo.
It has been long-debated the profile of α ϕϕ in accretion discs, especially about how its sign changes with latitude.In hydrodynamical cases, stratified rotating turbulence possesses a net kinetic helicity in each hemisphere, making the kinetic contribution to α positive in the upper plane and negative in the lower plane.The picture changes in the presence of MRI, where local shearing-box (Gressel 2010;Dhang et al. 2023a) and global (Hogg & Reynolds 2018;Dhang et al. 2020) simulations report that α ϕϕ has opposite signs to its hydrodynamic values within one to two disc scale heights.This has been attributed to the combined effects of differential rotation and magnetic buoyancy (Brandenburg 1998), or helicity flux (Gopalakrishnan & Subramanian 2023).Within one disc scale height, α ϕϕ has the expected signs as in the hydrodynamic case, i.e., positive in the upper plane and negative in the lower plane.
We now construct an α ϕϕ profile based on the results of The initial vertical profile of the dynamo driver α and turbulent diffusion η coefficients.
Near the disc mid-plane, the sign reversal of α ϕϕ is captured in an ad hoc way through an additional z dependence to equation ( 12): Note that cs has both radial (power-law) and vertical (exponential) dependence, which self-consistently gives the full profile of α (0) ϕϕ .The initial vertical profile of α (0) ϕϕ and η are plotted in figure 1, to be compared with the test-field results of Gressel (2010).
ϕϕ drives a helical dynamo that feedbacks onto the turbulent EMF, and is nonlinearly quenched.In previous GRMHD simulations with sub-grid dynamo terms, it is typical to use the parametrized form where Bst is some target saturation value.In the work by Mattia & Fendt (2022), the effects of different quenching formulae on jet launching and production of flares in the inner disc region has been explored.In this work, we consider a dynamical quenching model constrained by magnetic helicity conservation, similar to those applied to stellar and galactic dynamo models.
In the dynamical quenching formalism (Blackman & Field 2002), αij receives a back-reacting magnetic contribution α m ij whose trace is τ ⟨j • b⟩ /µ0ρ, i.e., proportional to the smallscale current helicity density ⟨j • b⟩ where j is the smallscale current helicity.We assume that α m ϕϕ has the same zdependence as that of α (0) ϕϕ , and the full form of α ϕϕ is taken to be where α m = τ |⟨j • b⟩| /3µ0ρ.
To connect α m to the large-scale fields, we consider the small-scale magnetic helicity associated with the current helicity, where a is the small-scale vector potential.If no helicity flux is present, conservation of magnetic helicity is hold locally at scale ≃ h, assuming that the dynamo has started with negligible total magnetic helicity.Here A is the large-scale vector potential.The right side of equation ( 17) becomes nonzero if (i) the helicity density is advected by a mean flow, or (ii) a Vishniac-type flux arises due to the inhomogeneities of smallscale fields.We now consider such two possibilities.Large-scale outflows like disc winds preferentially carry away large-scale magnetic fields since small-scale fields cannot survive the shredding by the shear flow when they buoyantly arise (Blackman & Pessah 2009).This leads to a deficit on the right side of equation ( 17), which we parametrize as a certain fraction fw of the large-scale magnetic helicity, If Uw is the typical wind speed leaving the disc and Ω −1 is the typical time scale for helicity loss, we have the rough estimate that i.e., fw is roughly equal to the Mach number of the disc wind.
In general Uw/cs varies radially, and more importantly, depends on the local magnetic field strength.For strong poloidal magnetic fields whose Alfvén speed is larger than the local sound speed, disc winds are magneto-centrifugally driven, whereas for weaker poloidal magnetic fields, disc winds are driven by the magnetic pressure of the toroidal fields.Focusing on the dynamo processes within the disc, we adopt constant values of fw in this work.Helicity flux may also arise from triple correlations of smallscale fields, as considered by Gopalakrishnan & Subramanian (2023).The derived flux depends on the gradients of u 2 and b 2 along the direction of the mean vorticity, i.e., in the z direction.Using equation ( 18) in Gopalakrishnan & Subramanian (2023), we find that the contribution from such terms is ∆ ⟨a • b⟩ ≃ αSSh b 2 (0.08ζ − 0.04) in the upper disc plane during one Keplerian time scale.For later convenience we also write it as a fraction of the large-scale magnetic helicity, where we have used b 2 /B 2 = ( b 2 /ρu 2 )(ρu 2 /B 2 ) = αSSζβL/2, and is the ratio between the thermal pressure and the magnetic pressure from the large-scale fields.Note that equation ( 20) gives a positive f f and hence a positive contribution to ⟨a • b⟩, in contrast to the thin-disc estimation of Gopalakrishnan & Subramanian (2023), due to the difference between the models.
Collecting all the terms, we have Combining both kinetic and magnetic contributions, we have where Note that equation ( 24) is only valid within the disc where the gas is sufficiently turbulent.In the numerical solutions, we do not allow α(βL) to become negative, i.e., α ϕϕ does not pass zero.Changes of signs for α ϕϕ mostly happen when a helical magnetic loop buoyantly rises into the corona, and then "untwists" due to the inverse transfer of magnetic helicity.
Given the full expression of α ϕϕ , the dynamo saturation level can be estimated from the dynamo number, i.e., the ratio between the products of the constructive and the destructive coefficients in the dynamo equation.As order-ofmagnitude estimates, we simply use α ϕϕ ≃ α(βL), ignoring its z-dependence, and f = 1 for the moment.The dynamo number for this αΩ dynamo is At dynamo saturation we have D(βL,st) = 1, yielding Equation ( 27) gives βL,st ≃ 74 if αSS = 0.1, and βL,st ≃ 11 if αSS = 0.3.These values are only ∼ 2 times larger than those derived from the exact dispersion relation.
A few remarks are made regarding connections to previous works: (i) Pudritz (1981) derived a critical value βL,P81 = 4/3cuα (where cu ≲ 1/2) based on energy balance between the magnetic and kinetic turbulent energy.Our estimation from the mean-field dynamo theory, equation ( 27), is ∼ 3 times larger than βL,P81 if αSS = 0.1.However, as we shall see in section 4, the actual saturation value of βL is much larger than βL,st, possibly due to the more complicated profile of α ϕϕ in the disc than that in the simplified scenario considered above.(ii) The energy equipartition between the large-scale field and the turbulent kinetic energy gives the saturated magnetic field is a few times smaller than its energy-equipartition value.(iii) The quenching formula (24) combined with equation ( 27) yields a faster quenching than the more often used parametrized quenching ( 14).For equation (24), the quenched value of α at dynamo saturation is times its kinematic value, whereas the typically used form gives i.e., only 50% quenching when βL = βL,st.Hence the latter formula will lead to a lower saturated value of βL than βL,st.

The incoherent α model
In this and the next subsections, we introduce the two nonhelical large-scale dynamo mechanisms that originate from forced shearing-box simulations (Brandenburg 2005;Yousef et al. 2008), and couple them to the thin-disc model.In contrast to the conventional α mechanism, these two mechanisms do not require a net kinetic helicity in the flow, but only the combination of turbulence and a shearing flow.
The IAE relies on the fluctuations of the kinetic helicity and therefore a fluctuating α term with zero mean.The possibility of dynamo driven by fluctuating helicity was first raise by Kraichnan (1976) and Moffatt (1978) in non-shearing turbulence, and numerically by Vishniac & Brandenburg (1997).More recently, Sridhar & Singh (2014) and Jingade et al. (2018) showed that in the presence of shear, such helicity fluctuations can be recast into an effective anisotropic turbulent magnetic diffusivity tensor which is capable of a dynamo action.
In dynamos driven by IAE, the turbulent EMF takes the form where the pseudo-scalar α is inhomogeneous in space and fluctuates in time.We denote the values of α in the unquenched regime by α0.
The α fluctuations are assumed to happen on meso-scales, i.e., at spatio-temporal scales larger than those of the turbulent flow but still smaller than the mean-field scale.We parameterise such scales using lα = ml for the length scale and τα = mΩ −1 for the time scale, where m > 1 is a model parameter.During the kinematic regime, m is solely a property of the turbulent flow, but unfortunately it is so far unclear how its value can be calculated theoretically.An attempt to determine it in simulations is on-going (Zhou and Jingade, in prep.), and preliminary results indicate that m is ≲ O(5) and scales weakly with the shear rate S as (Sτ ) 0.3 .In our thindisc model, since S = 3Ω/2 and hence Sτ = 3/2 is a constant, we shall take a constant and uniform m throughout the disc.
Another necessary ingredient of the model is what probability distribution function (PDF) the α fluctuations follow.Jingade & Singh (2021) have used a uniform PDF in their renovating-flow model, where α0 ∈ [−αmax, αmax] and αmax = u/3 is the fully helical value.Under such conditions they found that m ≥ 3 is needed for a large-scale dynamo.However, the actual PDF of α fluctuations in realistic turbulence is not clear.In forced shearing turbulence, Brandenburg et al. ( 2008) found α0 to follow a Gaussian distribution, α0 ∼ N (0, 0.2ηk2) = N (0, 0.2α 1/2 SS cs), which has a variance 25 times smaller than that from a uniform PDF, and hence less capable of driving a dynamo.On the other hand, the PDF of α fluctuations in MRI turbulence has not been determined in simulations yet.In this work, we explore both uniform and Gaussian PDFs for α0 to check their dynamo capability in a disc geometry.
Following Kraichnan (1976), the effective EMF for the incoherent α dynamo is where σ 2 = α 2 /αSSc 2 s is the normalized variance of the α fluctuations determined by the PDF.The dynamo growth rate in the kinematic regime is with σ 2 0 = α 2 0 /αSSc 2 s .The nonlinear saturation of the incoherent α dynamo is so far unclear.Numerical simulations of unstratified or stratified MRI turbulence suggest that the large-scale magnetic field strength typically saturate at O(1) times the energyequipartition value (Brandenburg et al. 1995;Davis et al. 2010;Shi et al. 2016;Guilet et al. 2022), βL,eq = 2/αSS [see equation ( 28)].To achieve this, we may use a parameterised quenching formula and correspondingly σ = σ0/(1 + qβL,eq/βL).The factor q is to be determined by the condition γIA(βL = βL,eq) = 0, giving q = √ mσ0 − 1. Collecting all the terms we have To implement meso-scale fluctuations of α as a sub-grid model, we first divide the simulation domain into cells of scale lα = ml = mα 1/2 SS h, and within each cell the value of α0 is uniform but randomly drawn from a given PDF at time intervals roughly equal to τα = mΩ −11 .An illustration of the coherent cells of α is in figure 2.
Using the specific forms equation ( 38), the normalized dynamo growth rate is found to be We now consider a quenching prescription for the shearcurrent effect.Similar to the incoherent-α model, we propose the parametrized form η rϕ,ϕr = η rϕ,ϕr,0 1 + qβL,eq/βL , and q is to be determined by requiring that the dynamo growth rate vanishes when βL = βL,eq.At αSS = 0.3, this gives approximately Hence an ad hoc quenching prescription can be

NUMERICAL SETUP
In this section we introduce the setup in the publicly available Pencil Code (Pencil Code Collaboration et al. 2021) to find numerical solutions of the dynamo models derived in the last section.The physical quantities are independent of the azimuthal angle, but vectors still have finite ϕ components, i.e., the simulations are 2.5-dimensional.We do not resolve the turbulent-scale fields, but only incorporate them as sub-grid physics including turbulent diffusion of mass, velocity and magnetic fields, and the dynamo coefficients.In this sense, the simulations shall be understood as numerical solutions of mean-field models, rather than DNS of accretion discs.
We consider azimuthally averaged MHD equations with a central gravitational source.Furthermore, since viscous and resistive dissipations are just energy transfers between the large-and small-scale fields, we assume such processes do not change the mean-field entropy and thereby consider a Here χ is the turbulent mass diffusion coefficient taken to be equal to ν, and Sij = (∂iUj + ∂jUi)/2 − δij∇ • U /3 is the rate-of-strain tensor.
To ensure that shocks are properly resolved, we add a shock-capturing diffusion term on the right-hand side of equation (45), which is a bulk viscosity of the form Here max5 means the maximum value in the neighboring five mesh points, and ∆r and ∆θ are the mesh sizes in the radial and polar directions, respectively.Hence the shock diffusion is the strongest for the strongly converging flows.The simulations are carried out in spherical coordinates, with the simulation domain being r ∈ [r0, 100r0] and θ ∈ [0, π].We use a resolution of Nr × N θ = 384 × 256 and static mesh refinement in the radial direction, i.e., when moving towards the disc inner region, we double the resolution inside 16r0, and further double it inside 8r0.
A summary of the runs is in table 2. For all the runs, the magnetic field is only initiated after the accretion flow has reached a steady state at t = 4000Ω −1 0 .

Boundary conditions
The boundary conditions mostly follow those of Mattia & Fendt (2022), with a few adaptions for the Pencil Code since the latter uses a finite-difference scheme and evolves the vector potential rather than the magnetic field.
In the radial direction, Ur is extrapolated in a power law with index −3/2 at the inner boundary, and uses outflow boundary condition at the outer boundary.U θ uses free boundary condition with vanishing second-order derivative at the inner boundary, and is ∝ r 0 at the outer boundary.U ϕ is kept at 98% times the local Keplerian value at both radial boundaries.For the density field and all the three components of the magnetic vector potential, free boundary condition is used for both boundaries.
In the polar direction, the usual symmetric or antisymmetric boundary condition is used.

Initial conditions
The initial density follows where pρ = −15/8 for a standard Shakura-Sunyaev thin disc.
We also impose a density floor where ρ floor,0 = 10 −3 ρ0.Within the disc, the initial azimuthal velocity is 98% of the Keplerian value, while in the corona the gas has zero initial velocity.

Model parameters
For all the simulations in this work, we use cs0 = 0.1 and γ = 1.4 for a cold disc, and GM = ρ0 = 1.
The Shakura-Sunyaev prescription gives ν = αSScsh ∝ ϖ pρ(γ−1)+3/2 = ϖ −3/4 .However, a steeper profile, ν ∝ ϖ −1 , is actually used to ensure a stable accretion flow.Hence the turbulent viscosity is kept fixed in space and time and given by ν = αSScs0h0 ϖ ϖ0 On the other hand, the turbulent resistivity is determined by the local gas properties through In the corona region where |θ − π/2| ≥ 2ϵ, we use fixed values η = ν to prevent magnetic field lines from accumulating and eventually crashing the simulation.For all the dynamo models, we restrict the dynamo-driving coefficients (i.e., all the components of αij and the offdiagonal components of ηij) inside the disc by applying to them an additional spatial profile which rapidly damps the dynamo drivers beyond two disc scale heights.

RESULTS
In this section we report the numerical solutions for each dynamo models with different choices of parameters.We first give a cross comparison between the fiducial runs for each dynamo models, and then explore the αΩ and the IAE models in more detail.A summary of the resulting magnetic field configurations is listed in the last column of table 2. In figure 4, we show the space-time diagram of the αΩ model in its steady state.The azimuthal component B ϕ reproduces the butterfly diagrams seen in previous local (Gressel 2010) and global (Hogg & Reynolds 2018) simulations, where the dynamo waves migrates away from the disc midplane.However, the radial component Br in the high latitude region, z/h ≥ 2, is rather weak and does not follow the pattern seen in Gressel (2010).This might result from that we have restrict the dynamo terms inside the disc.A snapshot of the field configuration is shown in figure 5, where the coherent scale of B ϕ is ≃ h in the vertical direction, and ≃ 10h in the radial direction.The alternating field polarity is in general agreement with global simulations of disc dynamos by Hogg & Reynolds (2016).
Overall, the amplified large-scale magnetic field resembles a dipolar configuration with a cycle period of ∼ 40 orbital time scale at r0, in agreement with previous global simulation results (Hogg & Reynolds 2018).Note that this is a few times longer than the estimate from local periodic-box dispersion relations (Gressel & Pessah 2015) and also previous sub-grid GRMHD simulations without α quenching (Bugli et al. 2014), which is about ∼ 10 times the orbital time scale.
In figures 6 and 7, we show the space-time diagrams and a snapshot of field configuration for the fiducial incoherentα model, run IA1.Although the magnetic field is somewhat more coherent in space than that in run AO1, the incoherentα model has much more temporal fluctuations.In contrast, figures 8 and 9 present the results for the shear-current model run SC1, and show the most coherent field structure in both space and time.For both of the nonhelical models, they do not present regular dynamo cycles, in contrast to global simulations.We hence conclude that even though they may co-  exist with the helical α effect in discs, they are not likely the dominant driver.

Dynamical versus parametrized quenching for α
We now compare the fiducial αΩ model with run AO2 which employs the often used parametrized quenching prescription, equation ( 14).Both runs have no helicity flux, i.e., f = 1.
In figure 10 we compare the disc magnetization for the two models.As expected, the two models share similar kinematic growth rates, but run AO2 produces slightly stronger magnetization, particularly at the dynamo minima.In figure 11 we show the space-time diagrams for run AO2.The cycle period is about 60 times 2π/Ω0, longer than the dynamically quenched case.
Most surprisingly, the parametrized quenching formula produces a quadruple rather than dipole configuration.In shearing-box simulations, Brandenburg (1998)   that using vertical or perfect-conductor boundary conditions for the magnetic field in the vertical direction determines the parity and cycle period of the dominant large-scale mode in the saturated state.In our cases, we do not have closed boundaries in the vertical direction, but it is still possible that the quenching formula yields an effective boundary condition at the disc-corona surface.

Effects of helicity fluxes
The effects of helicity fluxes are explored in runs AO3 and AO4, both having the helicity flux term but the former uses fw = 0 while the latter uses fw = 1.The space-time diagrams for the   the flux terms (either Vishniac flux or wind loss) do not have a noticeable impact on the cycle period, but only slightly change dynamo wave patterns and a decrease in the field strength (a factor of ≲ 2.5).

Efficiency of the incoherent-α model
Finally we analyze the dynamo capability of the IAE when using different models of α fluctuations.Both the coherent scale and the PDF of α fluctuations turn out to be critical.More specifically, we explore the results with (i) different values of m, and (ii) different PDFs.Future MRI simulations are needed to determined what model (i.e.coherent scale and PDF) best fit realistic turbulence.
In figure 13 we plot the time evolution of the mean βL in the  2018), who concluded that m ≥ 3 is necessary for a dynamo in their renovating-flow models in a periodic box, rather than in a thin disc.
A uniform PDF produces strong and weak α fluctuations with equal probability, but it is perhaps more plausible that weak to medium strength fluctuations appear more often in realistic turbulence.In forced shearing-box turbulence, Brandenburg et al. (2008) found the α fluctuations to follow a Gaussian distribution ∼ N (0, 0.2αmax).However, in such a case we observe no dynamo, i.e., the initial field decays and reaches a steady state where βL ≳ 10 9 .

CONCLUSION AND DISCUSSIONS
Implementing sub-grid dynamo terms in low-to mediumresolution GRMHD simulations allow for affordable yet selfconsistent amplification of large-scale magnetic fields.In this work, we formulated both helical and nonhelical dynamo terms in a Shakura-Sunyaev disc, and the derived mean-field models are numerically solved.
For the helical model, we considered several new effects compared to previous implementations: (i) an αij profile inspired by local MRI simulations, (ii) the dynamical quenching prescription which takes into account magnetic helicity evolution, and (iii) helicity flux terms.We demonstrate that the resulting large-scale field has a dipole confiuration, and reproduces the dynamo waves and cycle periods in previous local or global DNS simulations.
For the nonhelical models, we show that they are possible in thin-disc geometries when using optimistic parameters.However, the resulting dynamo patterns are not consistent with DNS simulations, either with the wrong parity, or with too much fluctuations.We hence conclude that such dynamos may co-exist with the helical mechanism in accretion discs, but are not the dominant ones.
We do not observe outflows in our simulations, which could be due to (i) a cold corona which renders the critical radius of Parker wind out of the simulation domain, rcrit = GM/2c 2 s ≃ 10 5 , (ii) the artificially large magnetic diffusivity in the corona which diffuses out the magnetic fields too quickly before they could relax to meet the Blandford-Payne condition, and (iii) the lacking of general-relativistic boundary conditions at the inner disc boundary.A realistic corona is likely heated by the reconnection of magnetic fields, and hence the heating rate depends on βL of the coronal fields.We leave such models in future work.
Disc magnetic fields with alternating polarity may lead to magnetic reconnection near the horizon with the formation of plasmoids and sometimes episodic jets (Parfrey et al. 2015;Mahlmann et al. 2020;Nathanail et al. 2020;Chashkina et al. 2021;Nathanail et al. 2022;Jiang et al. 2023).Using sub-grid models, Stepanovs et al. (2014) demonstrated such possibilities by switching on and off the dynamo terms by hand.Fendt & Gaßmann (2018) also observed pulsed ejections that are caused by time-varying toroidal fields.They have noticed that the pulsing time scale is 3−4 orders of magnitude shorter than those in protostellar jets, but the dynamo cycle period is also rather short, comparable to the local Keplerian time scale.In our improved model, dynamo cycle periods that are ∼ O(50) Keplerian time scale can be self-consistently produced, and their consequences on jet launching shall be examined.
We have taken αSS as a free parameter throughout, but in a more realistic model it should depend on βL.Analytical calculations (Rüdiger & Pipin 2000) and numerical simulations (Blackman et al. 2008) found αSS ∝ β −1 L when the saturated value of βL is used, and αSS ∝ β −1/2 L when the initial value of βL is used (Begelman & Armitage 2023).In any case, αSS should increase with decreasing βL, because a stronger magnetic field taps energy into larger turbulent eddies, and causes stronger turbulent diffusion.In such cases, some interesting mean-field dynamo features could be (i) a longer kinematic phase because the initial weak field does not produce energetic turbulence and therefore dynamo effects are weak, and (ii) a weaker quenching effect, because a growing large-scale field leads to stronger turbulence which in turn strengthens dynamo effects, i.e., against dynamo quenching.The latter effect coupled with accretion dynamics could possibly lead to non-trivial dynamo cycles.In fact, unifying large-scale dynamo and accretion theories in a mean-field approach has been suggested and investigated by several authors (Pessah et al. 2006;Blackman 2012;Mondal & Bhat 2023).
We have considered dynamo models in fully ionized discs where MRI is responsible for the turbulence.For weakly ionized flows such as young circumstellar discs, different dynamo models are possible (Béthune & Latter 2022).
Finally, we take the chance to highlight that the large-scale dynamo drivers, including the α terms and the shear-current effect, are all mean-field effects, and should only vary on temporal and spatial scales that are larger than the local turbulence scales.In the current study, we have carried out meanfield simulations which do not resolve the turbulence, and the dynamo coefficients are related to the mean gas density and the mean sound speed.We propose that a similar implementation should be carried out in GRMHD simulations with sub-grid physics, where the dynamo coefficients should only depend on ϕ-and time-averaged quantities, rather than the instantaneous values which necessarily include turbulent fluctuations.

Figure 2 .
Figure 2.An example of the constructed meso-scale cells (m = 3) of α fluctuations assigned with random values.The thick red curve indicates the disc-corona boundary at π/2 − θ = 2ϵ.

Figure 3 .
Figure 3. Maximal and minimal values of β L in the inner disc region for the fiducial runs.

Figure 4 .
Figure 4. Space-time diagram for the three components of the magnetic field at r = r 0 for the fiducial αΩ dynamo, run AO1.

Figure 5 .
Figure 5.The snapshot at t = 10550Ω −1 0 for the poloidal fields (solid curves with arrows) and the toroidal field (color), for run AO1.

Figure 6 .
Figure 6.Space-time diagram for the three components of the magnetic field at r = r 0 for the fiducial IAE dynamo, run IA1.

Figure 7 .
Figure 7.The snapshot at t = 10550Ω −1 0 for the poloidal fields (solid curves with arrows) and the toroidal field (color), for run IA1.

Figure 8 .
Figure 8. Space-time diagram for the three components of the magnetic field at r = r 0 for the fiducial SCE dynamo, run SC1.

Figure 9 .
Figure 9.The snapshot at t = 10550Ω −1 0 for the poloidal fields (solid curves with arrows) and the toroidal field (color), for run SC1.

Figure 10 .Figure 11 .Figure 12 .
Figure 10.Average values of β L in the disc region for the αΩ models with different quenching prescriptions, with α SS = 0.3.

Figure 13 .
Figure 13.Varying m = lα/l in run IA1.(a) Time evolution of the mean β L in the inner disc.(b) The kinematic growth rates with varying m.

Table 1 .
Summary of dynamo models.

Table 2 .
Summary of the simulations.The prefix in the run names AO refers to αΩ models, IA refers to incoherent α models, and SC refers to shear-current models.The last column labels the resulting dynamo modes after saturation: D for a dipole mode, Q for a quadruple mode, and F for strongly fluctuating polarization.