Jets from accretion disk dynamos: consistent quenching modes for dynamo and resistivity

Astrophysical jets are launched from strongly magnetized systems that host an accretion disk surrounding a central object. The origin of the magnetic field, which is a key component of the launching process, is still an open question. Here we address the question of how the magnetic field required for jet launching is generated and maintained by a dynamo process. By carrying out non-ideal MHD simulations (PLUTO code), we investigate how the feedback of the generated magnetic field on the mean-field dynamo affects the disk and jet properties. We find that a stronger quenching of the dynamo leads to a saturation of the magnetic field at a lower disk magnetization. Nevertheless, we find that, while applying different dynamo feedback models, the overall jet properties remain unaffected. We then investigate a feedback model which encompasses a quenching of the magnetic diffusivity. Our modeling considers a more consistent approach for mean-field dynamo modeling simulations, as the magnetic quenching of turbulence should be considered for both, a turbulent dynamo and turbulent magnetic diffusivity. We find that, after the magnetic field is saturated, the Blandford-Payne mechanism can work efficiently, leading to more collimated jets, that move, however, with slower speed. We find strong intermittent periods of flaring and knot ejection for low Coriolis numbers. In particular, flux ropes are built up and advected towards the inner disk thereby cutting off of the inner disk wind, leading to magnetic field reversals, reconnection and the emergence of intermittent flares.


INTRODUCTION
Astrophysical jets, that are highly-collimated beams of high-velocity material, are ejected from a variety of sources, such as young stellar objects (YSOs), microquasars (MQ), or active galactic nuclei (AGNs). They thus span a wide range in terms of extension, timescale and energy output, but there is common understanding that jet sources also share common properties that are responsible for jet launching. That is the presence of a central object, which is surrounded by an accretion disk (see e.g. Frank et al. 2014;Hawley et al. 2015;Pudritz & Ray 2019). Another fundamental prerequisite for the formation of a jet is the existence of a strong, largescale magnetic field. All launching processes that have been proposed through the decades (see e.g. Blandford & Znajek 1977;Blandford & Payne 1982; Lynden-Bell mattia@mpia.de, fendt@mpia.de * Member of the International Max Planck Research School for Astronomy & Cosmic Physics at the University of Heidelberg 1996), rely on the presence of a strong magnetic field with a favorable topology. This holds for both the nonrelativistic and the relativistic limit. The jet launching process can be studied by magnetohydrodynamic (MHD) simulations (see e.g. Uchida & Shibata 1985;Casse & Keppens 2002;Fendt 2006;Zanni et al. 2007;, mostly assuming a large-scale, strong, initial magnetic field. While this approach allows to understand a number of properties of the jet that is formed, it does not address the question about the formation of such a magnetic field. Several studies have pointed out a strong correlation between the disk magnetic field and the jet properties (Tzeferacos et al. 2009;Murphy et al. 2010;Stepanovs & Fendt 2016). Therefore, the understanding of the origin of the disk magnetic field represents a key step for understanding the formation of astrophysical jets.
For different central objects, several options for the magnetic field origin may be considered. For instance, in case of a stellar accreting object, the magnetic field may be provided by the star itself. This scenario is not feasible for AGN, since the black hole is not able to generate a magnetic field. Another scenario is that the jet-launching magnetic field is just advected through the accretion process from the ambient medium. Still, for protostars, a strong evidence for the advection of magnetic flux has not been found yet (Pudritz & Ray 2019). Therefore, a particularly interesting scenario is that the magnetic field can be generated and amplified by dynamo process that is working in the accretion disk. Disks are turbulent and in rotation, providing the two main ingredients for a dynamo (see below).
In this paper, we will indeed focus on the disk dynamo mechanism, because it presents a very general mechanism that can be responsible for generating a large magnetic flux in jet-launching accretion disks in both the protostellar and the AGN context. In particular, if we consider AGN, the disk magnetic field can also provide the magnetic flux that is needed to launch jets via the Blandford-Znajek mechanism.
The physics of cosmic dynamo action has been extensively studied in the last decades, in particular by numerical simulations (see e.g. Brandenburg & Subramanian 2005; Rincon 2019 for reviews on dynamo theory and their astrophysical applications). Astrophysical accretion disk dynamos are believed to have a turbulent origin, most probably caused by the magneto-rotational instability, MRI . The dissipation of the magnetic field has also been modeled in the same fashion (see, e.g., Casse & Keppens 2002).
Accretion disk dynamos have been suggested already decades ago (Pudritz 1981a,b;Brandenburg et al. 1995) in order to describe the strong magnetic field in turbulent accretion disks around black holes. A comprehensive modeling of jet launching by a dynamo-generated magnetic field is a challenging task because of the different physical mechanisms that operate over a very wide range of spatial and temporal scales. The spatial scale stretches from the scale on which turbulence occurs, to the scale on which the jet is launched. On the small scales, resolution studies of the MRI have been performed by a number of groups (see e.g. Guan & Gammie 2009;Davis et al. 2010;Hawley et al. 2013;Ryan et al. 2017). Modeling the disk and jet regions at the turbulent scales would be computationally hardly feasible, requiring unrealistically high computational resources.
For all these reasons, mainly two different types of dynamo action have been investigated. On one hand, socalled direct simulations show the direct amplification of the magnetic field by the turbulent motion (see, e.g., Gressel 2010;Bai & Stone 2013;Gressel & Pessah 2015;Riols, & Latter 2018;Hogg & Reynolds 2018;Dhang et al. 2020). On the other hand, there is the so-called mean-field dynamo approach (see, e.g., Krause & Rädler 1980;Rüdiger et al. 1995;Campbell 1999;Rekowski et al. 2000; Bardou et al. 2001;Chabrier & Küker 2006) that relies on (semi-)analytical solutions of the turbulent motion, allowing us to perform the large-scale numerical simulations required to consider the scales on which accretion disks and jets evolve.
In the present work we follow the mean-field approach. Whether a large-scale magnetic field can be generated by a mean-field dynamo has been numerically investigated by von Rekowski et al. (2003); ; Fendt & Gaßmann (2018); Mattia & Fendt (2020a,b). These results have demonstrated that a dynamo-amplified, strong magnetic field can evolve and is able to launch a jet or a disk outflows in general. In particular, it has been shown by von Rekowski & Brandenburg (2004); Dyda et al. (2018) how an accretion disk dynamo-generated magnetic field leads to an outflow also if combined with a central stellar magnetic field. Recently, the mean-field dynamo approach has been also extended to general relativistic MHD simulations in torii (Bucciantini & Del Zanna 2013;Bugli et al. 2014;Tomei et al. 2020Tomei et al. , 2021Del Zanna et al. 2022) or in thin accretion disks (Vourellis & Fendt 2021).
In this paper we want to further bridge the gap between analytical models and the time-dependent simulations. We focus on the non-linear evolution of the dynamo phenomena, i.e. how a strong disk magnetic field that has been dynamo-generated affects itself the dynamo action and the diffusivity within the accretion disk, which have both a turbulent origin.
Despite several options to parameterize a mean-field dynamo and turbulent diffusivity being developed in the literature, so far, based on analytical modeling of the disk turbulence (see e.g. Rüdiger, & Kichatinov 1993;Kitchatinov et al. 1994;Rüdiger et al. 1994), a consistent numerical treatment of the feedback has never been adopted in the context of the jets launching from a resistive and dynamo-active accretion disk. We thus compare different feedback strategies adopted in the context of accretion disk dynamos and we investigate how the different feedback models affect the jet launching process. Finally we investigate the effect of a -we think more consistent -diffusivity model which includes a quenching mode for strong magnetization. The latter is motivated by the fact that both the alpha dynamo and the magnetic diffusivity applied are believed to have the same source, namely the disk turbulence. Thus, a quenching of the turbulent dynamo by a strong field should also consider the quenching of turbulent dissipation at the same time. This has not yet been considered by thin-disk dynamo simulations launching jets.
The paper is organized as follows. In Section 2 we describe our numerical setup. In Section 3 we investigate the effects of different feedback dynamo models and their role in the jet launching process. In Section 4 we introduce a consistent turbulent model of both meanfield dynamo and magnetic diffusivity for thin accretion disks and we apply it to the large scale disk-jet simulation. We summarize our paper in Section 5.

MHD equations
We solve the time-dependent, resistive MHD equations with the PLUTO code 1 (Mignone et al. 2007) version 4.3 in spherical coordinates (R, θ, φ) assuming axisymmetry. We refer to (r, z) as cylindrical coordinates. The code solves the integral form of the MHD conservation laws through a high order finite volume method. To the standard PLUTO MHD equations we have added a mean-field dynamo term.
The set of primitive variables consist of, respectively, the plasma density ρ, the flow velocity v, the gas pressure p and the magnetic field B.
The set of conservation laws consists of the conservation of mass, ∂ρ ∂t (2) total energy, and magnetic field, The total energy density is computed, from the primitive variables, as with polytropic index Γ = 5/3, while the electric field is defined by (see Krause & Rädler 1980) where the first is the ideal term, the second is the diffusive term and the third is the dynamo term. The gravitational potential Φ g is provided by a central object mass M, i.e. Φ g = −GM/R. Previous studies (see e.g. Casse & Ferreira 2000a;Casse, & Ferreira 2000b;Zanni et al. 2007;Tzeferacos et al. 2013) have shown that the non-ideal heating and cooling processes may play a relevant role in the context of jet formation. Here, for the sake of simplicity, the cooling term Λ cool is set to be equal to the non-ideal (i.e. diffusive and dynamo) contribution of the electric field to the energy equation, as in Sheikhnezami et al. (2012); : We choose to neglect the non-ideal heating process in accordance with previous studies of jet launching from resistive disks (Casse & Keppens 2002).
The tensors α and η describe the effects of the meanfield dynamo and the magnetic diffusivity on the magnetic field evolution (see Sections 2.5 and 2.6).

Numerical setup
Since the PLUTO code solves the MHD equation in their adimensional form, intrinsic physical scales are not involved. We apply code units, for which every scale (time, length, energy, etc.) is normalized to its initial value at the innermost disk radius R in at the midplane. Therefore, velocities and are normalized to v k,in , the Keplerian speed at R in . As a consequence, the time is normalized to t in = R in /v k,in . Thus, 2πt in corresponds to one disk revolution at R in . For the sake of simplicity, all times that are shown without a normalization constant, are measured in units of t in .
The computational domain ranges in the radial direction between R ∈ [1, 100]R in , and in the angular direction, between θ ∈ [10 −8 , π/2 − 10 −8 ] [0, π/2] (where θ = 0 corresponds to the jet axis and θ = π/2 corresponds to the disk midplane). In order to be able to cover a large extension of the numerical grid, we apply stretched cells in the radial direction (∆R = R∆θ), while the spacing in the angular direction is equidistant. We adopt a resolution of [N R ×N θ ] = [512×256] through the whole paper. This gives a resolution throughout the paper of 32 cells per geometrical disk height H. We found that a lower resolution may lead to numerical issues, especially in case of strong interplay between the magnetic field and the dynamo in the innermost accretion disk regions.
The spatial reconstruction is performed employing the piecewise parabolic method for spherical coordinates (Mignone 2014), while the time integration is achieved through a third-order Runge-Kutta scheme. The flux is computed through a Harten-Lax-van Leer Contact (HLLC) Riemann solver (Toro et al. 1994;Li 2005) The solenoidality of the magnetic field is preserved through an upwind constrained transport method (Londrillo & del Zanna 2004). In particular, as in Mattia & Fendt (2020a,b), we employ a combination of the UCT-HLL for the ideal and resistive components, while the dynamo is added through an arithmetic average (see Mignone & Del Zanna 2021 for a more detailed study on the upwind constrained transport schemes). In order to preserve stability, we adopt a Courant-Friedrichs-Levy time stepping with CF L = 0.4 < 1/ √ N dim .

Initial conditions
The initial conditions are identical to those applied in Mattia & Fendt (2020a,b). For the sake of clarity, we summarize the fundamental steps. The simulation starts with a very weakly magnetized disk. The initial magnetization, defined as the ratio between the thermal and the magnetic pressure both measured at the disk midplane is µ in = B 2 in /2p in = 10 −5 . Therefore, we recover the initial disk structure by assuming equilibrium between thermal pressure gradients, gravity and centrifugal forces Fendt & Gaßmann 2018), This equation can be solved by assuming that all the hydrodynamical variables X scale as power laws of the radius R, i.e. X = X 0 R β X F X (θ), where X 0 is the corresponding variable evaluated at the inner disk radius, while F X is the angular dependence. We assume a polytropic gas, p ∝ ρ Γ , and that all characteristic velocities scale as the Keplerian velocity, β v φ = −1/2. Under these assumptions, the radial dependence of the hydrodynamical variables is determined by β ρ = −3/2 and β p = −5/2 (Blandford & Payne 1982). In order to recover the explicit value of F X (θ) we need to set two parameters, both evaluated at the inner disk radius: the disk density ρ 0 and the thermal disk height, defined as the ratio between the isothermal sound speed and the Keplerian speed = c s /v φ | θ=π/2 . We set ρ 0 = 1 and = 0.1.
Once the density and pressure distributions are recovered, we use the combination of the radial and angular component of Eq. 8 in order to compute the angular dependence of the toroidal velocity (as in Zanni & Ferreira 2009, 2013 but neglecting the viscous terms).
Above the disk we define an hydrostatic corona, with ρ c,in = 10 −3 ρ 0 . At the transition between accretion disk and corona the disk pressure equals the coronal pressure, involving a jump in density.
As pointed by , we define as geometrical disk height the region where density and toroidal velocity decrease significantly. Since through all the paper we assume = 0.1, we adopt a linear approximation H = 2r which is able to reproduce with great accuracy the relation between the thermal and the geometrical disk height. The study of the influence that the disk height has on the dynamo-generated magnetic field will be the subject of our future works.
All the simulations are initialized with a weak purely radial magnetic field defined by the vector potential ) The strength of the initial poloidal magnetic field is determined by B p,in = √ 2p 0 µ in , where µ in = 10 −5 is the initial magnetization along the disk midplane. By construction, the initial toroidal magnetic field B φ vanishes.

