High-order schemes for non-ideal 3+1 GRMHD: a study of the kinematic dynamo process in accretion tori

# High-order schemes for non-ideal $3+1$ GRMHD: a study of the kinematic dynamo process in accretion tori

## Abstract

We present the first astrophysical application of the ECHO code in its recent version supplemented by a generalized Ohm law, namely a kinematic study of dynamo effects in thick accretion disks. High-order implicit-explicit Runge-Kutta time-stepping routines are implemented and validated within General Relativistic MagnetoHydroDynamics (GRMHD). The scheme is applied to a differentially rotating torus orbiting a Kerr black hole, where the mean-field dynamo process leads to strong amplification of seed magnetic fields. We show that the interplay between the toroidal and poloidal components occurs qualitatively in the same fashion as in the Sun, butterfly diagrams are reproduced, and a typical time-scale for the field evolution is found, depending on the dynamo and resistivity numbers, which could explain periodicities as observed in several accreting systems.

\resetcounters

## 1 Introduction

Large-scale, ordered magnetic fields are believed to play a crucial role to activate the mechanisms responsible for the observed high-energy emission, from active galactic nuclei to gamma-ray bursts. However, it is still not clear which is the precise origin of such fields. While the process of collapse to the compact objects which are believed to power these sources is certainly able to amplify any existing frozen-in field, the question of its origin just shifts to the progenitors. Turbulent motions and small-scale instabilities, such as the Magneto-Rotational Instability (MRI) in accretion disks or the Tayler instability in proto-neutron stars, are indeed good candidates to enhance also the level of turbulent magnetic fields, but the resulting configurations will be highly tangled and thus not appropriate to generate the required large-scale fields.

A process naturally able to amplify ordered magnetic fields is that of dynamo. Here we shall consider the so-called mean-field dynamo processes (Moffatt, 1978), i.e. when small-scale turbulent (correlated) motions, invariably arising in astrophysical plasmas with very high fluid and magnetic Reynolds numbers, give rise to an effective ponderomotive force, along , in the Ohm law. This is capable in turn to increase any initial seed field when coupled to the Faraday evolution equation. The resistive-dynamo Ohm law can be written, for classical MHD, as (here , )

 E′≡E+\varv×B=ηJ+ξB, (1)

where is the coefficient of resistivity (due to collisions and to the mean-field effects), and is the mean-field dynamo coefficient (in the literature is most commonly employed) . The latter term is precisely responsible for the dynamo process, leading to the generation of exponentially growing modes in the kinematical phase, that is when the feedback of the increasing field on the other quantities can still be neglected. Differential rotation also plays a role in the dynamo process, and when both effects are at work we refer to an dynamo.

While in classical MHD the above relation simply defines the electric field, which is a derived quantity of and alone (since ), in the relativistic case the displacement current in the Maxwell equations cannot be neglected, and the evolution equation for the electric field must be solved as well. This difficulty was avoided in previous relativistic studies of kinetic dynamo in accretion disks (Khanna & Camenzind, 1996; Brandenburg, 1996). To our knowledge the first generalization to full General Relativity and its formulations for numerical relativity can be found in Bucciantini & Del Zanna (2013) (BDZ from now on), where also a systematic validation of the second-order numerical scheme, within the ECHO (Del Zanna et al., 2007) and X-ECHO (Bucciantini & Del Zanna, 2011) codes, is presented. Here we briefly summarize the formalism and we improve the numerical method employed to higher-order. Moreover, we present, as a first numerical application, a study of the kinematic dynamo action in thick accretion tori around Kerr black holes. More details and additional simulations will be presented in a forthcoming paper (Bugli et al., 2014).

## 2 Basic equations

The fully covariant formulation for a resistive plasma with dynamo action, first proposed in BDZ, is a relation for the quantities measured in the local frame comoving with the fluid 4-velocity , namely

 eμ=ηjμ+ξbμ, (2)

to be compared with Eq. (1), where , , and are the electric field, magnetic field, and current as measured by the observer moving with . When written in the formalism (Gourgoulhon, 2012), the Maxwell equations become

 γ−1/2∂t(γ1/2B)+∇×(αE+β×B)=0,(∇⋅B=0), (3)
 γ−1/2∂t(γ1/2E)+∇×(−αB+β×E)=−(αJ−qβ),(∇⋅E=q), (4)

