Finite-time Response of Dynamo Mean-field Effects in Magnetorotational Turbulence

Accretion disk turbulence along with its effect on large-scale magnetic fields plays an important role in understanding disk evolution in general, and the launching of astrophysical jets in particular. Motivated by enabling a comprehensive subgrid description for global long-term simulations of accretions disks, we aim to further characterize the transport coefficients emerging in local simulations of magnetorotational disk turbulence. For the current investigation, we leverage a time-dependent version of the test-field method, which is sensitive to the turbulent electromotive force (EMF) generated as a response to a set of pulsating background fields. We obtain Fourier spectra of the transport coefficients as a function of oscillation frequency. These are well approximated by a simple response function, describing a finite-time buildup of the EMF as a result of a time-variable mean magnetic field. For intermediate timescales (i.e., slightly above the orbital frequency), we observe a significant phase lag of the EMF compared to the causing field. Augmented with our previous result on a nonlocal closure relation in space, and incorporated into a suitable mean-field description that we briefly sketch out here, the new framework will allow us to drop the restrictive assumption of scale separation.


INTRODUCTION
Thirty years after its ultimate discovery by Balbus & Hawley (1991), the magnetorotational instability (MRI) is practically synonymous with accretion disc turbulence and is believed to be the key to understanding the structure and evolution of discs ranging from circumplanetary, to circumstellar (including those around black-holes or neutron stars, as well as the innermost and outer reaches of protoplanetary ones), all the way to active galactic nuclei.
Enhanced transport coefficients -stemming from correlated fluctuations in the MRI turbulence-have a profound impact on disc evolution.This can happen either directly, when disc accretion is enabled by angular momentum exchange through enhanced viscosity, or indirectly via a magnetocentrifugal disc outflow.Notably, the latter scenario requires large-scale ordered poloidal fields that are either created in situ by a disc dynamo (e.g., von Rekowski et al. 2003;Stepanovs et al. 2014), or -when an inherited large-scale field is invoked-are at least affected by enhanced field dissipation as a result of eddy diffusivity.Recent attempts of in-Corresponding author: Oliver Gressel ogressel@aip.decorporating sub-grid-scale physics into jet-launching simulations (e.g.Bucciantini & Del Zanna 2013;Fendt & Gaßmann 2018;Dyda et al. 2018;Mattia & Fendt 2020a,b;Vourellis & Fendt 2021) illustrate the need for comprehensive parametrizations that are ideally based on first-principles, resolved MRI simulations.
Depending on (i) the level of inherited / accumulated netvertical magnetic flux, and (ii) the relevance of the vertical disc structure, the MRI relies to a varying degree on the presence of an intrinsic dynamo of some sort to become a self sustained mechanism for powering disc accretion (see Rincon 2019, for an excellent review on this subject).A pronounced shortcoming of non-stratified box simulations is that they are very sensitive to the vertical aspect ratio (see Shi et al. 2016;Walker & Boldyrev 2017).
Notably, when including vertical stratification, both localbox (going back to Brandenburg et al. 1995) and global (as recent as Dhang et al. 2020) fully non-linear MRI simulations alike robustly develop near-periodic cycles in the (horizontally / azimuthally averaged) mean magnetic field with a characteristic propagation away from the disc midplaneproviding a natural explanation to sustaining the MRI via a large-scale dynamo action (e.g.Brandenburg 2005Brandenburg , 2008;;Blackman 2010).
The morphology of this so-called butterfly diagram -the hallmark of the mean-field dynamo-was previously found to depend somewhat on the amount of net-vertical magnetic flux (see, e.g., Gressel & Pessah 2015;Salvesen et al. 2016).We here nevertheless focus on the limit of negligible netvertical magnetic flux, which in a way is the crucial test for providing a robust accretion engine from MRI turbulence.Another important issue raised pertains to the onset of convective turnover (Bodo et al. 2012;Gressel 2013;Hirose et al. 2014), which was found to drastically affect the regularity of the magnetic-field cycles (see the discussion in Coleman et al. 2017).As with the net-vertical magnetic flux, we take a rather conservative stance and focus our investigation on the isothermal case, avoiding the complications associated with arguably more realistic thermodynamic representations.
A central question that remains unanswered is whether the dynamo wave can be reconciled with a conventional αΩ dynamo (i.e., driven via the interplay of helical turbulence with differential rotation), and/or whether its dynamics are enforced by the near-exact conservation of magnetic helicity at high magnetic Reynolds number (e.g.Vishniac 2009; Gressel 2010; Oishi & Mac Low 2011).While the cycle period as a function of shear rate can nicely be explained using the dispersion relation of a near-critical αΩ dynamo (Gressel & Pessah 2015), the propagation direction away from the midplane is still not well understood, possibly requiring a magnetic buoyancy contribution near the midplane (Brandenburg 1998).
Both the spatial non-locality (see Brandenburg & Sokoloff 2002) of the dynamo closure relation, the non-instantaneous aspects (i.e., so-called "memory effects", e.g., Hubbard & Brandenburg 2009), as well as their combined effect (see, e.g., Rheinhardt & Brandenburg 2012) have been demonstrated to influence the characteristics of the dynamo cycle.Another comprehensive example of how finite-time effects can influence dynamo-generated fields has been presented by Chamandy et al. (2013a,b) in the context of galactic magnetic fields.To complement our previous investigation of the scale-dependence / non-locality of the characteristic meanfield α effect in MRI turbulence (Gressel & Pessah 2015, sect. 3.4), we here investigate the potential role of finite-time effects in the mean-field closure relation.
Our paper is organized in the following manner: Section 2 briefly describes the numerical simulations and introduces the newly adopted non-instantaneous closure relation to the mean-field induction equation, as well as how it can be captured using the test-field method.We present the results obtained from a fiducial MRI shearing-box simulation in Section 3, and we discuss how these findings may be exploited in the future, in Section 4.

