Shedding light on the MRI driven dynamo in a stratified shearing box

We study the magneto-rotational instability (MRI) driven dynamo in a geometrically thin disc ($H/R\ll 1$) using stratified zero net flux (ZNF) shearing box simulations. We find that mean fields and EMFs oscillate with a primary frequency $f_{\rm dyn} = 0.017$ ($\approx 9$ orbital period), but also have higher harmonics at $3f_{\rm dyn}$. Correspondingly, the current helicity, has two frequencies $2f_{\rm dyn}$ and $4f_{\rm dyn}$ respectively, which appear to be the beat frequencies of mean fields and EMFs as expected from the magnetic helicity density evolution equation. Further, we adopt a novel inversion algorithm called the `Iterative Removal Of Sources' (IROS), to extract the turbulent dynamo coefficients in the mean-field closure using the mean magnetic fields and EMFs obtained from the shearing box simulation. We show that an $\alpha-$effect ($\alpha_{yy}$) is predominantly responsible for the creation of the poloidal field from the toroidal field, while shear generates back a toroidal field from the poloidal field; indicating that an $\alpha-\Omega$-type dynamo is operative in MRI-driven accretion discs. We also find that both strong outflow ($\bar{v}_z$) and turbulent pumping ($\gamma_z$ ) transport mean fields away from the mid-plane. Instead of turbulent diffusivity, they are the principal sink terms in the mean magnetic energy evolution equation. We find encouraging evidence that a generative helicity flux is responsible for the effective $\alpha$-effect. Finally, we point out potential limitations of horizontal ($x-y$) averaging in defining the `mean' on the extraction of dynamo coefficients and their physical interpretations.


INTRODUCTION
The problem of angular momentum transport is a key concept in a rotationally supported accretion disc (for a review, see Balbus & Hawley 1998).The current consensus is that a weak magnetic field instability, namely magneto-rotational instability (MRI; Velikhov 1959;Chandrasekhar 1960;Balbus & Hawley 1991;Balbus & Hawley 1992) is responsible for outward angular momentum transport and drives mass accretion in a sufficiently ionized accretion disc (e.g. as in X-ray binaries, inner part of AGN discs).Although linear MRI ensures outward angular momentum transport, it must be studied in the non-linear phase to account for different observable phenomena in accretion discs.
MRI in an accretion disc is either studied in a local set-up (shearing box; Balbus & Hawley 1992;Brandenburg et al. 1995;Hawley et al. 1995;Davis et al. 2010;Shi et al. 2010;Bodo et al. 2014;Bhat et al. 2016) or in a global simulation (Stone et al. 1999;Hawley 2001;Hawley et al. 2013;Beckwith et al. 2011;Parkin & Bicknell 2013;Hogg & Reynolds 2016;Dhang & Sharma 2019;Dhang et al. 2023).While a global approach is more desirable, it is computationally expensive.On the other hand, the shearing box approach offers an alternate path which is computationally less costly and can provide deep insights into the local processes in MRI-driven turbulence.
In the shearing-box approach (Goldreich & Lynden-Bell 1965), ⋆ E-mail:prasundhang@gmail.com we expand fluid equations to the lowest order of H/R, where H is the density scale height and R is the local radius.Therefore, this approach is valid only for geometrically thin discs with H/R ≪ 1. Depending on whether the vertical component of gravity (gz = −Ω 2 z) (producing a vertically stratified gas density) is considered in the momentum equation or not, shearing box simulations are of two types: stratified (gz ̸ = 0) and unstratified (gz = 0).Further, depending on whether the computational domain can contain net vertical magnetic flux, shearing box models can be classified into zero net flux (ZNF) and net flux (NF) models.Therefore, four possible combinations of the shearing-box model are: i) unstratified ZNF, ii) unstratified NF, iii) stratified ZNF and iv) stratified NF.This work considers a stratified ZNF shearing box model to explore the MRI dynamo in saturation.
Shearing box simulations provide a wide range of behaviour (e.g., convergence, turbulence characteristics etc.) depending on the shearing box model used (for details, we refer to readers to see Table 1 in Ryan et al. (2017)).However, it is to be noted that we will restrict our discussion to the isothermal (i.e.sound speed is constant) models where there is no explicit dissipation and the numerical algorithms provide the dissipation through truncation error at the grid scale.In the presence of an NF, unstratified shearing box simulations show a convergence (in terms of accretion stresses) and sustained turbulence (Hawley et al. 1995;Guan et al. 2009;Simon et al. 2009).On the other hand, stratified NF simulations present different accretion stresses depending on the net flux strength and sustained turbulence (Guan & Gammie 2011;Bai & Stone 2013).Unstratified ZNF models showed intriguing behaviour.Earlier isothermal unstratified ZNF studies (Fromang & Papaloizou 2007;Pessah et al. 2007) found decreased accretion stress and turbulence with increased resolution, implying non-convergence.However, later Shi et al. (2016) recovered convergence using a box with a larger vertical extent than the radial extent.In contrast, earlier stratified ZNF models (Davis et al. 2010) suggested that the models are converged till the resolution 128/H; however, recent studies (Bodo et al. 2011;Ryan et al. 2017) found the model loses convergent properties at higher resolution.
The convergence problem is closely related to the magnetic energy generation process in the MRI-driven flow.For the ZNF (absence of net flux) models, an MRI-driven dynamo must act to overcome the diffusion and sustain the zero net flux in the accretion flow.Earlier ZNF simulations in unstratified (Shi et al. 2016) and stratified (Davis et al. 2010;Bodo et al. 2014;Ryan et al. 2017) shearing boxes found MRI turbulence can self-generate large-scale magnetic fields attaining quasi-stationarity and sustaining turbulence.Riols et al. (2013) suggested that the non-linear MRI does not behave like a linear instability; rather, it provides a pathway for saturation via a subcritical dynamo process.This leads to the question of what kind of dynamo can be sustained in the MRI-driven accretion flow, small-scale or large-scale?The lack of convergence in ZNF models was attributed to the low numerical Prandtl number (Fromang & Papaloizou 2007; however, see Simon et al. 2009) and hence the inefficiency of smallscale dynamo to operate at small Prandtl number (Schekochihin et al. 2005;Bodo et al. 2011).However, it is unclear what happens when convergence is recovered in unstratified ZNF simulations with tall boxes (Shi et al. 2016).
Studying MRI dynamo is also important for understanding the generation of coherent large-scale magnetic fields determining the level of transport (Johansen et al. 2009;Bai & Stone 2013) and outflows from the accretion disc (von Rekowski et al. 2003;Stepanovs et al. 2014;Mattia & Fendt 2022).MRI, in principle, can generate magnetic fields coherent over several scale-heights (Dhang et al. 2023) and acts locally as a mean field in the absence of any external flux influencing convergence and the disc dynamics.
Generally, stratified models generate a more coherent large-scale field over the unstratified models (for a comparison, see Shi et al. 2016).Cyclic behaviour of azimuthally averaged magnetic fields (mean fields), popularly known as the butterfly diagram, is a typical feature observed in the stratified shearing box simulations (Brandenburg et al. 1995;Gressel 2010;Bodo et al. 2014;Ryan et al. 2017;Gressel & Pessah 2022).However, note that the presence of a strong magnetic net flux (Bai & Stone 2013;Salvesen et al. 2016), convection (Hirose et al. 2014;Coleman et al. 2017) etc. can alter the periodicity in the butterfly diagram.Although the cyclic behaviour of mean fields can be explained by invoking the interplay between shear and helicity (Brandenburg & Donner 1997;Gressel & Pessah 2015), some features, such as upward migration of the mean fields, still demand an explanation.
Several studies attempted to understand the underlying mechanisms of MRI dynamo using different approaches.While some of the studies (Lesur & Ogilvie 2008;Bai & Stone 2013;Shi et al. 2016;Begelman & Armitage 2023) invoked toy models to complete the generation cycles of radial and azimuthal fields, others (local: Brandenburg et al. 2008;Gressel 2010;Shi et al. 2016;Gressel & Pessah 2022;Mondal & Bhat 2023, global: Dhang et al. 2020) used mean-field theory to investigate the large-scale field generation in the MRI-driven turbulent accretion flow.Most of the studies characterising the turbulent dynamo coefficients in the regime of mean-field dynamo theory used state of the art 'Test Field' (TF) method (Gressel 2010;Gressel & Pessah 2015), while a few used direct methods such as linear regression (Shi et al. 2016), singular value decomposition (SVD; Dhang et al. 2020) to calculate dynamo coefficients in postprocess or statistical simulations to carry out combined study of the large-scale dynamo and angular-momentum transport in accretion discs (Mondal & Bhat 2023).In this work, we use a direct method, a variant of the cleaning algorithm ( Högbom CLEAN method;Högbom 1974), called 'Iterative Removal Of Sources' (IROS; Hammersley et al. 1992); mainly used in astronomical image construction to analyse MRI-dynamo in the mean-field dynamo paradigm.We modified the IROS method according to our convenience (for details, see section 2.3, also see Bendre et al. 2023) and used it to determine the dynamo coefficients by post-processing the data obtained from the stratified ZNF shearing box simulation.
The paper is organised as follows.In section 2, we describe details of shearing box simulations, basics of mean field closure used and techniques of the IROS method.Section 3 describes the evolution of MRI to a non-linear saturated state, spatio-temporal variations of mean magnetic fields, EMFs and periodicities present in different observables.The spatial profiles of calculated turbulent dynamo coefficients, the reliability of the calculation method (using both EMF reconstruction and a 1D dynamo model) and contributions of each coefficient to the mean magnetic energy equation are described in section 2.3.In section 5, we discuss the plausible reasons behind different periodicities present (in mean magnetic fields, EMFs and helicties), comparison of our work with the previous works, the possible importance of a generative helicity flux and limitations of the averaging scheme and mean-field closure used in decoupling contributions from different dynamo coefficients.Finally we summarized our key results in section 6.

