Toward 2D Dynamo Models Calibrated by Global 3D Relativistic Accretion Disk Simulations
Matthew D. Duez Courtney L. Cadenhead Zachariah B. Etienne Bernard Kelly Leonardo R. Werneck







             

§ ABSTRACT 
Two-dimensional models assuming axisymmetry are an economical way to explore the long-term evolution of black hole accretion disks, but they are only realistic if the feedback of the nonaxisymmetric turbulence on the mean momentum and magnetic fields is incorporated.  Dynamo terms added to the 2D induction equation should be calibrated to 3D MHD simulations.  For generality, the dynamo tensors should be calibrated as functions of local variables rather than explicit functions of spatial coordinates in a particular basis.  In this paper, we study the feedback of non-axisymmetric features on the 2D mean fields using a global 3D, relativistic, Cartesian simulation from the IllinoisGRMHD code.  We introduce new methods for estimating overall dynamo alpha and turbulent diffusivity effects as well as measures of the dominance of non-axisymmetric components of energies and fluxes within the disk interior.  We attempt closure models of the dynamo EMF using least squares fitting, considering both models where coefficient tensors are functions of space and more global, covariant models.  None of these models are judged satisfactory, but we are able to draw conclusions on what sorts of generalizations are and are not promising.



Toward 2D Dynamo Models Calibrated by Global 3D Relativistic Accretion Disk Simulations
    Leonardo R. Werneck
    September 9, 2024
=======================================================================================





§ INTRODUCTION



Untilted black hole accretion disks are three-dimensional, turbulent systems, but their fluid and field profiles can be viewed as combinations of an axisymmetric, slowly evolving “background” and non-axisymmetric “fluctuations.” Often, only the background is of interest, with fluctuations considered primarily for their feedback effects on the background. This perspective motivates the use of 2D (axisymmetric) simulations to model accretion flows.

Simply evolving the general relativistic magnetohydrodynamic (GRMHD) equations in 2D will almost certainly be inadequate; MHD flows are often not only quantitatively but even qualitatively different in 2D versus 3D, as illustrated by opposite turbulent energy cascades <cit.> and the 2D anti-dynamo theorem <cit.>.  What we aim for are 2D evolutions that “look like” 3D evolutions.  Note that even this is ambiguous.  Do we want 2D evolutions that resemble the azimuthal average of 3D evolutions, or do we want 2D evolutions that resemble representative 2D slices of 3D evolutions?  Both approaches might be interesting, but they will likely be very different.  Consider an evolution variable Ψ(r,θ,ϕ,t) with azimuthal Fourier decomposition ∑_mΨ̂_m(r,θ,t)e^imϕ.  Azimuthal averaging picks out the m=0 contribution, and if this is subdominant, the azimuthal average Ψ will be significantly smaller than a 2D slice of a 3D realization.  We might also expect most of the (r,θ) fluctuation to average out, so that azimuthal averaging might act as a low-pass filter, removing most eddies.  Given these differences, one must choose when designing a model whether the 2D fields represent azimuthal averages or a representative slice.  For this paper, we consider 2D models of azimuthally averaged 3D evolutions.

One might also worry that the distinction of background versus fluctuating turbulence is too sharp.  In terms of the Fourier decomposition of fields in the azimuthal direction, it seems safer to distinguish high-m modes from the background than, for example, m=1 or m=2 modes.  A possible criterion would be to consider an m 0 mode is “turbulence” only if the eddy turnover time is less than the orbital period.  For simplicity, consider unmagnetized incompressible turbulence with turbulent velocity ∼v_T at the largest length scale H (the disk scale height) and the usual Kolmogorov energy cascade. At scale ℓ, the velocity scale is v_ℓ and the timescale is τ_ℓ≈ℓ/v_ℓ.  As specific kinetic energy density flows to lower scales at the same rate ≈ v_ℓ^2/τ_ℓ for all ℓ, it follows that τ_ℓ≈ H^1/3ℓ^2/3v_T^-1.  Turbulent modes then have τ_ℓ<Ω^-1, where Ω is the orbital frequency.  Estimating m≈ R/ℓ at radius R, the eddy time is then sufficiently short if

    m > m_ crit≈ R H^1/2(Ω/v_T)^3/2 .

It might be necessary to evolve m<m_ crit (a compromise between 2D and 3D, assuming m_ crit turns out to be a small number), with the m>m_ crit contribution safely modeled by a turbulent/dynamo closure.


Filtering away the fluctuating component alters the evolution equations, as can be seen by azimuthally averaging the Newtonian induction equation:


    ∂_t𝐁   = ∇× (𝐁×𝐯) 
       = ∇× (𝐁×𝐯) 
       ≡∇× (𝐁×𝐯 + Δ 𝐄) .


To evolve these equations, we require a closure condition, providing the dynamo EMF Δ 𝐄 as a function of azimuthally averaged variables.  A commonly considered closure <cit.> is


    Δ E^i = α^i_j B^j + η^i_j J^j


where J^j is the current, B^j is the magnetic field, α^i_j is the dynamo alpha tensor, and η^i_j is the turbulent magnetic diffusivity tensor.

Dynamo corrections to 2D GRMHD simulations have been undertaken by a number of groups <cit.>, usually assuming an isotropic dynamo:  α^i_j=α_ dynδ^i_j, η^i_j=η_ dynδ^i_j.  The α_ dyn term enables the dynamo alpha effect:  toroidal field produces poloidal field.  The alpha effect enables 2D GRMHD simulations to maintain magnetic fields at the strength seen in 3D simulations.  (Note that using it for this reason implies a representative 2D slice 
rather than an azimuthal average interpretation.)  The η_ dyn term
is the turbulent magnetic diffusivity, which acts like a resistivity.  As α_ dyn is a pseudoscalar, it is expected to switch sign at the equator, so some such spatial dependence on α_ dyn is often stipulated. 

Dynamo terms are also sometimes added to 3D GRMHD simulations <cit.>.  In this case, the distinction between background and fluctuation (and hence the meaning of averaging) is done differently—the fluctuation is the MHD flow at spatial scales below the grid scale.  It is not clear a priori if the same Δ 𝐄 model should apply to the two cases.

Three-dimensional MHD simulations can be used to extract α^i_j and η^i_j.  Gressel and Pessah <cit.> calculated these using a test-field method. They studied shearing box simulations and averaged over boxes, so their distinction of mean versus fluctuating fields was not quite the same as the azimuthal average, but they were able to study dependence on shear and vertical field. Assuming α^i_j=α_ dynδ^i_j and η^i_j=0, Hogg and Reynolds <cit.> used global MHD simulations to extract α_ dyn in each hemisphere. Bendre  <cit.> introduced a method of extracting α^i_j and η^i_j by least-squares fitting via singular value decomposition (SVD), which they applied to radiatively inefficient accretion disks in Dhang  <cit.>. Both test-field and SVD methods find deviations in dynamo tensors from isotropy and report fairly consistent magnitudes. Dhang note the difficulty of extracting a clear dynamo signal in the disk region, as well as other indications that the field and flow are mostly non-axisymmetric.  The SVD method has also been used for solar dynamo and binary neutron star remnant calculations <cit.>.