METHODS
As in previous work, we solve the equations of isothermal, ideal magnetohydrodynamics (MHD) in a local shearing-box (e.g.Gressel & Ziegler 2007) frame of reference.Lacking explicit dissipation, the purist may call this an "implicit" large-eddy simulation (iLES).For practical purposes, we will nevertheless refer to these as direct numerical simulation (DNS).For brevity, we here only briefly recapitulate the essential properties of our numerical approach, and refer the reader to sections 2.1 et seqq. of Gressel & Pessah (2015) for a more in-depth discussion, motivating our particular choices.

Brief specification of the direct simulations
We here use local Cartesian coordinates, (x, y, z), but refer to some tensor coefficients in cylindrical components, (r, φ, z), for easier comparison with global models.Differential rotation is expressed via the parameter q ≡ d ln Ω/d ln r = −3/2 for a Keplerian rotation curve, and we use the "orbital advection" scheme of Stone & Gardiner (2010) to treat the background shear flow, v K ≡ q Ωx ŷ, with the benefit of a position-independent truncation error.The equations expressed in the local Eulerian velocity, v, are with the total pressure p ≡ p + B 2 /2, and the combined (i.e., tidal plus gravitational) effective potential defined in the locally co-rotating frame of reference at fixed angular frequency Ω ≡ Ωẑ.Horizontal boundary conditions are shear-periodic (see Gressel & Ziegler 2007, for details), and we apply standard outflow conditions in the vertical direction.We chose an intermediate box size of L x ×L y ×L z = H ×πH ×6H with a linear resolution of 32/H in all space dimensions -amounting to 32 × 100 × 192 cells in the radial (x), azimuthal (y), and vertical (z) coordinate directions, respectively.The initial plasma parameters in the disk midplane are β p = 800 and β p = 2.2 ×10 5 for the zeronet-flux contribution, and the additional net-vertical field, respectively.As previously, we include (i) an artificial mass diffusion term (see Gressel et al. 2011) to circumvent undue time-step constraints resulting from low-density regions in the upper disc corona, and (ii) replenish the mass lost via outflow through the vertical domain boundary to obtain an overall steady-state disc structure.

The non-instantaneous closure relation
Adopting the well established framework of mean-field magnetohydrodynamics (Krause & Raedler 1980), we seek a parametrization for the turbulent electromotive force, E ≡ v ×B with fluctuating magnetic and velocity fields defined as B ≡ B−B(z), and v ≡ v−v(z), respectively.1Here, and in the following, the overbar implies geometric averaging over horizontal slabs.This is the natural choice for the adopted box geometry and trivially satisfies the Reynolds rules required for a consistent mean-field description.By virtue of its definition, the EMF captures correlations in fluctuating velocity and magnetic field, whose non-zero mean appears as a source term on the right-hand-side of the (one-dimensional) mean-field induction equation By construction, this equation describes the long-term evolution of the (comparatively slowly changing) mean magnetic field under the effect of the underlying turbulence.To make progress over a direct simulation approach, the EMF is then typically expanded into a linear functional of the mean magnetic field and its gradients as where the indices i, j, l label the coordinates x, y and contraction over repeated indices is understood.Note that it is unnecessary to include the radial and azimuthal gradients in our case, which is due to the homogeneity of the turbulence in any given horizontal plane.
Under steady-state conditions, the second-rank tensors, α ij (z) and η ij (z) become time-independent and are expected to capture the statistical properties of the chaotic flow (see Krause & Raedler 1980).If the system at hand is sufficiently anisotropic (e.g., due to rotation) and inhomogeneous (e.g., due to vertical gravity/stratification), α ij (z) -as well as the off-diagonal elements of η ij (z)-are expected to be nonvanishing.Together, the tensors encapsulate the emergence of the mean EMF as a response to imposing an external mean magnetic field -or, in general, to the presence of a selfconsistently evolving mean field.The purpose of the present paper is to elucidate a possible finite-time character of this response.
Typically, the turbulent closure coefficients are thought to connect E(z, t) to the mean magnetic field, B(z, t), and its curl ε jzl ∂ z B l (z, t), in a local and instantaneous fashion.2How-ever, in contrast to this instantaneous characterization of the closure relation, the power-law nature of the turbulent cascade suggests that the space-time domain of dependence of E(z, t) is indeed finite -implying so-called "memory effects" (Hubbard & Brandenburg 2009), that is, a delayed (i.e., outof-phase) response to an applied mean field. 3Under the assumption of statistically stationary turbulence, a simple noninstantaneous closure relation (also see Gressel & Elstner 2020) can be formulated as a convolution integral in time of the form In the local box geometry, the integral kernels αij (z, t ) and ηij (z, t ) are functions of the vertical coordinate, z, only.Moreover, in its Fourier-space representation, the above relation can be expressed as a simple multiplication (see Hubbard & Brandenburg 2009, appendix A) with the Fourier transform, αij , of the kernel.That is (dropping the explicit z-dependence), we write attributing a spectral flavor to the mean-field coefficients.While all of these quantities are in general complex functions, we can obtain real values (of, e.g., E i ) in physical space and time by adding up the contributions from positive and negative frequencies.
Complementing the result of Gressel & Pessah (2015) on non-local, scale-dependent character of the mean-field effects in magnetorotational turbulence, we here aim to obtain frequency-dependent closure coefficients, corresponding to convolution kernels in the time domain.The frequency dependence can very naturally be obtained via the test-field (TF) method (Schrinner et al. 2005(Schrinner et al. , 2007) ) employing oscillating test fields as briefly outlined in the next section.

