Relativistic outflows from a GRMHD mean-field disk dynamo
DDraft version February 26, 2021
Typeset using L A TEX twocolumn style in AASTeX63
Relativistic outflows from a GRMHD mean-field disk dynamo
Christos VourellisChristian Fendt
Max Planck Institute for Astronomy, Heidelberg, Germany
ABSTRACTWe present simulations of thin accretion disks around black holes investigating a mean-field diskdynamo in our resistive GRMHD code (Vourellis et al. 2019) that is able to produce a large scalemagnetic flux. We consider a weak seed field in an initially thin disk, a background (turbulent)magnetic diffusivity and the dynamo action. A standard quenching mechanism is applied to mitigatethe otherwise exponential increase of the magnetic field. Comparison simulations of an initial Fishbone-Moncrief torus suggest that reconnection may provide another quenching mechanism. The dynamo-generated magnetic flux expands from the disk interior into the disk corona, becomes advected by diskaccretion, and fills the axial region of the domain. The dynamo leads to an initially rapid increasein magnetic energy and flux, while for later evolutionary stages the growth stabilizes. Accretiontowards the black hole depends strongly on the magnetic field structure that develops. The radial fieldcomponent supports extraction of angular momentum and thus accretion. It also sets the conditionsfor launching a disk wind, initially from inner disk area. When a strong field has engulfed the disk,strong winds are launched that are predominantly driven by the pressure gradient of the toroidal field.For rotating black holes we identify a Poynting flux-dominated jet, driven by the Blandford-Znajekmechanism. This axial Poynting flux is advected from the disk and therefore accumulates at theexpense of the flux carried by the disk wind, that is itself regenerated by the disk dynamo.
Keywords: accretion, accretion disks – dynamo – MHD – ISM: jets and outflows – black hole physics– galaxies: nuclei – galaxies: jets INTRODUCTIONOne of the fundamental prerequisites for the launch-ing of astrophysical jets is the support of a strong mag-netic field. A few mechanisms have been proposed thatcan explain the launching, acceleration, and collimationof these outflows, all of them further requiring the ex-istence of an accretion disk orbiting a central object(star or black hole). These processes focus either onjet launching from the surface of an accretion disk viamagneto-centrifugal acceleration (Blandford & Payne1982; Uchida & Shibata 1985; Pudritz & Norman 1986)or driven by a magnetic pressure gradient (Lynden-Bell1996), or from a black hole ergosphere (Blandford &Znajek 1977; Komissarov 2005), but commonly rely on asubstantial magnetic field strength and a favorable mag-netic field geometry. [email protected], [email protected]
Based on these theoretical ideas, a number of mag-netohydrodynamic (MHD) codes were developed, usingrelativistic gravity, in order to simulate the rotation ofthe accretion disk around a black hole and the subse-quent jet launching (Koide et al. 1999; Gammie et al.2003; De Villiers & Hawley 2003; Noble et al. 2006; DelZanna et al. 2007; Porth et al. 2017; Ripperda et al.2019). Jet launching simulations in the non-relativisticlimit have been pioneered by Casse & Keppens (2002)and then further developed (Zanni et al. 2007; Murphyet al. 2010; Sheikhnezami et al. 2012; Stepanovs & Fendt2014, 2016). The usual practice in these works is to ap-ply a strong magnetic field as part of the initial condi-tions. However, the strength and shape of the field mayexpress an unrealistic view of the actual physical case.Only after the disk-field system has evolved we can saythat we have approached a more authentic view of thesystem. Naturally, there is the downside that the originand the ab-initio evolution of the jet-launching magneticfield is practically unknown. a r X i v : . [ a s t r o - ph . H E ] F e b Vourellis & Fendt
While for protostellar jets and jets from compact starsthe jet launching magnetic field could in principle be in-duced by a strong stellar dipolar field (if not advectedfrom the interstellar medium), this option does not ex-ist for supermassive black holes. In AGN the jet-drivingmagnetic field thus must be generated by some kind of accretion disk dynamo mechanism , or be advected fromthe surrounding medium. The situation can be differ-ent for stellar mass black holes that form by mergersof magnetized neutron stars. Also, the environment ofstellar mass black holes resulting form core collapse su-pernova explosions could be highly magnetized (Cerdá-Durán et al. 2007; Obergaulinger et al. 2018; Matsumotoet al. 2020).Early studies have derived the structure of black holemagnetospheres in the stationary approach (Fendt 1997;Ghosh 2000; Fendt & Greiner 2001) (but see also Panet al. 2017; Mahlmann et al. 2018; Pu & Takahashi2020), as well as jet acceleration (Fendt & Greiner2001; Vlahakis & Königl 2001; Tomimatsu & Takahashi2003). Such steady-state modeling of black hole magne-tospheres could be very useful as a basis to study par-ticle acceleration and radiation (Takahashi et al. 2018;Rieger 2019). While steady state solutions are exact so-lutions of the physical equations, their stability remainsunclear. The same holds for the origin of the magneticfield structure considered.The accretion disk turbulence is strongly believed tobe generated by instabilities like the magneto-rotationalinstability (MRI) (Balbus & Hawley 1991, 1998). At thesame time, the MRI is known to amplify the magneticfield, giving rise to a turbulent dynamo effect. This istypically been investigated using high-resolution, idealMHD simulations (see e.g. Stone et al. 1996; Gressel &Pessah 2015) which are sometimes called direct dynamosimulations, as the dynamo effect is acting ab initio andwith no further prescriptions solely based the turbulentmotions of the gas as soon as a weak seed field is present.Under certain conditions, these small-scale fluctua-tions of velocity and magnetic field can be subject toa non-linear coupling, which then effectively representsa non-ideal MHD mechanism that amplifies the mag-netic field on larger scales. This mechanism is known as mean-field dynamo (Parker 1955; Steenbeck & Krause1969a,b) and enters the MHD induction equation as anon-ideal term. Similarly, the small-scale turbulence ef-fectively results in a turbulent resistivity (or diffusivity)in the induction equation. These small-scale turbulentmotions may be averaged resulting in a turbulent pat-tern on larger scales, typically leading to much strongercoefficients for the mean-field dynamo effect and the re-sistivity if compared to the small-scale values. The mean-field dynamo theory is a powerful tool toderive the large scale field structure in astrophysical ob-jects. In the limit of a kinematic dynamo the velocityfield is prescribed and remains constant while for largermagnetic field strength feedback from the field on theflow dynamics is expected and the non-linear dynamotheory must be considered.Overall, it is still not fully clear if such turbulent pro-cesses can result in the generation of a well ordered large-scale magnetic field that is needed to launch an outflow.Direct high-resolution simulations typically demonstratethe saturation of the turbulent dynamo by showing themagnetic energy evolution. However, they will hardlyprovide an estimate for the production of the large-scalemagnetic field. What is interesting for jet launching,however, is the generation of a large scale magnetic flux .This has been been achieved by simulations applying themean-field ansatz in non-relativistic studies (see below).The present paper aims at extending these studies tothe GRMHD context.Concerning the application of the mean-field dynamotheory for accretion disk and jet simulations, the lit-erature is rather sparse. Dynamos have in fact besuggested for generating the magnetic fields needed tolaunch jets early on. Pudritz (1981a,b) have applied the α Ω-dynamo theory for the case of thin accretion disks,demonstrating the growth of a magnetic field by meansof differential rotation of the disk (the Ω-effect), togetherwith the net effect of turbulent motions (the α -effect) ina vertically stratified turbulent disk.Simulations of magnetized shear flows by Bran-denburg et al. (1995) demonstrated how a dynamo-generated magnetic field can amplify the turbulence inthe flow, which, in turn, can amplify the magnetic fieldvia a dynamo process. Simulations of outflow launchingfrom a mean-field dynamo active accretion disk werefirst presented by von Rekowski et al. (2003) and vonRekowski & Brandenburg (2004). In Pariev et al. (2007)a possible origin of the dynamo mechanism in the caseof AGN accretion flows was suggested, triggered by apassing star that can heat up and perturb the magneticfield inside the accretion disk, altogether resulting in theinduction of a toroidal field into a poloidal magnetic flux.Stepanovs et al. (2014) and Fendt & Gaßmann (2018)presented simulations of jet launching from accretiondisk, where the jet driving magnetic field was self-generated by a mean-field disk dynamo. Different mag-netic field structures and thus jet parameters were ob-tained by exploring mean-field dynamos of differentstrength. Also, by switching off and on the α -dynamomechanism (externally) a periodic ejection of jets could GRMHD disk dynamo α , oscillating dynamomodes may occur resulting as well in pulsating ejectionsinto the jet (see also Dyda et al. 2018).We now turn to mean-field dynamo models in the rel-ativistic context. Early on the existence of a gravito-magnetic dynamo effect has been claimed (Khanna &Camenzind 1996a), which could, however, not realized insubsequent studies (Brandenburg 1996; Khanna & Ca-menzind 1996b). Fully dynamical GRMHD simulationsof mean-field dynamos have considered only recently.Applying GRMHD, Sądowski et al. (2015) parametrizedthe effect of a mean-field dynamo in combination withradiative transfer. However, their approach was in idealMHD regime neglecting magnetic diffusivity. Also a sim-plified dynamo model was used, where the dynamo ac-tion appears as small numerical correction for the evo-lution of poloidal magnetic field that drives it towards asaturated state. They applied their code to investigatethe evolution of thick disks at different accretion rates.The first fully covariant implementation of a dynamoclosure in a general relativistic context was presentedby Bucciantini & Del Zanna (2013) following the 3+1formalism. A mean-field dynamo was implemented inthe ECHO code (Del Zanna et al. 2007; Bucciantini &Del Zanna 2011), that was later applied to simulationsof kinematic dynamo action in tori around black holesdemonstrating the growth of the toroidal and poloidalfield components (Bugli et al. 2014). Recently, Tomeiet al. (2020) have extended these studies by simulatingfully dynamical cases including dynamo quenching.In our present paper, we aim for a broader descriptionof the mean-field dynamo in GRMHD simulations. Incomparison to the previously mentioned work we con-sider in particular the following points.(i) Compared to Bugli et al. (2014) in our paper wedo consider the feedback of the magnetic field on thedynamics of the system. A fully dynamical MHD studybrings a major advance over a kinematic study, and isparticularly essential concerning disk accretion and jetlaunching.(ii) Compared to Tomei et al. (2020) who consider thedynamical feedback by the field and also discuss the dis-tribution of the poloidal field energy, we will investigatein detail the geometry of the poloidal field that is gener-ated including the overall structure of the magnetic fieldlines. This type of analysis was not performed by Bugliet al. (2014).(iii) Compared to Tomei et al. (2020) we also showthe evolution of hydrodynamical quantities. This is in particular essential for a discussion of the disk accretionand the jet launching.We apply our resistive
GRMHD code rHARM3D (Qianet al. 2017, 2018; Vourellis et al. 2019) and run a seriesof fully dynamical simulations with an accretion diskmean-field dynamo. In particular, we aim to demon-strate the growth of the poloidal magnetic field, thus thegeneration of a large scale magnetic flux that is neededto launch winds or jets from the accretion disk.Our paper is structured as follows. In Section 2 wereview the basic GRMHD equations, discuss the mean-field dynamo including its implementation as well as itsquenching. In Section 3 we apply the dynamo effect forthe setup of a relativistic torus, which allows us to com-pare and test our results with the existing literature. InSection 4 we detail the initial conditions for our setup ofa thin disk mean-field dynamo. In Section 5 we presentsimulation results, discussing the evolution of the mag-netic field and the the accretion disk, the mass fluxes,and the launching of disk winds and relativistic jets. InSection 6 we summarize our work. THEORETICAL BACKGROUNDIn this section, we briefly review the basic equationsof resistive GRMHD as a basis for our calculations. Forthe metric we adopt the signature ( − , + , + , +) (Misneret al. 1973) and apply geometrized units, G = c = 1.Thus, length scales are expressed in units of the gravi-tational radius R g = GM/c , while time is measured inunits of the light travel time t g = GM/c . Vector quan-tities are written with bold letters while the vector andtensor components are indicated with their respectiveindices, with Greek letters running for 0,1,2,3 ( t, r, θ, φ )and Latin letters running for 1,2,3 ( r, θ, φ ).The usual 3 + 1 decomposition for GRMHD is used inorder to separate the time component from the spatialcomponents (3-dimensional manifolds). The space-timeis described by the metric g µν in Kerr-Schild coordinateswith g ≡ det ( g µν ). A zero angular momentum observerframe (ZAMO) exists in the spacelike manifolds, movingonly in time with the velocity n µ = ( − α, , , α = 1 / p − g tt is the lapse function. The gravitationalshift is β i = α g ti .In our simulations we apply our resistive coderHARM3D (Vourellis et al. 2019), now extended by in-cluding a turbulent mean-field dynamo, following thework of Bucciantini & Del Zanna (2013).2.1. Basic GRMHD equations
The Maxwell equations in covariant form ∇ ν F ∗ µν = 0 , ∇ ν F µν = J µ , (1) Vourellis & Fendt along with the conservation of the stress-energy tensor ∇ ν T µν = 0 , (2)describe the motion of a magnetized fluid in a generalrelativistic environment. J µ is the 4-current that satis-fies the electric charge conservation ∇ µ J µ = 0, T µν isthe stress-energy tensor and ∇ µ is the covariant deriva-tive. The electromagnetic field is described by the Fara-day and Maxwell tensors, F µν = u µ e ν − e µ u ν + (cid:15) µναβ u α b β ,F ∗ µν = − u µ b ν + b µ u ν + (cid:15) µναβ u α e β , (3)respectively, where u µ is the 4-velocity and e µ and b µ express the electric and magnetic field in the fluid restframe. The Levi-Civita symbol (cid:15) µνκλ is defined as (cid:15) αβγδ = √− g [ αβγδ ] , (cid:15) αβγδ = − √− g [ αβγδ ] , (4)where [ αβγδ ] are the conventional permutation sym-bols. The magnetic and electric field as measured by thezero angular momentum observer (ZAMO) are definedas B µ = − n ν F ∗ µν and E µ = n ν F µν , with B i = αF ∗ it , E i = − αF it , B = E = 0. The stress-energy tensor T µν can be split into a fluid and an electromagnetic compo-nent. The fluid component is T µν f luid = ( ρ + u + p ) u µ u ν + pg µν , (5)where ρ is the mass density, u is the internal energydensity and p is the thermal pressure. Pressure andinternal energy are connected through the equation ofstate for an ideal gas with u = p Γ − , (6)where Γ is the polytropic exponent. The electromag-netic component of the stress-energy tensor is T µν E M =( b + e ) (cid:18) u µ u ν + g µν (cid:19) − b µ b ν − e µ e ν − u α e β b γ (cid:0) u µ (cid:15) ναβγ + u ν (cid:15) µαβγ (cid:1) (7)Equations (1), (2) along with (6) are being solved byour code as a system of hyperbolic differential equa-tions. For a more detailed analysis of the equationsand their numerical implementation and solution we re-fer the reader to Bucciantini & Del Zanna (2013), Qianet al. (2017) and Vourellis et al. (2019).2.2. The mean-field dynamo
In ideal MHD, the electric field vanishes in the co-moving frame which is expressed as
EEE + vvv × BBB = 0 , (8) where vvv, EEE, BBB are the velocity and the electric and mag-netic field. In the more general case of non-ideal MHD,Ohm’s law is given by σE E E = σ ( EEE + vvv × BBB ) =
JJJ, (9)where E E E is the co-moving electric field, JJJ is the electriccurrent and σ is the electric conductivity.In the mean-field ansatz the turbulent fluctuations ofvelocity and magnetic field lead to an averaged, meanelectromotive force, δvvv × δBBB = α D BBB − β D JJJ, (10)that subsequently enters the induction equation. Here, α D describes the term that is responsible for the gener-ation of the magnetic field from the turbulent motions,the so-called α -effect, while the β d provides an addi-tional resistive term, that enhances the dissipation ofthe magnetic field by an enhanced turbulent diffusivityin addition to the ohmic resistivity (Parker 1955; Steen-beck & Krause 1969a,b).We can thus write the induction equation as ∂BBB∂t = ∇ × ( vvv × BBB ) + ∇ × αBBB − ∇ × ηJJJ, (11)where we now have defined α ≡ α D and the magneticdiffusivity η considering both the turbulent and ohmiccomponents. Renaming the coefficients also emphasizesthe fact that the mean-field dynamo theory starts fromthe ideal MHD when considering the correlated turbu-lent fluctuations, but after performing the averaging ofthe electromotive force term, an extra term for the non-ideal treatment of MHD appears.The first term on the r.h.s. of the induction equationrepresents a kinematic term that may induce a poloidalmagnetic field from a toroidal motion, the so-called Ω-effect. Combined with the action of the α -effect, thismay provide a dynamo cycle that induces a poloidalfield from a toroidal field ( α -effect) that then inducesa toroidal field from the poloidal component (Ω-effect),and so forth - the so-called α Ω-dynamo. Note that alsothe α -effect may induce a toroidal field from a poloidalcomponent. This is in particular important in stellardynamos with small shear in the toroidal motion, thusforming an α dynamo. However, for strong shear, forexample as in accretion disks following a Keplerian ro-tation or in relativistic tori, the Ω-effect will dominatethe α Ω-dynamo concerning the induction of the toroidalfield component.In a more general definition of the mean-field elec-tromotive force, the α D and the β D will represent ten-sors (Moffatt 1978). However, in most applications (see GRMHD disk dynamo α and η are spatially constant, theycan be moved in front of the curl operator in Equa-tion (11).We want to briefly connect to the direct dynamo simu-lations that were mentioned in the introduction. Thosesimulations, applying high resolution, actually may re-solve the turbulence pattern in the MHD flow and thusperform the averaging equation 10, thus providing an α D and β D , that can be than used for simulations inthe mean-field approach. As an example, we cite hereGressel (2010), who applied accretion disk shearing boxsimulations in order to derive profiles for the magneticdiffusivity and the dynamo coefficients. As complemen-tary to the mean-field approach, the direct simulationsserve as subgrid modeling , that can deliver the mean-fieldquantities from a subgrid scale.2.3. The GRMHD mean-field dynamo equations
The implementation of the mean-field dynamo mech-anism builds up on our previous works in which we haveintegrated resistivity in the form of turbulent diffusiv-ity in the original HARM code, establishing our resistiveGRMHD code rHARM3D (Qian et al. 2017, 2018; Vourelliset al. 2019).For our GRMHD treatment, we follow the closure re-lation introduced by Bucciantini & Del Zanna (2013).Here, the mean-field α -dynamo parameter is replacedby ξ = − α , just to avoid confusion with the gravita-tional lapse. Thus, in covariant form, in the fluid frameOhm’s law is written as e µ = ηj µ + ξb µ , (12)where j µ is the 4-vector of the electric current density.Now, the electric field can no longer be calculated by thecross product of fluid velocity and magnetic field andnew equations need to be formulated (see next Section).By setting ξ = 0 we get the resistive version of Ohm’slaw ( e µ = η j µ ). By setting η = 0 we get back to theideal MHD case e µ = 0.In order to derive the equations for the evolution of themagnetic and electric field we start by taking the tem-poral and spatial projections in Equation (1) separately.This provides the divergence condition ∂ j (cid:0) γ / B j (cid:1) = 0and an equation for the time evolution of the magnetic field, γ − / ∂ t (cid:16) γ / B i (cid:17) + (cid:15) ijk ∂ j ( E k α + (cid:15) knm B m β n ) = 0 , (13)along with Gauss’ law for the electric field ∂ j (cid:0) γ / E j (cid:1) = γ / q and the time evolution of the electric field, γ − / ∂ t (cid:16) γ / E i (cid:17) − (cid:15) ijk ∂ j ( α B k − (cid:15) knm E m β n ) == (cid:0) qβ i − αJ i (cid:1) . (14)where γ is the determinant of the spatial metric γ µν = g µν + n µ n ν . From Ohm’s law, Equation (12), we can get F µν u ν = ηI µ + η ( I ν u ν ) u µ + ξF ∗ µν u ν , (15)where I µ is 4-current as is seen by the ZAMO. Blackman& Field (1993) derived this equation only for the caseof a resistive fluid, however, the addition of a mean-field dynamo term is straightforward. The source term I ν u ν = − q is the electric charge density as measuredin the fluid frame (Komissarov 2007).By taking the temporal decomposition of Equation 12we get an expression for the electric charge in the formΓ E i v i = η ( q − Γ q ) + ξ Γ B i v i , (16)and from the spatial projection of Equation 12 we getan expression for the electric currentΓ (cid:2) E i + (cid:15) ijk v j B k − ( E k v k ) v i (cid:3) = η (cid:0) J i − qv i (cid:1) + ξ Γ (cid:2) B i − ( B k v k ) v i − (cid:15) ijk v j E k (cid:3) , (17)which we can now use to replace the source terms inEquation (14), resulting in the evolution equation forthe electric field γ − / ∂ t (cid:16) γ / E i (cid:17) − (cid:15) ijk ∂ j ( α B k − (cid:15) knm E m β n ) == − α Γ (cid:2) E i + (cid:15) ijk v j B k − ( E k v k ) v i (cid:3) /η + αξ Γ (cid:2) B i − (cid:15) ijk v j E k − ( B k v k ) v i (cid:3) /η − q ( αv i − β i ) . (18) Vourellis & Fendt
Discretizing the time evolution of the electric fieldleads after some “minor” algebraic calculations to E i (cid:20) ˜ η + Γ + ξ (Γ − η + Γ (cid:21) = − (cid:15) ijk ˜ v j B k + ˜ η (cid:20) Q i + ˜ v i ( Q k ˜ v k )Γ˜ η + 1 (cid:21) + ξ (cid:20) Γ B i − ˜ η ˜ v i Γ˜ η + 1 ( B k ˜ v k ) (cid:21) − ξ (cid:26) ˜ η(cid:15) ijk ˜ v j Q k + (Γ − B i − (˜ v j B j )˜ v i Γ + ˜ η (cid:27) − ξ Γ (cid:15) ijk ˜ v j B k Γ + ˜ η + ξ ˜ v i (cid:26) Γ˜ ηQ k ˜ v k + ξ B k ˜ v k (Γ˜ η + 1)(Γ + ˜ η ) (cid:27) , (19)where Q i = E i + ∆ t (cid:2) − (cid:0) αv i − β i (cid:1) q + (cid:15) ijk ∂ j (cid:0) α B k − (cid:15) klm β l E m (cid:1)(cid:3) ,q = γ − / ∂ k (cid:16) γ / E k (cid:17) , ˜ η = ηα ∆ t , ˜ v i = Γ v i . with E i representing the electric field as calculated inthe previous time step.Unfortunately, Equation (18) appears to be stiff. Forlow values of diffusivity, the terms on the r.h.s. evolvein different times than the timestep used for the timediscretization of of the electric field, resulting in unsta-ble solutions. The problem can be fixed by isolatingthe stiff terms and use an implicit method to calculatethem. The non-stiff part of the electric field, that isthe one which does not include diffusive terms, is cal-culated separately, and is denoted as Q i . For the stiff part we use an iterative method where for every timestep we calculate the electric field until its values haveconverged to the desired accuracy. A detailed analysisof the numerical implementation is presented in Qianet al. (2017).We clearly state that we have tested the implementa-tion of resistivity considering analytical, time-dependentsolutions of the induction equation (Qian et al. 2017;Vourellis et al. 2019) and find perfect agreement. Now,since the additional implementation of a mean-field dy-namo is achieved in exactly the same way as for the We note a correction here, namely the differently placed paren-theses for the term on the r.h.s of the 1 st line, and the differentsign for the term on the 4 th line in comparison to Bucciantini &Del Zanna (2013). resistivity, so our testing of resistivity supports the im-plementation of the dynamo term as well.2.4. Dynamo quenching
The exponential increase of the magnetic fieldstrength is a direct result of the mean-field dynamo pro-cess. However, a strong magnetic field will suppressthe MRI which is considered as one of the sources ofthe disk turbulence. Therefore, with the repression ofturbulence, also the dynamo will be suppressed, result-ing in a saturation level of the magnetic field generated.This has been demonstrated by direct simulations of theMRI (Stone et al. 1996; Gressel 2010; Gressel & Pes-sah 2015) and has been approved analytically (Ruedigeret al. 1995). The supression of MRI has recently beconsidered also in GRMHD (Bugli et al. 2018).In our approach, in absence of a more detailed turbu-lence model, the dynamo quenching must be enforced byinteractively lowering the value of the dynamo parame-ter, ξ , following certain physical criteria, motivated bymodels of turbulence.The usual procedure (that we refer to as "stan-dard quenching") prescribes a kind of equipartition fieldstrength beyond which the mean-field dynamo param-eter is reduced, until the dynamo becomes effectivelyinactive. For example, Bardou et al. (2001) have im-plemented a back-reaction on the α -dynamo parameterfrom the magnetic field by making this alpha α dependon the values of (see also von Rekowski et al. 2003;von Rekowski & Brandenburg 2004 ). A different ap-proach was undertaken by Fendt & Gaßmann (2018)who applied a disk diffusivity model where the strengthdynamo was quenched by increasing the magnetic dif-fusivity. This resulted an a smooth and long-term evo-lution for their jet launching simulations (up to several100000 disk rotations).However, the rate in which the magnetic diffusivityand the mean-field dynamo affect the magnetic field evo-lution can be very different.A leading parameter of mean-field dynamo action isthe dynamo number. For the case of a thin accretiondisk of scale height H , we apply the dynamo parameter ξ and magnetic diffusivity η for defining the total dynamonumber D ξ , following previous works in a non-relativisticenvironment utilizing a thin accretion disk (Bardou et al.2001; von Rekowski et al. 2003; Stepanovs et al. 2014), D ξ = |R ξ R Ω | = (cid:12)(cid:12)(cid:12)(cid:12) ξHη SH η (cid:12)(cid:12)(cid:12)(cid:12) (21)where the two fractions represent the Reynold’s num-bers R from the dynamo action turbulence (mean-field GRMHD disk dynamo α -dynamo) and from differential rotation (Ω-dynamo),respectively. The function S = r d Ω /dr measures theshear due to differential rotation and is calculated foreach cell separately.Depending on the values of the dynamo number, ξ will exponentially increase the magnetic field, while η would damp with a rate close to linear. Thus, if thedynamo number is too large, the change in diffusivityas suggested by Fendt & Gaßmann (2018) might not besufficient and a standard quenching model must also beemployed. Furthermore, by increasing the diffusivity,the numerical time step drastically decreases, increasingthe computational costs of the simulation.In our models, we follow the standard quenching ase.g. applied by Bardou et al. (2001). We define theplasma- β as the ratio between the gas and the magneticpressure and magnetization µ as its inverse β = 1 µ = p gas p magn . (22)We implement a quenching prescription that calculatesthe magnetization µ of the fluid and when it becomestoo large, the actual dynamo parameter ξ q is reduced, ξ q = ξ
11 + µ D /µ e q , (23)where ξ q is the quenched dynamo parameter, ξ its ini-tial value, µ D the disk magnetization averaged verti-cally over the disk and µ e q is a equipartition magne-tization, expressing the value of the disk magnetizationin which the mean-field dynamo is quenched by 50% (i.e. ξ q = ξ / θ direction so that thecalculation of the quenching is run by the same core inthe angular direction. For our 256 ×
256 grid we use 64cores with 4 cells each in the radial direction and 5 coreswith 51 cells each in the polar direction resulting in atotal of 320 cores. THE DYNAMO EFFECT IN A RELATIVISTICTORUSIn this section we investigate the mean-field dynamoeffect in a relativistic torus. This setup is essentially dif-ferent from that of a thin accretion disk in the sense thatthe differential rotation of a torus is stronger comparedto a disk, thus the Ω-effect of the mean-field dynamoplays a stronger role. This becomes clear if we comparethe dynamo numbers for the cases of the disk and thetorus. In a torus with a constant specific angular mo-mentum the angular velocity scales as ∝ r − while fora Keplerian rotating disk the scaling follows ∝ r − / .We assume a scale height H D = 0 . R for a thin disk astypical length scale for the disk dynamo action. Sim-ilarly, as typical length scale for the torus dynamo wetake the e − folding length scale H T of the density distri-bution. Considering the centre of the torus ( r = 15), thefirst contour is at z ≈ .
2, so H T /R = 0 . R T = R ∂ Ω T ∂R H η = 0 . η , (24) R D = R ∂ Ω D ∂R H η = 0 . η √ R. (25)Note that the dynamo number of the torus can be higherby an order of magnitude in the inner area, for largerradii, however, the disk dynamo number increases.We can then write R T / R D = 0 . . √ R ≈ R − / . (26)This relation changes of course with time as the disksevolve. Compared to these systems with strong shear,stellar dynamos are better described by a α mean-fielddynamo, and are found to be less stable and less sym-metric (Küker & Rüdiger 1999).Our simulation setup for this exemplary torus simu-lation simT is as follows. We start the simulation witha density and internal energy distribution following theEquation of State for an ideal gas, which reads as ρ = (cid:18) ( h − γ − Kγ (cid:19) / ( γ − , (27)where h ( r, θ ) is the specific enthalpy as it is calculatedin Fishbone & Moncrief (1976). The torus is rotatingaround a Schwarzschild black hole ( a = 0) with an inneredge at r in = 6 and its point of maximum density at r ρ max = 15. The internal energy is defined by the poly-tropic equation of state u = Kρ γ /γ −
1, with K = 10 − and γ = 4 / Vourellis & Fendt
Figure 1.
Initial conditions for our reference simulation.Shown is the initial distribution of the density, includingcontours separated by one e − folding scale height (logscale;upper left), the toroidal magnetic field B φ (linear scale; up-per right), the diffusivity η (linear scale; lower left), and thedynamo parameter ξ (linear scale; lower right). An initial toroidal magnetic seed field is defined insidethe torus, following the density distribution ( B φ ∝ ρ ) ofthe standard Fishbone-Moncrief torus applied in HARM(Gammie et al. 2003). This guarantees that the mag-netic field is confined inside the torus. The field direc-tion is positive in both hemispheres and it is normalizedconsidering a plasma- β = 10 .We prescribe a distribution for the magnetic diffusiv-ity that follows the density distribution of the torus witha maximum value of η = 10 − . The diffusivity basi-cally vanishes in the surrounding corona. Also the dy-namo action is limited to the area inside the torus. It isconstrained to an area smaller than distribution of theinitial magnetic field and the diffusivity in order to avoidfield generation in the corona of the torus. The dynamoparameter is constant in both hemispheres at a level of ξ = 10 − .The parameters of our torus simulation simT are sum-marized in Table 1. The initial conditions for this sim-ulation are shown in Figure 1. Figure 2.
Generation of a dipolar poloidal field. The colorgradient shows the density distribution in log scale and thewhite lines show the poloidal magnetic field. The initial con-dition starts with a toroidal field with high plasma- β ∼ .The poloidal field appears immediately in the area wheredynamo exists ( t = 1, upper left) and it evolves through t = 1000 , R ρ max = 15. Table 1.
Summary of parameters of the torus simulationrun. Listed is the ID of the run; the Kerr parameter a applied; the maximum magnetic diffusivity η ; the initialplasma- β and the dynamo parameter ξ .run a η β ξ simT − − We note that this setup is somewhat different from theone chosen by Tomei et al. (2020), as their initial torusis that of a so-called
Polish donut (see e.g. Abramowicz& Fragile 2013) and their seed field geometry is thatof poloidal loops inside the torus. With our approachwe prefer to connect our simulations to the literature ofHARM simulations which typically apply the Fishbone-Moncrief torus as initial condition, and While our initialsetup may look similar to Tomei et al. (2020) our resultscannot directly compared.
GRMHD disk dynamo
Field induction by a spatially constant dynamo
We now describe the field evolution of a torus inducedby a mean-field dynamo in more detail.When the simulations starts, the poloidal magneticfield appears immediately inside the torus in the areawhere B φ and dynamo coexist. As the simulation con-tinues, the field starts increasing in value due to the dy-namo while advecting towards the black hole followingpartially the accretion of the torus increasing the mag-netic field in the inner disk part and towards the blackhole. Simultaneously, diffusivity is trying to dampen themagnetic field. In the end, the R ξ determines where thefield is amplified or dampened. Since ξ is constant, theprofile of diffusivity determines the value of R ξ . As anexample of their values, in the center of the torus thenumbers are R ξ ∼
15 while its maximum values reaches R ξ ∼
150 in the boundaries of the dynamo distribution.In Figure 2 we see the snapshots of the evolution ofthe poloidal magnetic field. As mentioned before, in thebeginning, the field lines are restricted in the area where ξ = 0. However, they are eventually dragged along withthe material that is accreted towards the black hole. Attime t = 2000, as the dynamo has been working for ap-proximately 6 rotations of the torus center, we see a lowvelocity outflow being launched from the inner part ofthe torus. The magnetic field lines follow the low densityfluid showing the first indications of the development ofa jet. The launching point of the outflow is barely in-side of the innermost stable circular orbit (ISCO) andcan be attributed to the presence of the strong toroidalmagnetic field we see in the close atmosphere of the in-ner part of the torus. This strong toroidal field, whichhas been amplified by the Ω effect of the dynamo, canincrease the magnetic pressure and push material andpoloidal field outwards (Lynden-Bell 1996).The poloidal field continues to grow up to t ∼ β has reached values between 0 . − . . c (see Figure 3). The acceleration is sup-ported by magnetic pressure forces. This holds for thehigh speed floor material, but for the torus material thatmoves with a radial speed ’ .
01 (the light-red area inFigure 3, left) we expect that the driving is in addition also supported by gas pressure. Note also that along theaxis where the toroidal field vanishes we observe infalltowards the black hole.In Figure 4 we show the evolution of the average mag-netic energy for both the poloidal and the toroidal field.In both channels the energy increases over orders of mag-nitude until a saturation level appears around t = 1500.There is no full saturation as beyond t = 1500 the fieldenergy still increases by a factor of 100-1000, howeverthe strong increase of the dynamo action flattens sub-stantially. The flattening is interesting insofar as we donot apply dynamo quenching in this simulation. So thedynamo seems to self-quench. This has also been seenin the literature (see below). We hypothesize that this isdue to re-connection of the tangled magnetic field withtime, as Figures 2 and 3 show clear evidence for an in-creasingly turbulent state of the magnetic field. Due tothe application of resistivity, this field may physically re-connect, lowering the magnetic energy and heating theplasma.We think that the process works in two ways. First,the dynamo-generated magnetic field can transportedout of the dynamo-active area, however an area wherethere is still a non-negligible resistivity (basically thelight red area in Fig. 1 (upper right). Since the field istangled (see Figures 2) it will reconnect. However, thesame process may happen also inside the torus. Theresistivity is largest inside the torus (dark red areasin Fig. 1, right), thus, reconnection will be very effi-cient here. Considering this and also that the dynamo-generated field is strongest here (the ξ is large), recon-nection will work efficiently in lowering the net magneticfield energy.Interestingly, this points towards a new possible chan-nel for dynamo quenching by reconnection. Consideringmagnetic diffusivity, we may mostly think of diffusingaway magnetic flux and by this lowering the efficiencyof the dynamo. However, also reconnection is a result ofmagnetic diffusivity, that does naturally lowers the fieldenergy.3.2. The dynamo-generated torus magnetic field
The initial condition of the toroidal component of themagnetic field follows the density profile while keepinga positive sign in both hemispheres. A toroidal field isalso produced by the Ω effect by differential rotation ofthe induced poloidal field, increasing the initial toroidalfield.The induced toroidal field changes sign in the equato-rial plane, keeping its positive values in the upper hemi-sphere (where the radial field is negative). This results0
Vourellis & Fendt
Figure 3.
Distribution of the radial velocity at t = 2000 (left) and t = 3800 (middle) as well as plasma- β at t = 3800 (right)in the torus simulation. The color bars are chosen to highlight the 4 different velocity areas, such as the torus (white), thefunnel flow (dark red), the torus wind (light red), and the infall (dark blue), as well as the areas of high (green) and weak (red)magnetization. Figure 4.
Time evolution of the average magnetic energyacross our simulation grid for the poloidal and toroidal fieldcomponent. in a gradual change in the toroidal field across the torus,when the initial condition is superimposed by the newlygenerated field.At time t = 1000, the sign has changed already withinthe borders of the area where dynamo activity has pre-scribed. The change of sign of the toroidal field then con-tinues to gradually affect also the lower hemisphere. Attime t = 2000 this transformation has almost finished,resulting in a toroidal field distribution that changes signacross the equatorial plane.The dynamo-generated poloidal field is dominated bythe radial component. Note that the direction of theradial field component also affects the direction of thepoloidal field. Here we have the problem that thedynamo-generated field which is initially prescribed as aperfect dipole (with negative values in the upper hemi-sphere and positive values in the lower hemisphere), eventually turns into a form of stripes that show alter-nating positive and negative directions, which eventuallyresult in the formation of closed magnetic loops.This effect of inducing a layered structure of the ra-dial magnetic field components naturally also affects thedirection of the toroidal field, where we see a similar be-havior in later evolutionary stages. This effect becomesapparent from t = 500 − B r in the upper hemisphere, that wasinitially set to a negative radial field structure.Similarly, this is seen in all components of the mag-netic (and also in the electric field). From t = 1000 anadditional layer appears in the B r distribution, now witha negative value that is alternating with the previouslyinduced positive B r in the upper hemisphere.We observe this process for the whole duration of thesimulation run and eventually find a completely layeredmagnetic field structure of alternating field direction.We note that this layered structure shows up preferen-tially close to the equatorial plane of the torus. Here,we also find high dynamo numbers of the turbulent dy-namo exactly within previously mentioned layer. Inlater times this region of layered magnetic field directionspreads out over the whole torus and can eventually befound in the outflows as well. This dynamical processwas also reported in Bugli et al. (2014) and Tomei et al.(2020) despite a different prescription for the symmetryof the dynamo parameter.The time evolution and the spatial structure of themagnetic field evolving strongly reminds on dynamowaves that are generated during the kinematic phase of GRMHD disk dynamo Figure 5.
The radial magnetic field structure in simulation simT at time t = 2800. the mean field dynamo process. Dynamo waves are typ-ically discussed as oscillatory solutions of the kinematic dynamo problem. The existence of such dynamo waveshave been investigated early by Parker (1955) and manysubsequent studies (see e.g. Weiss et al. 1984; Hughes& Tobias 2010).Here, we solve for the dynamical dynamo problem, i.e.essentially considering the back reaction of the magneticfield generated on the gas. Nevertheless, in the earlystages of the simulations, we are still in the linear regimeof the dynamo action, and the field is not yet strongenough (high plasma- β ) to act on the gas. The situa-tion is thus comparable to the kinematic regime and theexcitation of dynamo waves may indeed be expected.More recent studies have detected also dynamo wavesin cases of high magnetic Reynolds numbers (Cattaneo& Tobias 2014; Nigro et al. 2017). This is essential toknow, as also in our simulations the magnetic Reynoldsnumbers R m ≡ ( vL ) /η ’ L ’
10, atypical velocity v = 0 .
1, and a maximum diffusivity of η = 0 .
001 (all in code units). We note here that much ofthe work cited here has been done for spherical stellardynamos. However, the case of of the relativistic torusis probably not too far from this, as compared to ourstudies of thin disks later on in this paper, the shear inthe torus is relatively weak.A characteristic feature of dynamo waves is that thetypical time scale of the fluctuations is similar to thediffusion time scale (see e.g. Giesecke et al. 2005 forspherical α -dynamos). For the simulations presentedin this section, this is indeed the case. With a maximum diffusivity ¯ η ’ − in the torus and the typical lengthscale in the torus ∆ R ’
1, we calculate a global diffu-sion time scale of τ diff = (∆ R ) / ¯ η ’ in code units.This is indeed similar to the time scale we measure forthe global variations of the dynamo-induced field struc-ture. Although the values considered are some roughestimates and vary drastically over the torus, we maysee this correlation as a hint for dynamo waves in oursimulations.Even at late stages of our simulation, the plasma-betain the torus is of the order of 100-1000, thus the dy-namo still acts in the kinematic regime. Nevertheless,the generated magnetic field has expanded and fills thespace between the rotational axis and the torus. Here,the plasma-beta is lower β ’ .
1, and the evolution isdominated by the magnetic field.Overall, we see indication that dynamo waves are ex-cited during the initial stages of the magnetic field evolu-tion on our torus model simulation. This is in particulardemonstrated in Figure 5 that shows the generation ofa layered structure of an anti-aligned field (here shownby the by the radial field component. This figure showstriking similarities to Figure 3 in Tomei et al. (2020),who, however, plot the magnitude of the field only, notits direction.3.3.
Comparison to literature studies
At this point it is interesting to compare to the exist-ing litrerature, in particular to the simulations of Buc-ciantini & Del Zanna (2013) on which the mean-fielddynamo closure of our code is based and the follow-uppapers (Bugli et al. 2014; Tomei et al. 2020) . As men-tioned above, while the geometry of the latter two ap-proaches look similar, the actual initial conditions arequite different. Our initial torus structure follows theFishbone-Moncrief solution (Fishbone & Moncrief 1976)that is applied in previous HARM simulations, whiletheir torus follows a different prescription. Overall, ourinitial torus extends to r ’
45, while Bugli et al. (2014)consider a quite smaller structure with r <
10. The
Pol-ish donut torus solution applied in Tomei et al. (2020)appears to be much larger, however of the full grid of r <
100 only the innermost radii r <
20 are shown.Also the initial field structure and the diffusivity dis-tributions are different.Given the different initial setup a quantitative com-parison to our results is not possible. However, we ob-serve similar structures evolving from our dynamo. Sim-ilar to Tomei et al. (2020) we observe the generation of Of course this value quite changes over the area of the torus Vourellis & Fendt bipolar magnetic "arms" around the inner edge of thetorus. Also, authors claim to detect dynamo waves be-ing dragged towards the black hole by accretion.We show the 2D magnetic field evolution at a some-what longer time step (in codes units t = 3500 vs. 1608in Tomei et al. 2020). Therefore we are able to followthe expansion of the generated magnetic flux away fromthe torus and towards the rotational axis. At t = 3500we indeed see a strongly magnetized axial region, thatmay allow for Blandford-Znajek jets in case of a rotatingblack hole.We note however that the dynamical evolution of thetorus and also of the dynamo action is quite faster intheir smaller sized torus setup. The 6 periods of torusorbits they show in their 2D images (and the 12 orbitsthey actually simulated) are calculated considering a ra-dius of the torus center of r c = 12, while the center ofour torus is further out r c = 15, resulting in about a fac-tor 2 longer rotational periods for our torus center andthus a similarly longer evolution time. The advectiontime of magnetic flux towards the axis is not affected bythat, however.Interestingly, Tomei et al. (2020) find a (second) satu-ration phase of field amplification even without quench-ing their mean-field dynamo, probably due to the ac-cretion dynamics together with an equipartition state ofthe fluid systemIn our simulations we also find a transition to a satura-tion state, more accurately, a strong break in the growthrate of both field components (see Figure 4 and our cor-responding discussion above). Tomei et al. (2020) arguethat the change is happening at the time when localequipartition is reached. We may confirm this as closeto the inner edges of the torus we also find equiparti-tion at this stage (see Figure 3 upper right for a latterevolutionary stage). However, we hypothesize that thechange in slope for the increase of magnetic energy canalso be due by reconnection of the tangled magnetic fieldthat is produced by the dynamo action (see discussionabove).We do not find the second saturation phase claimedby Tomei et al. (2020). The time scale of our simulationwould be sufficiently long, however, since our simulationsetup is different, in particular the dynamics of initialtorus, a proper comparison cannot be made here. Wenote again, that measured in orbital periods of the torus,the evolutionary time of their torus is somewhat longer,and we might therefore have not yet reached the secondsaturation phase (if this phase would be present at allin our setup). The poloidal field structure derived by Tomei et al.(2020) is visualized by showing the magnetic energy dis-tribution. Here, we also show poloidal field lines thatcan better trace the small-scale structure of the field, aswell as the field direction.To see the very structure of the vector field, meaningthe geometry of the magnetic field lines and the distribu-tion of the magnetic flux, is essential for understandingthe jet launching process . For example, for consideringmagneto-centrifugally driven disk winds, the field incli-nation is important (Blandford & Payne 1982).So far, for the test simulations provided in the liter-ature of GR-MHD dynamos (Bucciantini & Del Zanna2013; Bugli et al. 2014; Tomei et al. 2020), field lines forthe resulting magnetic field structure are only given forthe case of a neutron star dynamo (see Section 5.2.4 ofBucciantini & Del Zanna 2013). Thus we cannot com-pare our field structure in detail with the existing liter-ature. This would have been interesting as the essentialeffect of the dynamo action is in particular to generatepoloidal magnetic flux along the disk. INITIAL SETUP FOR A THIN DISK DYNAMOHere we describe the initial conditions for our dynamosimulations considering a thin accretion disk. The sim-ulation grid extends from just inside the event horizonto a outer radius of r = 80. The radial and polar di-mensions are based on the original grid of HARM
Thenumerical grid size used is 256 × ρ ( r, θ ) = (cid:20) Γ − r in r (cid:15) (cid:18) sin θ + (cid:15) ΓΓ − (cid:19)(cid:21) / (Γ − , (28)which is connected with the gas pressure by an idealequation of state P = K ρ γ . The aspect ratio of thedisk is ε ∼ .
1. The disk is surrounded by a corona withan initial density and pressure given by ρ cor ∝ r / (1 − Γ) , p cor = K cor ρ Γcor . (29)We choose K = 0 .
001 for the disk and for the corona K cor = 1 and force an initial pressure equilibrium be-tween he disk and the corona implying a density jumpbetween corona and disk surface. When the simulationsstarts, the initial corona starts collapsing, and is then Note that while the poloidal magnetic flux R ~B p d ~A may vanishon average, the magnetic energy ∝ B can still be substantial. GRMHD disk dynamo Figure 6.
Distribution of the magnetic diffusivity η , the dynamo parameter ξ , and the magnetic Reynolds number R ξ .Representations are in log scale. replaced by the floor values for density and pressure,described by a broken power (Vourellis et al. 2019).The disk is given an initial rotation profile followingPaczyńsky & Wiita (1980) as˜ u φ = r − / (cid:18) rr − R PW (cid:19) , (30)where R PW is a constant, here set equal to the gravita-tional radius R g .The initial seed magnetic field is purely radial andis prescribed only inside initial disk. Having a radialfield has the advantage of a vanishing shear betweenthe rotating disk and disk corona. We have also runsimulations starting with a purely toroidal initial field,however, resulting in a similar long-term magnetic fieldevolution. The initial magnetic field prescription in-volves purely radial magnetic field lines that convergein the black hole horizon, implemented numerically viathe magnetic vector potential, A φ ( r, θ ) ∝ exp " − (cid:18) .
05 tan θ (cid:19) . (31)The strength of the seed field is defined by the choice ofplasma- β = 10 .In the simulations we apply a profile of the mag-netic diffusivity that is somewhat modified compared toVourellis et al. (2019), but similar to the dynamo simula-tions by Stepanovs et al. (2014). This profile has a con-stant plateau across equatorial plane, but then quicklydrops in polar direction for θ < ◦ and θ > ◦ , η ( r, θ ) = η exp −
100 ( π/ − θ ) arctan( χ η · ε ) ! , (32) where η is the (initial) maximum value and χ η = 3characterizes the scale height of the diffusivity distri-bution. The profile of the mean-field dynamo followsStepanovs et al. (2014), ξ ( r, θ ) = ξ p r/r i n sin (cid:18) πχ ξ · ε · tan θ (cid:19) , (33)where χ ξ = 1 is a characteristic scale height for thedynamo. The latter is an important constraint, as bychoosing χ ξ < χ η we ensure that the dynamo action isalways enclosed in the resistive area. Otherwise we mayconsider arbitrarily large dynamo numbers (see below).Furthermore, the purely resistive outer layers of the diskwill allow for easier launching and mass loading of diskwinds (Vourellis et al. 2019). In Figure 6 we show thedistribution of the magnetic diffusivity η , the dynamoparameter ξ and Reynold’s number R ξ . THE ACCRETION DISK DYNAMOOur paper focuses on the generation of a strongpoloidal flux anchored to a thin accretion disk. As refer-ence simulation we consider the case of a Schwarzschildblack hole. A very low radial magnetic field with β = 10 is applied as a seed field for the dynamo action,following the (non-relativistic) approach by Stepanovset al. (2014); Fendt & Gaßmann (2018). Note that ourinitial disk structure is not in complete equilibrium withthe black hole (different from the non-relativistic sim-ulations). Thus, during the initial evolution the diskwill adjust on the relativistic potential, causing devia-tions on the seed field from the perfectly radial structure.These deviations, however, are minimal and significantly4 Vourellis & Fendt
Table 2.
Summary of thin disk simulation runs. Listed isthe ID of the run; the Kerr parameter a applied; the max-imum magnetic diffusivity η , typically located at the in-ner disk radius; the maximum (initial) dynamo parameter ξ ; and the quenching parameter β e q , corresponding to anequipartition field strength, beyond which dynamo action isstrongly quenched.run a η ξ β e q sim0 0 0.001 0.004 1000sim1 0 0.001 0.004 100sim2 0 0.001 0.004 10sim3 0 0.001 0.004 1sim0.1 0.9 0.001 0.004 1000sim1.1 0.9 0.001 0.004 100sim4 0 0.001 0 -sim5 0 0.001 0.004 no quenching smaller that the magnetic field that is generated by thedynamo later on. Simulations with a different initialfield structure (e.g. a toroidal field) lead to the samemagnetic structure in the saturation phase.Figure 7 shows exemplary evolution of the poloidalmagnetic field together with the evolution of the den-sity for reference simulation sim1 . When the simulationstarts, the initial radial structure of the magnetic fieldinside the disk is evolved by the dynamo action. Thegeometry of this dynamo-generated field does not followa clear dipolar structure.In addition to the dynamo action, the hydrodynamicturbulence that develops in the disk amplifies the stateof the magnetic field. We note here that this turbu-lence is a consequence of the existence of an anti-alignedpoloidal field structure that further develops when thepoloidal field is advected towards the horizon. The fieldgeometry close to the horizon (say for r <
5) resem-bles that of the Wald solution (Wald 1974; Komissarov2005). With our resistive approach, this leads to strongreconnection events that affect the mass loading of theoutflow as well as the accretion of material, and alsothe foot point of the magnetic field lines that guide theoutflow (see also Vourellis et al. 2019).We do not expect the MRI to be detected in our sim-ulations. In fact our resolution would be almost suffi-cient to resolve the some of the large wave modes. Sowhile the MRI might potentially be observed in an idealMHD study, for resistive MHD the excitation of the MRImodes is severely damped. Here, we refer to Fleminget al. (2000) and Longaretti & Lesur (2010) who inves-tigated in great detail the effect of resistivity on the evo-lution of the MRI. Qian et al. (2017), investigating the role of the MRI in a general relativistic torus in resistiveMHD came to the conclusion that the MRI is more andmore damped for increasing level of resistivity (althoughapplying a lower resolution).As soon as the dynamo-induced magnetic field isstrong enough to remove angular momentum from thedisk, accretion sets in. With accretion going on, thenewly generated field is advected towards the black hole.In the last stage of the simulation the environment closethe black hole is strongly magnetized by a poloidal fieldwith a dipolar configuration. However, at larger radiiand especially inside the disk, at the time scales investi-gated, the dynamo has not yet generated a clear struc-ture that provides a large scale magnetic flux. Note thatthese simulations are quite CPU-intensive, and the timescales we are reaching here of about t = 10000 corre-sponds to about 150 inner disk rotations. In contrast,the non-relativistic disk dynamo simulations (Stepanovset al. 2014; Fendt & Gaßmann 2018) reach up to several100.000 disk rotations.5.1. The evolution of the dynamo generated magneticfield structure
We now discuss the magnetic field structure that isgenerated by the disk dynamo (see Figure 7) and how itsgeometry and field strength depends on the parametersthat trigger the dynamo process, namely the strength ofthe dynamo ξ , the disk diffusivity η (in combination thedynamo number), and the quenching parameter µ eq .Figure 8 shows the evolution of the integrated poloidaland toroidal magnetic energy in three different regions(distinguished by their inner and outer radius). We em-phasize that for later times the dynamo induced poloidalfield remains substantially higher than the toroidal fieldcomponent. This again points towards an efficientlyworking α -dynamo, since only the α -process can amplifythe poloidal field component.For the innermost part of the accretion stream be-tween the black hole horizon and the initial inner diskradius (2 < r < t ’ t ’ GRMHD disk dynamo Figure 7.
Evolution of the disk density and the poloidal magnetic field (white lines) for simulation sim1 at times t =(0 , , , , , ,
0) point is the event horizon and the green line denotes the radiusof the ISCO. netic energy increase decreases even further to the pointit is almost constant. Still, we do not reach a constantvalue, as at t ’ < r <
30) but after t ’ t ’ Vourellis & Fendt
Figure 8.
Evolution of the poloidal (solid line) and toroidal(dotted line) magnetic energy integrated between three dif-ferent radii inside the accretion disk for simulation sim1 . Figure 9.
Evolution of mass flux rate measured at at r ∼ (3 , ,
11) (top) and normalized disk mass (bottom) forpolar angles 75 ◦ < θ < ◦ for simulation sim1 . The nega-tive values denote the accretion rate. The mass flux is alsonormalized with the initial disk mass. Figure 10 shows distribution of the plasma- β for thepoloidal and the toroidal magnetic field components attwo different time steps for simulation sim1 . Obviously,due to the dynamo activity, the plasma- β decreases with Figure 10.
Plasma- β at simulation times t = 4370 (top)and t = 12000 (bottom), considering the poloidal (left) andthe toroidal (right) magnetic field component, respectively,for simulation sim1 . time - first in inner disk area and later also in the maindisk body. Also the disk wind that is launched at latertime carries a low plasma beta. Overall, the toroidalmagnetic field component is dominating. This has beenobserved also in Vourellis et al. (2019) for a magneticfield that is initially prescribed and is consistent withthe literature of GRMHD disk simulations. For the timescales considered here, this tells us that the differentialrotation (Ω effect) still dominates the turbulent dynamo( α effect). At r ∼
15 along the equatorial plane thedynamo number characterizing differential rotation is R Ω ≈
59, while the maximum ξ dynamo number R ξ does not exceed 13 (see Eq. (21)).However, we see that the plasma- β measured in theinner part of the disk and in the accretion funnel doslightly increase during the later times of the simula-tion. This is a direct effect of the quenching mecha-nism we implemented. In addition, the magnetic fieldbecomes quite entangled in the later stages of the sim-ulation (see last panels of Figure 7). As a consequence, GRMHD disk dynamo Figure 11.
Time evolution of the absolute integrated mag-netic fluxes. Shown is the magnetic flux from the disk surface(top), thus integrated radially along the disk surface, andthe flux through a circular cross section over the polar angle0 ◦ < θ < ◦ at radius r ∼
11 and its corresponding areas inthe lower hemisphere (bottom). Simulations with differentlevels of dynamo quenching and Kerr parameter are consid-ered. there exist further effects that can lead in the dissipationof the magnetic field, in particular magnetic reconnec-tion. What can be clearly seen from the figure is thatexactly those areas of low plasma- β are also the areasfrom which strong outflows develop from the inner disk(see also Section 5.4 below).An essential condition for jet launching is a suffi-cient magnetic flux of the jet source, synonymous fora large-scale open magnetic field structure. As alreadybriefly mentioned, although we may measure a substan-tial (poloidal) magnetic energy ∝ B in the disk, thepoloidal magnetic flux R ~B p d ~A may vanish on average,if the field is strongly tangled.It is therefore interesting to see how the magnetic fluxevolves that is generated by our disk dynamo. This isshown in Figure 11 for simulations with different Kerrparameters and quenching thresholds. Here we choosetwo different representations of the absolute integratedmagnetic flux. One option is to integrate along a circleof constant radius (lower panel) and the other is to inte- grate along the disk surface (upper panel), using the B r and B θ components of the magnetic field, respectively.Note that radius of the circle of r = 11 is chosen suchthat it measures the magnetic flux close to the blackhole which potentially may launch a Blandford-Znajekjet. The angle for the integration is chosen so that weavoid to account for flux from inside the disk where thefield is more entangled.The increase in magnetic flux is similar to the timeevolution of the magnetic energy (see Figure 8). Differ-ences are due to the fact that when integrating the cu-mulative magnetic flux | R ~B p d ~A | , we effectively averageover fluxes in opposite directions, as discussed above.This clearly demonstrates that our disk dynamo alsogenerates a substantial magnetic flux when generatingmagnetic energy (see Figure 7 showing indeed large-scaleopen poloidal field lines). Interestingly, we find that alldynamo simulations evolve with an initially very similarbehaviour. At the later evolutionary stages they, how-ever, diverge, due to the different quenching thresholdthey obey. Simulation sim0 keeps an almost constantmagnetic flux during its final part, while for simulationswith a higher quenching threshold the flux increases witha steeper slope. Simulations sim2 and sim3 end earlierdue to the higher magnetization levels that is allowedby the quenching. Simulation sim0.1 that considers ahighly spinning black hole, also provides an ever increas-ing magnetic flux over the simulation time, however, thedynamo effect is somewhat delayed.For comparison, we have also performed test simula-tions sim4 without any dynamo action and sim5 withoutquenching. As expected, in simulation sim4 the initialradial field structure is conserved while slightly expand-ing in the surrounding disk corona. Minor weak out-flows can be detected from the disk surface that canbe attributed mostly to the local force balance betweenthe disk vertical pressure gradient (magnetic pressure isnegligible here) and the vertical gravitational force (in-duced by the black holes metric). In simulation sim5 themagnetic field develops as in simulation sim3 , however,the lack of a quenching mechanism generates a magneticfield that becomes extremely strong and that leads thecode to crash rather early. A similar behavior was alsoobserved by Tomei et al. (2020). Both simulations, sim4 and sim5 may serve as reference for further test simu-lations of our implementation of the dynamo and thedynamo quenching.5.2. Disk evolution - magnetic field and mass fluxes
We now investigate the evolution of the accretion disk,in particular the mass accretion rates, in more detail.Figure 9 (top) show the normalized mass flux through8
Vourellis & Fendt surfaces of constant radius r = (3 , , ◦ < θ < ◦ in order to ac-count for the initial disk structure and the evolution ofthe disk.Due to the weak seed field, the initial evolution is ba-sically hydrodynamic. While the initial radial structureof the magnetic field is ideal for angular momentum re-moval from the disk material, its strength is too low inorder to have a big impact.We find that right in the beginning of the simulation,the accretion rate shows large variations. This is dueto the lack of a perfect vertical hydrodynamic equilib-rium in the disk. This leads to distortions in the diskstructure, and to episodic accretion events of materialthat has moved inside the ISCO (note the initial posi-tion of the inner disk radius is at r = 10, well outsidethe ISCO). The inner part of the disk is quickly ad-vected and when it crosses the marginally stable orbitit disconnects from the disk and falls in the black holeleaving a “gap” behind it which quickly fills with ma-terial from the remaining inner part of the disk. Theinitial radial shape of the initial magnetic field supportsthe advection of disk material towards the black hole,however, as mentioned above, this field is not dynami-cally important. On the other hand, the motion parallelto the poloidal field also implies that no magnetic fluxis advected initially.The episodic accretion is repeated several times de-pending on the disk radius (see the accretion spikes inFigure 9, top). Further out along the disk, the accre-tion peaks correspond to disruptions that appear at thedisk surface displacing disk material (see Figure 7, topmiddle panel). This is a common consequence of thedynamics of the this disk in the relativistic environment(Vourellis et al. 2019).Around time t ’ t ∼ Figure 12.
The distribution of the dynamo parameter ξ and its variation due to quenching for simulation sim0.1 attimes t = 1000 , ary phase, however, showing a time variation in velocityand mass flux. Within the accretion disk the magneticfield is substantially entangled, forming loops that re-connect with each other. The disk starts expanding inthe vertical direction, but not yet developing any strongdisk wind.After t ∼ r >
10, which is located inside theaccretion disk. When the magnetic field is advected,the magnetic flux in the black hole magnetosphere in-creases along with the magnetization, but this area isneither dynamo-active nor affected by quenching. Thestrong field of the black hole magnetosphere results indivergences in the integration of the equations whicheventually lead to the termination of the simulation.When the disk field becomes advected, these field linesalso populate the area along the rotational axis resultingin a strong axial magnetic field that could potentiallylead into a Blandford-Znajek-driven outflow (in the caseof a rotating black hole, see Sect. 5.5).5.3.
Dynamo quenching
GRMHD disk dynamo sim1 has µ e q = 1 /β e q =0 .
01. This implies that as the plasma- β inside the diskdecreases from the initial value of 10 to an actual valueof 100, the quenching becomes stronger.Even for a plasma- β = 1000, in our parameter setupthe dynamo action is already quenched by ∼
11% and bythe time when plasma- β = 100 the dynamo parameterwill be half of its initial value. Note that the dynamoquenching works locally, defined by the local verticallyaveraged disk magnetization. Depending on this, theactual strength of the dynamo tensor is quenched.This is shown in Figure 12 where we show the distri-bution of the dynamo- ξ for two exemplary time stepsfor simulation run sim0.1 with a rapidly rotation blackhole. While at t = 1000 the dynamo still works atits initial strength, at t = 12 ,
000 quenching is clearlyseen at certain disk areas. The quenching is locally dif-ferent, implying that the dynamo works with differentstrength at different positions along the disk and, whileit smoothens the exponential amplification of the field,at the same time it restricts a stronger accretion rate.5.4.
Generating a disk wind
It is well known that a substantial disk magnetic fieldis essential for driving a strong disk wind (Blandford &Payne 1982; Casse & Keppens 2002; Ferreira 1997). Alsothe field inclination plays a role (Blandford & Payne1982). In contrast to simulations in the literature con-sidering jet launching from magnetized disk that ei-ther apply a pre-defined disk field (Zanni et al. 2007;Sheikhnezami et al. 2012; Stepanovs & Fendt 2016) orself-generate the disk field from a non-relativistic dy-namo (Stepanovs et al. 2014; Fendt & Gaßmann 2018),in our GRMHD dynamo simulations the resultingdiskmagnetic field does not have the smooth large-scalestructure that may support a Blandford-Payne diskwind due to the symmetry of the initial field.However, we can still detect the launching of strongoutflows. From our reference simulation sim1 (see Fig-ure 7) we can see the evolution of small scale outflowsfrom the disk surface which overall merge into a broaddisk wind. The initially weak winds are launched bya combination of the increased magnetic field pressure
Figure 13.
Radial velocity of the disk wind in the ref-erence simulation sim1 at times t = 8000 and 12000. Theblack dashed line follows the density contour ρ ≈ .
003 whichroughly indicates the size of disk and the transition to thecoronal region. gradient (mainly, but not only, from the toroidal com-ponent) and the re-arrangement of the accretion diskstructure in the gravitational field of the black hole (ad-vection of magnetic flux).In Figure 13 we show the distribution of the ra-dial velocity for the reference simulation at times t =8000 and 12000. We also display the density contour ρ ≈ .
003 that could be understood as a proxy for thedisk surface . A better choice would be transition frominflow to outflow, however, for most parts of the diskcorona close to the disk, there is large-sized area witha constantly positive outwards velocity. Instead we findpatches of both out-flowing and infalling material withvelocities somewhat below ( < . c ). This is the clearsignature of a turbulent outflow structure. We note thathere we are still dealing with the initial stages of outflowfrom the inner parts of the disk main body. Note thatthe magnetic field and the disk accretion is still evolv-ing, as time t = 10000 corresponds to only 50 inner diskorbits, much less than ones achieved in non-relativisticsimulations.An exception is a high-speed outflow that is rootedin the accretion funnel between the black hole and theISCO. Here the gas into which the magnetic field isfrozen-in orbits with high speed, resulting in radial out-flow velocities up to 0 . c .In the last evolutionary stages the picture changes,however. The disk wind is now dominated by out- As the disk evolves constantly, this choice is somewhat arbitrary.Still the location of the disk surface is essential to study of thedisk-outflow interrelation, in particular the mass loading of thewind. Vourellis & Fendt flowing patches moving with ’ . c . The infallingpatches have almost disappeared and we see a substan-tial disk wind.In the polar area above and below the black hole we seeconstant infall of material. Note this reference simula-tion is for a Schwarzschild black hole and no Blandford-Znajek-driving is possible (see below for rotating blackholes). This axial area stays almost free of magneticflux until the very final stages of the simulation. At thispoint in time, the accumulation of magnetic flux alsoleads to an increasing strength of the outflows from theinner part of the disk. These outflows are rooted in theaccretion disk and outside the accretion funnel as theyare seen in Figure 13.Concerning the disk dynamics, we observe a mixtureof positive and negative velocities, overall pointing to-wards a turbulent nature of the accretion disk. Strongaccretion is only detected inside the marginally stableorbit.Mow we want to quantify the fluxes carried by theoutflows. In Figure 14 we show the (normalized) ra-dial mass flux integrated along the polar angle at radius r ∼
32 for simulations sim0.1 and sim1.1 . These valuesexpress the rate at which material is flowing outwards(positive values) or flowing inwards (accretion, collapse)towards the black hole (negative values).We may compare three different regions along the po-lar angle. The first region is approximately 0 ◦ < θ < ◦ , which is the funnel region close to the rotationalaxis. The second region is within 25 ◦ < θ < ◦ whichrepresents the area where a high speed outflow may ap-pear. The third region is within 65 ◦ < θ < ◦ , which isthe area of a potential massive disk wind. The respec-tive areas of the lower hemisphere are also added to theintegration of the mass fluxes.The figure only shows the evolution of the mass fluxfor the late times of the simulation, this is the time whena disk wind is established. When comparing the numeri-cal values, we see that simulation sim0.1 shows strongerwinds with an average mass flux in the second and thirdregions of 5 . × − and 1 . × − . In contrast, in sim-ulation sim1.1 we find mass fluxes that are only about50% of that, namely 2 . × − , and 8 . × − , respec-tively.Note that the only difference in these simulations isthe threshold for the quenching (see Table 1). For sim1.1 the dynamo is quenched only at a higher mag-netization level, therefore producing a correspondinglystronger flux (compare also to Figure 11). We learn fromthis comparison that an accretion disk that carries astronger magnetic flux, will eventually generate an out- Figure 14.
Integrated normalized radial mass flux throughsurfaces of constant radius at r ∼
32 for simulations sim0.1 (top) and sim1.1 (bottom). The integration is done alongthree different circular sections, between 0 ◦ < θ < ◦ , 25 ◦ <θ < ◦ , and 65 ◦ < θ < ◦ including the symmetric sectionsin the lower hemisphere. Note the upper panels in the twogroups of panels has tick marks almost an order of magnitudebelow the others. flow of lower mass flux. We emphasize, however, thatsimulation sim1.1 terminates much earlier than simu-lation sim0.1 due to strong magnetic field close to theblack hole . We suspect that if simulation sim1.1 wouldnot have crashed, a correspondingly stronger wind wouldhave been formed.We close this section by noting that wind is strongenough to affect the disk mass evolution, thus changing This was possible due to the higher magnetization threshold wehave chosen for sim1.1 . GRMHD disk dynamo Figure 15.
Comparison of the angular distribution of nor-malized mass flux (red), Poynting flux per solid angle (green)and Lorentz factor (blue) for simulation sim0.1 at t = 7000at radius r = 8 and r = 44 for both hemispheres. Negativemass flux indicates accretion towards the black hole. TheBlandford-Znajek driven jet funnel is clearly distinguishedby the peaks in Lorentz factor and electromagnetic energyflux. The mass flux increases with radius demonstrating adisk wind that is dominated by disk material. the slope of the disk mass evolution towards a largermass loss, i.e. a more rapidly decreasing disk mass (seeFigure 9, bottom).5.5. Generating a Blandford-Znajek jet
Simulations with rotating black holes show the launch-ing of strong outflows from the black hole magneto-sphere, especially in later stages when the magnetic fieldhas engulfed the axial area. This is a clear signatureof Blandford-Znajek jets driven by the black hole ergo-sphere. In Figure 15 we show the angular profile of massflux, Poynting flux and Lorentz factor at radii r = 8 , t = 7000 for simulation sim0.1 .For the inner radius the mass flux has negative valuesaround the equatorial plane which change into positiveones for slightly higher and lower angles. Accretion to-wards the black hole is determined by the negative massflux, while the mass outflow in the radial direction is at-tributed to a low velocity wind from the inner disk andis expressed by the positive values. For angles closer tothe polar axis the mass flux decreases significantly. Inthis area the mass flux is determined by the floor den- sity chosen for the simulation, a feature that is commonto all GRMHD simulations in the literature. Note how-ever, that in our case only the axial outflow is affectedand not the disk wind that is loaded with disk mate-rial and that carries a gas density is substantially higherthan the floor values. The Poynting flux remains largelynegative for polar angles close to the equatorial plane.It is positive for angles closer to the axis, indicative ofa Poynting-dominated Blandford-Znajek jet. Similarly,this is accompanied by increasing Lorentz factor towardsthe axis, again implying the launching of a relativisticjet.The picture becomes even more clear for the respec-tive profiles derived for the larger radius. The Poynt-ing flux is predominately positive in all directions withthe higher values detected for the relativistic jet, simi-lar for the Lorentz factor. Note that this Poynting fluxis essentially generated by the disk dynamo and eitherlaunched directly from the disk surface, or, after ad-vection of magnetic flux along the accretion disk, beingfurther amplified by the black hole rotation.The mass flux, however, gives a somewhat differentpicture. It shows variations across the equatorial planewith quite large positive and negative deviations froman average value (about more than 4 times lower thanthe variations in the inner disk). At this point in time,at radius r = 44 the disk has completed almost 4 orbitsaround the black hole. This has to be compared to al-most 45 disk orbits at r = 8. The much slower evolutionof the outer parts of the disk, in combination with theentangled magnetic field results in a much more turbu-lent accretion pattern of the outer disk parts. Recon-nection also plays a role as the entangled magnetic fieldmay be anti-aligned (see also Vourellis et al. (2019) fora discussion of resistive effects for GRMHD jet launch-ing). Note however, that the evolutionary time of theoutflow is quite faster than the disk evolution, meaningthat the observed outflow dynamics provide an instanttracer of the launching conditions at the foot point ofthe outflow.In Figure 16 we show the evolution of integratedPoynting flux (integrated along circles of constant ra-dius, r = 32) for simulations sim0.1 and sim1.1 hostedby a rotating black hole with a = 0 .
9. We concentrate ontwo regions. One is the axial funnel area of 0 ◦ < θ < ◦ where the relativistic jet launched by the black hole mag-netosphere is propagating. The other is the disk windarea of 25 ◦ < θ < ◦ where massive outflows from thedisk are being launched. The values of the symmetricareas in the lower hemisphere are taken into account aswell. We focus on the later stages of the simulationswhen the Poynting flux increases substantially.2 Vourellis & Fendt
Figure 16.
Integrated Poynting flux (top) through surfaces of constant radius r ∼
32 for simulations sim0.1 (left) and sim1.1 (right), respectively. The integration is done along two segments, between 0 ◦ < θ < ◦ (funnel) and 25 ◦ < θ < ◦ (diskwind), and their corresponding areas in the lower hemisphere, respectively. The corresponding evolution of the absolute dynamoparameter ξ (bottom) , averaged in space for the whole disk and the inner area ( r < The phase of strong Poynting flux begins approxi-mately at the same time for both simulations. How-ever, the higher threshold in plasma- β for the dynamoquenching allows simulation sim0.1 to establish a highPoynting flux for a much longer time. In both simula-tions the Poynting flux measured in the jet funnel andin disk wind start with similar strength. Magnetic fluxis advected and when the area close the black hole andthe axial region is more and more encompassed by thedynamo-generated magnetic field, we see that the fluxwithin the funnel increases while the flux of the diskwind decreases (see the diverging lines in Fig. 16, toppanels). Our interpretation is that magnetic flux thatis first anchored in the disk, is advected to the funnel,thus increasing the funnel flux and reducing the diskflux. Especially for simulation sim1.1 this reconfigura-tion process is more obvious. In Figure 16 (bottom)we show the evolution of the averaged absolute value ofthe disk dynamo parameter ξ for the same two simula-tions. For simulation sim1.1 , for which the quenchingof the dynamo starts around β eq = 100, we notice theflat profile for the dynamo coincides with the absence ofPoynting flux until the sudden increase of the Poyntingflux. Later, when the flux from the disk wind decreases( t ’ SUMMARYIn this paper, we have extended our resistive GRMHDcode rHARM3D (Qian et al. 2017; Vourellis et al. 2019)implementing a mean- field dynamo in order to study thegeneration of a large-scale magnetic field and its effect onthe launching of outflows. We have focused on a setupconsidering a thin accretion disk around a black holesrunning a number of simulation with various thresholdsfor the dynamo quenching and also applying differentKerr parameters for the black hole rotation.
GRMHD disk dynamo ξ follows a sinusoidal profile in vertical direc-tion from the the equatorial plane and a 1 / √ r profilealong radius.In the following we summarize our results.(1) Our implementation of the mean-field dynamo isbased on the work of Bucciantini & Del Zanna (2013)and is an extension of the purely resistive version of ourcode (Vourellis et al. 2019). The dynamo parameter ξ (corresponding to the common dynamo-alpha in the lit-erature) is inserted into Ohm’s law as a source of themagnetic field, and via the Maxwell’s equations into anequation for the evolution of the electric field. We alsoapply a quenching mechanism to mitigate the exponen-tial increase of the magnetic field.(2) We first set up a torus initial setup for the simu-lations. When comparing our results to those of Bugliet al. (2014); Tomei et al. (2020), we find in general goodagreement in the sense that the dynamo produces simi-lar field structures within the inner torus. In addition tothe published literature we (i) continue the torus mag-netic field evolution to longer times scales, thereby be-ing able to follow the advection of magnetic flux towardsthe black hole rotational axis, and (ii) also provide themagnetic field lines of the field generated by the torusdynamo.(3) The magnetic field evolution of the torus shows asaturation without applying a classical dynamo quench-ing for the α -effect in the code as was reported by Tomeiet al. (2020). We hypothesize that in this particularcase the saturation of the dynamo process was estab-lished by reconnection of the dynamo-produced tangledmagnetic field. We propose this as another channel ofdynamo quenching that works self-consistently in a tur-bulent state of MHD, as long as physical resistivity isconsidered.(4) Overall, the dynamo works as expected from non-relativistic simulations with the magnetic field linesemerging from the disk interior and expanding into thedisk corona as the disk magnetic energy increases rapidlyalready at the beginning of the simulation. Dynamo ac-tion slows down mainly due to the quenching functionwe apply.The poloidal magnetic field dominates over thetoroidal component especially in the outer parts of thedisk where the difference can be an order of magnitude. The outer parts of the disk remain weakly magnetized(at the time steps investigated), but the disk magnetiza-tion still slowly increases till the final simulation times,with the toroidal field component again dominating.(5) We have also examined the evolution of the large-scale magnetic flux from the black hole-disk system forsimulations applying different quenching levels for thedynamo as well as for rotating and non-rotating blackholes, as the flux is an essential quantity for jet launch-ing. The time evolution of the magnetic flux basicallyfollows the evolution of magnetic energy and it is natu-rally affected by the quenching threshold. Even thoughthe first stage is almost identical to all simulations, theones with lower quenching threshold in plasma- β in-crease their magnetic field faster in the later stages tothe point that they also terminate faster. Black holerotation appears to induce a delay in the increase ofmagnetic flux and also leads to the higher levels of mag-netization.(6) The hydrodynamic evolution of the disk changeswith the growth of the magnetic field. Initially, the ra-dial structure of the seed field allows for episodic ac-cretion until the strong, dynamo generated field hasevolved. The existence of a dynamo-generated, strongvertical field component allows for the gradual appear-ance of disk winds, which in turn lead to a smoothingof the otherwise quite peaky accretion rate. At the finalstages of the simulations, the interplay between dynamoactivity and dynamo quenching, together with the diffu-sion of magnetic flux, allows to develop a stronger diskwind together with a relativistic axial jet.(7) We have identified small-scale disk outflows thatare launched during the earlier stages of the simula-tions from the inner part of the accretion disk and whichare interrelated with the presence of a strong magneticfield. These small-scale outflows maybe considered asdisk flares. At the later stages of the simulation strongdisk winds are ejected along the poloidal magnetic field,being mainly supported by the magnetic pressure of thetoroidal field component. We do not find indication forBlandford-Payne driven disk winds in our setup and forthe time scales considered, as the disk winds is kinemat-ically dominated (as also in the found in the simulationsby Vourellis et al. (2019).(8) In our simulations considering rotating blackholes, we observe an additional outflow structure thatis a highly relativistic jet, apparently driven by theBlandford-Znajek mechanism. This jet carries a largeelectromagnetic energy flux when it leaves the ergo-sphere into the axial regions above the poles of the blackhole (jet funnel). By measuring the Poynting flux both4 Vourellis & Fendt in the jet funnel and in the disk wind, we find that theflux in the jet funnel increases over time in expense ofthe magnetic flux in the disk wind. While the dynamo isquenched during later simulation times, and, thus, thedisk magnetic flux becomes saturated, advection of mag-netic flux along the disk towards the black hole leads tothe high levels of Poynting flux we observe in the funnel.In summary, we have applied, for the first time, themean-field dynamo theory in the GRMHD context ofthin accretion disks around black holes. We have in-vestigated the evolution of the dynamo-generated mag-netic field, the accretion disk and the outflows that arelaunched with the help of the evolved magnetic field. Inparticular we discuss the structure and evolution of thepoloidal magnetic flux.We find that the dynamo-generated magnetic fielddoes not follow a clear dipolar structure, in differenceto non-relativistic studies in the literature. Future workis needed to understand the mean-field dynamo action inGRMHD in more detail, in particular simulations run-ning for longer time scales. The lack of a well-aligned, large scale dipolar field does not allow to launch strongdisk winds via magneto-centrifugal driving. However,when the dynamo-generated magnetic field reaches acritical level, it is well capable to launch magnetic pres-sure driven disk winds and also highly relativistic jetsfrom the black hole ergosphere.ACKNOWLEDGMENTSThe authors are grateful to Scott Noble for the pos-sibility to use the original HARM3D code for furtherdevelopment (resistivity, dynamo). Preparatory workby Qian Qian was extremely helpful for completing thiswork. C.V. thanks the International Max Planck Re-search School for Astronomy and Cosmic Physics at theUniversity of Heidelberg (IMPRS-HD) for funding. C.V.also wants to thank Hans-Walter Rix for financial sup-port. We acknowledge a detailed report with many help-ful and interesting suggestions by an unknown referee.REFERENCES
Abramowicz, M. A., & Fragile, P. C. 2013, Living Reviewsin Relativity, 16, 1, doi: 10.12942/lrr-2013-1Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214,doi: 10.1086/170270—. 1998, Reviews of Modern Physics, 70, 1,doi: 10.1103/RevModPhys.70.1Bardou, A., von Rekowski, B., Dobler, W., Brand enburg,A., & Shukurov, A. 2001, Astronomy and Astrophysics,370, 635, doi: 10.1051/0004-6361:20010267Blackman, E. G., & Field, G. B. 1993, PhRvL, 71, 3481,doi: 10.1103/PhysRevLett.71.3481Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883,doi: 10.1093/mnras/199.4.883Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433,doi: 10.1093/mnras/179.3.433Brandenburg, A. 1996, ApJL, 465, L115,doi: 10.1086/310142Brandenburg, A., Nordlund, A., Stein, R. F., & Torkelsson,U. 1995, ApJ, 446, 741, doi: 10.1086/175831Bucciantini, N., & Del Zanna, L. 2011, Astronomy andAstrophysics, 528, A101,doi: 10.1051/0004-6361/201015945—. 2013, MNRAS, 428, 71, doi: 10.1093/mnras/sts005Bugli, M., Del Zanna, L., & Bucciantini, N. 2014, MNRAS,440, L41, doi: 10.1093/mnrasl/slu017Bugli, M., Guilet, J., Müller, E., et al. 2018, MNRAS, 475,108, doi: 10.1093/mnras/stx3158 Casse, F., & Keppens, R. 2002, ApJ, 581, 988,doi: 10.1086/344340Cattaneo, F., & Tobias, S. M. 2014, ApJ, 789, 70,doi: 10.1088/0004-637X/789/1/70Cerdá-Durán, P., Font, J. A., & Dimmelmeier, H. 2007,A&A, 474, 169, doi: 10.1051/0004-6361:20077432De Villiers, J.-P., & Hawley, J. F. 2003, ApJ, 589, 458,doi: 10.1086/373949Del Zanna, L., Zanotti, O., Bucciantini, N., & Londrillo, P.2007, A&A, 473, 11, doi: 10.1051/0004-6361:20077093Dyda, S., Lovelace, R. V. E., Ustyugova, G. V., Koldoba,A. V., & Wasserman, I. 2018, MNRAS, 477, 127,doi: 10.1093/mnras/sty614Fendt, C. 1997, A&A, 319, 1025Fendt, C., & Gaßmann, D. 2018, The AstrophysicalJournal, 855, 130, doi: 10.3847/1538-4357/aab14cFendt, C., & Greiner, J. 2001, A&A, 369, 308,doi: 10.1051/0004-6361:20010108Ferreira, J. 1997, A&A, 319, 340Fishbone, L. G., & Moncrief, V. 1976, ApJ, 207, 962,doi: 10.1086/154565Fleming, T. P., Stone, J. M., & Hawley, J. F. 2000, ApJ,530, 464, doi: 10.1086/308338Gammie, C. F., McKinney, J. C., & Tóth, G. 2003, ApJ,589, 444, doi: 10.1086/374594Ghosh, P. 2000, MNRAS, 315, 89,doi: 10.1046/j.1365-8711.2000.03410.x