where the electromagnetic fields and are those measured by the Eulerian observer, is the lapse function, the spatial shift vector, the 3-metric with determinant , is the charge density and the usual conduction current. In the language, Ohm’s law (2) translates into

 Γ[E+v×B−(E⋅v)v]=η(J−qv)+ξΓ[B−v×E−(B⋅v)v], (5)

here is the Eulerian velocity and the usual Lorentz factor, from which one can finally compute , so that the new equation for is

 γ−1/2∂t(γ1/2E)+∇×(−αB+β×E)+(αv−β)∇⋅E= −αΓη−1{[E+v×B−(E⋅v)v]−ξ[B−v×E−(B⋅v)v]}, (6)

reducing to for and to Eq. (1) in the non-relativistic case.

## 3 Numerical settings and convergence tests

In typical astrophysical situations the resistivity coefficient is very small and thus Eq. (6) is a stiff equation: terms can evolve on timescales , where the latter is the hyperbolic (fluid) timescale. We have then to solve, even for the simple kinematic case, a system of hyperbolic partial differential equations combined with stiff relaxation equations, that can be rewritten in the form

 X=γ1/2E:∂tX = QX(X,Y)+RX(X,Y), (7) Y=γ1/2B:∂tY = QY(X,Y), (8)

where are the non-stiff terms, whereas are the stiff terms requiring some form of implicit treatment in order to preserve numerical stability. To solve this system, we here extend the numerical methods in BDZ (up to second order in time and space) to higher orders. For spatial integration we use the full set of high-order shock-capturing methods implemented in the ECHO code, whereas for time integration we use the IMplicit-EXplicit (IMEX) Runge-Kutta methods developed by Pareschi & Russo (2005) and first introduced to resistive (special) relativistic MHD by Palenzuela et al. (2009).

In Fig. 1 we show a couple of 1D numerical tests in Minkowski spacetime, already described in BDZ. The first one is the resistive diffusion of a smooth current sheet, followed from to . Here , , and to reduce compressible effects (this test is made by coupling the Maxwell equations to the full set of GRMHD, we do not consider the kinematic case). In this limit the system remains non-relativistic and the analytical solution is

 By(x,t)=B0erf(x2√ηt), (9)

like for the classical diffusion equation. In the top row, left panel, we show the analytical solution at the initial and final times, and the numerical one at , in the case . In the right panel we show the errors on (comparing with a run at extremely high resolution) as a function of the grid points number . When we use the third-order SSP3(4,3,3) IMEX method, an overall second or third order space-time accuracy is correctly achieved, depending on the spatial reconstruction method employed. The best performing scheme is REC-MPE5 combined to DER-E6 (Del Zanna et al., 2007).

The second test concerns exponentially growing dynamo modes in a static, magnetized background. The analytical solution, for a given wave vector , is

 By(x,t)=B0exp(γt)cos(kx),γ=(2η)−1[√1+4ηk(ξ−ηk)−1], (10)

where is the dynamo growth factor. In Fig. 1 we show in the bottom row, left panel, the theoretical and measured amplification of for various combinations of (diamonds for , triangles for , squares for ), , and . The correspondence is always matched. The case of and should not have growing modes, due to the stabilizing effect of resistivity, but at late times () small-scale fluctuations due to truncation errors triggers the fastest possible growing modes anyway. In the right panel we show the errors on for various choices of the IMEX schemes (and for REC-MPE5, DER-E6). The nominal second or third order is achieved.

## 4 Kinematic dynamo in accretion tori

We choose as the background equilibrium the differentially rotating thick disk in Kerr metric and Boyer-Lindquist coordinates described in Del Zanna et al. (2007). We select the maximally rotating case with , resulting in an orbital period of (in code units) at the center. The domain is (logarithmically stretched) and , with a numerical resolution as low as due to the long-term runs needed to capture the full dynamics. The IMEX-RK is SSP3(4,3,3), thus third-order in time, and we use REC-MPE5 and DER-E6 for spatial reconstruction and derivation (thus fifth order in space for smooth profiles).