METHOD
The current work involves performing shearing box simulations of MRI-driven accretion flow, along with extracting dynamo coefficients using the mean-field dynamo model.In this section, we discuss details of the shearing-box simulation set-up, an introduction to the mean-field dynamo model and the IROS method used to determine turbulent dynamo coefficients using the simulated data.

Shearing-box simulation
We perform stratified zero net flux (ZNF) shearing box simulations to study the MRI driven dynamo in a geometrically thin disc (H/R ≪ 1).To do that, we solve ideal MHD equations in a Keplerian shearing box given by using the PLUTO code (Mignone et al. 2007) with x, y, z as the radial, azimuthal and vertical directions respectively.Here, ρ, P, v and B denote density, thermal pressure, velocity and magnetic fields, respectively.The terms gs = Ω 2 (2qxx − z ẑ) and 2Ωẑ × ρv represent the tidal expansion of the effective gravity and the Coriolis force respectively with Ω denoting orbital frequency.We use an isothermal equation of state Therefore, we do not need to solve the energy equation.Additionally, we use constrained transport (Gardiner & Stone 2005) to maintain divergence free condition for magnetic fields.We use the HLLD solver (Miyoshi & Kusano 2005) with second-order slope-limited reconstruction.Second-order Runge-Kutta (RK2) is used for time integration with the CFL number 0.33.Also note that despite our shearing-box model lacking explicit dissipation, we refer to it as the direct numerical simulation (DNS).
We initialize an unmagnetized equilibrium solution with density and velocity given by where q = 1.5 and ρ0 is the mid-plane (z = 0) density and is the thermal scale height.We set ρ0 = cs = Ω = 1, so that H = 1.Unless stated otherwise, all the length and time scales are expressed in units of H and Ω −1 respectively.We initialize a ZNF magnetic field given by with β0 = 10 4 defining the strength of the field and Lx, Ly, Lz denoting the size of the shearing-box.Our computational domain extends from −Lx/2 < x < Lx/2, −Ly/2 < y < Ly/2 and −Lz/2 < z < Lz/2.It has been found in earlier studies that shearing box results depend on the domain size; larger boxes tend to capture dynamo better than their smaller counterparts as well as smaller boxes demonstrate a transition to anomalous behaviour (e.g.see Simon et al. 2012;Shi et al. 2016).
To avoid these discrepancies, we choose a shearing box of size Lx × Ly ×Lz = 3H ×12H ×8H with a grid resolution Nx ×Ny ×Nz = 96 × 192 × 256 giving rise to a resolution of 32/H in the vertical direction.However, we must admit that there exists an issue with the convergence in stratified ZNF models as discussed in the section 1.We reserve the dependence of MRI dynamo on numerical resolution as a topic of future research investigation.
We use periodic and shearing-periodic (Hawley et al. 1995) boundary conditions in the y and x− boundaries, respectively.Outflow boundary conditions are implemented in the vertical (z) boundaries.A gradient-free condition is maintained for scalars and tangential components of vector fields at the boundaries.At the same time, vz ⩾ 0 for z > 0 and vz ⩽ 0 for z < 0 is set to restrict mass inflow into the domain at vertical boundaries.The z-component of the magnetic field is set by the divergence-free condition of the magnetic field.
Turbulent dynamo coefficient estimation involves analysis of time series of mean magnetic fields and EMFs obtained in shearing box simulation.Therefore, we dump the data quite frequently with data dumping interval ∆t = 0.2 Ω −1 and run it till t = 300 Ω −1 to have enough number data points in the time series of mean magnetic fields and EMFs.