The recent study by Jacquemin-Ide  <cit.>, using GRMHD disk simulations to study dynamo action from the nonlinear evolution of the magnetorotational instability, also presents some relevant findings.  They reported that m=1 to m=3 modes are crucial in dynamo generation of poloidal field (perhaps suggesting a suitable m_ crit). 
They also emphasize the importance of magnetic flux advection from the outer to the inner disk.  This might be a less prominent effect in the small disks (resembling short gamma ray burst setups) often studied in compact binary merger contexts.  Finally, they do estimate the viability of an α_ dyn and η_ dyn dynamo closure using correlation coefficients between Δ E_ϕ and B_ϕ or J_r.  The average Δ E_ϕ-B_ϕ correlation oscillates about zero and does not show a significant time-average; the Δ E_ϕ-J_r correlation does have a time average clearly distinct from zero, but is strongest with a time offset between cause and effect, suggesting it is mediated by additional dynamics.

The number of studies that extract dynamo coefficients from global disk simulations, particularly relativistic simulations, remains small given the parameter space of disk states, and not all existing results have high statistical significance.  Furthermore, dynamo coefficient extraction is of two types:  either a single average coefficient is extracted for the entire disk (or a hemisphere thereof), or dynamo coefficients are extracted at each (r,θ) point in space.  These spatial dependencies are presumably proxies for dependence on the local physical quantities and their derivatives.  If these dependencies could be made explicit, the resulting model for Δ E^i with only a few global fitting constants (as opposed to fitting scalar functions of space) could be expressed in covariant tensor form, making it immediately generalizable to general spatial coordinate systems and more easily generalizable to other accretion systems (varying disk thickness, magnetic flux, and rotation law).

In this paper, we present a new analysis of the azimuthally averaged evolution of a magnetorotationally turbulent disk around a Kerr black hole.  We use data from the Cartesian numerical relativity code IllinoisGRMHD <cit.>.  We present a detailed analysis of the degree of non-axisymmetry, introducing measures of its dominance in magnetic energy, kinetic energy, and angular momentum transport.  We find that most of the magnetic field and non-azimuthal velocity field are averaged out by azimuthal averaging, and momentum transport is predominantly non-axisymmetric. Thus, azimuthally averaged 3D ideal GRMHD resembles viscous hydrodynamics more than 2D ideal GRMHD.  We also propose ways of estimating an average α_ dyn and η_ dyn.  An elegant definition of the former comes from the Lorentz invariants of the azimuthally averaged field tensor.  For the latter, we take advantage of the fact that, at least for the first 10000 M, 2D ideal GRMHD overpredicts the magnetic energy compared to azimuthally averaged 3D GRMHD, so we estimate η_ dyn by the resistivity needed to achieve this level of field suppression.  For this purpose, we add a new phenomenological resistivity to the HARM code and compare with the 3D results.

Next, we attempt to extract α^i_j and η^i_j using SVD least-squares fitting.  We consider models for which coefficients are functions of (r,θ), functions of θ, or functions of local scalars and pseudoscalars (with spacetime dependency only from those scalars and pseudoscalars).  We introduce norms for the error in the best-fit model and variance in coefficient extraction to easily assess the reliability of models.  None of the models meet all of our pre-set standards.

Since these models assuming eq:NewtonianDynamo are not judged satisfactory, we assess certain alternative classes of models which maintain the same basic structure of the dynamo EMF. 
 First, we consider whether a successful global model could be created if the RMS deviation from axisymmetry of the velocity or magnetic field were known.  To be usable as a closure condition, an evolution equation for one of these variables that tracks the true evolution well enough would have to be introduced, similar to the evolution of the mean turbulent kinetic energy in k-ϵ models <cit.> (which might be generalized to MHD <cit.>).  In fact, such scalars do not provide acceptable dynamo models, so the motivation to devise such evolution equations does not arise. Second, we ask whether the residual EMF subtracting off the EMF computed from m=0 to m=3 components of 𝐯 and 𝐁 is more amenable to dynamo models of the form eq:NewtonianDynamo.  It is not.  We conclude that more general models for Δ E^i are likely required.

This paper is organized as follows.  In Sec. <ref>, we review the GRMHD and dynamo equations.  In Sec. <ref>, we provide details on the 3D simulation, in particular analyzing the contributions of non-axisymmetric fields.  In Sec. <ref>, we describe how to produce estimates of the isotropic α_ dyn, η_ dyn.  In Sec. <ref>, we describe how we produce fits to α^ij, η^ij and present quality of fits.  We summarize and present conclusions in Sec. <ref>.

We use units G=c=1 throughout.  Latin letters from the beginning of the alphabet (a–f) are spacetime indices, running 0–3.  Indices i, j, k are spatial, running 1–3.  When writing equations in geometric form, vectors and forms are written in boldface, with 2-forms getting a “2” superscript prefix (e.g., 2𝐁).



§ AN ANALYTIC DYNAMO MODEL



For use in a relativistic code, the dynamo closure condition must at least be spatially covariant; this will account for the unpredictable evolution of the coordinate system or the deliberate use of curvilinear coordinates or coordinates with designed radial or angular concentrations.  Azimuthal averaging introduces a preferred direction, the azimuthal Killing vector 𝐞_ϕ≡∂/∂ϕ, which generally will leave an imprint in the α and η tensors, so that 3D covariant equations for these tensors will be expected to include 𝐞_ϕ.  Whether the equation must be 4D covariant, i.e., Lorentz covariant, is more debatable.  It could be argued that the averaging procedure itself, which is defined to be azimuthal average at a fixed time, breaks Lorentz invariance.  Indeed, subgrid and effective viscosity prescriptions that are purely spatial, which have been widely and successfully used in numerical relativity <cit.>, can have similar justifications.  On the other hand, in one case, we found it important to keep terms first-order in velocity to properly recover the Newtonian limit <cit.>. 

Assume a foliation of the spacetime with the usual normal to the slice ñ = -N 𝐝𝐭 and 3-metric γ_ab = g_ab + n_a n_b.  The field tensor associated with the azimuthally averaged fields is


    F^ab   = n^a E^b - n^b E^a + ϵ^abcB_c  ,
    ⋆ F^ab   = B^a n^b - B^b n^a + ϵ^abcE_c ,


where E· n = B· n = 0, and ϵ^abc=n_dϵ^abcd.  (Here ϵ^abcd=|g|^-1/2[abcd] is the Levi-Civita tensor and [abcd] is the totally antisymmetric Levi-Civita symbol.)  For the rest of this section, we will suppress lines above averaged quantities, assuming that all quantities are azimuthal averages. 

The fluid frame is defined by the azimuthally averaged 4-velocity u^a, which can be decomposed as ,
where W is the Lorentz factor and 𝒱, defined so 𝒱^an_a=0, is the Eulerian velocity, not the transport velocity.  For the rest of this section, we will suppress lines above averaged quantities, assuming that all quantities are azimuthal averages.

The electric field uE^a in the fluid frame,
which is also the Lorentz force, is


    uE^a≡ F^ab u_b = W E^a  + n^a E^b 𝒱_b + ϵ^abc𝒱_b B_c .


In ideal MHD, this is zero, and the components parallel and perpendicular to 𝐧 must independently vanish, so


    MHDE^a = -1/Wϵ^abc𝒱_bB_c .


One can also compute the electric field using the transport velocity 3-vector v^i≡ u^i/u^t (cf. <cit.>):


    E_i = -1/Nϵ_ijk(v^j+β^j)B^k ,


where N is the lapse and β^i is the shift.


In the fluid frame, the magnetic field field is


    uB^a = -⋆F^ab u_b
        = W B^a + B^b𝒱_b n^a - ϵ^abc𝒱_b E_c .


A covariant expression that recovers the Newtonian limit (Eq. <ref>) is


    F^abu_b = -α^a_b uB^b
        + η^a_b ϵ^cbdeu_c∇_d uB_e .


Note that, in the last term, we can switch freely between ∇ and ∂ since d and e are antisymmetrized.

Another possibility would have been to use the actual current

    F^abu_b = -α^a_b uB^b
        + η^a_b uJ^b ,


which is different from the previous equation because of the displacement current ℒ_𝐧𝐄 term in Ampere's law.  However, the identification of the diffusivity term with ∇×𝐁 should be considered more fundamental than its identity with 𝐉, because it is motivated by an expansion in terms of one over the length scale of the magnetic field.  Hereafter, we will simply denote the curl of 𝐁 as 𝐉.

Unfortunately, eq:uB itself has the electric field in it.  Luckily, this term is proportional to velocity and subdominant in most regions. For a dynamo equation with only an isotropic α tensor (η^a_b=0), Most has shown that one can solve for E^a analytically without needing to assume small velocity <cit.>.  Instead, we will substitute the ideal MHD electric field (<ref>) in the η term, so that the weak dynamo case is properly recovered.  Then


    uB^a   = W B^a + B^b 𝒱_b n^a - ϵ^abc𝒱_b (-W^-1)ϵ_cef𝒱^eB^f 
       = B^b𝒱_b n^a + 1/W[B^a + B^b𝒱_b 𝒱^a] ,


where we have used W^2=1+𝒱^2.
As for the curl term in (<ref>)


    ϵ^abcdu_a∂_c uB_d = -W ϵ^bcd∂_c uB_d + ϵ^abcd𝒱_a ∂_c uB_d .


The first is just a spatial curl, albeit of the boosted B, which is easy to compute.  The second term is, in general, difficult to calculate because it will include time derivatives which are less easily available to an evolution code or post-processing analysis.  However, in most cases, this term will be small.  Suppose we are willing to limit ourselves to cases in which 𝒱^i and n^i are small enough (implying that the shift vector is small, which is usually the case when not very close to a black hole for commonly chosen gauge conditions) that we can keep only terms at lowest order in these.  Then 𝒱_0=B_0=0, and for spatial components of the free index b the index c must be the one that is the time component.  In the weak-dynamo limit, the time derivative of B^a is proportional to 𝒱^i, and for slowly-evolving metrics, ∂_tB_a is proportional to ∂_t B^b.  So this term, which looks first-order in 𝒱, is actually second-order, and will be dropped with slightly easier conscience.  We also drop the term B^b𝒱_b n^a in uB^a (see Eq. <ref>) as second-order.

For the component of Eq. (<ref>) perpendicular to 𝐧, and using (<ref>), we are left with


    W E^a =    -ϵ^abc𝒱_b B_c
        - α^a_b W^-1 [B^b + B^c 𝒱_c 𝒱^b] 
        - η^a_b W ϵ^bed∂_e [(B_d + B^c𝒱_c 𝒱_d)W^-1] .




§ NON-AXISYMMETRIC FEATURES OF THE 3D SIMULATION





 §.§ Methods of 3D simulation


The simulation we used was produced as part of the Event Horizon GRMHD code comparison project <cit.>.  The fixed spacetime is a Kerr black hole with dimensionless spin parameter a/M = 0.9375.  The gas is an ideal gas with adiabatic index Γ=4/3.  The initial state of the gas has constant specific angular momentum, defined as ℓ = u^tu_ϕ.  In the absence of magnetic forces, the gas forms an equilibrium torus with an inner radius of 6M and maximum density at 12M.  We seed a magnetic field at t=0, with the maximum of the magnetic pressure P_B^ max and the maximum of the gas pressure P_ gas^ max related as P_B^ max/P_ gas^ max=100.  The field is calculated from a toroidal vector potential A_ϕ∝ max(ρ/ρ_ max-0.2,0).  This initial state produces a “Standard and Normal Evolution” (“SANE”) accretion flow, i.e., one for which magnetic flux on the horizon does not significantly affect accretion. 

The IllinoisGRMHD <cit.> simulation uses a Cartesian FMR grid (utilizing the Cactus/Carpet infrastructure), in which the coarsest grid is a cube with a half-sidelength of 1750M. There are 7 different resolutions used, including 6 levels of refinement atop this coarse grid, and the finest-grid cube has a resolution of Δ x=Δ y=Δ z ≈ M/4.388571. The finest grid is composed of four cubes, each with a half-sidelength of 27.34375M, and is larger than the black hole, which has a radius of 1.348M. IllinoisGRMHD evolves a vector potential at cell edges to exactly preserve a finite difference version of the constraint ∇· 𝐁=0 <cit.>. The magnetized fluid is evolved using a high-resolution shock-capturing scheme, with PPM reconstruction and an HLLE approximate Riemann solver.

The disk was evolved for a time of 10000M, long enough to see many orbits of fully-developed turbulence.  Volume data of all evolution variables were stored on a  spherical-polar grid.  Data was output at a high cadence of 0.57M, but we find it sufficient to use only every Δ t=100M for time averages below.  When performing fits for dynamo coefficients, we use timesteps separated by Δ t=300M, where the increased spacing—roughly corresponding to the orbital period at the initial location of maximum density—is chosen to reduce correlation between data at adjacent steps while retaining sufficient data (31 steps), starting at t=1000M, after the magnetic field has become strong and the disk transition from the kinematic to the dynamical regime.

The evolution time is admittedly shorter than usual for dynamo studies.  Ordinarily, one would want to observe multiple cycles of the expected dynamo wave.  Assuming an αΩ dynamo with α_ dyn∼ 10^-4, a characteristic vertical scale H∼ 10M, and a characteristic orbital frequency Ω∼ 10^-2M^-1, we would expect a dynamo period ∼√(H/(Ωα_ dyn))∼ 10^4M <cit.>, meaning the full evolution time is not more than a cycle, and we will not be able to produce “butterfly” diagrams of the average toroidal field at different latitudes showing multiple dynamo cycles, as has been done in some longer disk simulations (e.g. <cit.>). 
We will, however, evolve comfortably long enough to observe Δ 𝐄 acting.



 §.§ Overview of 3D data




We define the following averages


    X(r,θ,t) ≡1/2π∫_0^2π  dϕ X(r,θ,ϕ,t) , 
    ⟨ X⟩(r) ≡∫_T^2Tdt∫_π/3^2π/3  dθ∫_0^2πdϕ√(γ) X(r,θ,ϕ,t)/∫_T^2Tdt∫_π/3^2π/3 dθ∫_0^2πdϕ√(γ) ,


where T=5000M.

To minimize effects of the metric on the spatial dependencies of vector and tensor components, we will use the normed components X^î≡√(γ_ii)X^i and X_î≡√(γ^ii)X_i when generating plots and fitting formulas. (Note that we drop the Einstein summation convention in the above definitions, so there is no sum over i.)

For a quantity X(x_A) which is a nonlinear function of the primitive variables x_A=(B^i,v^i,ρ,P), define the residual as


    Δ X = X - X(x_A) .


Suppose instead of retaining only the average of the spherical-polar components of 𝐯 and 𝐁, one were to retain the azimuthal Fourier decomposition, truncated at a particular m=m_t.  Call Ψ_m_t the approximation of Ψ with all azimuthal Fourier modes of m<m_t.  Then define the residual EMF Δ_m𝐄≡𝐁×𝐯 - 𝐁_m×𝐯_m.  Thus, Δ𝐄 = Δ_0𝐄.





In Fig. <ref>, we plot azimuthal and time averages of Δ E^î, B_î, and J_î.  A linear relationship between Δ E^î and B_î or J_î is not immediately suggested to the eye.  Of the three vectors, the B_î components are spatially smoothest, while J_î show small-scale behavior.  B_î are strongest near the poles, J_î in the disk, and Δ E^î in between.

In Fig. <ref>, we plot azimuthal and time averages of the components of the transport velocity v^î and the EMF residual beyond m=3, denoted Δ_3E^î in accord with the notation introduced above.  The structure of Δ_3E^î is similar to Δ E^î, suggesting that the former (representing higher-m modes) significantly contributes to the latter.  Velocity is dominated by the disk rotation and the polar outflow; little of the disk turbulence survives the averaging process.  We also plot three scalar quantities which might be expected to correlate with Δ_3E^î (but do not), whose discussion we defer to Sec. <ref>.




 §.§ Non-axisymmetric energies and fluxes


The azimuthal structure and the effect of azimuthal averaging are illustrated in Fig. <ref>, in which B^r and E_r are plotted as functions of θ and ϕ at one radius and time.  The azimuthal average is shown as a bar on the right.  We see that, at least for these quantities, azimuthal mean structure dominates only near the poles.  At intermediate latitudes (i.e. in and near the disk), there is small-scale structure (high m contribution), which averages out to value much smaller than the maximum magnitudes. 

As in <cit.>, we compute the Maxwell stress averaged over time and angle ⟨ w^rϕ⟩ as a function of r from the full magnetic field in 3D in two ways:  by integrating b^r̂b^ϕ̂ in Kerr-Schild coordinates and by integrating in a local frame comoving with the fluid with appropriate comoving volume element. (See <cit.> for details.)  We use the frame of the azimuthally averaged velocity for easier comparison with the stress from the azimuthal mean field.  [Of course, the azimuthally averaged velocity is itself a nonlocally-defined field, but once computed it has a local value at each point in 3D or 2D (r, θ) space.]  Like in <cit.>, we find that the two prescriptions give almost identical results.  As is commonly done, we define a Shakura-Sunyaev α_ SS by


    α_ SS≡⟨ w^rϕ⟩/⟨ P_G+P_B⟩ ,

where P_G and P_B are the gas and magnetic pressure, respectively.

The result, included in Fig. <ref>, closely resembles the left panel of Fig. 19 in <cit.>, especially the dip at smaller r in the lower resolution results, which is a new check on the consistency of IllinoisGRMHD with the codes with spherical-polar grids.  We also compute the Maxwell stress from the azimuthally averaged field 𝐁, using the same MHD formulas and in the azimuthally averaged comoving frame.  This is found to provide a much lower stress except very near the black hole, indicating that most of the angular momentum transport in the bulk of the disk accomplished by the magnetic field is done so by the non-axisymmetric field (the fluctuations, not the mean), a point also mentioned by Dhang  <cit.>.  Of course, using an ideal MHD formula for the stress tensor is not quite appropriate for the azimuthally averaged field tensor, which includes the effect of Δ 𝐄, but this residual EMF is itself a measure of the prevalence of nonaxisymmetry in the field and so does not affect the conclusion.

Because of the difference between the local 3D and azimuthally averaged velocity fields (𝐯 vs 𝐯), there is also a mean Reynolds stress which can transport angular momentum.  This is calculated in the frame comoving with the local azimuthally averaged velocity field:  t^r̂ϕ̂ = ρ h u^r̂u^ϕ̂.  (Note that the azimuthal mean velocity does not contribute, since if the field were axisymmetric, u^î would be zero in the frame comoving with the azimuthal average velocity.)  This Reynolds component is lower than the full (primarily non-axisymmetric) Maxwell component but still larger than the Maxwell stress of the mean field.

In addition to the shear-type stresses quantified by α_ SS, there can be a pressure-like contribution to the stress tensor from the non-axisymmeric Reynolds and Maxwell stresses.  We define the turbulent pressure as


    P_t ≡1/3Δ T^ab(g_ab+u_au_b) .


This is a very rough definition—applied instead to the stress-energy tensor of ideal MHD, it would give P + b^2/6, where b^2=b_ab^a and 𝐛 is the magnetic field in the 𝐮 frame.  We find this pressure to be close to the magnetic pressure (of the full field), which even after saturation is about a factor of 10 lower than the gas pressure.  (See Fig. <ref>, discussed below, which plots the associated energy densites.)  Thus, although turbulence dominates the angular momentum transport, it is very subdominant in pressure effects (e.g., determining the disk height).

We also compute time and angle-averaged energy densities in the disk.  For the full magnetic energy density, this is b^2/2, computed using the 3D local 𝐁 and 𝐮.  For the magnetic energy of mean field, we use the same formula (justified as above for stresses) using azimuthally averaged quantities.  The internal energy density of the gas is P/(Γ-1).  The kinetic energy of a fluid with 4-velocity 𝐮_1 relative to a frame 𝐮_2 can be defined as follows.  Start from the Reynolds term in the fluid stress tensor T_R^ab(𝐮_1)≡ρ h u_1^au_1^b.  The desired energy is in the 𝐮_2 direction, so we must project T_R^abu_2a.  For an energy density, we must also specify a 4-velocity indicating the direction of flux of this energy we are considering; this is the frame in which the spatial integral converting energy density to energy would be performed.  The natural choices are 𝐮_1 (for comoving frame) and 𝐧 (for regular spatial integration on the slice), and we choose the latter.  The Reynolds energy density is then u_2an_bT_R^ab(𝐮_1)=-(𝐮_1·𝐮_2)Wρ h.  We then subtract off the rest contribution when 𝐮_1=𝐮_2 to get

    e_K = -[(u_1· u_2) + 1]Wρ h .


There are three relevant 3-velocity fields in this problem: 1) 𝐯, the 3D local velocity; 2) its azimuthal mean 𝐯; and 3) the rotational velocity, which is 𝐯 with radial and meridional components removed.  From each 3-velocity field, a 4-velocity field can be constructed by imposing the normalization 𝐮·𝐮=-1 (meaning that the t component of the azimuthal mean 4-velocity is not exactly the azimuthal mean of the t component of the local 4-velocity).  We thus compute a turbulent kinetic energy density, comparing the 3D local velocity to the azimuthal average, and a rotational kinetic energy, comparing the azimuthal average velocity to the rotational velocity.  Energies are plotted in Fig. <ref>.  We omit the rotational energy, which dominates all others, as would be expected for an accretion disk.  The figure shows that the magnetic energy saturates more than an order of magnitude below the internal energy. However, the magnetic energy of the mean field saturates at a much lower value than this total magnetic energy.  This indicates that the magnetic energy is mostly non-axisymmetric (“turbulent”), and we infer that the energy in the turbulent magnetic field is close to the total magnetic energy.  The kinetic energy from the difference of azimuthal mean vs actual velocities (the turbulent kinetic energy) dominates over the non-rotational kinetic energy from the mean field, indicating that inside the disk eddy motion dominates over mean poloidal motions.  The total magnetic energy is seen to be very close to the turbulent kinetic energy, indicating that the kinetic and magnetic turbulent energies come into equipartition with each other.










§ ISOTROPIC DYNAMO MODELS





 §.§ Estimating an isotropic alpha effect


If one knows that the dynamo EMF is of an isotropic alpha form, extracting the pseudoscalar coefficient α_ dyn becomes straightforward.  Starting from the Newtonian formula 𝐄=𝐁×𝐯 + α_ dyn𝐁, one can take the dot product of both sides with respect to 𝐁 to get α_ dyn = E· B/B^2.  Remembering that B^2≫ E^2 (see Fig. <ref>), we can recognize this as the ratio of the two Lorentz invariants of the electromagnetic field.  We may, then, straightforwardly generalize the Newtonian formula as


    α_ dyn≡1/2⋆F_abF^ab/F_abF^ab .


a pseudoscalar, as expected.  Note that these are the invariants of the azimuthally averaged field tensors, not averages of the invariants of the field in 3D.

Because α_ dyn is expected to change sign across the equator, it would be inappropriate to compute an average over angles at a given radius, or an average over the entire grid (unless one were to presume this parity and perform a sign flip in one hemisphere); it is better to average over time and radius to extract an average α_ dyn as a function of θ.  This is plotted in Fig. <ref>.



A hemispherical average pseudoscalar α_ dyn had previously been extracted from a 3D disk simulation by Hogg and Reynolds <cit.>.  They find α_ dyn in the range 1–2× 10^-4.  In magnitude, this is quite consistent with Fig. <ref>.  There has been some disagreement over the sign of α_ dyn, with some finding it positive in the upper hemisphere and some finding it negative (always with the other hemisphere being opposite sign).  Our α_ dyn is negative in the Northern hemisphere, consistent with <cit.>.



 §.§ Calibrating isotropic diffusivity effect with 2D MHD


To estimate an effective η_ dyn, we use the fact that the energy of 𝐁 inside the disk in the 3D simulation is smaller not only than the magnetic energy in the disk of the non-averaged field in 3D, but also than the magnetic energy for a 2D MHD simulation of the same disk.  One could say that the difference in energy is transferred from the mean field to the non-axisymmetric field.  Since this field is not captured by the fields tracked in a 2D simulation, it is in a sense equivalent to internal energy.  The nonaxisymmetric effect that accomplishes this effective dissipative transfer thus acts as a diffusivity, and one can estimate an effective η_ dyn by finding what level of resistivity is needed to reduce the disk magnetic energy to the same degree.

We wish to implement a version of eq:relativistic_dynamo with α^i_j=α_ dynδ^i_j and η^i_j=η_ dynδ^i_j in a relativistic 2D MHD code.  Such codes usually evolve the magnetic field 2-form (or, equivalently, the densitized magnetic field vector) ^2B_ij≡ϵ_ijkB^k ≡ [ijk]B̃^k.

Including an effective resistivity will introduce second spatial derivatives in the induction equation.  We wish to obtain a spatially covariant equation while avoiding Christoffel symbols to the extent possible.  Thus, we continue with form notation where possible, although the spatial metric necessarily will be invoked in Hodge duals and in creating the B-field 1-form ^1𝐁=B_a dx^a.  We introduce an EMF 1-form 𝐑 to include all deviations from ideal MHD.  Continuing in form notation, the dynamo-modified induction equation is


    ∂_t2𝐁=-ℒ_𝐯2𝐁 + d𝐑 .


For the standard isotropic dynamo, 𝐑 is


    𝐑^ dyn=α_ dyn1𝐁 + η_ dyn⋆ d 1𝐁 ,


where in fact we use u𝐁 as defined in Eq. <ref> for the magnetic field.  In components,


    R^ dyn_k    =   α_ dynB_k + η_ dyn1/2ϵ_ijkA^ij , 
    
      A^ij   =   γ^iℓγ^jm(∇_ℓ B_m - ∇_m B_ℓ) 
       =   γ^iℓγ^jm(∂_ℓ B_m - ∂_m B_ℓ) .


(I.e., by using exterior derivatives, we avoid covariant derivatives.)  In order to avoid the timestep limitations of explicitly evolving a parabolic system, we promote R_k to an evolution variable driven to its dynamo form:


    ∂_t R_k = τ_ drive^-1(R^ dyn_k - R_k) .


We set τ_ drive^-1=τ̂_dΩ_K, where Ω_K is the Keplerian angular frequency; we report results for τ̂_d=0.5 but have checked that the behavior of all quantities is insensitive to it within a wide range.  Smaller τ_ drive will force R_k to track R^ dyn_k more closely, but we wish to keep η_ dyn/τ_ drive < 1 to avoid acausal propagation speeds.

The coefficient η_ dyn is set to c_sℓ, where ℓ is an effective mixing length; it is natural to suppose that it is similar in magnitude to the mixing length used to set the alpha effective viscosity.  This similarity is confirmed for shearing box studies (more precisely, that the turbulent magnetic Prandtl number is roughly unity <cit.>), but it should be remembered that the distinction between mean and turbulent field is different in these studies, so their applicability to azimuthally averaged dynamics cannot be certain.  We set


    ℓ = α_η c_sΩ_K^-1f_ρ f_β f_ BH ,


where α_η is a free constant (the symbol “α” chosen for its analogy to α_ SS, not to α_ dyn), and the factors f_ρ=ρ/(ρ+10^-4ρ_ init,max), f_β=P_ gas/(P_B+P_ gas), f_ BH= max[0, 1-(r/3M)^-2] suppress resistivity in the low-density region, the magnetosphere, and near the horizon, respectively.

We implement the above in the harmpi code, a version of the GRMHD code HARMPI <cit.> publicly available at github <cit.>.  Note that, unlike truly 4D covariant resistive MHD codes, we do not evolve the actual electric field, although making 𝐑 an evolution variable arguably accomplishes something similar.  Our less-precise implementation does have the advantage that no problem of stiffness arises in the limit of low η_ dyn.



In Fig. <ref>, we show the magnetic energy in the disk as a function of time for several simulations.  For the 3D simulation, we plot both total magnetic energy and the (much smaller) magnetic energy of the mean field.  We also show harmpi in 2D with ideal MHD, which produces a magnetic energy far too high during the evolution period to match the mean of the 3D field.  Better agreement can be achieved by adding an isotropic dynamo η term with α_η=0.018.



In Fig. <ref>, we plot magnetic flux Φ≡∫ |B^r| dA on the horizon.  For electromagnetic energy extraction and jet formation, it is important that suppressing the disk magnetic energy does not also suppress the horizon flux.  Our implementation of η does not add resistivity outside the high-density region or near the horizon, and we do see that the α_η=0.018 run maintains a horizon flux that is fairly constant and similar to 2D.  It remains to be seen how it would fare for longer timescales (not pursued for this project for lack of a 3D comparison) and whether the polar field, which eventually diminishes in 2D, needs to be replenished by a slight alpha effect. Furthermore, it is a separate challenge to show that, for a more extended magnetized disk (e.g., the more common MAD test disks <cit.>), magnetic flux could, over an extended period of time, accurately survive while advecting onto the horizon.



§ NON-ISOTROPIC DYNAMO MODELS



We now attempt to extract α^i_j and η^i_j by linear least-squares fitting.  Suppose one attempts a fit using N data points, each data point corresponding to one event—a particular spatial location at a particular time.  At each data point, there are three equations, one for each component Δ E^i.  For a given component i, combining all data points, we have the matrix equation that we wish to be satisfied


    𝐲_i ≅𝐀_i·𝐱_i ,


where 𝐲_i is length-N column vector made of all the Δ E^i, 𝐱_i is a 6-dimension column vector made of the α^i_j and η^i_j ∀ j, and 𝐀_i is a N× 6 matrix made of B^j and J^j ∀ j.  We use the symbol “≅” to mean that this is an equation deviation from which will be minimized (but, since there are more conditions than variables, will not be exactly satisfied).  The best fit model is a particular 𝐱_i=𝐱_i^ fit .

In choosing a residual function to be minimized, one must decide how data points are weighted.  One natural choice is to treat the absolute error at all points as equally important and minimize (𝐲_i -𝐀_i·𝐱_i^ fit)^T
(𝐲_i -𝐀_i·𝐱_i^ fit).  This will have the effect that the fit will care more to reduce relative error at points with strong field and high Δ E^i.  Alternatively, one might wish to weight relative errors at each point similarly.  One can do this by normalizing 𝐀 and 𝐲 at each spacetime event P used in the fit, dividing all components of each by ||y||_P, the L2 norm of the elements of 𝐲:

    ||y||_P=√(1/3∑_i=1^3y_i^2(P))

It is not obvious which sort of weighting is more desirable, so we have experimented with both.  For runs with weighting, we wish to avoid the possibility of fits dominated by attempts to reduce relative error in unimportant weakly-magnetized, very low-EMF regions.  Therefore, for these fits, we multiply 𝐀 and 𝐲 at each event P by the weighting factor


    w_P = max(||y||_P, 10^-3 y_ max)^-1 ,


where ||y||_P is the L2 norm ||y|| at P, and y_ max is the maximum value of ||y||_P across all P used for the given fit (i.e., the L∞ norm across events of the L2 norm at each event).



 §.§ Model validation


One measure of model accuracy is how well the model matches the original data.  We compute this accuracy by the following measure:


    χ_i^2 ≡(𝐲_i -𝐀_i·𝐱_i^ fit)^T
    (𝐲_i -𝐀_i·𝐱_i^ fit)/𝐲^T𝐲 .


This is related to the usual chi-squared statistic but with a different normalization.  More commonly, χ^2 is normalized by the expected error, so that its magnitude for a good fit is unity.  In this case, we normalize by the the data itself, so a very good fit would have χ^2≪ 1.  This is probably too much to hope for, so we would settle for χ^2<0.5 or so.

We also measure the variance, which provides the uncertainty in the best fit values 𝐱_i^ fit.  A variance estimate is provided by the SVD composition, valid for Gaussian statistics and uncorrelated data points.  In fact, even for the staggering in space and time we use to mitigate correlations, data points remain significantly correlated, so Bendre  <cit.> find (and we confirm) that the Gaussian variance significantly underestimates the uncertainty in 𝐱_i^ fit.

Instead, we estimate variance, similarly to Bendre , by computing 𝐱_i^ fit restricting the data used to different subsets of timesteps.  In our case, we compute 𝐱_i^ fit for even timesteps (𝐱_ even), for odd timesteps (𝐱_ odd), and for all timesteps (𝐱_ all).  Between two fits, 𝐱^A and 𝐱^B, define relative difference for each component i of 𝐱 as δ^A,B_i≡ (x^A_i - x^B_i)^2/(x^A_i + x^B_i)^2.  For each i, we get an overall δ_i by taking the maximum of δ^ even, odd_i, δ^ even, all_i, and δ^ odd, all_i.  For the fit to have some value, there must be at least one component of 𝐱^ fit that is consistent across fits (and this will probably be the one that contributes most to the fit), so for our overall variance measure, we minimize over components i:


    var(𝐱) = min_i[max(δ^ even, odd_i, δ^ even, all_i, δ^ odd, all_i)] .


Some fits involve more than one point in space, so for Models 2 and 4 below we also compute variances by skipping even and odd points in r and (if it varies within the fit) θ.  Then for each x_i component we maximize over relative differences for different time and space samples before extracting the best-constrained component.

This is a very generous measure, since we only require one dynamo coefficient to be well constrained.  For the model parameters to be meaningful, we require var([ 𝐱])<0.3.

Note that the more global the model, the more data points are relevant per fitting parameter.  This will tend to reduce variance, since one has very good statistics, but it is more difficult to have acceptable χ^2, since the model must work well for more points.  Conversely, for more local models, χ^2 will likely become smaller, but one should not be too encouraged by this, because one can fit anything as the number of conditions drops toward the number of fitting parameters, and the tendency to over-fit will often be revealed in a high variance.




 §.§ Pointwise models


We attempt two models for which α^i_j and η^i_j are taken to be functions of spatial coordinates but independent of time.  In Model 1, they are taken to be functions of r and θ, i.e., there is an independent fitting procedure applied at each point.  This is similar to the model construction of Dhang  <cit.>.  In Model 2, we take the dynamo coefficients to be functions of θ, i.e., they are presumed to be constant in both r and t.  Thus, compared to Model 1, each fit in Model 2 involves more data (corresponding to multiple points along a radial ray), which will tend to result in better variance but worse χ^2.  We perform fits both with and without weighting, meaning equations at data points are weighted either with the factor 1 or the factor w in eq:weighting.

Results for model validation are shown in Table <ref>.  For presenting average variance and χ^2, we divide the space into two regions:  “jet” for the region within π/3 of the axis, and “disk” for all other latitudes.  As we see, none of the models are very good by either of our validation measures.






 §.§ Global models



We now attempt to fit to a model for which α^ij and η^ij are not explicit functions of spatial coordinates.

For a covariant model, we would also have to build fitting formulas for these tensors out of physically relevant vectors and tensors, e.g., the metric γ^ij or the velocity v^i.  Since α^ij is axial, we require at least one pseudoscalar or pseudovector ingredient.  The construction of an azimuthally averaged system only makes sense if an approximate azimuthal Killing vector 𝐞_ϕ=∂/∂ϕ has been identified.
We can, thus, assume its availability and expect that it may appear in the formulas for dynamo tensors, both because of its physical relevance (as roughly the direction of mean velocity) and its relevance in constructing the 2D system (which breaks 3D rotational invariance).  We will use the normalized vector


    ê_ϕ≡ (g_ϕϕ)^-1/2∂/∂ϕ .


We get a pseudovector, corresponding to the direction of the symmetry axis, by taking a curl:


    ζ^i ≡ϵ^ijk∇_j(e_ϕ)_k = ϵ^ijk∂_j(e_ϕ)_k .


We will use the normalized version ζ̂≡ζ/|ζ|, which resembles the Cartesian z direction.

Some dynamo models also use ∇ρ, the gradient of the density, or they use a gradient of the turbulent velocity (assuming some knowledge of this) <cit.>.
For a thin accretion disk, the direction of density gradient would be close to sign(z)ζ̂ except on the equator.  To build the expected parity factor in α_ dyn (something like cosθ), we instead want a nearly radial third vector, which could be interpreted as the direction of gravity.  Since we have already assumed an approximate timelike Killing vector is available to us, 𝐞_𝐭, we might as well use it to define this direction


    𝐟   ≡    -∇(e_t· e_t) = -∇g_tt , 
    𝐟̂   ≡   𝐟/|f| .


 This then provides the pseudoscalar μ≡f̂·ζ̂, which is close to cosθ.  The resulting triad {ϕ̂, ω̂, 𝐟̂} provides a suitable vector basis inside the disk.  For global models including the polar region, it is not usable because on the poles ω̂ and 𝐟̂ are parallel or anti-parallel.  We can alternatively use the azimuthal Killing vector to define

    ϖ   ≡    -∇(e_ϕ· e_ϕ) = -∇g_ϕϕ , 
    ϖ̂   ≡   ϖ/|ϖ| .

Thus, 𝐟̂ resembles the spherical polar radial direction and is used to define the pseudoscalar μ, while ϖ̂ resembles cylindrical polar radial direction and is used as a basis vector.

We then have a triad of normalized but not orthogonal vectors in ê_ϕ, ζ̂, and ϖ̂ (and an associated dual cobasis); we denote components in this basis by hatted capital Latin letters, e.g. B_Î=B· e_Î.

The global dynamo model is a tensor equation for constructing the residual EMF Δ 𝐄 in terms of the magnetic field 𝐁, an axial vector, and some polar vector 𝐙, which we choose to be either 𝐉≡∇×𝐁 or B^2𝐯_T, where 𝐯_T is the transport velocity of the fluid.  For 𝐁, we report results using u𝐁 (and its curl for 𝐉), although we have also tried using the normal frame 𝐁 with similar results.

The tensor equation is


    Δ E^Î   = α^ÎK̂B_K̂ + η^ÎK̂Z_K̂ , 
    α^ÎK̂(𝐱,t)    = 
      a^ÎK̂μ^n_a(Î,Ĵ) ψ(𝐱,t) , 
    η^ÎK̂(𝐱,t)    = 
      e^ÎK̂μ^n_e(Î,Ĵ)ψ(𝐱,t) ,


where n_a(Î,Ĵ) and n_e(Î,Ĵ) are either 0 or 1, whichever is needed to obtain consistent parity for all terms.  

The scalar function ψ(𝐱,t) can be any scalar function constructed from the dynamical variables, with variables expected to be correlated with Δ 𝐄 being attractive choices.  Below, we consider several choices, including unity and the 2D spatial shear scalar σ=(σ_IJσ^IJ)^1/2.  We have also tried the gradient of the angular velocity |∇Ω|, the sound speed c_s, and the fraction of the pressure that is gas or magnetic (i.e. P_G/(P_G+P_B) and P_B/(P_G+P_B)).
We also try the scalar RMS deviation of 𝐁 and of 𝐯_T from their azimuthal average, e.g.,


    B_ rms = [|B-B|^2]^1/2 .