The spectral test-field method
The defining advantage of the TF method, compared with other methods of inference, is that it relies on analytically prescribed "test fields", that can be chosen to span a nondegenerate basis for determining all tensor coefficients in an unambiguous manner.This differs from direct inversion methods (see, e.g., discussion in Bendre et al. 2020), that are founded on the (potentially degenerate) mean fields, B(z, t), developing in the DNS.To invert equation ( 6), and solve for the tensorial closure coefficients, αij (ω) and ηij (ω), we apply the flavor of the method where the TFs, B (µ) (z), are quadruplets of trigonometric functions (also see, e.g., Brandenburg 2005;Sur et al. 2008;Brandenburg et al. 2008): For the purpose of determining the time-response, we here focus on a fixed vertical scale k z = k TF z = 2π/L z , with L z = 6 the vertical size of the box.This makes us sensitive to the coefficients representative of the largest scales available, and we refer the interested reader to section 3.4 of Gressel & Pessah (2015) for a complementary discussion about the scale dependence (also see Rheinhardt & Brandenburg 2012, for the general case of full spatio-temporal dependence).Having specified k TF z , we moreover use eleven spectral modes ω = ω TF = 1 /32, 1 /16, 1 /8, . . ., 16, 32 × 2π/P 0 , centered around P 0 = 1 /4 orbit for the temporal domain.We have arrived at this sampling interval by a combination of educated guessing and trial and error, and have found a posteriori that the relevant dynamic range appears to be covered.
In total, we are hence solving 4×11 = 44 additional induction equations, one for each of the TF fluctuations, B (r, t), alongside the DNS.In terms of the fluctuating velocity, u, these are Importantly, these equations are passive in that they do not influence the evolving magnetic fields in the original DNS.Since we are dealing with MRI turbulence -driven via an underlying genuinely magnetic instability -there likely are preexisting magnetic fluctuations, B 0 , that are statistically independent from the developing (horizontal) mean fields.Such fluctuations, if correlated with the turbulent velocity, may result in an additional EMF, namely E 0 ≡ u ×B 0 -which does, however, not enter our parametrization.To obtain the coefficients, we evaluate the corresponding mean electromotive force E (µ) ≡ v × B (µ) for each of the quadruplets from equation (7).A formal solution to equation ( 6) is then obtained as where the tensors ηij (k z , ω) (i.e., wrt. the current) and βijz (k z , ω) (i.e., wrt. the field gradients) are simply related via (see Hubbard & Brandenburg 2009) ηij = ε jkz βikz . (10) In order to arrive to a statistical sound basis and eliminate bychance fluctuations, equation ( 9) can simply be accumulated in time as needed.This appears to be particularly relevant for slowly oscillating TFs, motivating our preference of long simulations time over grid resolution.Note that in the timedependent case, that is, for ω TF 0, complex coefficients can arise, reflecting a frequency-dependent phase shift of the resulting EMF with respect to the originating oscillating TF.In practical terms, we replace the complex factor e iωt in equation (9) by either cos(ωt) or sin(ωt), in order to project the real / imaginary parts of the coefficients, respectively.In the appendix of Gressel & Elstner (2020), we have benchmarked the implementation of the described spectral TF method, using the simple test case of helical forcing in the strictly kinematic limit.

RESULTS
We present results from a single generic shearing-box simulation of the MRI in the presence of a weak (i.e., midplane β p = 2.2 ×10 5 ) net-vertical field, with L z = ±3 H and outflow boundary conditions and a moderate resolution of 32 grid cells per pressure scale-height.The simulation quickly reaches a steady state with a dimensionless accretion stress (i.e., Reynolds + Maxwell) of about 0.01, and we evolve the simulation for 555 orbits to obtain decent statistics.
In Figure 1, we plot vertical profiles of the dynamo α effect, where we show real (upper panel/red tones) and imaginary part (lower panel/blue tones) separately.The parametric  1, but for the symmetrized off-diagonal element α sym ≡ 1 2 (α φr + α rφ ) related to (differential) vertical turbulent pumping -which we note is distinct from the conventional diamagnetic pumping related to the anti-symmetrized γ ≡ 1 2 (α φr − α rφ ).
curves represent the variation with the imposed oscillation frequency, ω TF , of the TF inhomogeneity. 4ooking at the raw array of curves by eye makes it rather cumbersome to grasp anything but the most fundamental trends in the data.As a remedy, and to illustrate the basic features of the spectral response, we resort to sampling point values at the arbitrary location z = 2.67, roughly corresponding to the peak of the α profile.The real (red /' ') and imaginary (blue /' ') amplitudes thus obtained are shown in the inset of Fig. 1 along with best-fit response functions (solid/dashed, see Sect.3.1, below, for details).
We produce a corresponding plot -shown in Fig. 2-for the symmetric off-diagonal components, α sym ≡ 1 2 (α φr +α rφ ), of the dynamo tensor.One can see that apart from the overall amplitude, which is about a factor of two higher, these broadly match the characteristics of α φφ (z).
Unlike for the case of supernova-driven turbulence in the multi-phase interstellar medium (Gressel & Elstner 2020) -where the off-diagonal elements were found to be antisymmetric and distinct from the diagonal elements-the diagonal and off-diagonal components of the α tensor here (see Figs. 1 and 2, respectively) show a rather similar time response.The imaginary part displays a broad peak around ω = 2.At the same time, the real part has a moderate overshoot around ω = 1, before it reaches the asymptotic value for slowly-varying mean fields.
The general shape of the response can be understood by visualizing the overall character of (rotating) turbulence.Let us briefly recall what the α effect entails.It describes the emergence of a mean turbulent electromotive force as the direct consequence5 of the presence of a large-scale magnetic field.In the limit of high frequencies, this "presence" obviously looses its coherent character and magnetic fluctuations created by the u × B term in Eqn.(8) will tend to become uncorrelated with the velocity and their contribution to the EMF will consequently average out to zero -this is reflected in the vanishing amplitudes towards high frequencies.
Conversely, at low frequencies, we simply approach the previously reported (Gressel & Pessah 2015) amplitudes.In particular, the imaginary part of the effect also tends to zero in this limit, so that the α effect becomes a real number, implying that there is no longer a phase difference.The interesting regime falls in the region of intermediate frequencies, that roughly correspond to the eddy turnover time and/or rotational frequency of the turbulence.Here the finite-time character of the relation between the imposed large-scale field B(t) (as a "cause") and the turbulent E(t) (as a "response") becomes most obvious.In particular, owing to the noninstantaneous build-up of (correlated) magnetic fluctuations via the u ×B term, the imaginary part deviates from zero, implying a phase-lag between cause and effect.Since the α effect is thought to be related to the twisting-up of rising field loops by the rotation (see, e.g., Rüdiger & Kitchatinov 1993, for an analytic calculation of the expected rotational dependence), it appears natural that the effect peaks around ω Ω, where the pulsation of the assumed mean field matches the turnover of eddies that are affected by the Coriolis force.
In comparison with the α tensor, the diffusion coefficient, η T , (shown in Figure 3) shows a somewhat reduced coherence time.This broadly matches the expectation that the mixing aspect of the chaotic flow field depends somewhat less on the buildup of correlated motions and is hence preserved further into the limit of high frequencies.Moreover, the real part of η T (ω) remains strictly monotonic around ω = 1, hinting at a reduced influence of the orbital timescale on the mere random diffusion of field.