Mean field closure
Before describing the details of mean field dynamo theory and the closure used, we define what is meant by 'mean' and 'fluctuation' in our work.We define mean magnetic fields (B) as the x − y-averaged values as follows Ly /2 −Ly /2 B(x, y, z, t) dx dy. (10) Fluctuating magnetic fields are defined as Mean and fluctuations of the x− and z− components of the velocity are defined in the same way as those for magnetic fields, while the mean and fluctuation of y−component of velocity are defined as If we decompose the magnetic and velocity fields into mean and fluctuation and insert them into the magnetic field induction equation, we obtain the mean-field equation where we assume that microscopic diffusivity is vanishingly small (ideal MHD limit).Here mean EMF appears as a source term in equation 13.The crux of the mean-field dynamo theory is how to express mean EMF in terms of the mean magnetic fields.In general, the usual mean-field closure (Raedler 1980;Brandenburg & Subramanian 2005;Shukurov & Subramanian 2021) is given by, where we neglect higher than the first-order spatial derivatives and time derivatives of mean magnetic fields and αij, ηij are the turbulent dynamo coefficients which characterize the dynamo; and Jj = ϵ jzl ∂z Bl (z) is the current.Further, while calculating turbulent dynamo coefficients using direct methods (e.g.SVD (Bendre et al. 2020;Dhang et al. 2020), linear regression (Squire & Bhattacharjee 2016;Shi et al. 2016), it is also assumed that αij, ηij are constant in time.However, we find that in our simulation of MRI-driven accretion flow, the current helicity, which is potentially a primary component determining the αij, shows a reasonably periodic change over time with a period half the dynamo-period (for details, we refer the reader to section 3.3).This time-dependent feature of current helicity leads to considering a heuristic mean field closure defined as to capture the time dependence in αij.Here α 0 ij and α 1 ij are the time-independent (i.e.DC component) and time-dependent parts of αij respectively and Ω dyn = 2πf dyn = 2π/T dyn , with f dyn and T dyn are the dynamo frequency and period respectively.Further, one expects that ηij to be dominated by a DC component, because ηij-s are generally determined by the turbulent intensity of the flow, not by helicities.Thus for simplicity, we adopt a time-independent ηij.

Dynamo coefficient extraction method -IROS
We solve equation 16 in a least-square sense to extract the turbulent dynamo coefficients (α 0 ij , α 1 ij and ηij) using mean magnetic fields Bi(z, t) and EMFs Ēi(z, t) (with i ∈ {x, y}) obtained from shearing-box simulations described in section 2.1.Further we assume that these dynamo coefficients stay statistically unchanged during the quasi-stationary phase of evolution, i.e. the coefficients are independent of time.Hence, all the dynamo coefficients are only dependent on the vertical coordinate z.
As a first step of coefficient determination in this under-determined system, we construct the time series of length N , of mean EMFs Ēi(z, t1 . . .tN ), mean magnetic fields Bi(z, t1 . . .tN ) and currents Ji(z, t1 . . .tN ) from the DNS (i ∈ {x, y}).With these time series, we rewrite equation 16 at any particular z = z ′ as where the matrices y, A and x are defined as, Here the terms Ci(z ′ , t l ) = Bi(z ′ , t l ) cos (2Ω dyn t l + ϕ) (∀i ∈ {x, y}) which take into account the time dependent part of αij.For simplicity, we assume ϕ to be zero.
Our task is to determine the dynamo coefficients (x) by pseudoinverting equation 17.This task is complicated firstly by the fact that both components of mean-field and current have additive correlated noise and secondly by the fact that the y component of the meanfield is typically much stronger compared to the x component, due to the rotational shear (and by consequence the x component of current is much stronger than its y component).Typical schemes of the least square minimisation in such cases tend to underestimate the dynamo coefficients that are associated with the x component of mean-field (i.e. the coefficients α 0 ix and α 1 ix ), and those with the y component of mean current (i.e. the coefficients ηiy).To circumvent these issues, we rely upon the IROS method (Iterative removal of sources) (Hammersley et al. 1992) that we have recently adapted for such inversions in the dynamo context (Bendre et al. 2023).This method is based on Högbom clean algorithm (Högbom 1974) used in Radio Astronomy to construct an image by convolving multiple beams, iteratively locating and subtracting out the strongest source to model the rest of the dirty image.It is particularly useful when the relative contribution of some of the beams to the final image happens to be negligible.Such a situation is analogous to have only a few of the columns of A (the beams) largely contribute to y (an image).A brief outline of the method is as follows.
Firstly, at any particular z = z ′ we set all the dynamo coefficients, ) and ηij(z ′ ) to zero, (i.e we set X(z ′ ) = 0).Then to compute these coefficients, we iteratively estimate their magnitudes as follows.To derive the zeroth order estimates of these coefficients we fit every i th column of Y(z ′ , t) denoted as Yi(z ′ , t)), against the individual columns of A(z ′ ) (denoted as A k (z ′ )) separately as lines.Slopes and chi-square errors (χ 2 ik (z ′ )) associated with each fit are recorded.The individual chi-square errors are defined as Then the best fitted dynamo coefficient (the one which has the least chi-square error) is updated by adding to it, its zeroth order estimate multiplied by a small factor (ϵ < 1), called the loop-gain, while other coefficients are kept constant.For example, if the chi-squared error associated with the line fit Ex(z ′ , t1 . . .tN ) versus By(z ′ , t1 . . .tN xy ) is updated by adding to it a factor of mϵ.Subsequently, the contribution to the EMF associated with the best fitted coefficient, also multiplied by the ϵ is subtracted from the corresponding EMF component.For instance, using the same example, a factor of ϵ αxy(z ′ ) By(z ′ , t1 . . .tN ) is subtracted from Ex(z ′ , t1 . . .tN ).This residual EMF is then used as an actual EMF component to further compute higher order estimates of dynamo coefficients, and this process is repeated a suitable number of times until either all the dynamo coefficients converge to their respective constant values or all four chi-squared errors get smaller than a certain predefined threshold.All the aforementioned steps are then repeated at every z = z ′ .
We apply this method with ϵ = 0.1 for five hundred refinement loops to the time series of EMFs, mean-fields and currents obtained from the DNS data.While constructing these time series (from t = 100Ω −1 to 300Ω −1 ) with data dumping interval ∆t dump = 0.2 Ω −1 we make sure that they correspond to the quasi-stationary phase of the magnetic field evolution.
IROS method does not provide an estimate of errors on the calculated coefficients directly.We, therefore, calculate a statistical error of the dynamo coefficient by considering the five different realizations of time series.We construct five different time series of mean fields, currents and EMFs by skipping four data points in the time series.Specifically, the time series (t1, t2, . . .tN ) (of all components of mean-field, current and EMF) are split into (t1, t6 . ..), (t2, t7 . ..), (t3, t8 . ..), (t4, t9 . ..) and (t5, t10 . ..).We use these time series to calculate five sets of dynamo coefficients and calculate their standard deviations to represent the errors on the calculated coefficients.

RESULTS: SATURATION OF MRI, MEAN FIELDS AND EMF-S
We now turn to the results of our shearing box simulation of MRI in a geometrically thin disc, investigate its dynamo action in addition to discussing several important properties which illuminates the nature of the MRI dynamo.Most of our analysis of magnetic field generation will focus on the saturated state of MRI, when the disc is in the quasistationary phase.

Saturation of MRI
First, consider the time evolution of accretion stresses and magnetic energies.This will also allow us to determine the quasi-stationary phase of the MRI-driven turbulence.The top panel of Fig. 1 shows the time history of accretion stresses (Reynolds and Maxwell).Nor- malized Reynolds and Maxwell stresses are defined as where the averages are done over the whole volume.Here, .and ⟨.⟩z indicate the average over x − y and z respectively.Reynolds stress is due to the correlation of fluctuating velocity fields, while Maxwell stress is composed of a correlation between the fluctuating components as well as that between the mean components of the magnetic fields.Both the stresses grow exponentially during the linear regime of MRI, and eventually saturate around an average value when MRI enters into the non-linear regime.In our simulation, we find the volume and time-averaged (within the interval t = (100 − 300) Ω −1 ) values of Reynolds and Maxwell stresses are to be W Rey,av = 0.0048 and W Max,av = 0.0167 respectively.The ratio of Maxwell to Reynolds stress is W Max,av /W Rey,av = 3.5, close to 4, as predicted by Pessah et al. 2006 for q = 1.5 and similar to what is found in earlier numerical simulations (Nauman & Blackman 2015;Gressel & Pessah 2015).
The bottom panel of Fig. 1 shows how the volume-averaged mean (⟨ B2 ⟩z) and fluctuating (⟨B ′2 ⟩z) magnetic energies evolve over time.Like accretion stresses, magnetic energies oscillate about an average value in the quasi-stationary phase after the initial exponentially growing phase.It is also worth noting that the mean part of the magnetic field shows a larger time variation than the fluctuating part of the magnetic field.We point out an important point that the fluctuating magnetic field is stronger than the mean magnetic field, and the implication of this will be discussed in the latter part of the paper.
We see in Fig. 1 that the accretion stresses and magnetic energies start saturating around t = 40 Ω −1 .However, to remain safer, we consider the simulation in the time range t = (100 − 300) Ω −1 for dynamo coefficient calculation in the quasi-stationary state.

Evolution of mean fields and EMFs
The most preliminary diagnostic of the dynamo is to look at the spatio-temporal variation of the mean magnetic fields, popularly known as the butterfly diagram (e.g.see the review by Brandenburg & Subramanian 2005).Fig. 2 shows the butterfly diagrams for mean magnetic fields Bx and By along with the mean EMFs Ēx and Ēy.Here we note that the mean EMF acts as a source term in the mean magnetic field energy evolution equation.In particular, Ēy is responsible for the generation of poloidal field (here Bx) from a toroidal one due to an α-effect, which itself naturally emerges by a combined action of stratification and rotation (Krause & Raedler 1980) in our stratified shearing box simulation.At an early stage of evolution (around t ≈ 2 orbital period), both mean fields and EMFs show lateral stretches with changing the sign in the vertical direction, which is clearly due to channel modes of MRI (Balbus & Hawley 1992;Balbus & Hawley 1998).During saturation, both mean fields and EMFs show a coherent vertical structure which changes signs in time with a definite period.We find magnetic field components By and EMF Ēx show a very coherent spatio-temporal variation with a time period of ≈ 9 orbital period (2π/Ω), similar to the earlier studies of MRI dynamo (Brandenburg et al. 1995;Davis et al. 2010;Gressel 2010;Gressel & Pessah 2015;Ryan et al. 2017).This periodicity is semi-transparent in the butterfly diagram of Bx, while this is hardly apparent for Ēy.However, we note that periodicities exist in all components of mean fields and EMFs as will become clear below (see Fig. 4).

Evolution of kinetic and current helicities
The generation of large-scale magnetic fields by a dynamo action is often associated with helicity in the fluid velocity field.Assuming isotropic homogeneous turbulence, Krause & Raedler (1980) suggested a kinetic α-effect defined by responsible for magnetic field generation; where τc is the correlation time, and K hel = v ′ .∇× v ′ is the kinetic helicity.It is suggested that α kin accounting for the effects of the helical velocity field, takes the role of driver, while αmag (Pouquet et al. 1976) defined by is the non-linear response arising due to the Lorentz force feedback, gradually increasing and ultimately quenching the kineticα (Blackman & Brandenburg 2002;Subramanian 2002).Here, is the current helicity.Ideally, the effective α− effect, responsible for poloidal field generation, is supposed to be α dyn = α kin + αmag.Fig. 3 shows the spatio-temporal variation of α kin and αmag.We assume correlation time τc to be same for both α-s and τc = Ω −1 .The αmag changes sign with a time-period ≈ 5 orbital period (2π/Ω), roughly half of the dynamo period, with which the mean fields and EMFs change sign, while α kin does not show any explicit periodicity.We will postpone a detailed discussion on the periodicity of helicities to section 3.5 where we discuss periodicities associated with all important variables.

Co-existence of small and large scale dynamos
Both kinetic and magnetic-α-s are small close to the mid-plane as shown in Fig. 3, while this is not true of the random kinetic and magnetic energies ( e.g.see Fig. 6 where we illustrate vertical profiles of rms value of random fluid velocity and Alfven velocity ).At the same time, the amplitudes of the helicities increase away from the mid-plane.These features suggest that both small-scale dynamo (when magnetic field grows because of the random stretching and twisting of the magnetic fields due to turbulent fluid motion) and large-scale dynamo (when magnetic field grows under the action of helicities) co-exist in MRI-driven dynamo (Blackman & Tan 2004;Gressel 2010).The MRI-driven small-scale dynamo dominates magnetic field generation close to the disc mid-plane where stratification is unimportant and helicities are small.In contrast, at larger heights where stratification becomes important, and helicities are large, a helicity-driven large-scale dynamo governs the magnetic field generation (Dhang & Sharma 2019;Dhang et al. 2020).However, it is to be noted that αmag is larger than α kin by one order of magnitude, and hence it is very likely that the effective-α will be predominantly due to αmag.

Power spectra of mean fields, EMFs and helicities
The butterfly diagrams shown in the previous sections depict the apparent periodicities of mean fields, EMFs and helicities.We look at the power spectrum defined by where q(z, t) is any generic quantity to investigate the periodicities in greater detail.Here, the spatial average is done over different heights, namely z = 0 − H, z = H − 2H and z = 2H − 3.5H to study the variation of periodicities at different scale heights.Fig. 4 shows the power spectra of mean fields Bx, By (top panels), mean EMFs Ēx, Ēy (middle panels) and helicities K hel , C hel (bottom panels).It is noticeable that power spectra for mean fields and spectra peaks at the primary frequency f dyn = 0.017 Ω (equivalent to ≈ 9 orbital period), which was also visible in the butterfly diagrams.In addition to the primary frequency, the power spectra also show the presence of higher harmonics (at 3f dyn ), which went unnoticed in  the earlier works of MRI dynamo.Similarly, power spectra of current helicity C hel also show the presence of higher harmonics (at 4 f dyn ) in addition to the primary frequency at 2f dyn .We also note that dynamo frequency remains almost constant at different heights.However, kinetic helicity does not show any periodicity.Presence of a strong time variation in αmag and its dominance over α kin necessarily leads to the expectation that turbulent dynamo coefficients (α− coefficients) should harbour a time-dependent part (α 1 ij ) along with the traditional time-independent part (α 0 ij ) as discussed in section 2.3.

RESULTS: DYNAMO COEFFICIENTS FROM IROS
We obtained mean fields ( Bx, By), EMFs ( Ēx, Ēy) from the shearingbox simulation and use a modified version of IROS method (see section 2.3) to calculate time-independent and time-dependent turbulent dynamo coefficients characterizing the MRI dynamo.However, we find the x−y-averaging cannot remove all the signatures of the smallscale dynamo.The small-scale dynamo is expected to have a shorter correlation time of order few Ω −1 and contribute noise at the higher frequency end compared to the large-scale dynamo.Therefore, we further smooth the mean fields and EMFs using a low-pass Fourier filter and remove contributions from the frequencies f > fc.We consider three cases: (i) fc = 0.05 Ω (≈ 3f dyn ), (ii) fc = 0.12 Ω (≈ 6f dyn ) and (iii) fc → ∞ (unfiltered) to assess the effects of the small-scale dynamo on the dynamo coefficient extraction.

Time independent dynamo coefficients
Fig. 5 shows the vertical profiles of time-independent dynamo coefficients α 0 ij and ηij for different values of fc.Four panels at the top illustrate the vertical profiles of coefficients (α 0 xx , α 0 xy , ηxx, ηxy) associated with the x-component of EMF Ēx, while four panels at the bottom show profiles of those (α 0 yx , α 0 yy , ηyx, ηyy) associated with the y-component of EMF Ēy.
The 'coefficient of most interest' out of the calculated ones is α 0 yy , which plays a vital role in producing the poloidal field (here Bx) out of the toroidal field ( By) (also see section 4.4) via an αeffect, associated with the helicities (see section 3.3 ).The coefficient α 0 yy shows an anti-symmetric behaviour about the z = 0 plane, with a negative (positive) sign in the upper (lower)-half plane (for |z| < 2).For |z| > 2, the sign of α 0 yy tends to be positive (negative) in the upper (lower)-half plane.Earlier studies of MRI dynamo in local (Brandenburg 2008;Gressel 2010;Gressel & Pessah 2015) and global (Dhang et al. 2020) frameworks also found a similar trend in α 0 yy .However, it is to be noted that our study suggests a stronger negative α 0 yy in the upper-half plane compared to that in the earlier studies.The negative sign in the upper half plane is attributed to the buoyant rise of magnetic flux tubes under the combined action of magnetic buoyancy and shear (Brandenburg & Schmitt 1998;Brandenburg & Subramanian 2005; see also Tharakkal et al. (2023)).Brandenburg & Schmitt (1998) also suggested that negative αyy is responsible for the upward propagation direction of dynamo waves seen in the butterfly diagrams of MRI-driven dynamo simulations (e.g.see Fig. 2).Another different way of looking at the origin of the effective α is to link it to the helicity flux as envisaged by Vishniac 2015 and Gopalakrishnan & Subramanian 2023.We discuss this possibility in section 5.3.
The off-diagonal terms of the α-coefficients are related to turbulent pumping.This effect is responsible for transporting large-scale magnetic fields from the turbulent region to the laminar region.We found α 0 xy and α 0 yx to be antisymmetric and α 0 xy > α 0 yx unlike the previous studies (Brandenburg 2008;Gressel & Pessah 2015) which found α 0 yx ≈ α 0 xy .This resulted in a strong turbulent pumping γz = (α 0 yx − α 0 xy )/2, transporting large-scale magnetic fields from the disc to the corona as shown in the top panel of Fig. 6.We also compare the relative importance of turbulent pumping (γz) and wind (vz) in advecting the magnetic field upward (in the upper half-plane) at different heights.Vertical profiles of γz and vz in the top panel of Fig. 6 shows that at low heights (|z| < 2.5), turbulent pumping is the dominant effect over the wind, while the effects of wind become comparable or larger than the pumping term at large scale-heights (see also Fig. 11).
The theory of isotropic kinematically forced turbulence predicts that γz is supposed to be in the direction of negative gradient of turbulent intensity (v ′2 ) (Krause & Raedler 1980), that is, in the negative z-direction (in the upper-half plane) in our simulation.This is opposite to what has been found in Fig. 6.However, it is to be noted that MRI turbulence in a stratified medium is neither isotropic nor homogenous.Minimal τ -approximation (MTA) suggests that in a stratification and rotation-induced anisotropic turbulent medium, which includes the quasi-linear back reaction due to Lorentz forces, where τ is the correlation time and it is assumed that ρ = 1 (see equation (10.59) in Brandenburg & Subramanian 2005).The last term in equation 24 vanishes because all the variables are functions of z alone.Therefore, equation 24 and the bottom panel of Fig. 6 illustrating the vertical profiles of v ′2 and v ′2 A imply that sign of turbulent pumping obtained from MTA supports that obtained from extracted dynamo coefficients.We found turbulent diffusion tensor ηij to be anisotropic with ηxx > ηyy and having a significant contribution from the offdiagonal components ηxy and ηyx.Different values of diagonal components of ηij imply that mean field components Bx and By are affected differently by the vertical diffusion (also see section 4.4).It is worth mentioning that ηyy ≈ 0 for the fc = 0.05 case, while it is slightly negative for the other two cases.This is somewhat different from the earlier studies (Gressel 2010;Gressel & Pessah 2015), which calculated dynamo coefficients using the TF method and found ηxx ≈ ηyy > 0. Out of the two off-diagonal terms of the diffusion tensor, ηyx is of particular interest.It is suggested that a negative value of ηyx can generate poloidal fields by the shear-current effect (Squire & Bhattacharjee (2016)).However, we find ηyx to be always positive, nullifying the presence of a shear-current effect in our stratified MRI simulation.
Finally, we discuss the effects of filtering the time series of mean magnetic fields and EMFs on the extracted dynamo coefficients.Fig. 5 illustrates how the dynamo coefficients vary if we filter out the contribution above a cut-off frequency fc with (i) fc = 0.05 Ω (≈ 3f dyn ), (ii) fc = 0.12 Ω (≈ 6f dyn ) and (iii) fc → ∞ (unfiltered).Broadly speaking, while the coefficients (α 0 xx , α 0 yx , ηxy, ηyy ) associated with Bx and its derivative in the mean-field closure (equation 16) show larger variations at higher heights with fc, those (α 0 xy , α 0 yy , ηxx, ηyx) associated with By and its derivative are less affected by the filtering process.Especially, ηyy tends to be more positive with fc = 0.05, which is more desirable.To summarize, filtering out the time series of mean magnetic fields and EMFs helps to remove the signature of the small-scale dynamo and to obtain noise free coefficients.