In an actual 2D simulation, these functions would not be available.  Rather, new evolution variables would have to be introduced with evolution equations designed to suitably approximate them.  However, designing such dynamics is only a worthwhile exercise if it is found that the actual RMS measures of non-axisymmetry provide useful ψ(𝐱,t) variables.

We also consider two alternatives for what constitutes the EMF residual:  Δ𝐄 and Δ_m𝐄.

We produce two types of global models. For Model 3, we use data across all r, θ, and t to produce a single fit for the numbers α^ÎK̂, η^ÎK̂.  For Model 4, we use the time average of the data and fit for the numbers α^ÎK̂, η^ÎK̂ using points across all r, θ.

We introduce a couple of variations.  First, it is possible that the dynamics is very different in the jet and disk regions, so we create fits using data only within certain latitude regions.  For this purpose, “disk” is defined as the region between π/3 < θ < 2π/3 and “jet” is defined as the regions θ<π/5 and θ>4π/5, so the intermediate coronal region is excluded from both.  This division is slightly different from the “disk” and ”jet” designations for reporting Models 1 and 2 quality measures in Table <ref>, but the separation here plays a different role.  It is used to determine which data is used for creating a given fit, rather than just which region is being averaged over for reporting a given norm after fitting is completed.

Also, one could argue that a dynamo model needn't capture small-scale features of Δ E^Î, so in some cases we try smoothing data before fitting.  This effectively reduces the amount of independent data and can be expected to give higher χ^2.  We convolve the data with a smoothing kernel, here chosen to be a triangle function in each direction centered on the point.  Specifically, if the half-length of the triangle is N_ sm points, then variable f at point i, j (indexing location in r and θ, respectively) is smoothed as follows:

    f_i,j   →   ∑_m=-N_ sm^N_ sm∑_n=-N_ sm^N_ smw_mn/𝒵 f_i+m,j+n , 
    
      w_mn   =   (1-|m|/N_ sm+1)(1-|n|/N_ sm+1) , 
    𝒵   =   ∑_m=-N_ sm^N_ sm∑_n=-N_ sm^N_ sm w_mn .
   
Thus, data at each point is replaced by a weighted sum of its neighbors, and the amount of smoothing is determined by the half-length of the triangle, N_ sm.

Quality of fit measures are reported in Table <ref>. 
 Due to the large number of events used to fit a small number of global parameters, best-fit values to at least some dynamo coefficients can now be determined precisely. However, χ^2 reveals that none of these best fits are good matches for the data, at least not good enough according to our criteria for viability.  We see many cases with χ^2 close to 1, which occurs because the functions 𝐁 and 𝐙 do not well overlap with 𝐄 (see Fig. <ref> and <ref>), and none of our chosen correlates ψ successfully corrects for this, so the best fit is obtained by making the model values 𝐀·𝐱^ fit small compared to the data 𝐲.






§ SUMMARY AND CONCLUSIONS



Modeling non-axisymmetric effects in accretion disks is challenging because, for some quantities of interest, the non-axisymmetric/fluctuating part dominates over the  axisymmetric/mean part.  Inside the disk, the mean field is often itself more like residual noise than a dominant field component, and it is not surprising that non-axisymmetric effects cannot be easily correlated to its local values.  One crucial feedback of fluctuations onto the mean structure is the angular momentum transport, often captured via viscous hydrodynamics.  Mean fields remain important in the polar region, and for the long-term maintenance and evolution of this region, and for the late-time oscillating toroidal mean field, 2D models will need to model the mean field while incorporating non-axisymmetric effects.  Using results for a 3D Cartesian FMR GRMHD code, we have quantified the effective dynamo α_ dyn and η_ dyn.  We then tested whether this model is a good fit to the residual EMF, considering a wide range of pointwise and global models.  We have considered dynamo closures of growing complexity.  1) We considered models with Δ 𝐄 a function of 𝐁 and 𝐉, perhaps with added proportionality to some other function of variables available to the 2D mean system.  2) We generalized this to use proportionality to the RMS magnetic or kinetic energy, assuming that some evolution system involving only mean fields could be devised to roughly mimic these.  3) We generalized further and considered that only the high-m residual Δ_3 𝐄 need be fit by a closure involving mean fields.  None of these provided satisfactory models (although, of course, evolving m=1 to m=3 modes of all variables would capture a significant fraction of the non-axisymmetric effects in itself even without a closure for Δ_3 𝐄, at the cost of being a low-azimuthal resolution 3D model rather than a 2D model).  In particular, designing evolution equations for B_ RMS or v_ RMS does not seem promising because we have failed to find a way to use them effectively.  More general and sophisticated models could still be attempted, e.g., providing evolution equations for 
Δ 𝐄 itself, analogous to turbulence models that evolve the Reynolds stress.

On the other hand, this is a stricter test than what is required of viscous hydrodynamics, in that it is usually not considered necessary to show that the residual stress tensor is well-fitted by a viscous stress tensor with some kinematic viscosity that is either a constant or is some specified function of azimuthal mean variables.  Rather, we insist that the viscosity capture one particular effect and calibrate a coefficient setting ν so that one effect has the appropriate magnitude.  Similarly, if one identifies the key features needed for long-term evolution, e.g., reducing poloidal field strength in the disk while maintaining it in the poles, simple models can be calibrated to do this. 

The main interest of this study is to have explored some new ways of quantifying dynamo effects in 3D simulations.  The main limitation of this study is the narrowness of the data to which it has been applied.  More could be learned by carrying out a similar analysis on a run that continues much longer for multiple dynamo cycles.  It is also necessary to run with different accretion disk-black hole configurations:  disks around lower-spinning black holes, different seed fields (e.g. toroidal), larger disks, disks that become magnetically arrested within a reasonable evolution time, radiatively efficient disks.



We acknowledge S. Bose, M. Forbes, F. Foucart, F. H. Nouri, and E. Most for useful conversations and encouragement, and we express deep gratitude to A. Tchekhovskoy for making the HARMPI code publicly available.  M.D. gratefully acknowledges support from the NSF through grants PHY-2110287 and PHY-2407726 and support from NASA through grant 80NSSC22K0719. Z.B.E. gratefully acknowledges support from NSF awards AST-2227080, OAC-2227105, PHY-2110352, and PHY-2409654, as well as NASA awards ISFM-80NSSC21K1179 and TCAN-80NSSC24K0100. 
B.J.K. gratefully acknowledges support from NASA LISA Preparatory Science award 80NSSC24K0360; this material is also based upon work supported by NASA under award number 80GSFC21M0002.
Computational resources for performing the GRMHD simulations were provided by the WVU Research Computing Thorny Flat HPC cluster, which is funded in part by NSF OAC-1726534.
Post-simulation analysis was performed in part on the Pleiades cluster at the Ames Research Center, with support provided by the NASA High-End Computing (HEC) Program.