Boundary conditions
The physical evolution is heavily determined by the choice of boundary conditions. A "wrong" choice of boundary conditions may easily lead to a non consistent or unphysical scenarios (see e.g. a recent review on these issues by Boneva et al. 2021).
The boundary conditions we adopt here are a slightly revised version of the approach in ; Fendt & Gaßmann (2018), and are reported in Table 1 for convenience. Standard symmetry conditions are applied along the disk mid-plane and along the rotational axis.
Equatorial symmetry (i.e., in this case, symmetry for the component B θ and anti-symmetry for the components B R and B φ ) conditions does have certain consequences. Concerning the magnetic field it imposes that the toroidal magnetic field does vanish along the equatorial plane. Thus the poloidal field intersects the midplane perpendicularly. While with our approach we constrain the disk mid-plane to the equatorial plane, we expect similar effects also from a bipolar approach, unless a large-scale bipolar asymmetry will build up. In Fendt & Gaßmann (2018) the bipolar setup lead to some disk warping, however, on the large scales the field remained symmetric. This implies that across the disk mid-plane (in that case different from the equatorial plane of the grid) the toroidal field changes sign, leading to a vanishing toroidal field over there -similar to our present setup. For the choice of the sign of α (in both papers) a large scale flux emerged.
The situation may be different for mean-field dynamo of opposite sign. This may be relevant in particular for bipolar simulations considering both hemispheres. Rekowski et al. (2000); Bardou et al. (2001) have shown that a change of sign for the α (applied to both hemispheres leads to a quadrupolar large scale magnetic field geometry. These results, obtained by a semianalytical treatment has been confirmed by numerical simulations of a scalar mean-field dynamo, finding either dipolar (negative α, upper hemisphere) or quadrupolar fields (positive α, upper hemisphere) (Fendt & Gaßmann 2018). This effect can, obviously, not be investigated if only one hemisphere is considered. However, such an investigation would be beyond the aim of the present paper, which emphasis is to apply more consistent dynamo and quenching models. A different polarity of the dynamo would also be in contrast with the analytical model of Rüdiger et al. (1995), and we thus choose to not investigate further the polarity of α.
Other reasons to constrain ourselves to one hemisphere is just simplicity. We intend to investigate the pure dynamo effect that should not be affected by the dynamics of the disk mid-plane.
In fact, our prescription for the α and η distribution is fixed in space (and centered on the mid-plane) 3 . To allow for disk warping would require to determine the physical disk mid-plane for each time step and update for the α and η distribution consequently. This has not been done previously and it is beyond the scope of the paper.
Both the inner and outer radial boundaries are separated into a disk (θ > π/2 − 3 ) and a coronal (θ < π/2 − 3 ) region. The extent of the inner disk boundary is somewhat broader than the initial disk height. The reason is that the disk, especially in case of strong accretion, may slightly inflate, but all material delivered by disk accretion must be able to be absorbed by the boundary.
Along both the disk and the coronal inner boundary we prescribe v θ = 0. For the disk part, density, pressure and toroidal velocity are extrapolated from the domain quantities by a power law. The radial velocity is also computed by extrapolating a power law, however we allow only for negative radial velocity.
The density and pressure for the inner coronal boundary condition are set to their initial value, but multiplied by a factor ρ c,in (e.g. p c,in for the pressure, see Eq. 9), and describing a fixed radial inflow of v R = 0.1. Together, this allows for a more stable evolution between the disk and coronal boundary (compared to Mattia & Fendt 2020a,b) as strong density gradients may appear at the interface.
We prescribe a poloidal magnetic field boundary condition along the boundary, that is defined by satisfying the divergence-free condition. For that, we need to prescribe the θ-component (only). At the inner disk boundary we prescribe a fixed inclination of the poloidal magnetic field, with an angle where ϕ is the angle between the magnetic field and the initial disk surface. At the inner coronal boundary we require the conservation of the magnetic flux. The B φ component is recovered by extrapolating a power law in the radial direction ∝ R −1 .
For the outer boundary conditions, the density, pressure and toroidal magnetic field are computed by imposing a power law extrapolation, while the three velocity components follow a zero-gradient prescription. However, in addition we require that there is no radial inflow into the coronal region and no outflow from the disk region. Finally, the component B θ is computed assuming a power-law ∝ R −1 , while the radial magnetic field is computed satisfying the solenoidality condition.

The dynamo model
As we have shown recently (Mattia & Fendt 2020b), the non-isotropic disk dynamo model of Rüdiger et al. (1995) proved to have some considerable advantages, e.g. the reduced number of parameters required to describe the dynamo tensor and the greater stability (compared to a scalar dynamo model) when the initial magnetic field has a non-zero vertical component. Therefore, in this paper, we apply the dynamo tensor derived by Rüdiger, & Kichatinov (1993); Rüdiger et al. (1995) within the thin disk approximation, thus with negligible Table 1. Shown is the choice of how the variables in the outermost domain cell are extrapolated across the respective ghost cells of the boundary. Outflow denotes standard, zero-gradient PLUTO outflow boundary conditions. Alternatively, divergence-free conditions are involved, a prescribed field inclination, or flux conservation.
where the symbol • corresponds to the element-wise product of two vectors, and with the adiabatic sound speed c s at the disk midplane. The vectorᾱ 0 determines the strength of the dynamo tensor, and F α (z) describes the vertical profile of the alpha-effect. Naturally, we also need to assume a sufficiently high disk ionization. The vectorq α is a generic dynamo-quenching function. The specific form of the dynamo quenching will be described in Section 2.7. As pointed by Rüdiger et al. (1995), the effects of the rotation on the turbulence can be described by the Coriolis number, where Ω is the frequency of revolution and τ c is the turbulence correlation time, which can be recovered only by direct simulations (see e.g. Gressel 2010; Nauman & Blackman 2015; Gressel & Pessah 2015. The exact connection between the local shearing box simulations and the large-scale mean-field dynamo is still unclear, since the values found by the local approach (Ω * 1) are 10 smaller than the ones required in order to recover the amplitude of the meanfield dynamo. Future multi-scale dynamo simulations will hopefully solve this problem.
In order to investigate the effect of strong and weak dynamos, we select three values of the Coriolis number.
With Ω * = 10 we investigate the strong dynamo regime, while Ω * = 5 refers to the moderate dynamo regime, and Ω * = 1 is the weak dynamo regime. For comparison, all of our simulations are performed in all these three regimes.
The explicit form ofᾱ 0 is in cylindrical coordinates (Rüdiger et al. 1995;Mattia & Fendt 2020b): In the thin disk approximation we can safely assume α 0,R α 0,r , however, for the sake of consistency, all the dynamo vector components are transformed into spherical coordinates. The profile function F α (z) confines the dynamo within the accretion disk A sinusoidal function (as in Bardou et al. 2001) is adopted instead of a linear function (as in Rüdiger et al. 1995;Rekowski et al. 2000) in order to avoid sharp discontinuities between the disk and the coronal area.

The diffusivity model
The magnetic diffusivity tensor is assumed to have a diagonal structure.
As in Zanni et al. (2007); Sheikhnezami et al. (2012), we adopt an α-prescription where v A is the Alfvén speed at the disk midplane, H denotes the geometrical disk height, and F η (z) describes the vertical profile of the magnetic diffusivity. As pointed out by , all the diffusivity models that we investigate can be represented as where c s is the adiabatic sound speed at the disk midplane. The diffusivity profile is anisotropic, following (Ferreira & Pelletier 1995) The quantity α ss represents the dimensionless parameter measuring the strength of turbulence (Shakura & Sunyaev 1973). Implicitly, the magnetic diffusivity is assumed to have a turbulent nature, caused by the MRI . We apply where µ| π/2 is the magnetization derived along the disk midplane. However, the magneto-rotational instability can be triggered by both the poloidal and the toroidal field (Fromang 2013). As the latter vanishes along the mid-plane by construction, we adopt the following prescription, if not specified otherwise, where the magnetization µ D is defined by the average total magnetic pressure over the geometrical disk height and by the mid-plane gas pressure. This approach has also been adopted by , although a different diffusivity model was applied.
For the strength and the anisotropy of the diffusivity tensor we follow Kitchatinov et al. (1994); Rüdiger et al. (1995); Rekowski et al. (2000); Mattia & Fendt (2020b), This is the same model as in Kitchatinov et al. (1994); Mattia & Fendt (2020b), although the terms are arranged differently 4 for the sake of clarity. Quite different models for the anisotropic diffusivity have been employed in the last decades (see e.g. Casse & Ferreira 2000a;Ferreira & Casse 2013). In our approach, the strength of anisotropy is not an independent quantity, but depends directly on the Coriolis number. In the low dynamo regime, Ω * → 0, the diffusivity is quasi-isotropic. The anisotropy becomes larger when the Coriolis number reaches higher values.
We apply a profile function that allows a smooth decrease of the diffusivity outside the disk region.

Dynamo quenching models
The dynamo action can be understood as a macroscopic effect of the local magneto-rotational instability, 4 following Mattia & Fendt 2020b which results in an additional hyperbolic term in the induction equation. As the presence of a strong large-scale magnetic field is able to suppress the MRI, the same will happen to the dynamo process in the accretion disk as soon as the dynamo-amplified magnetic field becomes strong enough. Dynamo action will then be quenched.

Diffusive dynamo quenching
The Diffusive Dynamo Quenching (DDQ) has been proposed by  in the context of jet launching simulations from resistive accretion disk and by  in order to saturate the dynamo amplification of the magnetic field. With this approach, no direct quenching on the dynamo or the magnetic diffusivity is applied.
Instead, the infinite exponential increase of the magnetic field is prevented by a strong increase in the magnetic diffusivity. In contrast to the "standard" diffusivity models (e.g. Casse & Keppens 2002;Zanni et al. 2007) applied in non-ideal simulations of jet launching regions, the quantity α ss is defined as follows, The main advantage with this choice of quenching mode is that it avoids the accretion instability (Campbell 2009), which may suppress the jet launching process and, in addition, is prone to numerical issues. On the other hand, the strong dependence of α ss on the magnetization may lead to un-physically high values of the magnetic diffusivity.

Standard dynamo quenching
So far there is no general consensus about how to calculate the critical magnetization value for the quenching. Since the turbulence that is causing the turbulent dynamo effect is supposed to be a consequence of the MRI, the saturation of the dynamo action should most probably depend on the relative magnetic field strength at the disk mid-plane.
Moreover, a quenching based on the disk magnetization (see e.g. Vourellis & Fendt 2021) is in agreement with the fact that the MRI is excited by both the poloidal and the toroidal magnetic field 5 . Thus, we start with the most simple approach for an isotropic quenching model (henceforth Standard Dynamo Quenching, SDQ) (Ivanova & Ruzmaikin 1977;Brandenburg & Subramanian 2005;Moss et al. 2016) that is basically de-  Rüdiger, & Kichatinov (1993) for the the different dynamo components as a function of the disk magnetization pending depending on the disk magnetization, We point out that such a global -thus non-localquenching prescription is not easy to fully parallelize in the code. As investigated by Vourellis & Fendt (2021), a weak parallelization in the θ−direction (i.e. a parallelization with only a few number of cores in the θ−direction) overcomes this problem with only little additional computational costs. For a local quenching (using the local magnetization value on the grid cell, see Tomei et al. 2020), that is straightforward to parallelize, the main idea of turbulence generation in the disk mid-plane in combination with the generation of a large scale magnetic flux gets somehow lost. Also, the dynamo process itself becomes unstable when every grid cell applies a different strength of the dynamo -again a conflict with the aim of generating a large-scale magnetic flux.

Non-isotropic dynamo quenching
As for the strength of the mean field dynamo, also the feedback of the magnetic field on the dynamo cannot always be approximated with a single scalar function. Thus, similar to the definition of an an-isotropic dynamo tensor, the quenching of the dynamo effect is also tensorial, thus acting in different strength on the dynamo tensorial components.
Here, we consider a feedback model (henceforth Nonisotropic Dynamo Quenching, NDQ) that follows from an analytical study of turbulence (Rüdiger, & Kichatinov 1993), which has elaborated different quenching functions for different components of the α−tensor. It does not only depend on the strength but also on the orientation of the dynamo-amplified magnetic field. For our purpose, the suppression of the mean-field dynamo is first computed in cylindrical coordinates (Rüdiger, & Kichatinov 1993) and then converted to spherical coordinates, where we have defined (26) The quenching parameter β is defined as β = µ D /µ 0 . The dependence of the quenching components on the magnetization is shown in Fig. 1.
Essentially, the quenching depends on the disk magnetization as ∝ µ −3/2 for the radial and toroidal component, while follows ∝ µ −1/2 for the θ-component. Therefore, we can expect a more sudden and rapid saturation of the magnetic field at lower magnetization. We point out that 0 ≤ q αz ≤ 1, regardless of the values of µ D or B z,D .

Dynamo number and turbulence parameter
In order to investigate a set of different simulations, the comparison of dimensionless parameters which do not have a strict dependence on the initial parameter space, plays a key role for understanding the physical evolution. In the context of the mean-field dynamo simulations the dynamo number D plays the essential role when it comes to understand the efficiency of the dynamo process, The dynamo number is the product of the azimuthal Reynolds number which defines the balance between magnetic diffusion and Ω-effect, and the magnetic Reynolds number which defines the balance between magnetic diffusion and α−effect. Both processes contribute to the amplification of the magnetic field. While the component α φ is related to the amplification of the poloidal components, the Ω-effect amplifies the toroidal magnetic field component. Note, however, that our approach considers as well an α-effect for the toroidal field component, thus considering formally an α 2 Ω dynamo. Through the different tensorial components both the toroidal and poloidal components of the magnetic field are amplified. However, in our case, the α-effect for the toroidal field is always less efficient than the Ω-effect, due to the strong shear that is present in a Keplerian accretion disk (in comparison to e.g. a star). We therefore can neglect this effect for our further discussion, and consider only the α φ for the magnetic Reynolds number, and, thus, also for the Dynamo number.
High dynamo numbers imply the possibility of an efficient dynamo process, low numbers vice versa. The critical dynamo number, separating the two regimes, depends on the details of the model setup (see e.g. Stepinski & Levy 1988, 1990Torkelsson & Brandenburg 1994) and has to be found by applying a series of parameter runs. For example, Brandenburg & Subramanian (2005) find a critical dynamo number D 10, below which the magnetic field cannot be amplified. By connecting galactic dynamo simulations with accretion disk simulations, Boneva et al. (2021) The size of the system is denoted by H, here represented by the disk height. By construction, the dynamo components vanish at z = H. We therefore compute the quantity α φ at z = H/2.
Another key parameter in disk simulations as well as in jet launching simulations is the turbulence parameter α ss , which parametrizes the strength of disk turbulence, i.e. the disk turbulent viscosity (Shakura & Sunyaev 1973). On one hand, it represents a direct link between the disk magnetization and the disk diffusivity (see Eq. 20), on the other hand it can be recovered both from observations and direct simulations. As pointed by King et al. (2007), observational evidences show that a range 0.1 < α ss < 0.4 is required to provide a good description of the behaviour of fully ionized, thin accretion disks. Nevertheless, numerical simulations of direct turbulence recover values which are an order of magnitude below the observational value.

DYNAMO FEEDBACK MODELS
Quenching the dynamo tensor prevents an infinite field amplification. Different quenching methods lead to a different saturation of the magnetic field. As mentioned above, quenching of the turbulent dynamo is a physical consequence of the process that produces turbulence. Ideally, quenching models are derived from first principles of turbulent plasmas.
Very general correlations have been found between the accretion disk magnetization and the jet speed or the jet collimation (Fendt 2006;Pudritz et al. 2006;Stepanovs & Fendt 2016), demonstrating that a high disk magnetization is tightly correlated with a high jet velocity.
These correlations can be extended, then linking the strength of the dynamo with the jet speed, as a stronger dynamo implies a stronger field amplification (Fendt & Gaßmann 2018;Mattia & Fendt 2020a,b). For this reason, we do investigate the interplay between the amplified magnetic field and the dynamo, and how it affects the jet launching process.
In Figure 2 we show the density distribution of the disk-jet structure, together with the magnetic field geometry, for different feedback models (from left to right, respectively, the diffusive quenching, the standard quenching and the non-isotropic quenching) and different Coriolis numbers (from top to bottom, respectively, Ω * = 1, 5, 10). Overall, we see that the magnetic field, near the rotation axis, the magnetic field amplified by the dynamo, has evolved into a large-scale open geometry. The magnetic field structure, together with its amplification, leads to a highly collimated outflow.
However, for lower Coriolis numbers we notice a more turbulent outflow from the inner disk region. Such a magnetic field distribution suggests that the outflow is driven by the toroidal magnetic pressure gradient rather than by magneto-centrifugal forces. These simulations show a major magnetic loop, whose distance from the inner disk depends on the Coriolis number and on the feedback model. In addition, another loop may emerge, then indicating the presence of a dynamo inefficient zone (Mattia & Fendt 2020a).

Magnetic field amplification
The primary effect of the dynamo tensor is the amplification of the magnetic field. However, below a critical value of the Coriolis number, even in presence of a dynamo effect, the magnetic field is not amplified or is only weakly amplified. As a direct consequence, for example a fast and collimated outflow cannot be launched.
We identify a critical value for the Coriolis number as the one where the dynamo timescale is longer than diffusion timescale (Dyda et al. 2018;Fendt & Gaßmann

2018),
Note that quenching methods may act in quite a different way, as the initial strength of diffusivity for the diffusive quenching is ∼ 2 − 3 orders of magnitudes weaker than the one from the standard quenching. This difference strongly reflects on the existence of a critical Coriolis number.
When applying the diffusive quenching method, we have recovered a critical value (i.e. the value below which the magnetic field is not immediately amplified) of the initial Coriolis number Ω * C 0.15 (Mattia & Fendt 2020b). On the other hand, when applying the standard quenching method, or the non-isotropic quenching methods we have developed, a critical value of the Coriolis number (in order to determine whether the initial dynamo can amplify the magnetic field) Ω * C 2 is found.
This difference can be seen in Fig. 3, where we show the time evolution of the poloidal magnetic energy from radius R = 10 to the end of the domain (R out = 100), while applying a Coriolis number of Ω * = 1 (dotted lines). As expected, the poloidal magnetic energy of the diffusive quenching method, which is shown in Fig. 3, is amplified stronger and more rapidly than for the cases of the other quenching models.
Moreover, the field amplification is preceded by a short decrease. The reason behind this is that, at t = 0, the diffusive timescale is shorter than the dynamo timescale. Thus, the magnetic field is diffused away, leading to a decrease in the magnetization and, therefore, of the magnetic diffusivity. Once the magnetic diffusivity has decreased, the dynamo timescale becomes again shorter than the diffusive time scale. When the Coriolis number is higher than its critical value (i.e. in all the cases considered but the ones defined by the dotted red and green lines in Fig. 3), the amplification of the magnetic field occurs almost instantly. Thus, when the magnetic field increases, also the dynamo quenching increases, suppressing the dynamo action and lowering the amplification of the magnetic field. In order to disentangle the impact of the quenching methods we applied Coriolis numbers of Ω * = 5 and Ω * = 10. We find that the field amplification is faster and also stronger when the diffusive quenching model is applied (see Fig. 3).
We notice that the magnetic energy that originates from dynamo action by employing the diffusive quenching method and a Coriolis number Ω * = 1 is approximately the same of the one obtained by the standard quenching method and Ω * = 10 until t = 3000. Since the diffusive dynamo quenching model does not involve a suppression of the dynamo tensor, the magnetic field saturates as soon as the diffusivity is strong enough to counterbalance the dynamo effect. On the other hand, the standard quenching method features both an increase of the magnetic diffusivity and a decrease of the dynamo. Thus, the same Coriolis number will lead to a different strength of the disk poloidal field depending on the feedback model. However, as shown in Fig. 4, the presence or the absence of the dynamo inefficient zones from the early stages (Mattia & Fendt 2020a) plays a key role in the jet structure and evolution. We find that a low Coriolis number leads to the formation of dynamo inefficient zones regardless of the quenching model, in agreement with Mattia & Fendt (2020b). On the other hand, we also find that dynamo inefficient zones are present for large Ω * , which we did not find in our previous works. The formation of these zones may be connected with the increase of grid resolution that is coupled with the HLLC Riemann solver we now apply, which together provides a better resolution of the disk substructures (probably smeared out by the more diffusive HLL solver).
Finally, the results obtained by applying the nonisotropic dynamo quenching model show no difference from the standard dynamo quenching model for low Coriolis number. However, for higher Coriolis numbers, the different suppression of the dynamo in the non-isotropic dynamo quenching model (α φ ∝ µ −3/2 D , while α φ ∝ µ −1 D in the standard dynamo quenching model) leads to a different saturation of the magnetic field. More specifically, we find that the magnetic field, in the non-isotropic dynamo quenching model, is saturated towards a lower disk magnetization than the one obtained by the standard dynamo quenching model.

Dynamo number and turbulence parameter
The dynamo number is traditionally used to indicate the strength and efficiency of dynamo activity. A high dynamo number indicates an efficient dynamo, thus leading to strong field amplification. Vice versa, a low dynamo number indicates that dynamo cannot act efficiently anymore, and the magnetic field generated has reached its saturation value -either established by a strong magnetic diffusivity, thus by diffusive quenching (diffusive dynamo quenching), or by suppressing the dynamo activity itself, thus by direct quenching of the dynamo-alpha (standard dynamo quenching, nonisotropic dynamo quenching). The critical dynamo number (see Sec. 2.8) differentiates the two regimes.
Our simulations (see Fig. 5) confirm earlier predictions of Brandenburg & Subramanian (2005); Boneva et al. (2021), that is that (i) field amplification does not occur for D 10 (red in Fig. 5), while (ii) amplification occurs when the dynamo number supersedes its critical value since the magnetic diffusivity decreases (blue in Fig. 5).
Moreover, we find that the critical value of the dynamo number does not depend on the feedback model applied or the Coriolis number that is given. This suggests, essentially, that the amplification and the saturation of the magnetic field is a very general property of the mean-field dynamo approach, and does not depend on certain modeling details. Note, that the exact value of the critical dynamo number can be influenced by the numerical resolution applied and the numerical algorithms adopted (see, e.g. Stepinski & Levy 1988, 1990Torkelsson & Brandenburg 1994).
The main difference between the quenching models we apply, is the absence of dynamo inefficient zones for low values of the Coriolis number and the feedback models which imply a suppression of the dynamo. A possible explanation is that a magnetic field reversal can be maintained only if the dynamo does not vanish. If the magnetic field is not constantly amplified by the meanfield dynamo (because of the quenching), the field reversal zones are able to reconnect and be diffused away. This is not possible if the dynamo is not suppressed. However, since the standard dynamo quenching model and the non-isotropic dynamo quenching model lead to a suppression of the dynamo (and not in the diffusivity), the dynamo inefficient zones are more likely to be suppressed. We also point out that the presence of dynamo-inefficient zones, which are not strictly related to the component α φ , is still possible. In this regard, On the other hand, the turbulence parameter α ss (Shakura & Sunyaev 1973) does not explicitly depend on the dynamo tensor, but only on the magnetic diffusivity. Therefore it can be applied as a useful tool in order to understand the evolution of the magnetic diffusivity once the magnetic field amplification took place. As we can see from Fig. 6, the feedback models which include the suppression of the dynamo and the diffusive dynamo quenching model show several differences with regard of how the turbulence parameter depends on the Coriolis number once the magnetic field is saturated. In particular, we find that the results applying the diffusive dynamo quenching model show a unique dependence on the Coriolis number.
As pointed out in the previous section, for the diffusive dynamo quenching model, the magnetic diffusion is the only process that is able to saturate the meanfield dynamo. Because of the different amplification of the magnetic field for different Coriolis numbers, different in both the strength of the magnetic field and the timescale of amplification, the disk magnetization saturates to different levels (see Fig. 3). Because of the strong, i.e. quadratic, dependence of the diffusivity η on the disk magnetization µ, in the diffusive dynamo quenching model this dependence is reflected in the fact that also the α ss is found to depend on µ.

Accretion and ejection
The different saturation levels for the magnetic field strength play a key role in the dynamical evolution of the accretion disk and the subsequent jet launching process. This is demonstrated in Fig. 7 showing the accretion and the ejection rate of the launching area. The accretion rate is computed by integrating the net radial mass flux through the disk at R = 7, with the disk opening angle defined as θ S ≡ arctan(2H/r). Similarly, the ejection rate is computed by integrating between R in and R = 7 along the disk surface,Ṁ The interrelation between the accretion rate and the Coriolis number confirms, regardless of the feedback mode, previous results obtained by Fendt & Gaßmann (2018). In particular, we find that a stronger dynamo leads to a stronger accretion rate. We investigated the impact of the magnetic field topology on the accretion process (i.e. the presence or the absence of dynamo inefficient zones) in detail previously (Mattia & Fendt 2020a). This influence is also confirmed by comparing the standard dynamo quenching model for Ω * = 10 with the diffusive dynamo quenching model for Ω * = 1. Here, despite a similar amplification of the poloidal magnetic field (see Fig. 3), the presence (DDQ model, Ω * = 1) or absence (SDQ model, Ω * = 10) of the dynamo inefficient zones, plays a key role in the accretion of material.
On the other hand, the mass ejection, acting on much shorter timescales, shows less pronounced differences for the variety of quenching methods or the different values of the Coriolis number. However, we observe that lower Ω * generally lead to a higher ejection-to-accretion ratio, which is in good agreement with the findings of Mattia & Fendt (2020b). In case of a strong dynamo, the isotropic quenching models show that < 50% of the accreted material is ejected. This is in good agreement with previous resistive jet launching simulations that do not consider dynamo action, but start from a prescribed large-scale magnetic field (see e.g. Zanni et al. 2007;Sheikhnezami et al. 2012). For lower Coriolis numbers, accretion requires more time to be established because of the slower amplification of the magnetic field. On the other hand, the non-isotropic quenching model shows that almost all the matter accreted is ejected. This is a direct consequence of the feedback model which leads to a saturation of the magnetic field at lower magnetizations (se Fig. 3). As a consequence, the magnetic diffusivity, which depends on the disk magnetization, is less amplified than in the isotropic feedback models, leading to a lower accretion. For this reason, during certain periods of time, the ejection-accretion rate (defined asṀ ej /Ṁ accr ) may actually exceed unity. This may imply that disk volume regions of very low mass or density may be present for some time until they are replenished from the mass reservoir at larger disk radii.
We observe a particularly interesting period at time ≈ 3500 for the standard quenching model (STQ) with high Coriolis number Ω * = 10. A sudden drop in the accretion and, consequently, in the ejection rate appears. When looking at the strength of the mean-field dynamo φ−component as a function of radius, we notice that at t = 3500 it becomes significantly stronger than immediately before or after.
We find that the reason behind this sudden change is a small decrease in the toroidal magnetic field strength, which seems to amplify the dynamo because of the smaller quenching term (which is related to the magnetic field strength, see Fig. 8). In particular, a weaker magnetic field (and therefore a weaker magnetization) leads to a lower quenching and therefore a stronger amplification of the magnetic field. In addition, low magnetic field strength triggers a lower magnetic diffusivity, which, in turn, curbs the accretion process, because it implies a lower diffusivity. Once the magnetic field is amplified again by the dynamo, the system goes back to a more stable configuration.

Jet properties
All our models, which consider a feedback of the magnetic field on the dynamo, lead to a quasi-steady, saturated state. Therefore, our dynamo simulations should retrieve, qualitatively, the essential jet properties found by simulations that do not invoke a dynamo process (see e.g. Tzeferacos et al. 2009;Murphy et al. 2010;Figure 8. Evolution of the dynamo at each radius for the Ω * = 10 and standard dynamo quenching model simulation. The angle θ = 85.5 • corresponds to the thermal disk height (and to half of the geometrical disk height). At this angle the angular profile of the dynamo reaches its maximum absolute value.
Stepanovs & Fendt 2016), regardless of the method for the dynamo feedback.
In Mattia & Fendt (2020b) we discovered a unique numerical correlation between the Coriolis number Ω * (and therefore the dynamo strength) and the asymptotic jet speed. In this paper, we now make a further step in this regard and investigate this interrelation for different feedback models. For this purpose, we select certain magnetic flux surfaces (thus contours of the vector potential, respectively poloidal magnetic field lines), and compute the local disk magnetization and the poloidal velocity of the corresponding outflow along that field line.
We do this for each of the code outputs (uniformly distributed and separated by ∆t = 100), starting from t = 700, i.e. the time when a jetted outflow is already formed and has reached the outer boundary, until the final time step of each simulation. The results are shown in Fig. 9, where the two panels present, respectively, the values obtained for radii 1.5 < R < 5 and 5 < R < 10.
As we can see, the radial distance makes a difference concerning the smoothness of the interrelation. For radii 5 < R < 10 the interrelation looks very well defined , while for smaller radii this correlation is partially broken. At small radii we notice the presence of vertical 'columns', i.e. zones with similar magnetization that exhibit a variety in jet speed. This is mainly due to the time evolution of these parameters: both the disk magnetization as well as the jet speed may vary in time for the same radial distance, during the same simulation. Such variations occur more frequently in the inner disk region because of different timescales involved.
Our understanding is the following. Small temporal variations of the dynamo leads to suppression and amplification of the magnetic field in the inner disk region (as seen on the left region of Fig. 8). As a consequence, small differences in the magnetic field can lead to different outflow speed in the inner launching region. Such differences are leading to internal shock within the jet, affecting the outflow. The small variations in the inner disk magnetic field are also able to affect the ratio between the poloidal and the toroidal magnetic field, which can be related to the jet speed and collimation. This is a result of the magneto-centrifugal acceleration involved.
On the other hand, at larger radii the system has reached saturation towards a steady state, leading to a more narrow interrelation. At even larger radii, R > 10, the low magnetization and the slow rotation lead to a weak disk magnetization and, therefore, to a slower outflow speed. We point out that at these large radii, the disk has accomplished only few revolutions, and the whole inflow-outflow structure has not yet settled into a quasi-steady state.
However, as a key result, despite showing differences in both the disk magnetization and the jet speed, the different feedback models we have examined show in general a unique trend. They all follow a very similar relation between the two quantities, suggesting that this relation jet speed versus disk magnetization does not depend on the dynamo process, the diffusivity model, or the quenching method. It simply confirms the general relation between these leading inflow-outflow parameters that have been discovered previously (Stepanovs & Fendt 2016;Mattia & Fendt 2020b).
Still, the feedback of the magnetic field on the dynamo action plays a key role for the saturated disk magnetization. This holds in particular, because of the more efficient suppression of the dynamo, the non-isotropic dynamo quenching model reaches the saturation of the magnetic field already at a lower magnetization levels, about one order of magnitude below.

A CONSISTENT QUENCHING MODE FOR DIFFUSIVITY
In the last section we have considered the effects of applying different dynamo feedback modes, meaning how to realize the physical effect of quenching the dynamo activity by a strong magnetic field. Here, we want to go one step further towards a self-consistent modeling of mean-field dynamos. That is to consider the backreaction of the magnetic field on the magnetic diffusivity. Because of their common origin -the turbulence of the disk material -both the quenching of the mean-field dynamo and the magnetic diffusivity should be treated in a similar way. A strong global magnetic field suppresses the turbulence and, thus, both the turbulent dynamo effect and and the turbulent magnetic diffusivity.
In this section we will put this on more physical grounds, considering a quenching model that follows from analytical mean-field theory and that incorporates effects on both the mean-field dynamo and the magnetic diffusivity, following the prescriptions of ; Rüdiger et al. (1994) (henceforth Consistent Turbulence Quenching, CTQ).

Quenching of turbulent magnetic diffusivity
For the feedback of the magnetic field on the dynamo -the dynamo quenching -we follow Rüdiger, & Kichatinov (1993), as described in Section 2.7.3. As in Section 2.6 the an-isotropic components of the magnetic diffusivity can be computed directly in spherical coordinates.
Here we apply the quenching model following Kitchatinov et al. (1994); Rüdiger et al. (1994), with The dependence of the diffusivity quenching components on the magnetization is shown in Fig. 10. We point out that the dependence of magnetic diffusivity on the   Kitchatinov et al. (1994); Rüdiger et al. (1994) the η−diffusivity as function of the disk magnetization disk magnetization is also determined by the disk turbulence parameter α ss . Here, we model this applying α ss ∝ √ µ D (as in Eq. 20).
Mean-field dynamo models applying a diffusivity quenching (as in Rüdiger et al. 1994) have been applied in the context of galactic dynamos (Schultz et al. 1994;Elstner et al. 1996), although the quenching model was never coupled with a non-isotropic diffusivity. Here, because of the rapid disk rotation and the strong magnetization needed for jet launching, we have included both an-isotropic effects. In the limits of slow rotation and weak magnetization, the diffusivity tensor becomes isotropic. When either a rapid disk rotation or a strong magnetization becomes relevant, the isotropy of the diffusivity tensor is broken. This is the first time, that such modeling with a higher degree of more self-consistency, has been applied in the context of jet launching simulations from accretion disks.

A reference Simulation
In order to investigate the effects of the feedback of the magnetic field on the diffusivity, we have chosen to focus on the case Ω * = 1 as a reference simulation.
In Fig. 11, we show the time evolution of the density and poloidal magnetic field of this simulation. We see that the saturation of the magnetic field occurs on a short timescale, i.e. already until t = 1000, corresponding to 50 revolutions of the inner disk. At this point in time, the magnetic energy is amplified by an order of magnitude, while the magnetic field lines are already opened up to a radius R = 70 in the outflow region.
Since our approach considers an α 2 -Ω dynamo, the toroidal magnetic field is amplified by the disk rotation and by the radial dynamo component (Mattia & Fendt 2020a). Although the dynamo components are quenched by the magnetic field, the Ω-effect still acts on the toroidal magnetic field component. This, combined with the shorter timescale of the Ω-effect (compared to the amplification of the poloidal field through the dynamo action), leads to a faster and stronger amplification of the toroidal magnetic field.
Note that our initial condition is that of a purely radial field. However, as in Fendt & Gaßmann (2018), we point out that the dynamo-amplified magnetic field should not depend on the initial conditions.
The formation of a magnetic loops (rooted at foot points of different radius in the accretion disk) strongly indicates, that, at least in the early evolutionary stages, the launching mechanism for this initial outflow is that of a tower-jet, thus a magnetic pressure driven-outflow Figure 11. Mass density (colors) and magnetic field lines of the reference simulation (Ω * = 1, CTQ feedback mode) at times t = 0, 250, 500, 750, 1000. Figure 12. Evolution of the poloidal magnetic field disk energy for different disk portions. Solid lines show the poloidal magnetic energy, while dashed lines show the total magnetic energy (poloidal + toroidal). The radii that are labeled denote the lower integration boundary, while the upper integration boundary is at the end of the domain, R = 100. (Lynden-Bell & Boily 1994;Lynden-Bell 1996). This is further supported by the inclination of the magnetic field towards the disk surface, as being not favorable for a magneto-centrifugal driving of the outflow. As the system evolves, the magnetic loops diffuse outwards, more and more poloidal field lines break up and the magnetic field geometry reaches the inclination required for a Blandford-Payne-like outflow. The system evolves further until a quasi-steady state is reached. At this point, the system consists of a highly magnetized accretion disk and a super-Alfvénic disk wind, which evolves into a high-speed outflow. The Alfvén surface is located at ≈ 5 thermal disk scale heights above the disk surface.

Amplification of the magnetic field
The amplification of the poloidal magnetic energy within the accretion disk is shown in Fig. 12. The amplification of the magnetic field is stronger in the innermost accretion disk regions, because of the combination of the Ω-and the α−effects, in agreement with Fendt & Gaßmann (2018) The poloidal magnetic energy is amplified by about 2 orders of magnitude and occurs, mostly, within t = 3000. By comparing the red line of Fig. 12 with the dotted green line of Fig. 3, we notice that the feedback on the diffusivity has a very minor impact on the disk poloidal magnetic field until t = 4000. However, the oscillations in the magnetic energy at t ∼ 5000, t = 6000, t = 8000 and t ∼ 12000 are a consequence of this novel diffusivity feedback model (see Section 4.2.3).
We observe a different evolution compared to our previous models without feedback on diffusivity (Mattia & Fendt 2020b). The magnetic energy does not undergo any intermittent decrease before the magnetic field reaches a quasi-steady state. We think that this difference can be explained by a combination of effects that rely on the chosen feedback model. At first, the quenching of the dynamo leads to a lower disk magnetic diffusivity, just because of the lower magnetization level at which the magnetic field saturates. Then, because of the lower magnetic diffusivity, the disk mass loss is less compared to a diffusive quenching model (where the dynamo tensor is not explicit suppressed), therefore the sound speed (which affects both dynamo and the diffusivity) shows no decrease. Moreover, the magnetic diffusivity is suppressed to an even lower level, just because of the feedback model (quenching of turbulent magnetic diffusivity). Figure 13. Mid-plane quantities at t = 50000 for the reference simulation. The solid lines represent the radial profiles of some selected MHD quantities at the disk mid-plane, while the dotted and dashed lines represent, respectively, their corresponding power-law approximation and their solution assuming self-similarity. On the y-axis are the corresponding variables, normalized in code units.
At t 20000 the magnetic field has saturated out to a radius R 50. The field still continues to be amplified in the outer disk region, as this part of the disk has performed so far only few rotations and is not yet in dynamic equilibrium. Also, since the α−effect is less efficient here, it takes more time to saturate the dynamo. Nevertheless, the evolution of the outer disk can weakly affect the disk evolution, most likely by triggering episodic ejections (see Section 4.2.3).

Disk structure
Because of the magnetic field amplification and the evolution of the magnetic field structure due to dynamo activity, also the disk structure evolves through time. This is partly due to the change of forces acting on the disk material, but also due to the change of angular momentum balance and the subsequent re-distribution of disk material.
At time t = 50000 the magnetic field has saturated in most parts of the disk (of the size we investigate). We may therefore investigate the disk structure by looking at the profiles of its mid-plane quantities. In Figure 13 we display the radial profile of some leading MHD quantities measured at the disk midplane and compare them with an idealized radial self-similar solution of the steady-state MHD equations (see Blandford & Payne 1982). We also compare with the power law approximation (obtained by performing a linear fit on the logarithmic radial profiles) for each quantity.
We see that the disk kinematics remains unaffected by the dynamo action as the disk rotation remains Keplerian, i.e. β v φ = −1/2. The density profile power-law index changes, however, from β ρ = −3/2 to β ρ = −5/4, as the mass is mostly accreted from the inner disk, and the whole disk looses only very little mass.
We find that both the sound speed and the Alfvén speed show very small deviations from the self-similar solution. The power-law indexes obtained are, respectively, β cs ≈ −4/9 and β v A ≈ −5/9. The radial dependence of the magnetization that is obviously strongly affected by dynamo action, can be recovered by computing the ratio of the Alfvén speed vs sound speed, which corresponds to computing the difference in the respective power-law indexes, We find a very good agreement with the results obtained by ; , suggesting that the properties of the saturated state do not depend (or depend only weakly) on the quenching model for dynamo action and diffusivity.

Intermittent ejection
As pointed before 4.2.2, also the system that undergoes consistent feedback concerning the dynamo action and magnetic diffusivity evolves to a quasi-steady state. However, we diskover an interesting, intermittent feature. Between t = 4000 and t = 12000 the magnetic energy features substantial oscillations, especially in the inner disk region (see Fig. 12). We now want investigate the origin of such processes and the consequences on the disk and jet structures. When the quenching models of Section 2.7 are adopted, the feedback applies only to the mean-field dynamo, finally leading to a saturation of the magnetic field. Note, however, that even if the dynamo is quenched, for a low diffusivity it can lead to a re-amplification of the magnetic field. This is what is happening if the diffusivity follows a consistent quenching model.
We point out that the amplification process of the magnetic field depends not only on the strength of the dynamo and the diffusivity tensors, but also on the radial and angular dependence of the tensor components. This is an essential property of the non-linear quenching models. The dependence of the magnetic field amplification on the radial and angular profiles for dynamo and diffusivity plays a key role in the dynamo and diffusivity models, in particular, when the magnetic field shows a reversal within the accretion disk. Under such circumstances, the magnetic diffusivity is sufficiently high in order to trigger reconnection processes, but, on the other hand, not strong enough to saturate the magnetic field amplification. As a result, flux ropes and current sheets are more prone to be formed in the accretion disk, Figure 14. Time evolution of the intermittent ejection for the reference simulation (Ω * = 1, CTQ feedback model). From left to right, snapshot of the poloidal velocity (colors), superimposed with magnetic field lines (defined as the contour lines of the vector potential), are shown at times t = 5500, 5880, 6080, 6500. On the bottom right part of each panel a 5 × 5 Rin cut from the inner most region is shown, indicating that the jet turbulence is ejected from the innermost disk region. before they are lifted above the disk surface, similar to what have been described in Yuan et al. (2009).
These flux ropes are able to reconnect above the disk surface in the disk corona, since the given magnetic diffusivity profile. Note that also the polarity of the toroidal field is opposite to the one of the launching region of the accretion disk. As a consequence, these flux ropes, after undergoing magnetic reconnection, are advected towards the accreting object and are able to disrupt the jet launching.
These features are highly interesting and may have essential relevance for jet launching conditions, such that this intermittent behavior may be related to the generation of jet knots. The understanding of the physical processes behind the flaring activity is not straightforward as they result from the temporal evolution of the highly non-linear resistive MHD equations.
Our understanding is as follows. The magnetic field geometry emerging from the dynamo activity and the subsequent opening of flux loops enhances the field reversal process typical of the mean-field dynamo action. While the field structure in these areas are prone to reconnection, the low field strength over there also implies a lack of magnetic pressure support. Further, since the reconnection area is not rooted in the disk via the magnetic field (no lever arm), it is only slowly rotating. In combination -lack of pressure support and centrifugal forces -gravity wins and leads to a collapse of this area towards the central object, thereby cutting off the inner disk wind.
Lateron, because of the opposite polarity, the magnetic field in the inner disk decreases, leading to a restoration of the dynamo. Then the magnetic field is reamplified, leading to a strong ejection and acceleration of the disk matter. As the field is amplified, the quenching on the mean-field dynamo saturates the magnetic field, which goes back to both strength and topology that it had before this episodic fast ejection.
The time period between to consecutive flares increases after each flare, since the reconnection process occurs further and further from the launching region. (the first flare appears at t = 5000, the second at t = 6000, the third at t = 7500 and the last at t = 10000). We expect less reconnecting plasma and thus less variability in the jet once the disk has reached a saturated state.

A parameter Study
In order to investigate the similarities and the differences between the consistent turbulence quenching and the dynamo quenching methods, we have performed simulation runs applying different Coriolis numbers Ω * .
The results are shown in Fig. 15. In the left panel we see the time evolution of the poloidal disk magnetic energy, that is essentially the dynamo-generated field amplification by the α-effect. We notice that for a higher Coriolis number, the differences between the non-isotropic dynamo quenching model and the consistent turbulence quenching model are only little. This finding is interesting as it may sound counter-intuitive -wouldn't one expect that a stronger magnetization is leading to a stronger quenching on the magnetic diffusivity, and therefore to more differences with the model without the consistent quenching?
However, note that a stronger dynamo implies also a stronger feedback on the dynamo (as shown in Section 3). Because of the feedback (quenching) on the diffusivity, in combination the α− and the Ω-effect lead to a stronger amplification of the toroidal magnetic field 6 . Since the dominating launching mechanism at the early evolutionary stages of jet launching is the toroidal pressure-dominated launching (tower jet), the quenching on the dynamo and the diffusivity is mainly triggered by the toroidal field.
As a result, the quenching of the diffusivity plays a minor role in the amplification and saturation of the poloidal field. Therefore, despite the quenching of magnetic diffusivity, which would be thought to lead to a higher dynamo number, we find that these numbers are almost identical as shown in Fig. 15 (right panel) for the case Ω * = 10 at t = 10000. This is consistent with the usual understanding that the dynamo number is a useful tool in order to characterize the onset of field amplification for a dynamo (in the disk).
On the other hand, although working in the accretion disk, the quenching on the diffusivity plays a key role also for the jet dynamics. The interrelation between the disk magnetization and the jet speed for the new quenching setup shows substantial differences from those by applying only a quenching on the dynamo (see Fig. 16). While we find a clear correlation between the disk magnetization and the jet speed in Section 3.4 above, the simulations with the consistent feedback thus the quenching of diffusivity show almost no correlation between these two quantities. On the other hand, the simulations with the consistent turbulence quenching and high Coriolis number show an outflow with an extremely high degree of collimation (see Fig. 17). This lack of correlation is also present at larger radii.
We recognize that at t = 10000 the magnetic loops, which are the main driver for the magnetic tower at initial evolutionary stages, are diffused outwards, and the outflow that is launched form the inner disk is now driven by the magneto-centrifugal mechanism, i.e. 6 which occurs on a timescale shorter than the amplification of the poloidal field launched sub-Alfvénic and subsequently supersedes in the Alfvén and the fast magnetosonic speed. However, because we also have strong amplification of the toroidal field, the collimation process is very rapid (in terms of spatial scales). The Alfvén surface is rather close to the disk surface where the jet is launched. As a consequence, magneto-centrifugal acceleration can happen only along a short distance and the final speed remains rather low. This is in particular true for the simulations with higher Coriolis number (and therefore stronger amplification of the magnetic field). The Alfvén surface moves further down to the disk surface. This is a very interesting effect. These jets are less efficient concerning the Blandford-Payne acceleration mechanism (when compared with jet launched by dynamo disks with different feedback models, as in, e.g., Fendt & Gaßmann 2018 or the previous section of this paper), but very efficient concerning the Blandford-Payne collimation mechanism. Note that the latter is indeed a consistent result as the material along collimated field lines cannot be accelerated magneto-centrifugally anymore.
In summary, we find that our consistent feedback model that quenches also the disk diffusivity strongly impacts the launching process. This holds in particular for a weak dynamo (low Coriolis number Ω * = 1), because of the magnetic field reversals that are induced in the jet launching region and also the intermittent disruption of the jet followed by the production of a flare. However, even for high Coriolis numbers, when the differences are less pronounced, the jet structure is severely affected by this new feedback model, as acceleration decreases and collimation increases. One may hypothesize about a "critical" or "optimal" Coriolis number for the dynamo that can produce the fastest jets. However, this issue need further detailed analysis, which is beyond the scope of this paper.
So far we have investigate only thin Keplerian disks (i.e. neglecting the non-diagonal components of the dynamo tensor). The differences we find by applying a small Coriolis number (representing a small effect of the rotation on the turbulence) suggest that we may expect even more structured and probably unstable jets that are produced by thicker accretion disks.

CONCLUSIONS
In this paper we presented MHD mean-field dynamo simulations in the context of large-scale jet launching. Our simulations were performed in axisymmetry, evolving all three vector components for the magnetic field and velocity. We employed the resistive MHD code PLUTO 4.3 (Mignone et al. 2007), which we have ex- Figure 15. Magnetic field evolution for the NDQ and CTQ quenching models. In the left panel is shown the Time evolution of the disk poloidal magnetic energy. In the right panel only the cases with Ω * = 10 are considered. Shown are the respective radial profiles for dynamo numbers, the φ components of magnetic diffusivity at midplane, and the mean-field dynamo α at half disk-height , all at t = 10000. tended by implementing the mean-field dynamo action as detailed in Mattia & Fendt (2020a).
In extension of our previous works on mean-field dynamo-driven jets Mattia & Fendt 2020b), in this paper we have essentially investigated the feedback of the dynamo-generated magnetic field on the dynamo activity and the disk diffusivity. This is a further step towards a consistent numerical modeling of mean-field disk dynamos.
In summary we have applied (i) different quenching models for the mean-field dynamo, and (ii) an analytically derived formalism for the mean-field dynamo and turbulent diffusivity, that consistently incorporates the suppression of turbulence by a strong and ordered magnetic field.
The following summarizes our approaches and results. 1) We have numerically investigated how different dynamo quenching models affect the magnetic field evolution and thus the jet launching process. More specifically, we compared the most common quenching strate-gies (Brandenburg & Subramanian 2005;) with the analytical model of Rüdiger, & Kichatinov (1993). The latter model has the advantage to be a more consistent approach, based on an analytical model of turbulent dynamo theory.
2) Essentially, we find that a stronger feedback by the magnetic field on the dynamo leads to a saturation of the magnetic field that is generated by the dynamo at lower disk magnetization. On the other hand, the so-called standard quenching or the diffusive quenching lead to a stronger saturated magnetic field. Nevertheless, the launching process and the jet structure that emerges are affected by the possible evolution of the dynamo inefficient zones, and that even for similar values of the disk magnetization.
3) The diffusive quenching model typically leads to a very stable disk-jet connection, just because of the continuous production of large-scale magnetic flux. However, the strong coupling in the model between the disk magnetization and the magnetic diffusivity, in combination with the absence of a quenching term of the dynamo tensor, may potentially lead to unphysically high values of diffusivity. This problem has been solved by applying a more consistent quenching model on the dynamo and a diffusivity model where α ss ∝ √ µ D .

4)
In agreement with previous studies, when we apply a feedback only on the dynamo, we recover an interrelation between the disk magnetization and the outflow speed regardless of the Coriolis number (i.e. the strength of the dynamo) or the feedback model. While the feedback model plays a key role for in the saturation of the magnetic field, the relations between the inflowoutflow parameters seems to be independent of the dynamo/diffusivity model applied. This interrelation holds only when no quenching of the magnetic diffusivity is applied (see point 8 of this section). 5) We applied and investigated the effects of a more consistent quenching model which encompasses the suppression of the turbulence by a strong ordered magnetic field, for both the dynamo tensor (Rüdiger, & Kichatinov 1993) and the magnetic diffusivity Rüdiger et al. 1994). Such an approach has never been applied in the context of jet launching by dynamo-generated magnetic fields. 6) We found that, in the early evolutionary stages, the jet is driven by the magnetic pressure. Once the magnetic field has saturated in the inner disk region, and the magnetic loops are opened up and their central part has diffused outwards, the magneto-centrifugal Blandford-Payne outflow is produced. 7) We found that, by applying a consistent turbulence quenching model, reconnection processes lead to the formation of flux ropes (with opposite magnetic field polarity respect to the disk magnetic field) that are accreted and disrupt the jet. As a consequence, the magnetic field is re-amplified in the launching accretion disk region, leading to a very fast intermittent ejection. 8) When applying the quenching model that consistently quenches dynamo and diffusivity, for higher Coriolis numbers, we do not find the established interrelation between jet speed and disk magnetization. Instead, the high Coriolis number is associated with a more collimated jet. The strong toroidal field induced leads to rather short acceleration distances, such that these jets gain only little speed, but a high degree of collimation.
In this paper we have concentrated on thin Keplerian disk, an assumption that also constraints the strength of certain dynamo tensorial components (then strictly related to the non-diagonal components of the mean-field dynamo). However, we expect that the variety of disks found around astrophysical objects, the disk thickness may play a key role in the launching process, and may produce a variety in the magnetic field geometry and, subsequently, in the jets produced. Our future simulations will investigate the dynamo action on thick disks and will be presented in a forth-coming paper.