Time dependent dynamo coefficients
We discussed the time-dependent nature of αmag in the previous sections.Effective α-effect is expected to be determined by αmag, especially at the larger scale heights where it is of larger amplitude.While α tensor is expected to have the time-dependent part, η-tensor is supposed to have only the time-independent part, as it only depends on the turbulent intensity (see section 2.3).Fig. 7 shows the vertical profiles of time-dependent α-tensor components for the fiducial fc = 0.05 Ω case.For comparison, we also plot vertical profiles of the time-independent α−s in Fig. 7.We find that the coefficients αxx and αyx associated with Bx in the mean-field closure (equation 16) have stronger time-dependence compared to those coefficients (αxy and αyy) associated with By.Overall, the amplitudes of α 1 ij are much smaller than the α 0 ij implying that the time-independent α−s are predominantly governing the dynamo action.Additionally, we observed that (not shown in Fig. 7) α 1 ij -s in the fiducial case (fc = 0.05) are relatively smaller than the other two cases (fc = 0.12 and Unfiltered).

Verification of method
To verify the reliability of the determined dynamo coefficients we reconstruct the EMFs using the calculated coefficients and run a 1D dynamo model.  .Vertical profiles of time-independent turbulent dynamo coefficients (α 0 ij , η 0 ij ) in MRI simulation calculated using IROS method.A low-pass Fourier filter with a cut-off frequency fc removes the contribution from the small-scale dynamo.We used two values of fc: fc = 0.05 Ω and fc = 0.12 Ω.The results are compared to the case when IROS is applied to the unfiltered data obtained from DNS. Shaded regions associated with each line in the plot represent ±1σ statistical error on calculated coefficient as described in section 2.3.

Reconstruction of EMFs
Fig. 8 shows butterfly diagrams of the EMFs ( Ēx,f , Ēy,f ) used to determine the turbulent dynamo coefficients and the EMFs ( Ēx,r, Ēy,r) reconstructed using calculated coefficients and mean fields for fc = 0.05.Here it is to be noted that Ēx,f , Ēy,f are the smoothed EMFs obtained by filtering (by using a low-pass filter) EMFs Ēx, Ēy from DNS respectively.We can see a close match between the broad features, such as the dynamo cycle period, in the smoothed and reconstructed EMFs, implying the goodness of fit.
Further, we investigate the residual of the filtered and reconstructed EMFs, defined by for Ēx over Ēy is expected as Ēx obtained from DNS shows a more regular, coherent space-time variation when compared to Ēy.

1D dynamo model
We additionally run a 1D dynamo model using the calculated dynamo coefficients and mean velocity field vz.In particular we solve equation 13, or in component form for Bx and By with α 0 ij and ηij obtained using the IROS method.We note that Bz = 0 as a consequence of the zero net flux (ZNF) assumption in our model.The initial profiles of Bx and By at are taken directly from the DNS, at time t = 100 Ω −1 roughly consistent with the beginning of the quasi-stationary phase in the DNS.Vertical Top panel: profiles of turbulent pumping (γz) and mean vertical outflow (vz).They act in the same direction, transporting mean fields vertically outward.Bottom panel: vertical profiles of average fluctuating velocity (v ′2 ) and fluctuating Alfven speed A suggest similar sign of γz as calculated using IROS.
profile of vz is taken as a constant throughout the evolution and is also extracted from the direct simulations by averaging it over time throughout the quasi-stationary phase, over which it roughly stays constant.Additionally, for the profiles of dynamo coefficients α 0 ij (z) and ηij(z), we first smooth them with a box filter and also cut them off above and below three scale heights, and use them in the 1D dynamo model.We do this mainly to avoid the numerical instability at boundaries noting that these profiles are sharply flayed outside of that range.Note that only the time independent parts of the dynamo coefficients are used in the mean field equations, since the contributions of α 1 ij are negligible compared to the time-independent part.
Furthermore, it must be noted that there is a contribution to the diffusion from the mesh grids.We do a rough estimation of numerical diffusion as follows η0 = v ′ rms ∆x, where we consider the smallest one among the relevant velocities (v ′ rms , cs, vA) in the problem.Therefore, we add a correction term η0 ≈ 10 −3 (with ∆x = 1/32 and v ′ rms = 0.1) to the diagonal components of diffusivity tensor ηij to consider the contribution from the mesh to the magnetic field diffusion.This also helps us to stabilize the 1D dynamo solution.
With this setup, we solve the system of equations (equation 26) with a finite difference method over a staggered grid of resolution ∆z = 1/32, same as the z resolution of DNS.The outcome of this analysis is presented in Fig. 10, where the top and bottom panels show the butterfly diagrams of Bx and By obtained using the 1D dynamo model respectively.We find both x and y-components of mean fields flip sign regularly with a cycle of ≈ 9 orbital period, similar to what is found in DNS (see Fig. 2).Thus, applying calculated coefficients to the 1D dynamo model successfully reproduces broad features of spatiotemporal variations mean magnetic fields.