Characteristic response function
Pertaining to the non-instantaneous closure relation, and for the case of simple helically-forced turbulence, Hubbard & Brandenburg (2009) have demonstrated a frequency- dependence of the form of an "oscillating decay", where Θ(t) simply denotes the Heaviside step function, enforcing causality by suppressing dependence on future times.Translated into Fourier space, the spectral shape function becomes α with independent coefficients A 0 , τ c , and ω 0 -and with corresponding expressions for the other two coefficients of interest.For the purpose of curve-fitting the frequency response, we write equation ( 12) separated into real and imaginary part as with separate sets of fit parameters A 0 , τ c , and ω 0 for the three coefficients αφφ (ω), αsym (ω), and ηT (ω), respectively.We note that -while we express the complex dependence in terms of two separate functional shapes via equations ( 13) and ( 14)-we fit the real and imaginary branches simultaneously in practical terms.The real part (solid lines) and imaginary parts (dashed lines) of the fitted curves are overlaid in the insets of Figures 1-3, and represent the data rather well.We, moreover, A 0 τ c ω 0 ω 0 τ c dynamo α φφ 0.88 3.92 0.20 0.78 off-diagonal α sym 1.86 4.06 0.20 0.80 diffusivity η T 0.56 3.05 0.17 0.52 Table 1.Fit coefficients for the complex response function defined via equations ( 13) and ( 14), along with the derived parameter ω 0 τ c 1.  1 above.The fit amplitude, A 0 , simply follows the original curves from the previous figures and is omitted here.
report best-fit values for coefficients (sampled at z = 2.67) in Table 1, where we also provide the dimensionless number, ω 0 τ c , which we naively expect to be on the order of unity.As a consequence of the vertical stratification in density, one may conjecture that the largest eddy-size depends on the height in the disc.This in turn should be reflected in the τ c and ω 0 fit coefficients that we determine.To test this assumption, we plot the two numbers in the upper and lower panels of Figure 4, respectively.It can again be seen that the dynamo α term and the off-diagonal elements are very similar, and at the same time, distinct from the diffusive coefficient.While τ c remains fairly constant for the latter, the former two show a pronounced increase (by a factor of three) of the turbulent correlation time towards the disc midplane.In contrast to this, their oscillatory parameter, ω 0 , remains fairly even in that limit but instead shows a moderate peak around z = 1.75 H.As consistent with the monotonic profile seen in Fig. 3, the ω 0 parameter is reduced in the η T coefficient, and even drops to quite small values towards the disc midplane.
What precisely causes the observed trends is unclear at this point, and it is important to keep in mind that MRI turbulence is critically affected by magnetic forces so that intuition from hydrodynamic turbulence may be of limited value.Irrespective of this, comparing mean-field models with and without variation in τ c will allow to establish whether the seen variations do have an impact on the appearance of the butterfly diagram.

