Dynamo action in thick disks around Kerr black holes: high-order resistive GRMHD simulations

We present the first kinematic study of an $\alpha\Omega$-dynamo in the General Relativistic Magneto-HydroDynamics (GRMHD) regime, applied to thick disks orbiting around Kerr black holes and using a fully covariant mean field dynamo closure for the Ohm law. We show that the $\alpha\Omega$-dynamo mechanism leads to a continuous exponential growth of the magnetic field within the disk and to the formation of dynamo waves drifting away or toward the equatorial plane. Since the evolution of the magnetic field occurs qualitatively in the same fashion as in the Sun, we present also butterfly diagrams that characterize our models and show the establishment of an additional timescale, which depends on the microscopic properties of the turbulent motions, possibly providing an alternative explanation to periodicities observed in many high-energy astrophysical sources where accretion onto a rotating black hole is believed to operate.


INTRODUCTION
Ordered, large-scale magnetic fields are believed to be a fundamental ingredient of the accretion processes that power many of the astrophysical sources of high-energy emission, like jets from Active Galactic Nuclei (McKinney & Blandford 2009) or Gamma Ray Bursts (Bucciantini et al. 2009;Rezzolla et al. 2011). This is particularly true for the Blandford-Znajek mechanism (Blandford & Znajek 1977) where the rotational energy is extracted from a rotating black hole through large-scale magnetic fields penetrating the ergosphere and twisted by the rotation of the surrounding spacetime. Unfortunately it is still not clear how such magnetic fields originate. The process of collapse to the compact objects which are believed to be the engines that power these sources could amplify any preexisting frozen-in field. However, recent MHD simulations of collapsing dense cores (Hennebelle & Fromang 2008;Santos-Lima, de Gouveia Dal Pino & Lazarian 2012) have shown how the magnetic braking due to a preexisting large-scale magnetic field can extract angular momentum from the cloud fast enough to prevent the formation of a disk. Turbulent reconnection can efficiently inhibit the magnetic braking by removing magnetic flux, and therefore a turbulent magnetic field allows the formation of a disk. However, turbulent motions and small-scale instabilities in accretion disks, such as the Magneto-Rotational Instability (MRI) (Balbus 2003), could enhance the level of turbulent magnetic fields, but the tangled final configuration would not be able to drive large scale outflows. matteo@mpa-garching.mpg.de A solution for the origin of the necessary large-scale magnetic field is a dynamo process. If one considers the presence of small-scale correlated fluctuations in the flow due to turbulence (that invariably arises in astrophysical plasmas with very high fluid and magnetic Reynolds numbers) a mean-field dynamo mechanism (Moffatt 1978) could lead to a large-scale effective electromotive force capable of generate and sustain a large-scale magnetic field. While this process has been vastly studied in classical MHD (mostly the solar dynamo), there are only a few applications in general relativity. In particular, Khanna & Camenzind (1996) have tested the role of a differentially rotating absolute space in exciting dynamo modes in accretion disks in Kerr metric in the simplified regime where the displacement current is neglected and for a fixed velocity field. This approximation was shown to be reasonable by Brandenburg (1996), but it cannot be achieved consistently in a covariant formalism, if causality is to be preserved. It might also not be appropriate when fast and variables motions, as expected in accretion, near the event horizon with velocities close to the speed of light are allowed to develop. The first formulation of a fully covariant closure of the GRMHD equations for a non-ideal plasma with a mean-field dynamo mechanism has been presented by Bucciantini & Del Zanna (2013), but so far it has never been applied to an actual astrophysical system. In this paper we present the first example (to our knowledge) of a resistive GRMHD simulation of a thick disk orbiting around a Kerr black hole with a mean-field dynamo mechanism. Given a small magnetic seed, we follow its evolution through the kinematic phase, assuming therefore that its feedback on the dynamics can still be neglected. We use the ECHO code (Del Zanna et al. 2007 (Pareschi & Russo 2005), to integrate the Maxwell equations in the 3+1 formalism (Gourgoulhon 2012) using the fully covariant mean-field dynamo closure for the Ohm's law given by Bucciantini & Del Zanna (2013).
The plan of this paper is as follow. In Section 2 we briefly describe the mean-field dynamo mechanism and its implementation in the ECHO code. Then in Section 3 we illustrate the setup of our simulations, and our choice for the dynamo properties. We finally show the results of our study and present our conclusions in section 4. In the following we set c = G = 1 and absorb all the factors √ 4π in the definition of the electromagnetic field.

DYNAMO MODEL AND EQUATIONS
For classical MHD, the introduction of a mean-field mechanism leads to the following form of the resistive-dynamo Ohm law: where E, B and J are respectively the electric field, magnetic field and current measured by an eulerian observer, E is the electric field measured in the frame comoving with the fluid-velocity v, η is the magnetic resistivity due mostly to mean-field effects (and in smaller part to collisions) and ξ is the mean-field dynamo coefficient (most commonly employed in the literature as α dyn ≡ −ξ). The latter term is a direct result of the assumption of correlated turbulent motions in the plasma, and provide the means for the generation of currents J ≡ ξ/ηB parallel to the magnetic field B, in contrast with the ideal case where J = (∇ × B) ⊥ B. Following Bucciantini & Del Zanna (2013), we adopt a fully covariant form of Ohm's law for a resistive plasma with a mean-field α-dynamo effect given by: where e µ , b µ and j µ are respectively the electric field, magnetic field and current measured by an observer comoving with the fluid 4-velocity u µ . The spatial projection of Eq.
(2) is then given by: and provides a mean to express the current J in terms of E, B and v. Here Γ ≡ (1 − v 2 ) −1/2 is the Lorentz factor and q = ∇ · E is the local charge density. Since only the electric and magnetic fields evolve through time, we just need to integrate the Maxwell equations, which in the 3+1 formalism are: with γ the determinant of the spatial metric tensor γij, α the lapse function and β the shift vector. Computing J from Eq. (3) and substituting it in Eq. (5) we obtain the final form of the equation for the evolution of the electric field: The above equation applies, rather than its ideal limit, when η is not completely negligible, as it occurs when in addition to the usual ohmic diffusivity we take into account the turbulent contribution. In this case, we do expect the terms ∝ η −1 to influence heavily the evolution of the electric field, and since they may evolve on a time-scale τη much shorter than the MHD time-scale some sort of implicit time integrator must be employed to guarantee stability to the integration of Eq. (6).

The background model
The assumption of kinematic dynamo implies that the background flow structure can be assumed as given. Our disk model corresponds to an unmagnetized thick torus with constant specific angular momentum L, in dynamical equilibrium, orbiting around a Kerr black hole (in Boyer-Lindquist coordinates) with specific angular momentum a = 0.99MBH, as described in Font & Daigne (2002). Distances and time are expressed in units of MBH.
The center of the torus is located at a radius rc = 5, while its inner edge is at rin = 3. Here the temperature is pc/ρc 7×10 −3 , where p and ρ are density and pressure. With these choices we determine ρ, v and p at every point of the disk. The disk is assumed to be surrounded by a low density (ρatm 10 −5 ρc ) hot (patm/ρatm 0.5) atmosphere in hydrostatic equilibrium and corotating with the disk itself (so that the atmosphere is also in keplerian-like rotation). As we will show, the choice of the initial seed magnetic field is unimportant (apart from the parity with respect to the equator) for the evolution of the dynamo, given that the fastest growing mode will always be selected. Model can be initialized either with a purely toroidal seed field (BT ) or a purely poloidal seed field (BP ). The initial field is always confined within the disk. The initial electric field is then set equal to the ideal value Given that in this letter we aim at presenting a first astrophysical application of a mean-field dynamo in GR to a situation of interest, our choice for the dynamo properties is mostly dictated by simplicity. We will therefore use the simplest possible choice for the parameter η and ξ, that still preserve some of their expected global properties. We are aware that a sub-scale model for the turbulence will be required to build a proper realistic model, however, we hope to show here the feasibility of this approach, and derive a rough understanding of its working.
Because of the kinematic approximation we adopted, we expect the behavior of the αΩ-dynamo to be affected by the values of η and ξ alone, since the other free parameters of the model affect only fixed quantities. Here we investigate the response of the system to a variation of this two parameters.
Let's consider the magnetic resistivity first. We assume the resistivity to be confined in the disk, where turbulent motions can act to dissipate the mean field. The resistivity in the disk scales as: where the central density ρc is also the maximum density, ρ * atm is the minimum one in the atmosphere, and η disk is the normalization value of the resistivity. The atmosphere is assumed to be an ideal conductor. The other option, corresponding to the opposed regime, is to assume it to be highly resistive, but we will not consider here this latter case. Since it is not possible to set η = 0 using the IMEX schemes, we have set ηatm = 10 −5 (which is also a lower bound to the resistivity in the entire domain). We have verified that at the resolution of our runs, such value gives results indistinguishable from the ideal case (this is the value of our numerical resistivity).
In the same way we assume mean-field dynamo to act only inside the disk. The fact that the dynamo parameter ξ never appears at the denominator in the resistive GRMHD equations allows us to set its value in the atmosphere to zero. Being ξ proportional to the plasma kinetic helicity, it has therefore to be odd in the axial coordinate z = r cos θ, so that ξ(r, θ) = −ξ(r, π − θ). In analogy with the resistivity, the simplest form we can choose is: in the atmosphere.
The profiles of resistivity and mean-field dynamo coefficient are shown in Fig. 1. A mean-field dynamo process taking action in an accretion disk leads to a αΩ-dynamo, where the toroidal magnetic field is enhanced by the differential rotation (the Ω-effect) and the poloidal component increases through the mean field dynamo (known in literature as the α-effect). In the kinematic case, the dynamics of such a system is characterized by the following dynamo numbers: which estimate the efficiency of the dynamo processes (due respectively to the mean-field dynamo and the differential rotation) in enhancing the magnetic field amplitude against the dissipation due to a finite resistivity of the plasma. In this relations R represents a typical high-scale of the disk (for thick disks we used the radius of the center), while ∆Ω is a typical shear (we have chosen the difference in angular velocity between the center and the inner edge of the disk).

Numerical methods and settings
Our domain extends in the range r = [r+ + 1.5, 25], θ = [π/4 − 0.2, 3π/4+0.2], where r+ is the radius of the event horizon. All of our simulations were performed on a numerical grid of 256 × 256 points. The grid is logarithmically stretched in the radial direction in order to have a larger resolution near the black hole (222 grid points within r = 10). The domain has been shaped to fully contain the disk. We have verified that the magnetic field never migrates beyond the polar regions that we have excised. We have also verified that at this resolution the results are converged.
Developing what has been done by Bucciantini & Del Zanna (2013), we adopted an IMEX Runge-Kutta scheme (Pareschi & Russo 2005;Palenzuela et al. 2009) to perform the time integration of Eq. (4) and Eq. (5). In particular, extending previous results, we have implemented the SSP3(4,3,3) scheme, that guarantees third order accuracy in time. This has been coupled, for the hyperbolic part, to a HLL solver with a 5 th order spatial reconstruction routine (MP5), which allow us to reach a global 3 rd order accuracy and to get convergent results already at modest resolutions (for more comprehensive numerical tests see Del Zanna, Bugli & Bucciantini 2014).

RESULTS AND DISCUSSION
In Table 1 we show the various runs of our study. They differ by the value of η, ξ and the initial profile of the magnetic field. In all of the models where the initial magnetic field is toroidal (symmetric to the equator), the α-effect generates from the beginning a poloidal antisymmetric component: therefore the magnetic field has quadrupolar symmetry. On the other hand, starting with a poloidal field (antisymmetric to the equator) a toroidal antisymmetric com- ponent arises due to the Ω-effect, thus the field has dipolar symmetry.
In either cases the magnetic field evolves to an eigenstate of the system, that is reached after a relaxation-time that ranges, depending on the particular model, roughly between 3 Pc and 30 Pc, with Pc 76.2GMBH /c 3 the orbital period of the disk center. As shown in Fig. 3, the eigenstate is characterized by a dynamo wave that propagates away from the equatorial plane for ξ > 0. For the models with ξ < 0 the propagation occurs towards the equatorial plane, qualitatively in agreement with the solar dynamo case and consistent with our definition of ξ ≡ −α dyn (Moffatt 1978). The magnetic field appears to originate approximately where the dynamo parameter is stronger, and then drifts from that location. As the field migrates toward the atmosphere, or toward the equator (depending on the sign of ξ), resistive effects become proportionally stronger and the field is dissipated.
Besides this drifting, the amplitude on the poloidal and toroidal components of the magnetic field grow exponentially with the same growth rate. This because, in the kinematic regime we have assumed, there is no quenching effect. The growth rate increases for larger values of ξ. The ratio between the maximum of the poloidal component and the maximum of the toroidal one ranges from 0.04 to 0.36 and shows an oscillating behavior indicative of a phase difference between these two components.
In order to quantitatively characterize the migration of the magnetic field within the disk, we produced for each model a butterfly diagram considering the value of the toroidal field along the trajectories of its maxima at different times. From the diagrams for Models 1, 2 and 7, shown in Fig. 4, is now evident the periodicity of the eigenstate and its spatial scale. For larger values of ξ we retrieve smaller periods and smaller spatial scales, while decreasing the action of the mean field dynamo both periods and spatial scales increase.
It is interesting to note that the particular symmetry of the initial magnetic fields does not affect the dynamical characteristics of the eigenstate selected by the problem: comparing for example the results for Models 1 and 7 (which differ only by the initial field) we can see how the growth rate, the period, the spatial scale and the ratio between poloidal and toroidal component are roughly the same. Moreover, observing the butterfly diagrams is evident how the initial symmetry of the magnetic field at the passage at the equator (quadrupolar for the first six models, dipolar for the rest of them) is preserved during all the evolution of the system. This means that the state reached by the system is degenerate in the parity at the equator of the magnetic field.
The turbulence proprieties, that determine the characteristics of the mean-field dynamo, can in principle set a time-scale (the dynamo wave period) which might be unrelated to the large-scale dynamics. Of course the same turbulence that sets the dynamo coefficient ξ, as well as the mean field resistivity η, should be related with other properties of the disk like the typical viscosity, and might be influenced by the cooling properties and the ionization state of the plasma. Understanding these possible connections, and their observable consequences, is however beyond the scope of this paper. We hope to extend this initial introductory work to a more realistic regime, where a more meaningful choice for η and ξ, based on a turbulent model, could allow us to properly evaluate the importance of mean-field effects (Brandenburg & Sokoloff 2002). This could help our understanding of some of the unexplained time-variability observed in accreting systems (Gilfanov 2010), that cannot be naively associated with orbital periodicities.
The relation between mean-field dynamo and MRI in thick disks is also a fundamental aspect to be investigated. So far, this has been studied numerically only in the classical MHD limit in the shearing box local approximation (for a review Blaes 2013). When vertical stratification is included, magnetic buoyancy couples with the turbulent dynamo effects (Brandenburg et al. 1995) and quasiperiodical field amplification, migration, and reversal can be observed (Davis, Stone & Pessah 2010;Flock et al. 2012). It would be interesting to reproduce these effects of self-generated turbulence in the GRMHD case as well.
Further developments might also include: the quenching of the α effect; the dynamical feedback of the magnetic field on the evolution of the hydrodynamical quantities, requiring the integration   Fig. 3) as a function of time. The magnetic field has been normalized to its maximum value in the whole domain at every time, which is measured in units of the orbital period of the center of the disk Pc. The parameter s ≡ z/ sin χ 0 represents the position on the trajectories. of the full set of GRMHD equations. Moreover, in this study we have fixed not only the disk structure but also the mass and the spin of the black hole. The role of these parameters needs also to be investigated, to understand if there are possible observable quantities that could help us to constrain them. This could be of potential interest especially in view of the large improvements in precision of the measurements of the mass and the spin of a black hole of the last few years (Risaliti et al. 2013).