Mean magnetic energy equations
It is challenging to calculate dynamo coefficients uniquely in the presence of both shear and rotation (Brandenburg et al. 2008) as there are many unknowns (see also discussion in section 5.4).Therefore, it is worth seeing how different terms involving turbulent dynamo coefficients contribute to the mean magnetic energy equation to make physical sense.The mean magnetic energy evolution equation is obtained by taking the dot product of the mean-field equation (equation 13) with the mean magnetic field B and given by ∂ ∂t ∂ ∂t where Fig. 11 shows the space-time plots of different terms involving mean flow (vz) and turbulent dynamo coefficients (αij, ηij) in the mean magnetic energy evolution equations.The top six panels in Fig. 11 describe the terms in the x-component of the magnetic energy equation (equation 27), while the bottom seven panels illustrate terms in the y-component of the magnetic energy equation ( 28) at different heights and times.
Fig. 11 provides a fairly complicated picture to account for the generation-diffusion scenario of the mean magnetic fields.Broadly speaking, the poloidal field ( Bx) is predominantly generated by an α-effect (the term Tα yy in Fig. 11).However, there is a significant contribution from αyx (the term Tα yx in Fig. 11) in generating Bx in larger scale-heights.Toroidal field generation is mainly due to the presence of shear, here differential rotation, (TS in Fig. 11) which converts poloidal fields to the toroidal fields.However, it is worth noting that there is a minute contribution from the αxx, generating a toroidal field out of the poloidal field by an α-effect, (as in an α 2 -Ω dynamo).The dominance of α-effect in generating a poloidal field and that of Ω-effect (shear) in generating a toroidal field imply the presence of an α − Ω type dynamo in MRI-driven geometrically thin accretion disc.This is similar to what has been found in the study of the dynamo in an MRI-driven geometrically thick accretion disc Additionally, we observed the presence of higher harmonics at 3f dyn , which went unnoticed in earlier MRI simulations (see section 3.5).Unlike the mean fields and EMFs, current helicity shows periodicities at different frequencies 2f dyn and 4f dyn , respectively.The presence of the frequencies in the mean EMFs, mean fields, and current helicities can be understood better if we follow the magnetic helicity density evolution equation (e.g.see Blackman & Field (2000); Subramanian & Brandenburg (2006); Kleeorin & Rogachevskii (2022); Gopalakrishnan & Subramanian (2023)), where h b = ⟨A ′ .B ′ ⟩ is magnetic helicity density, A ′ is the fluctuating vector potential, and FH is the helicity flux.Roughly speaking, magnetic helicity is related to current helicity (and αmag, see equation 22) by some length scale and therefore, we can investigate equation 30 to shed light on the time variation of current helicity.
The component of the EMF along the mean magnetic field generates mean magnetic and associated current helicities.Now to consider the effect of the DC term in αyy, we assume that this is the dominant term in generating the poloidal field, which is a valid approximation, as we noted in section 4.4 (also see Fig. 11).Then, Ē. B ≈ α 0 yy B2 y , which is a source term in equation 30.Now, e.g. for simplicity, if we assume, By ∼ sin (2πf dyn t), then magnetic and current helicities which are ∝ B2 y , will have primary frequency of 2f dyn .This explains the generation of magnetic and current helicities at a primary frequency, twice that of By, i.e. 2 f dyn .This current helicity can now add to the α-effect, which combined with the mean field in the dynamo equation 13, can lead to secondary EMF and mean fields components oscillating at 3f dyn which in turn sources helicity components at 4f dyn and so on.These primary and secondary frequency components, limited by noise, are indeed seen from the analysis of our simulations.

Dynamo coefficients, comparison with earlier studies
Earlier studies calculating turbulent dynamo coefficients using the simulation data and the mean field closure (equation 16) in the local (Brandenburg et al. 1995;Brandenburg 2008;Gressel 2010;Gressel & Pessah 2015;Shi et al. 2016) and global (Flock et al. 2012;Hogg & Reynolds 2018;Dhang & Sharma 2019;Dhang et al. 2020) (Dhang et al. 2020) studies used direct methods to quantify dynamo coefficients.However, it is important to note that while several authors used a linear regression method assuming few constraints on the diffusion coefficients (namely, ηxx = ηyy), we use the IROS method without any constraints on the coefficients.Like most of the earlier local and global studies, we find a neg-ative αyy close to the mid-plane in the upper half-plane.However, direct methods seem to capture negative signs better than TF; which can be realised by comparing αyy profiles in our work (also in Shi et al. (2016)) and in Gressel (2010).Additionally, we find stronger turbulent pumping (compared to that in the TF method), transporting large-scale magnetic fields from the disc to the corona, similar to that found in global MRI-dynamo studies (Dhang et al. 2020).Additionally, for the first time, we ventured to calculate the timedependent part of αij inspired by the periodic behaviour of αmag.However, we found that the amplitudes of the time-dependent part of α−s (α 1 ij ) are much smaller than that of the time-independent α−s (α 0 ij ).Therefore, we suspect that the time-independent α−s are predominantly governing the dynamo action.
Diffusivity coefficients ηij in our work are found to be quite different from that in the earlier local studies (Brandenburg 2008;Gressel 2010;Gressel & Pessah 2015;Shi et al. 2016), with ηxx ̸ = ηyy and ηyy ≈ 0. Several earlier studies (Shi et al. 2016;Zier & Springel 2022) found ηyx < 0) in their unstratified and stratified MRI simulations after imposing a few constraints on the coefficients (e.g.ηyy = ηxx, ηxy = 0 etc.) and they proposed shear-current effect (Raedler 1980;Rogachevskii & Kleeorin 2004;Squire & Bhattacharjee 2016) generating poloidal fields in addition to α-effect.Recently, Mondal & Bhat (2023) carried out statistical simulations of MRI in an unstratified ZNF shearing box and found ηyx < 0 proposing 'rotation-shear-current effect' and the 'rotation-shear-vorticity effect' responsible for generating the radial and vertical magnetic fields, respectively.However, like some other studies (TF: Brandenburg 2008;Gressel 2010;Gressel & Pessah 2015), SPH:Wissing et al. 2022), we find ηyx ⩾ 0, unless we impose a constraint on ηyy being a positive fraction of ηxx.If we assume ηyy = fη ηxx while calculating the coefficients, we find negativity of ηyx is an increasing function of the factor fη (see Fig. A1 and Appendix).However, we find that the quality of fit is compromised slightly and histograms of the residual of filtered (input) and reconstructed EMFs get broader (with higher standard deviation) with the assumption ηyy = fη ηxx.We refer the reader to see Appendix for details.

Helicity flux and the DC α−effect
The coefficient α 0 yy represents the α−effect responsible for poloidal magnetic field generation out of the toroidal field.We found an anti-symmetric profile of α 0 yy about the disc-midplane similar to the earlier studies (e.g.see Brandenburg 2008;Gressel 2010).However, it is to be noted that understanding of the physical mechanism determining the vertical profile of α 0 yy is incomplete.E.g.Brandenburg & Schmitt (1998) proposed a buoyancy-driven dynamo to explain the negative sign of αyy in the upper half plane.Here, we propose a different way of looking at the origin of α−effect by connecting it to a generative helicity flux.
In order to understand the DC value (time-independent) of the α-effect, we take the time average of equation 30.The term ∂h b /∂t averages to zero, and one gets the well-known constraint (Blackman 2016 where ⟨⟩ indicates a time average.This shows that in the absence of helicity fluxes, the average EMF parallel to the mean field, responsible for the generation of poloidal from the toroidal mean field, is resistively (or catastrophically) quenched.Of the several helicity fluxes discussed in the literature, the generative helicity fluxes as envisaged in Vishniac (2015)   and we have taken q = 3/2.Adopting estimates for the correlation time τ ∼ Ω −1 , correlation length λ ∼ H/2 and using the vertical profiles of various physical variables from the simulation, we calculate the vertical profile of α 0 yy hc due to the generative helicity flux.This is shown as a solid line in Fig. 12 and for comparison, we also show 10α 0 yy (for fc = 0.05 case) from the IROS inversion.It is encouraging that the (α 0 yy ) hc predicted by the generative helicity flux is negative in the upper-half plane of the disc and has a qualitatively similar vertical profile as that determined from IROS inversion.The amplitude, however, is larger, which perhaps indicates the importance also of the neglected diffusive and advective helicity fluxes which act as sink terms in equation 31.

Vanishing ηyy, missing information?
In section 4.4 we pointed out that wind carries away the mean magnetic field and acts as the effective sink of its energy.However, the poloidal field is also expected to be diffused by ηyy, and a positive ηyy is required for diffusion.Instead, we find a vanishingly small (in some regions even negative) ηyy, which leads us to two possible thoughts: either it is impossible to recover ηyy in the direct methods, or there is incompleteness in the closure we used to retrieve the coefficients.Here, we discuss both possibilities.
It is clear from equation 16 that the turbulent diffusion coefficients are associated with the currents, which are calculated by taking the zderivative of mean magnetic field components.Calculating derivative makes the currents noisy, especially Jy, as it involves a derivative of Bx, which is fairly incoherent over space and time, as can be seen from the butterfly diagram of Bx (Fig. 2).Additionally, also note that the y-component of EMF is also noisy.Thus, the coefficients associated with Jy and Ēy turned out to be error-prone and difficult to calculate.This pattern has been noticed by earlier works (Squire & Bhattacharjee 2016), which used direct methods other than IROS, used in the current work.
In general, mean EMF can be expressed in terms of symmetric, anti-symmetric tensors and mean fields as follows, where we neglect the higher than first-order spatial derivatives and time derivatives of mean fields (Raedler 1980;Brandenburg & Subramanian 2005;Schrinner et al. 2007;Simard et al. 2016).The co-efficients α and γ represent the symmetric and anti-symmetric parts of αij tensor in equation 13.The coefficients αij is related to αxx and αyy, while γ represents turbulent pumping.The term η is a ranktwo tensor representing diffusivity.The coefficient δ is interpreted as a magnetic field generating term (Raedler 1980).The κ-term is a third-rank tensor having a complicated influence on mean fields.( Therefore it is evident from equation 34 that it is impossible to decouple a few coefficients (coefficients in the last three identities) as there are more unknown coefficients than independent variables ( B, Ē) and the actual diffusion coefficients (ηij) might be different from the calculated ones (ηij).

SUMMARY
We carried out stratified zero net flux (ZNF) shearing-box simulations of MRI using ideal MHD approximation.We characterised the MRI-driven dynamo using the language of mean field dynamo theory.The turbulent dynamo coefficients in the mean-field closure are calculated using the mean magnetic fields and EMFs obtained from the shearing box simulation.For this purpose, we used a cleaning (or inversion) algorithm, namely IROS, adapted to extract the dynamo coefficients.We verified the reliability of extracted coefficients by reconstructing the EMFs and reproducing the cyclic pattern in mean magnetic fields by running a 1D dynamo model.Here we list the key findings of our work: • We find mean fields and EMFs oscillate with a primary frequency f dyn = 0.017 Ω (≈ 9 orbital period).Additionally, they have higher harmonics at 3f dyn .Current helicity αmag has two frequencies 2f dyn and 4f dyn respectively.These frequencies can be understood from mean-field dynamo effective dispersion relation and helicity density evolution equation, respectively (for details, see section 5.1).
• Our unbiased inversion and subsequent analysis show that an α−effect (αyy) is predominantly responsible for the generation of poloidal field (here Bx) from the toroidal field ( By).The differential rotation creates a toroidal field from the poloidal field completing the cycle; indicating that an α − Ω-type dynamo is operative in MRI-driven accretion disc.
• We find encouraging evidence that the effective DC α−effect can be due to a generative helicity flux (section 5.3).
• We find strong wind (vz) and turbulent pumping (γz) carry out mean fields away from the mid-plane.Interestingly, they act as the principal sink terms in the mean magnetic energy evolution equation instead of the turbulent diffusivity terms.
• The unbiased inversion finds an almost vanishing ηyy, while ηxx and ηyx are positive.However ηyx and ηyy are strongly correlated; if one imposes an arbitrary prior that ηyy = fηηxx, then one finds for increasing fη, an increasingly negative ηyx which has been interpreted as evidence of shear-current effect for generating poloidal fields (see Appendix A).
• We point out that defining mean fields by planar averaging can necessarily introduce degeneracy in determining all the turbulent dynamo coefficients uniquely.This may have important consequences for the physical interpretation of the dynamo coefficients (see section 5.4).

Figure 2 .Figure 3 .
Figure 2. Spatio-temporal variation of mean magnetic fields, Bx (top left panel), By (bottom left panel) and mean EMFs Ēx (top right panel) and Ēy (bottom right panel).Mean magnetic field component By and y-component of EMF Ēx show a coherent change in space and time (with a time period ≈ 9 orbital period (2π/Ω)), while the spatio-temporal patterns in Bx and Ēy are less coherent.
Figure5.Vertical profiles of time-independent turbulent dynamo coefficients (α 0 ij , η 0 ij ) in MRI simulation calculated using IROS method.A low-pass Fourier filter with a cut-off frequency fc removes the contribution from the small-scale dynamo.We used two values of fc: fc = 0.05 Ω and fc = 0.12 Ω.The results are compared to the case when IROS is applied to the unfiltered data obtained from DNS. Shaded regions associated with each line in the plot represent ±1σ statistical error on calculated coefficient as described in section 2.3.
Fig.9shows the histograms of the normalised residuals δ Ēx/| Ēx| and δ Ēy/| Ēy| calculated within the region of different heights, namely between 0 − H, H − 2H and 2H − 3H, for the fc = 0.05 Ω case.All the histograms peak about the region close to zero.However, a Gaussian fit of the histograms shows that the mean of the distribution always deviates from zero.Additionally, a careful comparison of histograms of δ Ēx/| Ēx| and δ Ēy/| Ēy| tells that fit is better for Ēx than that for Ēy, especially at larger scale-heights.Better quality fit

Figure 7 .Figure 8 .
Figure7.Vertical profiles of time-dependent turbulent dynamo coefficients (α 1 ij ) in MRI simulation calculated using IROS method for fc = 0.05.Shaded regions associated with each line in the plot represent ±1σ statistical error on calculated coefficient as described in section 2.3.The time-dependent α−s are compared to time-independent α− for the same value of cut-off frequency fc = 0.05 Ω.Smaller amplitudes of α 1 ij compared to that of α 0 ij implies that that the time-independent α−s are predominantly governing the dynamo action.

Figure 11 .
Figure 11.Contributions of different terms involving mean flow (vz) and turbulent dynamo coefficients (α ij , η ij ) to the x-(top six panels) and y-(bottom seven panels) components of mean magnetic energy evolution equation (equations 27 and 28).Each term in the poloidal and toroidal magnetic energy equations are multiplied by the factors 10 6 and 10 4 respectively.Poloidal field ( Bx) generation is primarily attributed to an α-effect (the term Tα yy ), while shear (the term T S ) dominates the toroidal field generation; thus implying and α − Ω type of dynamo.Winds carry mean fields out of the computational box and contribute largely as the sink term in the mean magnetic energy evolution equation.
If we define mean fields and EMFs as the x − y-averaged quantities, then mean field closure reduces to equation 15.The symmetrised coefficients in equation 33 and non-symmetrised coefficients in equation 15 are related as αxx = αxx,