DISCUSSION & CONCLUSIONS
The presented results clearly display the finite-time character of the mean-field dynamo effect emerging in stratified zero-net-flux MRI turbulence.The approximate functional form presented in the preceding section will enable to incorporate the effects into a more comprehensive mean-field description of the evolution of large-scale magnetic fields in accretion discs.
As a first step in that direction, one may neglect the ω 0 contribution, related to oscillatory behavior at intermediate frequencies.In this case, equation ( 12) simplifies to with corresponding expression for η T .With the further approximation τ (α) c = τ (η) c (cf.Table 1 for judging to what degree this is justified), the characteristic time, τ c , enters as a relaxation time and implies a time-dependent (i.e., noninstantaneous) EMF response.In contrast to the algebraic relation from equation (4), this needs to be modeled as an extra PDE (also see Rheinhardt & Brandenburg 2012) complementing the mean-field induction equation (3), and where α (0) ij (z) and η (0) ij (z) represent the above "A 0 " coefficient (now measured for each tensor element individually).A very basic complete closure model would arise in combination with the non-locality derived in Gressel & Pessah (2015).Expressed as a characteristic length, l c , relating to the finite "domain of dependence", the simplest form for the right-hand-side would then become where the l c appears as a smoothing term.This, however, still neglects the possible advection (i.e., via a term u • ∇ E i ) of the EMF with the disc outflow, u z , as well as potential effects related to ∇•u = ∂ z u z .These may act to (de-)compress the EMF -similar to the (−∇ • u) B contribution to the ∇ × (u × B) term in the induction equation itself.Moreover, if one were to restore the effect related to ω 0 , one would obtain a wave-like second-order time derivative of the EMF on the left-hand-side, as well as time derivative terms related to B appearing on the right-hand side (M.Rheinhardt, private communication).A straightforward Cranck-Nicolson discretization of equation ( 17) has already been implemented into a simple dynamo code.In view of the mentioned complications, we however defer a detailed mean-field treatment along these lines to a later point in time.
Coming back to the inference of the dynamo coefficients from DNS, an often mentioned shortcoming of the TF method pertains to the absence of magnetic fluctuations stemming directly from the simulation.Operating in the socalled "quasi-kinematic" realm (also see discussion in Gressel & Pessah 2015), the QK-TFM is agnostic to the presence of a possible E 0 , and the velocity v(r, t) is the only manifest link in equation ( 8) to the physical evolution traced by the DNS.While this may indeed be seen as a shortcoming of the current approach, we highlight that it merely implies that the detected mean-field effects likely are not exhaustive, but simply restricted to the chosen ansatz.A promising avenue to accounting for these additional contributions has first been laid out by Rheinhardt & Brandenburg (2010) for a simplified set of equations (i.e., lacking the pressure and self-advection terms in the momentum equation).More recently, a workable solution has been found also for the complete set of MHD equations (Käpylä, Rheinhardt & Brandenburg 2021).It appears natural to test this approach for MRI turbulence as well.

Figure 1 .
Figure1.Profiles of the real part (upper panel /' ') and imaginary part (lower panel / ' ') of the dynamo α φφ (z) coefficient as a function of height in the disc.Each curve represents one specific value of ω TF , and the solution sampled at z = 2.67 is plotted in the inset as a function of angular frequency, ω TF , of the pulsating TFs, spanning three decades in dynamic range.

Figure 4 .
Figure 4. Height dependence of the two fit parameters from equation (12) describing the temporal behavior of the response.Circles indicate the values shown in Table1above.The fit amplitude, A 0 , simply follows the original curves from the previous figures and is omitted here.