In the kinematic case the dynamics is basically determined by the two dynamo numbers, corresponding to the respective importance of the and effects with respect to diffusion. We take and , where is the radial extension of the thick disk (and vertical too), the difference in the angular velocity between the disk center and the inner radius. The choice of these two numbers determine the values of and at the disk center. The coefficients are then modulated within the disk as the square root of the mass density and proportional to it, respectively, while we take and in the low-density atmosphere (co-rotating to minimize shear flow numerical dissipation). The dynamo coefficient is further multiplied by because we want to obtain a combination of the dynamo effect which is symmetrical with respect to the equator.

The initial seed field () can be taken either as purely toroidal or purely poloidal, in the latter case it is derived from a potential proportional to the square of the local pressure. In Fig. 2 we show a simulation with (maximum values at the disk center), corresponding to and , with an initial toroidal field. This is shown as a function of time (in units of ) and we can clearly see the periodical formation and migration of structures from the center to higher latitudes, always confined within the disk. In Fig. 3 we show the maximum value of the toroidal and poloidal orthonormal components of the field in the disk and their ratio, as a function of time. We observe the exponential increase of both components and the saturation of the ratio around a value , while the oscillating behavior is simply due to a difference in phase between the dynamo modes for the two components.

The periodical generation and migration of structures resembles very much what happens in the Sun, where the is believed to be responsible for the 11-years sun-spots cycle. This is usually tracked in the so-called butterfly diagrams, with periodical field reversal and migration towards the solar equator. In Fig. 4 we report the position of the peaks along an appropriate new coordinate as a function of time. In the left panel the run previously described is shown, with migration from the equator to higher latitudes ( in the northern hemisphere). In the right panel the sign of has been reversed (and the field initialized to a purely poloidal one), to better reproduce the solar situation (we recall that in solar dynamo applications).

To conclude, in this paper we have shown the first implementation of high-order IMEX-RK methods for the resistive-dynamo non-ideal GRMHD version of the ECHO code. Preliminary runs of the kinematic dynamo effect in thick accretion tori around maximally rotating Kerr black holes have ben performed as an astrophysical test case. Similarly to the solar cycle, we have found periodical generation and migration of magnetic structures, with exponential growth of the maximum value of the field, expected to be quenched by non-linear feedback when the full set of GRMHD will be solved. The (half) cycle of the dynamo process is seen to be around orbital periods, but this value can change according to the parameters chosen. In any case the turbulence properties can set an additional timescale, through the mean-field dynamo effect, unrelated to the larger-scale quantities. This property may help to explain periodical modulations in the accreting process and in thus the source luminosity itself (e.g. Gilfanov, 2010).

### References

1. Brandenburg, A. 1996, \apjl, 465, L115
2. Bucciantini, N., & Del Zanna, L. 2011, \aap, 528, A101. 1010.3532
3. — 2013, \mnras, 428, 71. 1205.2951
4. Bugli, M., Del Zanna, L., & Bucciantini, N. 2014, MNRAS, in press
5. Del Zanna, L., Zanotti, O., Bucciantini, N., & Londrillo, P. 2007, \aap, 473, 11. 0704.3206
6. Gilfanov, M. 2010, in Lecture Notes in Physics, Berlin Springer Verlag, edited by T. Belloni, vol. 794 of Lecture Notes in Physics, Berlin Springer Verlag, 17. 0909.2567
7. Gourgoulhon, E. 2012, 3+1 Formalism in General Relativity, vol. 846 of Lecture Notes in Physics, Berlin Springer Verlag (Springer)
8. Khanna, R., & Camenzind, M. 1996, \aap, 307, 665
9. Moffatt, H. K. 1978, Magnetic field generation in electrically conducting fluids (Cambridge University Press)
10. Palenzuela, C., Lehner, L., Reula, O., & Rezzolla, L. 2009, \mnras, 394, 1727. 0810.1838
11. Pareschi, L., & Russo, G. 2005, Journal of Scientific Computing, 25, 129
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters