text:
# A Numerical Method for General Relativistic Magnetohydrodynamics

## 1 Introduction

Accretion into black holes is believed to account for a wide variety of astrophysical phenomena, from solar-mass black holes in X-ray binaries, to supermassive black holes in active galactic nuclei. The basic theory of black hole accretion was laid out more than thirty years ago in a number of now-classic papers (e.g. Novikov & Thorne 1972; Shakura & Sunyaev 1973; Lynden-Bell and Pringle 1974). Since then, theoretical progress has been steady, but perhaps at a slower rate than might have been hoped. The problem is complex, involving time-dependent, three-dimensional dynamics of magnetized plasmas in the relativistic potential of Kerr black holes. Solutions are impossible to obtain analytically beyond simplified cases that rely upon time-stationarity and spatial symmetry. Because of this, numerical experiments must play an increasingly important role in driving theory. The importance of numerical simulations in the investigation of such complex physics has long been recognized, but until recently the necessary computational hardware was not available.

We are now in a position to pursue full 3D simulations with general relativistic MHD codes. Again, the basics were laid out thirty years ago by Wilson (1972) who pioneered both hydrodynamic and magnetohydrodynamic (MHD) two-dimensional simulations of black hole accretion (Wilson 1975; 1977; 1978). Since then there have been a significant number of efforts by a number of groups to develop accretion simulations for both Newtonian and relativistic gravitational fields, and for both hydrodynamics and magnetohydrodynamics. This paper details our current effort to develop a fully general relativistic MHD accretion code.

Why focus on general relativistic MHD? In fact, for many applications Newtonian gravity (or a pseudo-Newtonian mock-up of a black hole) is sufficient. In some important cases, however, it is not. The detailed physics of gas interacting with a black hole provides, in principle, a rare test of strong-field general relativity, but only if we understand the unique properties of accretion in the fully relativistic case. Distinctive and astrophysically interesting effects are expected from the Kerr metric, including those due to the lack of a stellar surface, the presence of an innermost stable circular orbit, and the dragging of inertial frames. The need to include MHD in an accretion simulation is also now quite clear. Magnetic fields play an essential role in the outward transport of angular momentum in accretion disks through the action of the magnetorotational instability (MRI; Balbus & Hawley 1998). This, in turn, means the simulations must be three-dimensional, for only in three dimensions can there be a self-sustaining (and accretion-sustaining) MHD dynamo.

In this paper we describe an algorithmic approach to the three-dimensional equations of general relativistic MHD that is based upon techniques first developed for axisymmetric hydrodynamics around black holes (Hawley, Smarr & Wilson 1984a, 1984b; hereafter HSWa and HSWb). The hydrodynamic portion of the code was recently tested and employed in three-dimensional simulations in the Kerr metric by De Villiers & Hawley (2002). Here we cast the equations of relativistic MHD in a form that is compatible with the HSW formulation. We will describe how those equations are evolved numerically. We will also discuss a suite of tests suitable for relativistic codes, and show how our code performs against those tests.

It is important to recognize that the development of numerical algorithms in any area of physics with few analytic solutions, and even fewer experimental observations, is likely to be difficult. Progress will be made only by exploring a variety of approaches. Such a situation necessarily precludes the development of only one numerical algorithm, or even a single code implementation of a promising algorithm. Some approaches will have strengths for certain applications, but it is unlikely that one numerical technique will be optimal for all problems. In any case, experience gained, both with what works and what does not, will be essential as we work toward developing a consensus within the community. With this in mind, our own efforts have concentrated on extending the established general relativistic hydrodynamic approach of HSW, since our main area of interest, that of accretion flows, has been well handled (in the hydrodynamic case) by this solver, and by similar algorithms designed to perform Newtonian MHD simulations. Our test suite has been developed in conjunction with Gammie, McKinney, & Tóth (2002), who use an alternative scheme for general relativistic MHD.

The plan of this paper is as follows. In §2 we write the equations of general relativistic MHD in the form that they will be solved in the code. In §3 we describe the numerical implementation of our algorithm. In §4 we describe the one-dimensional Minkowski (special relativistic) test problems. In §5 we describe the test problems in the Kerr and Schwarzschild metrics, and present our test results. The summary and conclusions are presented in §6.

## 2 Equations of General Relativistic Magneto-Hydrodynamics

We wish to study the evolution of a magnetized fluid in the background spacetime metric of a Kerr black hole. We work in Boyer-Lindquist coordinates, for which the line element has the form,

$$ds^2=g_{tt}dt^2+2g_{t\varphi }dtd\varphi +g_{rr}dr^2+g_{\theta \theta }d\theta ^2+g_{\varphi \varphi }d\varphi ^2.$$
(1)
In keeping with Misner, Thorne, & Wheeler (1973), we use the metric signature $`(,+,+,+)`$, along with geometrodynamic units where $`G=c=1`$; the black hole mass is unity, $`M=1`$. The determinant of the 4-metric is $`g`$, and $`\sqrt{g}=\alpha \sqrt{\gamma }`$ where $`\alpha `$ is the lapse function, $`\alpha =1/\sqrt{g^{tt}}`$, and $`\gamma `$ is the determinant of the spatial $`3`$-metric.

The state of the relativistic test fluid at each point in the spacetime is described by its density, $`\rho `$, internal energy, $`ϵ`$, $`4`$-velocity $`U^\mu `$, and isotropic pressure, $`P`$, which is related to the first two scalars via the equation of state of an ideal gas, $`P=\rho ϵ(\gamma 1)`$, where $`\gamma `$ is the adiabatic exponent.

The equations of general relativistic MHD are the law of baryon conservation,

$$_\mu (\rho U^\mu )=0,$$
(2)
where $`_\mu `$ is the covariant derivative, the conservation of stress-energy,

$$_\mu T^{\mu \nu }=0,$$
(3)
where $`T^{\mu \nu }`$ is the energy-momentum tensor for the fluid, and Maxwell’s equations,

$`_\mu F^{\mu \nu }`$ $`=`$ $`4\pi J^\nu ,`$ (4)
$`_\mu {}_{}{}^{}F_{}^{\mu \nu }`$ $`=`$ $`0,`$ (5)
where $`F^{\mu \nu }`$ is the electromagnetic field strength tensor, $`J^\mu =(\rho _c,J^i)`$ is the current $`4`$-vector, and $`\rho _c`$ the charge density. The dual tensor is defined as

$`{}_{}{}^{}F_{}^{\mu \nu }`$ $`=`$ $`{\displaystyle \frac{1}{2}}ϵ^{\mu \nu \delta \gamma }F_{\delta \gamma },`$ (6)
where $`ϵ^{\mu \nu \delta \gamma }=(1/\sqrt{g})[\mu \nu \delta \gamma ]`$ is the contra-variant form of the Levi-Civita tensor. Maxwell’s equations are supplemented by the equation of charge conservation $`_\mu J^\mu =0`$.

The energy-momentum tensor consists of perfect fluid and electromagnetic parts, $`T^{\mu \nu }T_{(fluid)}^{\mu \nu }+T_{(EM)}^{\mu \nu }`$, and has the form

$$T^{\mu \nu }=\rho hU^\mu U^\nu +Pg^{\mu \nu }+\frac{1}{4\pi }\left(F_{}^{\mu }{}_{\alpha }{}^{}F^{\nu \alpha }\frac{1}{4}F_{\alpha \beta }F^{\alpha \beta }g^{\mu \nu }\right).$$
(7)
Following Lichnerowicz (1967) and Anile (1989), we define the magnetic induction and the electric field in the rest frame of the fluid,

$`B^\mu `$ $`=`$ $`{}_{}{}^{}F_{}^{\mu \nu }U_\nu ,`$ (8)
$`E^\mu `$ $`=`$ $`F^{\mu \nu }U_\nu .`$ (9)
We adopt the ideal MHD limit and assume infinite conductivity (the flux-freezing condition), wherein the electric field in the fluid rest frame is zero, i.e., $`F^{\mu \nu }U_\nu =0`$.

We combine the definition of the magnetic induction (8) with (6) and the condition for infinite conductivity to obtain

$$F_{\mu \nu }=ϵ_{\alpha \beta \mu \nu }B^\alpha U^\beta ,$$
(10)
where $`ϵ_{\mu \nu \delta \gamma }=\sqrt{g}[\mu \nu \delta \gamma ]`$. The orthogonality condition

$$B^\mu U_\mu =0$$
(11)
follows directly from (8) and the anti-symmetry of $`F^{\mu \nu }`$.

Using these results, it is possible to rewrite the electro-magnetic portion of the energy-momentum tensor as

$$T_{(EM)}^{\mu \nu }=\left(\frac{1}{2}g^{\mu \nu }b^2+U^\mu U^\nu b^2b^\mu b^\nu \right),$$
(12)
where $`b^\mu B^\mu /\sqrt{4\pi }`$, and $`b^2=g^{\mu \nu }b_\mu b_\nu `$.

The induction equation (5) is solved using the equivalent form

$$_\delta F_{\alpha \beta }+_\alpha F_{\beta \delta }+_\beta F_{\delta \alpha }=0.$$
(13)
Working directly with $`F_{\mu \nu }`$ replaces covariant derivatives with simple coordinate derivatives. This is the basis for the Constrained Transport (CT) formulation of Evans & Hawley (1988; hereafter EH88) for the induction equation. In CT we make the identification

$$^r=F_{\varphi \theta },^\theta =F_{r\varphi },^\varphi =F_{\theta r},$$
(14)
where $`^j`$ are the CT magnetic field variables. With $`t`$ for one of the indices in (13) we get the evolution equations for $`^j`$, while permuting over the spatial indices yields $`_rF_{\varphi \theta }+_\theta F_{r\varphi }+_\varphi F_{\theta r}=0,`$ which is the familiar divergence-free condition. The CT magnetic field is considered the fundamental expression of the magnetic field for evolution in the induction equation. Equation (10) relates the CT field to the field $`B^\mu `$ (or $`b^\mu `$) which is fundamental in the definition of $`T_{(EM)}^{\mu \nu }`$.

The induction equation (5) can also be rewritten by substituting definitions,

$$_\alpha \left(U^\alpha b^\beta b^\alpha U^\beta \right)=0.$$
(15)
By expanding this equation using the product rule and applying the orthogonality condition (11), we obtain the identity

$$U_\nu b^\mu _\mu U^\nu =0.$$
(16)
Proofs of the results for (10), (12), and (15) can be found in Appendix A.

In preparation for discretizing these equations we list the specific variables that will be used in the code. First, the transport velocity (also known as the coordinate velocity) $`V^\mu `$ is defined

$$V^\mu =\frac{U^\mu }{U^t},$$
(17)
where $`U^t=W/\alpha `$, and $`W`$ is the relativistic gamma-factor, which can be expressed as

$$W=\left[g^{tt}\left(g_{tt}+2g_{t\varphi }V^\varphi +g_{rr}(V^r)^2+g_{\theta \theta }(V^\theta )^2+g_{\varphi \varphi }(V^\varphi )^2\right)\right]^{1/2},$$
(18)
a result which follows directly from the definition of transport velocity and the normalization condition on the $`4`$-velocity,

$$g_{\mu \nu }U^\mu U^\nu =1.$$
(19)
The magnetic field of the fluid is described by two sets of variables, the constrained transport (CT) magnetic field, $`^i`$, and magnetic field $`4`$-vector $`b^\mu `$. The latter is fundamental to the definition of the total four momentum,

$$S_\mu =(\rho h+b^2)WU_\mu ,$$
(20)
and the normalization condition

$$g^{\mu \nu }S_\mu S_\nu =\left[(\rho h+b^2)W\right]^2,$$
(21)
which is algebraically equivalent to (19). Finally, we define auxiliary density and energy functions $`D=\rho W`$ and $`E=Dϵ`$. The set of variables $`D`$, $`E`$, $`S_\mu `$, $`^i`$, $`V^i`$, and $`b_\mu `$ will be those on which the numerical scheme is built.

Since we are working in a coordinate basis, we also record here the identities for the divergence of a four-vector and a tensor

$$_\mu \left(fv^\mu \right)=\frac{1}{\sqrt{g}}_\mu \left(\sqrt{g}fv^\mu \right),$$
(22)
$$_\mu \left(\beta ^{\mu \nu }\right)=\frac{1}{\sqrt{g}}_\mu \left(\sqrt{g}\beta ^{\mu \nu }\right)+\mathrm{\Gamma }_{ϵ\mu }^\nu \beta ^{\mu ϵ}.$$
(23)
### 2.1 The Equations of Baryon, Energy, and Momentum Conservation

The equation of baryon conservation (2) is unaltered by the presence of magnetic fields, and can be expanded in terms of the code variables to read (keeping in mind that the metric is stationary),

$$_tD+\frac{1}{\sqrt{\gamma }}_j(D\sqrt{\gamma }V^j)=0.$$
(24)
The equation of energy conservation is derived by contracting (3) with $`U_\nu `$,

$$U_\nu _\mu T^{\mu \nu }=U_\nu _\mu \left\{\left(\rho h+b^2\right)U^\mu U^\nu +\left(P+\frac{b^2}{2}\right)g^{\mu \nu }b^\mu b^\nu \right\}=0.$$
(25)
By using identity (16) and the law of baryon conservation (2), we recover the local energy conservation law

$$_\mu \left(\rho ϵU^\mu \right)+P_\mu U^\mu =0,$$
(26)
which is unaffected by the presence of magnetic fields in the ideal MHD limit and corresponds to an equation of entropy conservation. Applying the definition for the auxiliary energy function $`E`$, the energy equation is rewritten as follows:

$$_t\left(E\right)+\frac{1}{\sqrt{\gamma }}_i\left(\sqrt{\gamma }EV^i\right)+P_t\left(W\right)+\frac{P}{\sqrt{\gamma }}_i\left(\sqrt{\gamma }WV^i\right)=0,$$
(27)
which is the same as in HSW.

The momentum conservation equations follow from

$$_\mu T_{}^{\mu }{}_{\nu }{}^{}=_\mu \left\{\left(\rho h+b^2\right)U^\mu U_\nu +\left(P+\frac{b^2}{2}\right)\delta _{}^{\mu }{}_{\nu }{}^{}b^\mu b_\nu \right\}=0.$$
(28)
Using the definition of momentum $`S_\nu `$ the first term in the preceding expression can be rewritten as $`S_\nu V^\mu /\alpha `$ and simplified to $`S_\nu S^\mu /\alpha S^t`$, and also using $`\mathrm{\Gamma }_{\mu \nu }^{}{}_{}{}^{\gamma }\mathrm{\Pi }_{\gamma }^{}{}_{}{}^{\mu }=\frac{1}{2}\mathrm{\Pi }^{\gamma \mu }_\nu g_{\mu \gamma }=\frac{1}{2}\mathrm{\Pi }_{\gamma \mu }_\nu g^{\mu \gamma }`$, which holds for any symmetric tensor, we rewrite the momentum equation as

$$\frac{1}{\alpha \sqrt{\gamma }}_\mu \sqrt{\gamma }S_\nu V^\mu +\frac{1}{2\alpha }\frac{S_\alpha S_\beta }{S^t}_\nu g^{\alpha \beta }+_\nu \left(P+\frac{b^2}{2}\right)\frac{1}{\alpha \sqrt{\gamma }}_\mu \alpha \sqrt{\gamma }b^\mu b_\nu \frac{1}{2}b_\alpha b_\beta _\nu g^{\alpha \beta }=0.$$
(29)
To obtain the final form of the equations, multiply (29) by the lapse $`\alpha `$, split the $`\mu `$ index into its space ($`i`$) and time ($`t`$) components, and restrict $`\nu `$ to the spatial indices ($`j`$) only:

$$_t\left(S_j\alpha b_jb^t\right)+\frac{1}{\sqrt{\gamma }}_i\sqrt{\gamma }\left(S_jV^i\alpha b_jb^i\right)+\frac{1}{2}\left(\frac{S_ϵS_\mu }{S^t}\alpha b_\mu b_ϵ\right)_jg^{\mu ϵ}+\alpha _j\left(P+\frac{b^2}{2}\right)=0.$$
(30)
The $`\nu `$ index can be restricted to the spatial indices because the equation that arises from $`\nu =t`$ for the time components of momentum and magnetic fields is redundant, corresponding to a total energy conservation equation. In our formalism, we solve equation (27) separately for the internal energy. The time component of $`S_\mu `$ is obtained from the normalization condition on the momentum (21), and on the relations between components of the magnetic field.

We note that there are alternate ways to write this equation that, while analytically equivalent, correspond to significant differences in numerical implementation. One is to evolve directly the total stress energy $`T^{\mu \nu }`$ and solve algebraically for the primitive variables. This approach is taken by Koide, Shibata, & Kudoh (1999), Komissarov (1999), and Gammie, McKinney & Tóth (2002). Yet another alternative is to write the divergence of the electromagnetic stress-energy tensor in the form $`J_\mu F^{\mu j}`$ (the familiar $`J\times B`$ force from Newtonian MHD) and regard it as a source term for the hydrodynamic momentum evolution; this approach was used by Wilson (1978).

### 2.2 The Induction Equation

The equations of energy and momentum conservation have been expressed in terms of the magnetic field $`4`$-vector $`b^\mu `$ since this yields expressions structurally compatible with the existing hydrodynamic equations. This has important advantages in the numerical implementation of these equations. However, the induction equation is evolved with the constrained transport (CT) approach of EH88 using the variables $`^i`$. The practical consequence is that we must translate between the two sets of variables.

Working with the second of Maxwell’s equations, (5), we have previously identified the space-space components of $`F_{\mu \nu }`$ with the CT field $`^j`$. It is easy to show that we obtain the same result using (15),

$$_\mu \left(U^\nu b^\mu b^\nu U^\mu \right)=_\mu \sqrt{\gamma }W\left(V^\nu b^\mu b^\nu V^\mu \right)=0.$$
(31)
By splitting the above equation into two pieces and defining the CT magnetic field as

$$^i=\sqrt{4\pi }\sqrt{\gamma }W(b^iV^ib^t),$$
(32)
it is possible to write

$`_j\left(^j\right)`$ $`=0`$ $`(\nu =0),`$ (33)
$`_t\left(^i\right)_j\left(V^i^j^iV^j\right)`$ $`=0`$ $`(\nu =i),`$ (34)
the starting point for the EH88 CT scheme. In addition, the time-space components of $`F_{\mu \nu }`$ can be obtained by applying the condition of infinite conductivity, $`F^{\mu \nu }U_\nu =0`$, to obtain

$$F_{tr}=V^\varphi ^\theta V^\theta ^\varphi ,F_{t\theta }=V^r^\varphi V^\varphi ^r,F_{t\varphi }=V^\theta ^rV^r^\theta .$$
(35)
These can be identified with $`^j`$, the CT electromotive forces (EMFs, EH88).

The divergence-free character of the magnetic field is maintained by evolving the induction equation using these CT variables. Hence the magnetic induction $`b^\mu `$ is derived from the CT magnetic field. We obtain the spatial components using (32). To obtain $`b^t`$, we use the definition of the magnetic induction (8),

$$b^\mu =\frac{1}{\sqrt{4\pi }}{}_{}{}^{}F_{}^{\mu \nu }U_\nu =\frac{1}{2\sqrt{4\pi }}ϵ^{\mu \nu \delta \gamma }U_\nu F_{\delta \gamma },$$
(36)
and we expand this definition and substitute for the Faraday tensor. This allows us to obtain the time-component of the magnetic field,

$`b^t`$ $`=`$ $`{\displaystyle \frac{W}{\sqrt{4\pi }\alpha ^2\sqrt{\gamma }}}\left({\displaystyle \frac{V^r^r}{g^{rr}}}+{\displaystyle \frac{V^\theta ^\theta }{g^{\theta \theta }}}+{\displaystyle \frac{g^{tt}V^\varphi g^{t\varphi }}{g^{tt}g^{\varphi \varphi }(g^{t\varphi })^2}}^\varphi \right).`$ (37)
This relationship can also be obtained from the orthogonality condition (11). Using this result and the relations in (32) we obtain, after some algebra, a useful expression for $`b^2`$,

$`b^2`$ $`=`$ $`{\displaystyle \frac{1}{4\pi \gamma W^2}}\left({\displaystyle \frac{_{}^{r}{}_{}{}^{2}}{g_{rr}}}+{\displaystyle \frac{_{}^{\theta }{}_{}{}^{2}}{g_{\theta \theta }}}+{\displaystyle \frac{g^{tt}_{}^{\varphi }{}_{}{}^{2}}{g^{tt}g^{\varphi \varphi }(g^{t\varphi })^2}}\right)`$
$`+{\displaystyle \frac{1}{4\pi \gamma \alpha ^2}}\left({\displaystyle \frac{V^r^r}{g^{rr}}}+{\displaystyle \frac{V^\theta ^\theta }{g^{\theta \theta }}}+{\displaystyle \frac{(g^{tt}V^\varphi g^{t\varphi })^\varphi }{g^{tt}g^{\varphi \varphi }(g^{t\varphi })^2}}\right)^2.`$
The advantage of this form is that it guarantees that $`b^2`$ is positive. The derivation is supplied in Appendix A.

## 3 Implementing GRMHD

The GRMHD code evolves time-explicit, operator-split, finite difference forms of equations (24), (27), (30), and (34). Variables are placed on a fixed spatial grid using Boyer-Lindquist coordinates. They are advanced in time over a timestep size $`\mathrm{\Delta }t`$ that remains fixed for the duration of the simulation, and is determined by the extremal light-crossing time for a grid zone, as described in HSWb. The evolution algorithm is a three-dimensional generalization of the solver described in HSWb extended to include the contribution of magnetic fields. The hydrodynamic portion of the solver was described in De Villiers and Hawley (2002).

One timestep consists of three sequential sub-steps: the induction step, the transport step, and finally the source step. The induction step updates the magnetic field vector $`^i`$ using the Constrained Transport framework of EH88. This is discussed in greater detail below. After the update, the CT field is transformed to the magnetic field in the fluid rest frame, $`b^i`$. Velocities $`V^i`$ are computed using the normalization condition during the source step.

Three versions of GRMHD have been used in development and testing: (1) a 1D version with Minkowski metric to do Alfvén and shock wave tests, described here in §4, (2) a 2D axisymmetric version in the Kerr metric to perform the tests described in §5, and (3) the full 3D production version of the code. The 3D version of the code uses message passing parallelism with domain decomposition, where the global grid is partitioned into subgrids, with each subgrid assigned to a processor. Data on each subgrid are evolved independently during a timestep and data on subgrid boundaries are exchanged when required through message-passing calls. This results in a highly scalable code that exhibits good speedup over the full range of practically realizable subgrids.

### 3.1 Transport and Source Steps

The transport step, as its name implies, handles the transport terms in the continuity, energy, and momentum equations,

$$\left[_t𝒳\right]_{transport}=\frac{1}{\sqrt{\gamma }}_i\left[\sqrt{\gamma }𝒳V^i\right],\left(𝒳D,E,S_j\right).$$
(39)
Other than the modified definition of the inertia, these expressions are completely analogous to those in HSWb.

The source step handles the remaining terms in the energy and momentum equations. The energy terms

$$\left[_tE\right]_{source}=P_tW\frac{P}{\sqrt{\gamma }}_i\left[\sqrt{\gamma }WV^i\right]$$
(40)
are evolved as described in HSWb. The momentum source terms contain new elements, namely

$$\left[_tS_j\right]_{source}=\alpha _t\left[b_jb^t\right]+\frac{_i\left[\sqrt{\gamma }\alpha b_jb^i\right]}{\sqrt{\gamma }}\left[\frac{S_ϵS_\mu }{S^t}\alpha b_\mu b_ϵ\right]\frac{_jg^{\mu ϵ}}{2}\alpha _j\left[P+\frac{b^2}{2}\right].$$
(41)
In the following expressions, superscript indices for the current time step are understood. The spatial differencing of some value $`f`$ is shown in compact form through the shift operator $`D_i(f)`$; any averaging of zone- and face-centered quantities that may be required to properly center calculations is done as described in HSWb.

The time derivative of the magnetic fields is evaluated from the newly-obtained values from the induction step, and stored values from the previous time step,

$$S_{j}^{}{}_{}{}^{n+1}=S_j+\alpha \left[\left(b_jb^t\right)^{(n+1)}\left(b_jb^t\right)\right].$$
(42)
The magnetic gradient term discretizes readily,

$$S_{j}^{}{}_{}{}^{n+1}=S_j+\frac{\mathrm{\Delta }t}{\sqrt{\gamma }}\frac{D_i\left(\sqrt{\gamma }\alpha b_jb^i\right)}{\mathrm{\Delta }x^i},$$
(43)
as does the pressure gradient term,

$$S_{j}^{}{}_{}{}^{n+1}=S_j\mathrm{\Delta }t\alpha \frac{D_j\left(P+b^2/2\right)}{\mathrm{\Delta }x^j}.$$
(44)
Finally, the accelerations due to gradients in the metric,

$$S_{j}^{}{}_{}{}^{n+1}=S_j\frac{\mathrm{\Delta }t}{2}\left[\frac{S_ϵS_\mu }{S^t}\alpha b_\mu b_ϵ\right]_jg^{\mu ϵ}.$$
(45)
are computed using $`S^t`$ and $`S_t`$ obtained from $`S_j`$ using momentum normalization (21); the metric derivatives, $`_jg^{\mu ϵ}`$, are evaluated analytically over the whole grid and stored in memory at the beginning of the simulation.

In the source step, after the pressure gradient and acceleration terms have been computed, the transport velocity $`V^i=S^i/S^t`$, and the Lorentz factor $`W=\alpha U^t`$ are obtained using velocity normalization. To guarantee that the normalization is well-behaved in evacuated zones, we work with the momenta $`S_i`$ rather than applying the normalization condition directly to the velocities. Although the two normalizations are analytically equivalent, use of the momentum protects against numerical rounding errors that cause spurious floating point exceptions when extremely tenuous material is moving at velocities close to that of light.

The source step also includes an artificial viscosity in both the energy and momentum equations to provide a mechanism to increase the entropy of the fluid in shocks. The artificial viscosity is unchanged from HSWb, except for modifying the definition of the inertia to include the magnetic energy.

### 3.2 The Induction Step: A General Relativistic Form of the MOCCT Algorithm

The CT framework of EH88 is based on two basic ideas: the use of $`F_{\mu \nu }`$ as the fundamental variable so that the induction equation reduces to simple partial derivatives, and a staggered-grid centering of terms so that the numerical divergence of the CT magnetic field is constrained to be zero. The CT magnetic field component $`^j`$ is located on the $`j`$ face of a zone (a cube in three dimensions), and the CT EMFs, $`^i`$, that are differenced to evolve $`^j`$ are located on the zone edges that define the $`j`$ face.

The CT formalism imposes no stringent conditions on exactly how the EMFs are computed. EH88 noted that some of the magnetic field variables appear in the discretrized equations for the EMFs in the guise of transport terms, requiring that they be upwinded for stability, while others appeared as shear terms, apparently requiring no special treatment. EH88 used simple space-centered differencing for those terms, but Stone & Norman (1992) found that this prescription is inadequate. The results from simple Alfvén pulse propagation tests where the Alfvénic velocity exceeds that of the background fluid found substantial undamped dispersion error, resulting in excessive zone-to-zone noise.

Stone & Norman (1992) proposed to calculate the CT EMFs by solving a one dimensional analytic linear Alfvén wave equation to obtain improved values of $`^j`$ and $`V^j`$ in the construction of the EMFs. This approach became known as the Method of Characteristics Constrained Transport (MOCCT). Since then several variants have been proposed (e.g. Hawley & Stone 1995; Clarke 1996) to further improve the computation of Newtonian MHD. While the traditional approach to MoC does not carry over readily to GR MHD, it is nevertheless possible to address the problem by using the same motivation that underlies MOCCT. Instead of using an analytic Alfvén wave solution we use a numerical solution of the simplified one-dimensional induction and momentum equations to obtain the improved EMFs.

Our implementation of CT for the GRMHD code is accomplished as follows. Following the example of Hawley & Stone (1995), the starting point for constructing, for instance, the $`z`$-EMF,

$$^z=V^x^yV^y^x,$$
(46)
lies in obtaining the discretized solutions to the 1D incompressible MHD equations. On a staggered mesh, $`^z`$ is located on zone edges in such a way that it is constructed from $`x`$\- and $`y`$-components relocated in the $`y`$\- and $`x`$-directions respectively. To take a specific example, we would like obtain an edge-centered estimate $`_{}^{y}{}_{}{}^{}`$ to construct the above EMF. We do this by solving the 1D linearized induction equation in the $`x`$-direction (i.e. $`_y_z0`$) for small transverse ($`y`$) perturbations. This simplified induction equation has the form,

$`_x^x`$ $`=`$ $`0,`$ (47)
$`_t^y`$ $`=`$ $`^x_x\left(V^y\right)V^x_x\left(^y\right).`$ (48)
By applying the usual differencing we obtain for the left-hand side,

$$_t^y=\frac{_{}^{y}{}_{}{}^{}\overline{^y}}{\mathrm{\Delta }t},$$
(49)
where the overbar denotes the original EH88 upwinding. We can now readily obtain estimates of the starred quantity,

$$_{}^{y}{}_{}{}^{}=\overline{^y}+\frac{1}{2}\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}\left[V^xD_x(^y)+^xD_x(V^y)\right],$$
(50)
where the factor of $`1/2`$ ensures proper time-centering of the starred quantity (the starred values are centered between the $`n`$ and $`n+1`$ time levels). A similar approach can now be used to generate $`_{}^{x}{}_{}{}^{}`$, and produce an improved estimate of $`^z`$. The procedure is just as easily applied to the other two EMFs to complete the process for all CT variables. Note that because of the fully covariant nature of the CT equations, this process applies to any coordinate system.

For consistency, we also need to evolve the linearized momentum equation for the transverse velocities. Care must be taken in deriving these equations since, in contrast to the CT induction equation, metric terms are explicitly involved in the expressions. Continuing with the example for $`^z`$, and regarding the $`x`$ direction to be either $`r`$ or $`\theta `$, and $`y`$ a direction transverse to this, we obtain

$$\left(_t+\xi V^x_x\right)V^y\kappa ^x_x^y0.$$
(51)
where,

$`\xi `$ $`=`$ $`{\displaystyle \frac{\rho hW^2b^2}{\rho hW^2+b^2}},`$ (52)
$`\kappa `$ $`=`$ $`{\displaystyle \frac{\alpha ^2}{\left(W\sqrt{4\pi \gamma }\right)^2\left(\rho hW^2+b^2\right)}}.`$ (53)
(The details are found in Appendix A.) Note that the above equations are similar to the Newtonian results, but with the relativistic factors $`\xi `$ and $`\kappa `$. By simply applying the usual differencing, however, we readily obtain expressions for an improved estimate of the transverse velocity,

$$_tV^y=\frac{V_{}^{y}{}_{}{}^{}\overline{V^y}}{\mathrm{\Delta }t}$$
(54)
where the overbar denotes the edge-centered average, from which,

$$V_{}^{y}{}_{}{}^{}=\overline{V^y}+\frac{1}{2}\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}\left[\xi V^xD_x(V^y)+\kappa ^xD_x^y\right].$$
(55)
In the Kerr metric, the calculation in the $`\varphi `$ direction is somewhat different due to the $`g^{t\varphi }`$ metric term. In this case we obtain

$$\left[_t+\left(\xi V^\varphi +\lambda \right)_\varphi \right]V^y\kappa ^\varphi _\varphi ^y0.$$
(56)
where $`\xi `$ and $`\kappa `$ are unchanged and

$$\lambda =\frac{2b^2g^{t\varphi }}{g^{tt}\left(\rho hW^2+b^2\right)}.$$
(57)
Combined with the starred quantities for the CT variables, we now have a complete set of variables from which to construct the EMFs.

## 4 1D Test Cases - Minkowski Spacetime

The first series of tests is designed to verify code performance in the special relativistic limit. For this purpose we use a one-dimensional version of the code and the Minkowski metric. The tests include Alfvén wave propagation, fast and slow magnetosonic shocks, rarefaction waves, and both relativistic and non-relativistic shock tubes. Many of these tests have appeared previously in the literature in one form or another. In particular, Komissarov (1999) (hereafter referred to as K99) presented a series of challenging test problems for a special relativistic MHD Godunov scheme, and Stone et al. (1992) presented a test suite for Newtonian MHD codes. Our presentation will conform to the formats of these published results to facilitate comparison.

### 4.1 Alfvén Wave Propagation

The propagation of a linear Alfvén wave is simple in conception, but, as demonstrated by Stone & Norman (1992), can be quite revealing as a numerical test. The special relativistic Alfvén wave includes the displacement current normally neglected in Newtonian MHD. This is what limits wave speeds to less than $`c`$.

We consider a fixed background magnetic field $`^x`$ in Minkowski spacetime with constant fluid velocity $`V^x`$ and a small transverse perturbation with velocity $`V^y`$ and CT magnetic field $`^y`$. The induction and momentum equations for the perturbed quantities have the form given in equations (48) and (51). Equation (51) can be combined with the induction equation (48) by taking partial derivatives \[$`\left(_t+V^x_x\right)`$ for the former and $`_x`$ for the latter\] to get an equation for $`V^y`$,

$$\left[_t^2+\left(1+\xi \right)V^x_t_x+\left(\xi (V^x)^2\kappa (^x)^2\right)_x^2\right]V^y=0$$
(58)
which describes the propagation of Alfvén waves with speeds

$$v_{A}^{}{}_{}{}^{(\pm )}=\frac{V^x\pm \eta \sqrt{\eta ^2+W^2}}{1+\eta ^2}\pm \sqrt{\frac{b^2}{\rho h+b^2}}(forV^x=0),$$
(59)
where $`\eta ^2=b^2/(\rho hW^2)`$. (The equation for $`^y`$ can be obtained in a similar manner.) Figure 1 plots $`v_{A}^{}{}_{}{}^{(\pm )}`$ as a function of the background velocity of the fluid ($`V^x`$) for four values of $`\beta `$ and illustrates several key properties of the wave speeds: (1) for $`V^x=0`$, waves move with equal and opposite speed (for all $`\beta `$); (2) $`v_{A}^{}{}_{}{}^{(\pm )}`$ is constrained to the upper bound of $`c1`$; (3) for a given $`\beta `$, there is a $`V_0^x`$ that yields a stationary wave, $`v_{A}^{}{}_{}{}^{()}=0`$; (4) for $`V^x>V_0^x`$, both waves move in the same direction, with $`v_{A}^{}{}_{}{}^{()}`$ lagging behind the main fluid flow, and $`v_{A}^{}{}_{}{}^{(+)}`$ leading; (5) as $`V^xc`$, all wave speeds converge, $`v_{A}^{}{}_{}{}^{(\pm )}c`$.

We describe four tests, listed in Table 1, carried out in a periodic box of width $`x_{max}x_{min}=3.0`$. In each case the fluid has constant density $`\rho =1`$ and energy density $`ϵ=10^2`$. Individual tests are distinguished by the amplitude of the $`x`$-component of the CT magnetic field $`^x`$, given here in terms of $`\beta =\sqrt{2P/b^2}`$, and the amplitude of the background $`x`$-component fluid velocity $`V^x`$. The waves are initialized with a square pulse in the transverse velocity $`V^y`$ with amplitude $`A_0=10^3`$, namely

$$V^y=\{\begin{array}{cc}A_0\hfill & \text{for }1.0<x<1.5\hfill \\ A_0\hfill & \text{for }1.5x<2.0\hfill \end{array}.$$
(60)
The transverse magnetic field is initially zero, $`^y=0`$.

Table 1 summarizes the results and Figure 2 overlays $`V^y`$ for each test at time $`t_{final}`$ at three grid resolutions, $`512`$, $`1024`$, and $`2048`$ zones. In test ALF2 the fast and slow pulses have wrapped around and re-entered the grid from the left, the fast pulse having done so twice. In ALF4, the fast pulse has also wrapped around while the slow pulse is stationary. (In this linear regime the pulses overlap cleanly provided that the spatial resolution is sufficiently high to minimize losses due to numerical diffusion.) Measuring from the center of each pulse (originally at $`x=1.5`$) we verify that the observed pulse speeds agree with (59). With a stationary background, ALF1, the pulses separate from the initial perturbation with equal and opposite velocities and also with equal amplitudes such that $`A_0=A^++A^{}`$. In a moving background, the splitting of the velocity pulses is no longer symmetrical, as is especially obvious in ALF2 (although $`A_0=A^++A^{}`$ still holds). The slow velocity pulse has a greater amplitude than the fast pulse. We find that the ratio of the velocity amplitudes scales inversely with their Lorentz factors,

$$\frac{A^+}{A^{}}=\frac{W(v_{A}^{}{}_{}{}^{})}{W(v_{A}^{}{}_{}{}^{+})}.$$
(61)
It is easy to verify that this holds for each model in Table 1. The pulses in the magnetic field variable, $`^y`$, on the other hand, are equal in magnitude and opposite in phase for all $`V^x`$. The numerical values given in Table 1 can be verified by direct computation using results presented in Appendix B.

### 4.2 Magnetosonic Shock and Rarefaction Tests

Shock waves and rarefactions test a code’s ability to respond to discontinuities and to maintain the jump conditions that correspond to the conservation properties of the flow. In this section we consider a number of one-dimensional shock tests, both relativistic and non-relativistic, that involve nonzero transverse magnetic fields.

We begin with relativistic magnetosonic shock wave tests; the initial left ($`x<0.0`$) and right ($`x0.0`$) states are given in Table 2. In each case we verified that these states satisfy the relativistic shock invariants (see Appendix C). In this and in subsequent tables we note those tests that were described in K99.<sup>1</sup><sup>1</sup>1K99 magnetic field variables $`_{K99}^i`$ are related to our CT-variables by a normalization factor, $`^i=\sqrt{4\pi }_{K99}^i`$. Also, a revised version of K99 Table 2 was given by Komissarov (2002).

The first test, the slow magnetosonic shock, is shown in Figure 3, which plots $`U^x`$, $`\rho `$, gas pressure $`P`$, and magnetic pressure $`P_{mag}`$ at $`t=2.0`$, by which point the shock front has moved from $`x=0.0`$ to $`x=0.5`$. These graphs are quite similar to the results of K99. There is a small step in density in the high-density medium to the right of the moving shock. When we repeat our tests at lower resolution, our step feature resembles the “wiggle” visible in K99. The origin of this small step may be due to a small discrepancies in the initial state (i.e. the data in K99, which we use here, differ slightly from what we compute using the method outlined in Appendix C).

Figure 4 shows plots of $`U^x`$, $`\rho `$, gas pressure $`P`$, and magnetic pressure $`P_{mag}`$ for the fast magnetosonic shocks at $`t=2.5`$. In the case of the stationary shock (I), we see that artificial viscosity has slightly altered the shock front, and that there is an overvalue in the left state $`U^x`$, and undervalue in $`\rho `$ and $`P_{gas}`$. In the case of moving shock (II), the shock front has traveled to $`x=0.5`$, as expected for the given shock speed and the value of $`t_{final}`$, and the left and right states agree with the analytic values. In the case of the $`W3.8`$ shock (III), a small mismatch has developed in the post-shock gas pressure near the shock edge. In addition, the speed of the shock slightly exceeds the expected value $`0.2`$. As discussed by HSWb and studied in detail by Norman & Winkler (1986), the artificial viscosity used in the code to increase the entropy of the fluid through a shock produces sufficient heating for $`W<2`$, but increasingly underestimates the heating as $`W`$ increases. The results of this shock test are consistent with these earlier findings.

We next consider a pair of challenging rarefaction tests designed by K99. The initial left and right states of these tests are given in Table 3. Figure 5 shows a plot of $`U^x`$, gas pressure $`P`$ and density $`\rho `$ for the switch-off fast rarefaction at $`t=1.0`$ and the switch-on slow rarefaction at $`t=2.0`$. Our results generally agree with those of K99. The fast rarefaction extends from $`x=0.5`$ to $`x=+0.6`$ and the slow rarefaction extends from $`x=0.5`$ to $`x=+0.9`$. There is high-frequency noise in the plot of $`U^x`$ and $`P`$ to the right of the slow rarefaction. We attribute this noise to an artifact of the routine that converts CT-variables to fluid-frame magnetic field. This test is particularly sensitive to the averaging techniques that are an unavoidable part of staggered-mesh discretizations. This noise is reduced at lower grid resolution, where increased numerical diffusion acts to damp out the oscillations.

Next we consider a series of shock tube tests that combine strong shocks and rarefactions of various types. We ran two test problems from K99, the Brio & Wu (1988) shock tube test as described in Stone et al. (1992), an alternate version of the Brio & Wu shock tube using the same initial state but with $`\gamma =4/3`$ instead of $`\gamma =2`$ to ensure that the sound speed is not relativistic, and a non-relativistic shock tube from Ryu & Jones (1995). In the non-relativistic shock tubes the values of pressure and density are rescaled to ensure that the sound speed is much less than $`c`$. Table 4 summarizes the various states for the shock tube tests, following the convention of Stone et al. (1992). Figures 6, 7, 8, and 9 show our numerical results for the evolution of the shock tubes at a grid resolution of $`2048`$ zones. Some ringing is present in the non-relativistic shock tubes, at the location of the compound wave and at the slow shock in the density variable, for instance, whereas this ringing is absent from the relativistic shock tubes. This appears to be a manifestation of decreased dispersion error at relativitistic velocities.

## 5 1D and 2D Test Cases - Kerr Spacetime

To validate the general relativistic aspects of the code we require code tests that involve magnetic fields and the Kerr geometry. Test problems with analytic solutions include one dimensional magnetized Bondi inflow, and a two-dimensional Bondi flow along a split monopole field; in these tests the field does not affect the dynamics of the problem. We also perform a one-dimensional magnetized inflow test problem where the field does act dynamically (Gammie 1999). Finally, we have considered two-dimensional magnetized constant-$`l`$ tori that are subject to the magnetorotational instability. While there is no analytic solution to this problem, the results can be compared to those obtained in non-relativistic simulations.

### 5.1 Magnetized Bondi Flow and Split Monopole Tests

The analytic expressions for the hydrodynamic Bondi solution of HSWa are unchanged in the presence of a radial magnetic field (see proof in Appendix A). However, this result does not necessarily hold for a numerical solution to the Bondi problem. A non-trivial test of the magnetic components of the code consists, therefore, in maintaining the equilibrium Bondi solution for any magnitude of radial magnetic field. The magnetized Bondi problem is set up in the Schwarzschild metric ($`a=0`$). The radial grid extends from just outside the event horizon at $`r=2.20M`$ to $`r=25.0M`$; the critical point is at $`r_{crit}=8.0`$. Grids of $`64`$, $`128`$, $`256`$, $`512`$, and $`1024`$ points are used, and the magnetic field at the critical point is set to $`\beta =100,10`$, and $`1`$, where the $`\beta `$ parameter sets the grid-averaged ratio of gas to magnetic pressure. The Bondi solution is generated in the initialization routine by a numerical root-finder on a logarithmically scaled grid of $`1024`$ points. This solution is then sampled at the required resolution for the particular test run (to ensure proper alignment of the coarser grids). To complete the test, the unmagnetic case is also presented. The solution is allowed to evolve to a time sufficient to allow numerical convergence on each grid ($`t_{final}=100M`$).

Figure 10 shows the results of spatial accuracy tests. The spatial accuracy of variable $`𝒳`$ is obtained from the $`L_1`$ norm over an interval $`r_{min}rr_{max}`$,

$$L_1(𝒳)=\frac{1}{r_{max}r_{min}}𝒳_{final}𝒳_0𝑑r,$$
(62)
where $`𝒳_{final}`$ denotes the converged solution, and $`𝒳_0`$ the original solution, as described above. Different regions of the flow are dominated by numerical errors either from a particular routine in the numerical solver or a particular code variable. Here, the bounds of integration were chosen to range from $`r_{min}=7.0`$ to $`r_{max}=25.0`$, which represents a region of the flow where pressure error dominates. Figure 10 is for the variables $`d`$, and is typical of error curves for all code variables. The code is clearly first-order accurate for these variables.

In the split monopole test, we initialize a 2D $`(r,\theta )`$ grid with the radial Bondi solution, but the radial magnetic field has opposite polarity above and below the equator, hence establishing a current sheet along the equator. This test is made into a numerical stability test by applying small enthalpy perturbations to the initial state (without this, the code is simply a set of 1D radial tests). The test is repeated with increasingly large values of magnetic field, until the code breaks. There is no known analytic solution for the evolution of a current sheet, but it can be used as a probe of the code’s stability in the presence of strong, fluctuating current sheets. The split monopole test was carried out with $`\beta =100,10,1`$ at the critical point, and a grid of $`180^2`$ scaled logarithmically in both the $`r`$ and $`\theta `$ coordinates, with an initial $`1\%`$ random enthalpy perturbation. The test was allowed to run to $`t=1250M`$ ($`10`$ free-fall times). For $`\beta =100,10`$, the code variables show no visible change over the course of the run; the enthalpy perturbations do not grow, nor do they trigger a numerical instability. For $`\beta =1`$, the late-time solution shows evidence of reconnection at the equator near the horizon, and erratic behavior. The effective $`\beta `$ where reconnection occurs is of order $`10^1`$. Given the moderate resolution, the results are consistent with the spatial accuracy tests in that the region near the horizon requires high resolution in order to minimize numerical error.

### 5.2 Magnetized Gammie Inflow

Gammie (1999) describes a simplified 1D model for inflow from the inner edge of an accretion disk in the Kerr geometry that serves as an excellent test problem for relativistic MHD. This solution is similar to the Weber & Davis (1967) rotating magnetized wind which was used as a test problem Stone et al. (1992). The model assumes that the inflow is cold ($`h1`$), time-steady, purely radial, confined to the equatorial plane, and initiated from the edge of an accretion disk located at the marginally stable orbit ($`r_{mso}`$). At $`r_{mso}`$ the solution is assumed to satisfy the conditions $`U^r(r_{mso})=0`$ and $`\mathrm{\Omega }_FV^\varphi (r_{mso})=1/(r_{mso}^{}{}_{}{}^{3/2}+a)`$. The flow is characterized by a fast magnetosonic critical point. The numerical procedure to find a particular solution to this inflow is described in Appendix D.

Our test consists in replicating the solution described by Gammie et al. (2002) for a flow near a Kerr black hole with $`a=0.5`$ ($`r_{horizon}=1.866`$, $`r_{mso}=4.233`$, $`\mathrm{\Omega }_F=0.108588`$), with $`^r=0.5`$, $`F_M=1`$, $`F_L=2.815344`$, $`F_E=0.9083782`$. The critical point for this test case is $`(r_{crit},U_{crit}^r)=(3.616655,0.04054696)`$.

Figure 11 shows a plot of the solution for this test problem on a logarithmically scaled grid of $`1024`$ points. The inner edge of the grid is at $`r_{min}=2.00`$ and the outer edge of the grid is at $`r_{max}=4.04`$. In Boyer-Lindquist coordinates, the event horizon cannot be part of the computational domain because of the coordinate singularity. In the Gammie inflow test, the marginally stable orbit, which is the formal outer boundary of the problem, is also excluded since the density is divergent as $`r`$ approaches $`r_{mso}`$ from below. The plot overlays the late-time solution (after $`1000`$ time steps) on the initial state, but the two curves are indistinguishable at the chosen plot scale. Figure 12 shows the results of convergence tests carried out analogously to those for the Bondi problem, with the range of integration going from $`r_{min}=2.7`$ to $`r_{max}=4.04`$. Here we have an example of a highly-resolved near-horizon region, and the convergence rate is close to second order.

### 5.3 Constant-$`l`$ Disks

The final set of tests involve constant-$`l`$ disks ($`l=U_\varphi /U_t`$, the specific angular momentum), solutions of the axisymmetric GR hydrodynamic equations described in HSWa. Here, we add weak poloidal magnetic field loops that overlay the hydrodynamic solution; this will trigger the MRI. We quantify the strength of the magnetic field through the $`\beta `$-parameter, the ratio of the volume-averaged gas pressure to magnetic pressure.

The initial magnetic field is obtained from the definition of $`F_{\mu \nu }`$ in terms of the $`4`$-vector potential, $`A_\mu `$,

$$F_{\mu \nu }=_\mu A_\nu _\nu A_\mu .$$
(63)
Restricting the field to poloidal loops is done by using $`A_\mu =(A_t,0,0,A_\varphi )`$, where

$$A_\varphi =\{\begin{array}{cc}k(\rho \rho _{min})\hfill & \text{for }\rho \rho _{min}\hfill \\ 0\hfill & \text{for }\rho <\rho _{min}\hfill \end{array}.$$
(64)
Using (63), it follows that $`^r=_\theta A_\varphi `$ and $`^\theta =_rA_\varphi `$. This choice of $`A_\varphi `$ produces poloidal field loops that coincide with isodensity contours. The constant $`k`$ is set by the input parameter $`\beta `$. The constant $`\rho _{min}`$ sets a suitable minimum density within the disk, and it is chosen to keep the initial magnetic field away from the outer edge of the disk.

An axisymmetric (2D) magnetized constant-$`l`$ torus was set up in both the Schwarzschild and Kerr metrics. For the Schwarzschild case, the disk has a specific angular momentum $`l=4.5`$, an initial pressure maximum at $`r=15.3M`$, and an orbital period at the pressure maximum of $`T_{orb}=376M`$. For the Kerr case, the black hole is maximally rotating, $`a=1`$, and the disk has a prograde specific angular momentum $`l=4.3`$, an initial pressure maximum at $`r=15.4M`$, and an orbital period at the pressure maximum of $`T_{orb}=386M`$. This particular choice of parameters yields a disk that is similar to the Schwarzschild case, although the equipotentials which define the overall disk structure are notably different, as can be seen in the left column of Figure 13. For both tests the average field strength is $`\beta =100`$, the initial state is perturbed with random $`1\%`$ enthalpy fluctuations. The simulations are run for $`10`$ orbits. The MRI develops after a few orbits (center column, Fig. 13), and is soon fully developed (right column, Fig. 13). By the tenth orbit, the MRI-induced turbulence has settled down considerably, as is to be expected since the imposition of axisymmetry precludes the development of the azimuthal modes which sustain the MRI. The observed qualitative differences between the Schwarzschild and Kerr cases can be attributed to the different shapes of the potential cusp near the horizon.

The density plots bear a strong resemblance to the previous axisymmetric MRI studies of Hawley (2000) in a pseudo-Newtonian potential; this indicates that the numerical solver is well behaved in the moderately strong gravitational fields that exist near the pressure maximum, where the implied code comparison is taking place. Since MRI studies are the main area of application of this code, it is reassuring that we can both trigger the MRI in a general relativistic treatment of MHD, and that we qualitatively reproduce established results. Gammie et al. (2002) obtain similar results with their MRI torus test.

## 6 Discussion

In this paper we document the development of a General Relativistic MHD code that is a direct extension of the general relativistic hydrodynamic solver of Hawley, Smarr and Wilson (1984b). As before, the equations of relativistic MHD are written in a nonconservative form, supplemented here by the relativistic induction equation. The magnetic field is added to the definition of the inertia for the four momentum and the magnetic field adds several source terms to its evolution.

The equations are evolved on a staggered grid using the time-explicit model of “transport $`+`$ source.” The induction equation is evolved using the Constrained Transport structure of EH88. We implement a version of CT that is inspired by the Method of Characteristics technique employed in the popular ZEUS code (Stone & Norman 1992), yet is simpler to implement in a general relativistic framework.

We have performed a number of 1D and 2D tests to validate the code and to note its limitations. The tests include a suite of non-relativistic and relativistic MHD shock tube tests, the fast and slow rarefaction and slow magnetosonic shock tests of Komissarov (1999), and simple Alfvén pulse propagation. The formulation for these test and their results are documented to establish a basis for comparison against other existing and future codes.

The main shortcoming of our code revealed by these tests is due to the already-known limitations of the artificial viscosity algorithm of the underlying hydrodynamic solver. Improvements to this algorithm are an ongoing area of research, but is important to stress that these limitations do not pose an undue restriction on our main area of study: dynamical, MHD-driven accretion flows around spinning black holes.

Although the set of general relativistic black hole accretion test problems is limited, we have also verified that the code reproduces the standard results for Bondi flow in the presence of magnetic fields, and have observed the triggering and evolution of the MRI in thick accretion tori for both the Kerr and Schwarzschild black holes. The MRI evolution can be compared qualitatively with pseudo-Newtonian simulations (e.g. Hawley 2000). One new test problem with an analytic solution is the magnetized Gammie (1999) inflow test. This solution examines the behavior of a dynamically important magnetic field in the strong-field region of the Kerr metric. The results of these tests provide evidence that this new solver is correctly reproducing magnetized flows outside the event horizon of a Kerr black hole. This new GR MHD code should provide a valuable tool with which to study accretion flows in the Kerr metric.

This work was supported by NSF grant AST-0070979 and NASA grant NAG5-9266. We thank Charles Gammie with whom we collaborated to develop and apply a suite of code tests for general relativistic MHD.

## Appendix A Proof or Various Results

Some of the following proofs use two identities for contractions of the Levi-Civita tensor,

$$ϵ^{\alpha \beta \delta \gamma }ϵ_{\alpha \lambda \mu \nu }=\delta _{}^{\beta \delta \gamma }{}_{\lambda \mu \nu }{}^{},$$
(A1)
$$ϵ^{\alpha \beta \delta \gamma }ϵ_{\alpha \beta \mu \nu }=2\delta _{}^{\delta \gamma }{}_{\mu \nu }{}^{},$$
(A2)
and identities for the $`\delta `$-symbol,

$$\delta _{}^{\alpha \beta }{}_{\sigma \rho }{}^{}=\delta _\sigma ^\alpha \delta _\rho ^\beta \delta _\rho ^\alpha \delta _\sigma ^\beta ,$$
(A3)
$$\delta _{}^{\alpha \beta \delta }{}_{\lambda \mu \nu }{}^{}=\delta _\lambda ^\alpha \delta _\mu ^\beta \delta _\nu ^\delta \delta _\mu ^\alpha \delta _\lambda ^\beta \delta _\nu ^\delta +\delta _\mu ^\alpha \delta _\nu ^\beta \delta _\lambda ^\delta \delta _\nu ^\alpha \delta _\mu ^\beta \delta _\lambda ^\delta +\delta _\nu ^\alpha \delta _\lambda ^\beta \delta _\mu ^\delta \delta _\lambda ^\alpha \delta _\nu ^\beta \delta _\mu ^\delta .$$
(A4)
To prove (10), it is easiest to work backwards from the result. Expand the right-hand side of (10) using (8) and (6); the resulting expression simplifies through identities (A1) and (A4), followed by the velocity normalization condition, infinite conductivity, and antisymmetry of $`F^{\mu \nu }`$:

$`ϵ_{\alpha \beta \mu \nu }B^\alpha U^\beta `$ $`=`$ $`{\displaystyle \frac{1}{2}}ϵ^{\alpha \delta ϵ\gamma }ϵ_{\alpha \beta \mu \nu }U_\delta U^\beta F_{ϵ\delta }`$
$`=`$ $`{\displaystyle \frac{1}{2}}\delta _{}^{\delta ϵ\gamma }{}_{\beta \mu \nu }{}^{}U_\delta U^\beta F_{ϵ\delta }`$
$`=`$ $`{\displaystyle \frac{1}{2}}\left(U_\beta U^\beta F_{\mu \nu }U_\mu U^\beta F_{\beta \nu }+U_\mu U^\beta F_{\nu \beta }U_\nu U^\beta F_{\mu \beta }+U_\nu U^\beta F_{\beta \nu }U_\beta U^\beta F_{\nu \mu }\right)`$
$`=`$ $`F_{\mu \nu }.`$
To prove (12), expand the EM component of tensor (7) using (10), and simplify using identities (A1)-(A4), followed by the velocity normalization condition, and orthogonality condition (11):

$`T_{(EM)}^{\mu \nu }`$ $`=`$ $`{\displaystyle \frac{1}{4\pi }}\left(g^{\mu \gamma }F_{\gamma \alpha }F^{\nu \alpha }{\displaystyle \frac{1}{4}}F_{\alpha \beta }F^{\alpha \beta }g^{\mu \nu }\right)`$
$`=`$ $`{\displaystyle \frac{1}{4\pi }}\left[g^{\mu \gamma }\left(ϵ_{\delta ϵ\gamma \alpha }B^\delta U^ϵ\right)\left(ϵ^{\rho \sigma \nu \alpha }B_\rho U_\sigma \right){\displaystyle \frac{1}{4}}\left(ϵ_{\delta ϵ\alpha \beta }B^\delta U^ϵ\right)\left(ϵ^{\rho \sigma \alpha \beta }B_\rho U_\sigma \right)g^{\mu \nu }\right]`$
$`=`$ $`{\displaystyle \frac{1}{4\pi }}\left(g^{\mu \gamma }\delta _{\delta ϵ\gamma }^{\rho \sigma \nu }B^\delta U^ϵB_\rho U_\sigma {\displaystyle \frac{g^{\mu \nu }}{2}}\delta _{\delta ϵ}^{\rho \sigma }B^\delta U^ϵB_\rho U_\sigma \right)`$
$`=`$ $`{\displaystyle \frac{1}{4\pi }}\left({\displaystyle \frac{g^{\mu \nu }B^\alpha B_\alpha }{2}}+B^\mu B^\nu B^\alpha B_\alpha U^\mu U^\nu \right)`$
$`=`$ $`{\displaystyle \frac{1}{2}}g^{\mu \nu }b^2+U^\mu U^\nu b^2b^\mu b^\nu .`$
To prove (15), the second of Maxwell’s equations, $`_\mu {}_{}{}^{}F_{}^{\mu \nu }=0`$, can be rewritten by substituting definitions, and using identities (A2) and (A3):

$`_\alpha {}_{}{}^{}F_{}^{\alpha \beta }`$ $`=`$ $`{\displaystyle \frac{1}{2}}_\alpha ϵ^{\alpha \beta \delta \gamma }F_{\delta \gamma }`$
$`=`$ $`{\displaystyle \frac{1}{2}}_\alpha ϵ^{\alpha \beta \delta \gamma }ϵ_{\rho \sigma \delta \gamma }b^\rho U^\sigma `$
$`=`$ $`_\alpha \delta _{\rho \sigma }^{\alpha \beta }b^\rho U^\sigma ,`$
which simplifies, using (A3), to $`_\alpha \left(U^\alpha b^\beta b^\alpha U^\beta \right)=0`$.

To prove (2.2), we begin with an identity derived from the expression for $`W`$ (18),

$$W^2\alpha ^2=g^{tt}W^2=(g_{tt}+2g_{t\varphi }V^\varphi +g_{rr}(V^r)^2+g_{\theta \theta }(V^\theta )^2+g_{\varphi \varphi }(V^\varphi )^2)^1,$$
(A5)
or

$$g_{tt}+g_{rr}(V^r)^2+g_{\theta \theta }(V^\theta )^2+g_{\varphi \varphi }(V^\varphi )^2=\frac{\alpha ^2}{W^2}2g_{t\varphi }V^\varphi .$$
(A6)
Now expand $`b^2`$,

$$b^2=g_{tt}(b^t)^2+g_{rr}(b^r)^2+g_{\theta \theta }(b^\theta )^2+g_{\varphi \varphi }(b^\varphi )^2+2g_{t\varphi }b^tb^\varphi ,$$
(A7)
and substitute

$$b^i=\frac{^i}{\sqrt{4\pi }W\sqrt{\gamma }}+b^tV^i$$
(A8)
for the spatial components of the magnetic field. Making use of the above identity, it follows that

$`b^2`$ $`=`$ $`\left[{\displaystyle \frac{\alpha ^2}{W^2}}+2g_{t\varphi }V^\varphi \right](b^t)^2+2g_{t\varphi }\left[{\displaystyle \frac{^\varphi }{\sqrt{4\pi }W\sqrt{\gamma }}}+b^tV^\varphi \right]b^t`$
$`+`$ $`{\displaystyle \frac{2}{\sqrt{4\pi }W\sqrt{\gamma }}}\left[g_{rr}^rV^r+g_{\theta \theta }^\theta V^\theta +g_{\varphi \varphi }^\varphi V^\varphi \right]b^t`$
$`+`$ $`{\displaystyle \frac{1}{4\pi W^2\gamma }}\left[g_{rr}(^r)^2+g_{\theta \theta }(^\theta )^2+g_{\varphi \varphi }(^\varphi )^2\right].`$
Now substitute the expression for $`b^t`$ (37) and simplify, noting that

$$g_{\varphi \varphi }V^\varphi +g_{t\varphi }=\frac{g^{tt}V^\varphi g^{t\varphi }}{g^{tt}g^{\varphi \varphi }(g^{t\varphi })^2}$$
(A10)
to get, after a few steps, expression (2.2).

To derive (51) we first obtain the components of the fluid-frame magnetic field, which are, to leading order,

$`b^t`$ $``$ $`g_{xx}{\displaystyle \frac{WV^x^x}{\alpha ^2\sqrt{4\pi \gamma }}},`$
$`b^x`$ $``$ $`{\displaystyle \frac{W^x}{\sqrt{4\pi \gamma }}},`$
$`b^y`$ $``$ $`{\displaystyle \frac{^y}{\sqrt{4\pi \gamma }W}}+g_{xx}{\displaystyle \frac{WV^xV^y^x}{\alpha ^2\sqrt{4\pi \gamma }}},`$
$`b^2`$ $``$ $`g_{xx}\left({\displaystyle \frac{^x}{\sqrt{4\pi \gamma }}}\right)^2,`$
and also note the approximate value for $`W`$,

$$W^21g_{xx}(V^x)^2\alpha ^2.$$
(A11)
We then expand the momentum equation and substitute the expressions for $`b_y`$ and $`b^t`$. The left-hand side of the momentum equation reads

$`_t\left(S_y\alpha b_yb^t\right)`$ $`=`$ $`_t\left(g_{yy}(\rho h+b^2)W^2\alpha ^1V^yg_{yy}\alpha b_yb^t\right)`$
$``$ $`_t\left(g_{yy}\rho h{\displaystyle \frac{W^2}{\alpha }}V^y+g_{yy}b^2{\displaystyle \frac{W^2}{\alpha }}V^y\left[1{\displaystyle \frac{g_{xx}(V^x)^2}{\alpha ^2}}\right]{\displaystyle \frac{g_{xx}g_{yy}^x^yV^x}{\alpha \left(\sqrt{4\pi \gamma }\right)^2}}\right)`$
$``$ $`{\displaystyle \frac{g_{yy}(\rho hW^2+b^2)}{\alpha }}_t\left(V^y\right){\displaystyle \frac{g_{xx}g_{yy}^xV^x}{\alpha \left(\sqrt{4\pi \gamma }\right)^2}}_t\left(^y\right).`$
Recall that $`x`$ stands for either $`r`$ or $`\theta `$, and $`y`$ stands for a transverse direction, either $`r`$, $`\theta `$, or $`\varphi `$.

As for the left hand side, we neglect the purely geometric terms and retain only

$`{\displaystyle \frac{1}{\sqrt{\gamma }}}_x\sqrt{\gamma }\left(S_yV^x\alpha b_yb^x\right)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{\gamma }}}_x\sqrt{\gamma }\left(g_{yy}(\rho h+b^2){\displaystyle \frac{W^2}{\alpha }}V^yV^x\alpha g_{yy}b^yb^x\right)`$
$``$ $`{\displaystyle \frac{1}{\sqrt{\gamma }}}_x\sqrt{\gamma }\left(g_{yy}\rho h{\displaystyle \frac{W^2}{\alpha }}V^yV^x\alpha g_{yy}{\displaystyle \frac{^x^y}{\left(\sqrt{4\pi \gamma }\right)^2}}\right)`$
$``$ $`{\displaystyle \frac{g_{yy}\rho hW^2V^x}{\alpha }}_x\left(V^y\right)+{\displaystyle \frac{\alpha g_{yy}^x}{\left(\sqrt{4\pi \gamma }\right)^2}}_x\left(^y\right).`$
Combine these results and replace $`_t\left(^y\right)`$ with the results for the induction equation,

$`_t\left(V^y\right)+{\displaystyle \frac{\rho hW^2V^x}{(\rho hW^2+b^2)}}_x\left(V^y\right)`$ $``$ $`{\displaystyle \frac{g_{xx}^xV^x}{4\pi \gamma (\rho hW^2+b^2)}}_t\left(^y\right)`$
$`+`$ $`{\displaystyle \frac{\alpha ^2^x}{4\pi \gamma (\rho hW^2+b^2)}}_x\left(^y\right)`$
$``$ $`{\displaystyle \frac{g_{xx}^xV^x}{4\pi \gamma (\rho hW^2+b^2)}}\left(^x_x\left(V^y\right)V^x_x\left(^y\right)\right)`$
$`+`$ $`{\displaystyle \frac{\alpha ^2^x}{4\pi \gamma (\rho hW^2+b^2)}}_x\left(^y\right).`$
By regrouping these terms and substituting the approximate expressions for $`W`$ and $`b^2`$ given above, equation (51) follows immediately. A similar procedure yields expression (56) for transverse velocities when $`x\varphi `$.

To prove that the Bondi solution is unchanged by the presence of a radial magnetic field, we revisit the calculations presented in HSWa. We are looking for a time-independent purely radial solution to the equations of GRMHD in the presence of a constant CT radial magnetic field $`^r=F_{\theta \varphi }=B`$ in the Schwarzschild metric. We will work with the more general formalism of equations (2) and (3).

The equation of continuity (2) is not dependent on magnetic fields, so

$$_\mu \left(\rho U^\mu \right)=0_r\left(\rho U^r\right)=0,$$
(A12)
leads to the same conserved quantity as in the hydrodynamic case,

$$\sqrt{g}\rho U^r=C_1.$$
(A13)
The energy-flux conservation equation is the $`t`$-component of (3),

$$_\mu T_{}^{\mu }{}_{t}{}^{}=0_rT_{}^{r}{}_{t}{}^{}=0,$$
(A14)
and contains an explicit dependence on the magnetic field. However, by expanding the $`rt`$ component of the energy momentum tensor,

$$T_{}^{r}{}_{t}{}^{}=\left(\rho h+b^2\right)U^rU_tb^rb_t,$$
(A15)
and using equations (32) and (37) to obtain the fluid-frame magnetic field,

$$b^r=\frac{^r}{\sqrt{4\pi }}\frac{W}{\sqrt{\gamma }},b^t=g_{rr}b^rV^r,b^2=g_{rr}\frac{(b^r)^2}{W^2},$$
(A16)
we obtain after simple substitutions

$$T_{}^{r}{}_{t}{}^{}=\left(\rho h+b^2\right)U^rU_tb^2U^rU_t.$$
(A17)
Clearly, the magnetic contributions cancel (this is simply restating the general result that there is no net force due to magnetic fields parallel to the direction of motion). We therefore recover the hydrodynamic expression,

$$_rT_{}^{r}{}_{t}{}^{}=\left(\rho h\right)U^rU_t=0,$$
(A18)
which leads to the same conserved quantities as in the hydrodynamic case,

$$\sqrt{g}\rho U^rhU_t=C_2,$$
(A19)
and

$$hU_t=C_3.$$
(A20)
## Appendix B Amplitude of Alfvén Pulses

To prove (61), we used the method of characteristics. We write the momentum and induction equations as a matrix equation,

$$A_xU+B_tU=0$$
(B1)
where $`U(x,t)=(V^y(x,t),^y(x,t))^T`$ and

$$A=\left[\begin{array}{cc}\frac{V^x}{(1+\eta ^2)}& \frac{\eta ^2}{^x(1+\eta ^2)}\\ ^x& V^x\end{array}\right],B=\left[\begin{array}{cc}1& \chi \\ 0& 1\end{array}\right],\chi =\frac{\eta ^2V^x}{^x(1+\eta ^2)}.$$
(B2)
subject to the initial state $`U_0(x,0)=(f(x),0)^T`$ We obtain the eigenvalues from the roots of $`F(\lambda )=detA\lambda B`$, which are simply the expressions for $`v_{A}^{}{}_{}{}^{(\pm )}`$, (59).

To solve this system, we seek a matrix $`P`$ such that $`PBU=z`$ where $`z_i(x,t)=g_1(x\lambda _it)`$ and $`g(x,0)=PBU_0(x,0)`$. The rows of matrix $`P`$ satisfy $`A^T\lambda _iB^TP_i=0`$. It is straightforward to show that the matrix $`P`$ has the form,

$$P=\left[\begin{array}{cc}\zeta & 1\\ \zeta & 1\end{array}\right],\zeta =\frac{^x(1+\eta ^2)}{\eta \sqrt{\eta ^2+W^2}}.$$
(B3)
We obtain the results

$`\zeta V^y+(1+\zeta \chi )^y`$ $`=`$ $`\zeta f(xv_{A}^{}{}_{}{}^{}t),`$ (B4)
$`\zeta V^y+(1\zeta \chi )^y`$ $`=`$ $`\zeta f(xv_{A}^{}{}_{}{}^{+}t),`$ (B5)
which yield

$`V^y(x,t)`$ $`=`$ $`\left({\displaystyle \frac{1\zeta \chi }{2}}\right)f(xv_{A}^{}{}_{}{}^{}t)+\left({\displaystyle \frac{1+\zeta \chi }{2}}\right)f(xv_{A}^{}{}_{}{}^{+}t),`$ (B6)
$`^y(x,t)`$ $`=`$ $`{\displaystyle \frac{\zeta }{2}}\left(f(xv_{A}^{}{}_{}{}^{}t)f(xv_{A}^{}{}_{}{}^{+}t)\right).`$ (B7)
As discussed in the text, the velocity amplitudes are clearly asymmetrical, except for the case where $`V^x=0`$, and the magnetic field amplitudes are always equal, but opposite in sign. It is easy to show that the ratio of the velocity amplitudes is indeed the ratio of the gamma factors, $`W(v_{A}^{}{}_{}{}^{(\pm )})`$. It is also possible to verify by direct computation that the amplitudes of the pulses given in the text agree with the above formulas.

## Appendix C MHD Shocks

Relativistic MHD shocks have been discussed by Anile (1989) and Komissarov (1999). Their work builds upon the monograph of Lichnerowicz (1967) who systematically developed the relations for GR MHD shock invariants. We summarize here the relevant results from Lichnerowicz, translating symbols to our own notation. We use these results to construct our shock tests.

MHD shocks are discontinuities in one or more of the MHD variables whose propagation through spacetime lies on a hypersurface $`\mathrm{\Sigma }`$. It is useful in the analysis of these shocks to construct a $`4`$-vector $`n^\mu `$ normal to $`\mathrm{\Sigma }`$. Shocks can be characterized by invariant scalar and $`4`$-vector quantities—the jump conditions—which are derived readily from the equations of GRMHD. MHD shocks can be assigned to one of two categories, tangential and non-tangential, and this classification hinges on whether the invariants have non-zero projections onto the normal $`n^\mu `$. Slow and fast magnetosonic shocks are examples of non-tangential shocks, so it is this type of shock that is of interest here.

The analysis of non-tangential shocks was shown by Lichnerowicz to reduce to the study of five scalar quantities, $`\rho `$, $`h`$, $`\eta =b^\mu n_\mu `$, $`𝒜=\rho U^\mu n_\mu `$, and $`b^2`$ which are related through five positive-definite scalar invariants (derivable from the jump conditions),

$`\left[𝒜\right]`$ $`=`$ $`0,`$ (C1)
$`\left[\right]`$ $`=`$ $`\left[h\eta \right]=0,`$ (C2)
$`\left[\right]`$ $`=`$ $`\left[{\displaystyle \frac{\eta ^2}{𝒜^2}}{\displaystyle \frac{b^2}{\rho ^2}}\right]=0,`$ (C3)
$`\left[\right]`$ $`=`$ $`\left[\alpha +{\displaystyle \frac{q}{𝒜^2}}\right]=0,`$ (C4)
$`\left[𝒦\right]`$ $`=`$ $`\left[\left(\rho ^2+𝒜^2\right){\displaystyle \frac{h^2}{\rho ^2}}+k^2\left(2{\displaystyle \frac{h}{\rho }}\right)\right]=0,`$ (C5)
where $`\left[𝒳\right]=𝒳_+𝒳_{}`$ indicates the change in a quantity $`𝒳`$ across the shock, and

$`q`$ $`=`$ $`P+{\displaystyle \frac{b^2}{2}},`$ (C6)
$`k^2`$ $`=`$ $`b^2𝒜^2,`$ (C7)
$`\alpha `$ $`=`$ $`{\displaystyle \frac{h}{\rho }}.`$ (C8)
Magnetosonic shocks are characterized by the parameter $`\alpha `$. The determination as to whether a magnetosonic wave is fast or slow hinges on the sign of the pre- and post-shock values of parameter $`\alpha `$. Physically realizable shocks are characterized by $`\alpha _{(+)}\alpha _{()}>0`$. A slow shock has $`\alpha _{(+)}<\alpha _{()}<0`$ and a fast shock has $`0<\alpha _{(+)}<\alpha _{()}`$.

To construct an initial state, we solve the scalar invariants using the Mathematica FindRoot function. First, we postulate a left state in the rest frame of the shock, where the normal vector $`n^\mu `$ has an especially simple expression ($`n^\mu =(0,1,0,0)`$). Next, we set up the scalar invariants as a system to be solved by the FindRoot function, supplying an initial guess for the right state. This is an iterative procedure, requiring adjustments to the initial guess. Once a solution is found, it is possible to construct a moving shock by repeating the above process with a small non-zero shock speed, for which the normal vector now has the form $`n^\mu =W(v_{sh})(v_{sh},1,0,0)`$. Once a solution is found the process is repeated, gradually increasing the shock speed and using the previous converged right state as a starting point for the solution at the new shock speed. This is a pragmatic solution to a problem that has a high degree of complexity, and it allows us to set up a family of closely related shock solutions to aid in testing.

Given that the FindRoot function has converged on a right state, it is possible to compute the parameter $`\alpha `$, and classify the shock. It is also prudent to verify that the left and right states are true solutions by verifying that they satisfy the shock invariants by direct computation. Once this is done, the left and right states can be transferred directly to a numerical grid.

## Appendix D Determining Conserved Quantities for Magnetized Gammie Inflow

The magnetized Gammie inflow is characterized by three conserved quantities,

$`F_M`$ $`=`$ $`2\pi r^2U^r\rho ,`$ (D1)
$`F_L`$ $`=`$ $`2\pi r^2T_{}^{r}{}_{\varphi }{}^{}=2\pi r^2\left[\rho U^rU_\varphi {\displaystyle \frac{g^{rr}g^{\theta \theta }}{4\pi }}F_{r\theta }F_{\theta \varphi }\right],`$ (D2)
$`F_E`$ $`=`$ $`2\pi r^2T_{}^{r}{}_{t}{}^{}=2\pi r^2\left[\rho U^rU_t+{\displaystyle \frac{g^{rr}g^{\theta \theta }}{4\pi }}F_{r\theta }F_{t\theta }\right],`$ (D3)
which are readily derived from the equations of baryon and momentum conservation. Note that Gammie writes $`𝒟r^2g^{rr}g^{\theta \theta }`$. From the induction equation,

$`_r^r`$ $`=`$ $`0,`$ (D4)
$`_r\left(V^\varphi ^rV^r^\varphi \right)`$ $`=`$ $`0,`$
we obtain $`^r=F_{\theta \varphi }=const`$, and $`F_{t\theta }=V^\varphi ^rV^r^\varphi =\mathrm{\Omega }_FF_{\theta \varphi }`$ which follows from the boundary condition. Infinite conductivity ($`U^\mu F_{\mu \nu }=0`$) yields

$$F_{r\theta }=\frac{F_{\theta \varphi }}{U^r}\left(U^\varphi \mathrm{\Omega }_FU^t\right).$$
(D5)
Finally, the velocity normalization condition yields

$$1=g^{tt}(U_t)^2+2g^{t\varphi }U_tU_\varphi +g^{\varphi \varphi }(U_\varphi )^2+(U^r)^2/g^{rr}$$
(D6)
where the mix of upper and lower indices was made to coincide with code variables.

The critical point is a saddle point in the two-dimensional parameter space of the energy function $`F_EF_E(r,U^r;F_L)`$. It is apparent that $`F_E`$ as given above is not a simple function of these parameters, and it is necessary to substitute the other expressions in order to obtain the desired form for $`F_E`$. Since the resulting expression for $`F_E`$ is quite lengthy, we will only sketch the procedure. The process of substitution and elimination described here was carried out using Mathematica.

First, the conserved quantity $`F_M`$ is used to obtain $`\rho (r,U^r)=F_M/(2\pi r^2U^r)`$. Next, the conserved quantity $`F_L`$ is used to isolate $`U_\varphi `$,

$$U_\varphi =\left[\frac{F_L}{2\pi r^2U^r}+\frac{𝒟}{4\pi r^2}F_{\theta \varphi }F_{r\theta }\right]\frac{1}{\rho U^r}.$$
(D7)
Now substitute the expressions for $`\rho (r,U^r)`$ and $`F_{r\theta }`$ expressed in terms of $`U_\varphi `$ and $`U_t`$, to yield an expression of the form $`U_\varphi =A+BU_t`$ where $`A`$ and $`B`$ are expressions involving the constants $`F_L`$, $`F_M`$, $`F_{\theta \varphi }`$, and $`\mathrm{\Omega }_F`$, as well as $`U^r`$ and metric terms. Now substitute this expression for $`U_\varphi `$ in the velocity normalization condition and “solve” the resulting quadratic for $`U_t`$, denote this solution $`(U_t)^\pm `$ to emphasize that there are two distinct roots. Finally, take the conserved quantity $`F_E`$ and substitute the expression for $`F_{r\theta }`$ (expressed in terms of $`U_\varphi `$ and $`U_t`$) and then perform a final substitution for $`U_\varphi `$ and $`U_t`$ by $`(U_t)^\pm `$. Simplify the resulting expression, which has the desired functional form: $`F_E(r,U^r;F_L)`$.

The critical point must satisfy the conditions

$`_rF_E(r_{crit},U_{crit}^r;F_L)`$ $`=`$ $`0,`$
$`_{U^r}F_E(r_{crit},U_{crit}^r;F_L)`$ $`=`$ $`0,`$
and the conserved quantity $`F_E`$ must also satisfy the constraint at the outer boundary

$`F_E(r_{crit},U_{crit}^r;F_L)`$ $`=`$ $`F_E(r_{mso},0;F_L).`$
These three conditions allow us to solve for $`r_{crit}`$, $`U_{crit}^r`$, and $`F_L`$.

In order to find the critical point, we will eventually use the FindRoot function of Mathematica, but we first need to obtain a rough location for the saddle point. To do this, we produce a coarse contour plot over the range $`r_{horizon}<r<r_{mso}`$ and a large range of $`U^r`$, say $`0.5<U^r<0`$ for a given choice of root of $`(U_t)^\pm `$. If the contour plot does not show evidence of a saddle point, we choose the other root. Figure 14 illustrates this for our specific test problem. The saddle point is shown to lie near $`(r_{crit},U_{crit}^r)(3.65,0.04)`$. This is then used as an initial guess for the Mathematica FindRoot function, from which the parameters used in the test were obtained.

To generate an initial state within our code, we apply a numerical root-finder to the energy function $`F_E(r,U^r;F_L)`$ (coded using the Mathematica FortranForm command), and generate a numerical solution for $`U^r`$ (using the parameters in the text). Once $`U^r`$ is found, we construct $`U_t`$ and $`U_\varphi `$ using the algebraic expressions obtained from Mathematica (also coded using FortranForm). Once these variables have been initialized the remainder of the code variables can be obtained in a straightforward manner.


format: text

text:
# Comparison of Analytical Methods of 𝐸1 Strength Calculations in Middle and Heavy Nuclei

## 1 Introduction

Gamma-emission is one of the most universal channel of the nuclear deexcitation processes which can accompany any nuclear reaction. Gamma-decay as well as photoabsorption can be described through the use of the radiative strength (RS) functions . The electron- positron decay is also depended on shape of these functions. There are radiative strength functions of two types. The photoexcitation strength function is connected with cross-section of the gamma-ray absorption. The average radiative width $`\mathrm{\Gamma }_\lambda `$ is determined by a gamma- decay (downward) strength function $`\stackrel{}{f}_\lambda `$, $`\stackrel{}{f}_\lambda \mathrm{\Gamma }_\lambda (ϵ_\gamma )\rho (U)/3ϵ_\gamma ^3\rho (Uϵ_\gamma )`$, where $`\rho (U)`$ is total density of the initial states in heated nuclei at initial excitation energy $`U`$ and $`ϵ_\gamma `$ is the gamma-ray energy. Therefore the RS function is important constituent of the compound nucleus model calculations of capture cross sections, gamma-ray production spectra, isomeric state populations and of competition between gamma-ray and particle emission. The gamma-ray strengths are ingredient of time-consuming calculations of different nuclear processes, for these reasons simple closed-form expressions are preferred for strength functions. The approaches based on recent achievements in nuclear studies are also preferable to improve the reliability of the RS function expressions.

Dipole electric gamma-transitions are dominant, when they occur simultaneously with transitions of other multipolarities and types. According to the Brink hypothesis , the Lorentzian line shape with the energy- independent width (SLO model) is used widely for calculations of the dipole RS function. This approach is probably the most appropriate simple method for a description of the photoabsorption data in medium-weight and heavy nuclei . In the case of the gamma-emission, the SLO model strongly underestimates  the gamma-ray spectra at low energies $`ϵ_\gamma <1MeV`$. A global description of the gamma-spectra by the Lorentzian can be obtained in the range $`1<ϵ_\gamma <8MeV`$ only when the giant dipole resonance (GDR) parameters are inconsistent with that ones for photoabsorption data. On the whole, the SLO approach overestimates the integral experimental data like the capture cross sections and the average radiative widths in heavy nuclei (-).

The first model with correct description of the $`E1`$ strengths at the energies $`ϵ_\gamma `$ near zero was proposed in Ref.. Thereafter an enhanced generalized Lorentzian model (EGLO) was developed and analyzed in Refs. for a unified description of the low-energetic and integral data. The EGLO radiative strength function consists of two components for spherical nuclei: the Lorentzian with the energy and temperature dependent empirical width and additional term  corresponding to zero value of $`\gamma `$ \- ray energy. The EGLO method reproduces the experimental data on gamma-decay rather well in the mass region $`A=50÷200`$. It should be noted that the EGLO and SLO expressions for the gamma-decay strength function of heated nuclei are in fact the parameterizations of the experimental data. They are in contradiction with some aspect of the recent theoretical studies, specifically:

1) the shapes of the radiative strengths within these approaches do not consistent with general relation between the RS function for gamma-decay of heated nuclei and the imaginary part of the nuclear response function on electromagnetic field ;

2) the damping width of the EGLO model is proportional to collisional component of zero sound damping width in the infinite Fermi- liquid when the collisional (two-body) relaxation is taken into account only. However the important contribution to the total width in nuclei is also given by the fragmentation (one-body) width arising from nucleon motion in self-consistent mean field . This width is almost independent of the nuclear temperature and is not included in EGLO model. Note that the SLO-model width is independent of energy; that is, it has properties of the fragmentation width but it is not allowing for collisional damping.

Recently new closed-form models for the E1 strength were proposed in Refs. (-) which avoid some of these defects at least in approximate way.

Specifically, the approach considered in Refs.(-) is in line with detailed balance principle . This method was named previously as the thermodynamic pole approximation. It is denoted below as the modified Lorentzian (MLO) approach, because a resulting expression for the dipole RS function has a Lorentzian shape scaled by an enhanced factor. An expression for energy-dependent damping width of the MLO model includes also component corresponding to fragmentation contribution in some simplified way. Damping width appearing in the generalized Fermi liquid (GFL) model  of the E1 strength has fragmentation component too; it is taken as resulted from the dipole-quadrupole interaction. Below we investigate these new approaches to compare their with EGLO and SLO models.

## 2 Simple closed-form models for radiative <br>strength function calculations

The expression for dipole radiative strength function has the following form in SLO model

$$\stackrel{}{f}_{SLO}(ϵ_\gamma )=8,67410^8\sigma _r\mathrm{\Gamma }_r\frac{ϵ_\gamma \mathrm{\Gamma }_r}{(ϵ_\gamma ^2E_r^2)^2+(\mathrm{\Gamma }_rϵ_\gamma )^2},MeV^3,$$
(1)
where $`\sigma _r`$ (in mb), $`\mathrm{\Gamma }_r`$ and $`E_r`$ (in MeV) are GDR parameters.

The expression for gamma-decay dipole RS function within EGLO model  is the following,

$$\stackrel{}{f}_{EGLO}(ϵ_\gamma )=8.67410^8\sigma _r\mathrm{\Gamma }_r\times $$

$$\times [\frac{ϵ_\gamma \mathrm{\Gamma }_k(ϵ_\gamma ,T_f)}{(ϵ_\gamma ^2E_r^2)^2+(ϵ_\gamma \mathrm{\Gamma }_k(ϵ_\gamma ,T_f))^2}+0.7\mathrm{\Gamma }_k(ϵ_\gamma =\frac{0,T_f)}{E_r^3}],$$
(2)
where $`T_f`$ is the temperature of the final state. The energy-dependent width $`\mathrm{\Gamma }_k(ϵ_\gamma ,T_f)`$ is taken proportionally to the collisional damping width in Fermi-liquid scaled by an empirical  function $`K(ϵ_\gamma )`$:

$`\mathrm{\Gamma }_k(ϵ_\gamma ,T)=K(ϵ_\gamma )\mathrm{\Gamma }_r\left[ϵ_\gamma ^2+(2\pi T)^2\right]/E_r^2,`$ (3)
$`K(ϵ_\gamma )=\kappa +(1\kappa )(ϵ_\gamma ϵ_0)/(E_rϵ_0).`$
Dipole strength of the GFL model , $`\stackrel{}{f}_{E1}\stackrel{}{f}_{GFL}`$, can be given in the following form for spherical nuclei

$`\stackrel{}{f}_{GFL}(ϵ_\gamma )`$ $`=`$ $`8.67410^8\sigma _r\mathrm{\Gamma }_r{\displaystyle \frac{K_{GFL}ϵ_\gamma \mathrm{\Gamma }_m(ϵ_\gamma ,T_f)}{(ϵ_\gamma ^2E_r^2)^2+K_{GFL}(\mathrm{\Gamma }_m(ϵ_\gamma ,T_f)ϵ_\gamma )^2}},`$ (4)
$`K_{GFL}`$ $`=`$ $`\sqrt{E_r/E_0}=(1+F_1^{}/3)^{1/2}/(1+F_0^{})^{1/2}=0.63,`$
where $`F_0^{}`$ and $`F_1^{}`$ are the Landau parameters of the quasiparticle interaction in isovector channel of the Fermi system; $`F_0^{}=1.49`$ and $`F_1^{}=0.04`$; $`E_0`$ is the average energy of one-particle one-hole states forming GDR. The Eq.(4) is an extension of original expression of the GFL model  for wide range of the gamma-ray energies: term $`K_{GFL}(\mathrm{\Gamma }_mϵ_\gamma )^2`$ is added to the denominator of the Eq.(4) to avoid singularity of the GFL approach near GDR-energy in a way similar to the other models for E1 strength but with additional factor $`K_{GFL}`$ in order to keep standard relationship between value of the RS function at the GDR energy and peak value $`\sigma _r`$ of the photoabsorption cross-section.

The width $`\mathrm{\Gamma }_m`$ in (4) is taken as a sum of a collisional damping width, $`\mathrm{\Gamma }_C`$, and a term, $`\mathrm{\Gamma }_{dq}`$, which simulates the fragmentation width,

$$\mathrm{\Gamma }_m(ϵ_\gamma ,T)=\mathrm{\Gamma }_C(ϵ_\gamma ,T)+\mathrm{\Gamma }_{dq}(ϵ_\gamma ),\mathrm{\Gamma }_CC_{coll}\left(ϵ_\gamma ^2+4\pi ^2T^2\right),$$
(5)
$`C_{coll}=16m\sigma (np)/4\pi ^29\mathrm{}^2`$. Here, a magnitude of the in-medium cross section $`\sigma (np)`$ is taken proportional to the value of the free space cross section $`\sigma _f(np)=5fm^2`$ (near Fermi surface) with a factor $`F`$. The form of collisional component corresponds to damping width in infinite the Fermi-liquid and $`\mathrm{\Gamma }_{dq}`$ is considered as resulted from spreading GDR over surface quadrupole vibrations due to dipole-quadrupole interaction:

$$\mathrm{\Gamma }_{dq}(ϵ_\gamma )=C_{dq}ϵ_\gamma \overline{\beta }_2\sqrt{1+\frac{E_2}{ϵ_\gamma }}=C_{dq}\sqrt{ϵ_\gamma ^2\overline{\beta }_2^2+ϵ_\gamma s_2},s_2=E_2\overline{\beta }_2^2.$$
(6)
Here, $`C_{dq}=(5\mathrm{ln}2/\pi )^{1/2}=1.05`$; $`E_2`$ is energy of the first excited vibrational $`2^+`$ state; $`\overline{\beta }_2`$ is effective deformation parameter of nuclear surface. It is determined by reduced electric photoabsorption transition rate for the ground state to $`2^+`$ state transition . The global parameterization  is adopted for $`s_2`$: $`s_2E_2\overline{\beta }_2^2=217.16/A^2`$ , when experimental data for even-even nuclei and in the case of odd and odd-odd nuclei.

It should be noted that there is not consistency of GFL model with detailed balance principle in systems with constant temperature and the collisional components of the damping width of the GFL model can have negative values in some deformed nuclei if magnitude $`C_{coll}`$ in Eq.(5) is found from fitting width $`\mathrm{\Gamma }_m`$ to experimental value of GDR width.

Expression for dipole RS within MLO model, $`\stackrel{}{f}_{E1}\stackrel{}{f}_{MLO}`$, is obtained by calculating the average radiative width in nuclei with microcanonically distributed initial states (-,)). In spherical nuclei the dipole radiative strength function is proportional to the strength function for nuclear response on dipole field with frequency $`\omega =ϵ_\gamma /\mathrm{}`$. An analytical semiclassical expression for nuclear response function in cold and heated nuclei with excitation of GDR and theory-supported expressions for damping are used. As the result the radiative strength function within MLO model has the following general form (in MeV)

$$\stackrel{}{f}_{MLO}(ϵ_\gamma ,T_f)=8,67410^8\sigma _r\mathrm{\Gamma }_rL(ϵ_\gamma ,T_f)\frac{ϵ_\gamma \mathrm{\Gamma }(ϵ_\gamma ,T_f)}{(ϵ_\gamma ^2E_r^2)^2+(\mathrm{\Gamma }(ϵ_\gamma ,T_f)ϵ_\gamma )^2},$$
(7)
where scaling factor $`L(ϵ_\gamma ,T_f)=1/[1exp(ϵ_\gamma /T_f)]`$ is essential for low-energy radiation and leads to consistency of the Eq.(7) with detailed balance principle in systems with constant temperature .

Three variants of the modified Lorentzian model (MLO1, MLO2 and MLO3) are considered below. They include different expressions for damping width $`\mathrm{\Gamma }(ϵ_\gamma ,T_f)`$. The expression for MLO1 approach has the following form (-)

$$\mathrm{\Gamma }_{MLO1}(ϵ_\gamma ,T_f)\beta \gamma _c(ϵ_\gamma ,T_f)\frac{E_r^2+E_0^2}{(E_r^2E_0^2)^2+(\gamma _c(ϵ_\gamma ,T_f)ϵ_\gamma )^2},$$
(8)
$$\beta =\left(1+\frac{E_r^2}{E_0^2}\right)^2\frac{E_0^2}{2},\gamma _c(ϵ_\gamma ,T)=\frac{2\mathrm{}}{\tau _c(ϵ_\gamma ,T)}.$$

The energy $`E_0`$ of particle-hole state is taken as equal the harmonic oscillator energy $`\mathrm{}\omega _0=41/A^{1/3},MeV`$. In this model the relaxation time within doorway state mechanism  is used

$$\frac{\mathrm{}}{\tau _c(ϵ_\gamma ,T)}=b\left(ϵ_\gamma +U\right),b=\frac{E_r}{4\pi }\frac{F}{\alpha },\alpha =\frac{9\mathrm{}/16m}{\sigma ^{free}(np)},F=\frac{\sigma (np)}{\sigma ^{free}(np)}.$$
(9)
The damping widths in MLO2 and MLO3 models are taken in approximation of independent sources of dissipation  as a sum of the collisional damping width $`\mathrm{\Gamma }_C`$ and a term which simulate the fragmentation component of the width: $`\mathrm{\Gamma }(ϵ_\gamma ,T_f)=\mathrm{\Gamma }_C(ϵ_\gamma ,T_f)+\mathrm{\Gamma }_F(ϵ_\gamma )`$. The component $`\mathrm{\Gamma }_C`$ is inversely proportional to the collision relaxation time $`\tau `$ in isovector channel at dipole distortion of the Fermi surface . The magnitude of $`\mathrm{\Gamma }_F`$ is taken proportionally to the wall formula value  $`\mathrm{\Gamma }_w`$ with a scaling factor $`k_s`$: $`\mathrm{\Gamma }_F(ϵ_\gamma )=k_s(ϵ_\gamma )\mathrm{\Gamma }_w,\mathrm{\Gamma }_w=36.43A^{1/3}(MeV)`$. As a result, the expression for the damping width has the following form,

$$\mathrm{\Gamma }_{MLO2,3}(ϵ_\gamma ,T_f)=\mathrm{}/\tau _c(ϵ_\gamma ,T_f)+k_s\mathrm{\Gamma }_w.$$
(10)
The energy-dependent power approximation is adopted for simplicity for factor $`k_s`$ : $`k_s(ϵ_\gamma )=k_r+(k_0k_r)|(ϵ_\gamma E_r)/E_r|^{n_s}`$, when $`ϵ_\gamma <2E_r`$, and $`k_s(ϵ_\gamma )=k_0`$ if $`ϵ_\gamma 2E_r`$, where the quantities $`k_0k_s(ϵ_\gamma =0)`$, $`k_rk_s(ϵ_\gamma =E_r)`$ determine the contribution of the ”wall” component to the width at zero energy and GDR-energy, respectively. The value of the $`k_r`$ is obtained from fitting the GDR width $`\mathrm{\Gamma }_r`$ at zero temperature by the expression (10) with $`ϵ_\gamma =E_r`$. The quantities $`k_0`$ and $`n_s`$ are some parameters which are obtained below from condition of correct description of a general behaviour of the experimental gamma-decay strengths. The collisional damping width $`\mathrm{\Gamma }_C`$ is taken in form given by Eq.(9) in the case of MLO2 model. The expression for relaxation time according to the Fermi-liquid approach is used for MLO3 model:

$$\frac{\mathrm{}}{\tau _c(ϵ_\gamma ,T)}\frac{F}{\alpha }\left[\left(ϵ_\gamma /2\pi \right)^2+T^2\right].$$
(11)
## 3 Analysis of the dipole strength function calculations

The RS functions were calculated on following nuclei: $`{}_{}{}^{198}Au`$, $`{}_{}{}^{168}Er`$, $`{}_{}{}^{156}Gd`$, $`{}_{}{}^{158}Gd`$ $`{}_{}{}^{146}Nd`$, $`{}_{}{}^{144}Nd`$, $`{}_{}{}^{196}Pt`$, $`{}_{}{}^{148}Sm`$, $`{}_{}{}^{182}Ta`$,$`{}_{}{}^{174}Yb`$, $`{}_{}{}^{172}Yb`$, $`{}_{}{}^{90}Zr`$, $`{}_{}{}^{139}Ba`$, $`{}_{}{}^{137}Ba`$, $`{}_{}{}^{120}Sn`$, $`{}_{}{}^{106}Pd`$. The shape of the RSF within all models was investigated and compared with experimental data; the $`\chi ^2`$-criteria was used. In Fig.1 the shape of the RSF in $`{}_{}{}^{146}Nd`$ is given as an example. Here, only the MLO2 model is shown, because of similarity of the RSF for different MLO approaches. The following global set of parameters was used $`F=1,k_s=0.3,n_s=0.1.`$

This set can be recommended for calculations within MLO models. The optimal parameters for each nuclear can be used for more precise calculation. Calculation within GFL and MLO models are rather similar, and they are in close agreement with experimental data.

The GFL and MLO models can be also used to estimate the M1 strength $`\stackrel{}{f}_{M1}(E_\gamma )`$ in a wide interval of $`\gamma `$-ray energies on the base of the ratio $`R=\stackrel{}{f}_{E1}(B_n)/\stackrel{}{f}_{M1}(B_n)`$ for the E1 and M1 strength functions at neutron binding energy $`B_n`$. Then the M1 strength function is calculated by the following relationship

$$\stackrel{}{f}_{M1}(E_\gamma )=\frac{\stackrel{}{f}_{E1}(B_n)}{R}\frac{\varphi _{M1}(E_\gamma )}{\varphi _{M1}(B_n)},$$
(12)
where $`\varphi _{M1}(E_\gamma )`$ is a function describing shape of the dipole magnetic radiative strength function. The magnitude of the $`R`$ can be obtained from experimental data or systematics : $`R=\stackrel{}{f}_{E1}(B_n)/\stackrel{}{f}_{M1}(B_n)=0.0588A^{0.878},B_n7MeV.`$

Two models  for a function $`\varphi _{M1}`$ are usually used:

1) $`\varphi _{M1}(E_\gamma )=const`$ according to single- particle model;

2) $`\varphi _{M1}(E_\gamma )`$ is taken in form of the SLO model(1) corresponding to spin-flip giant resonance mode with the following global values  for energy and damping width (in $`MeV`$): $`E_r=41A^{1/3}`$, $`\mathrm{\Gamma }_r=4`$.

The radiative strength functions for $`E1+M1`$ transitions were calculated for the following nuclei: $`{}_{}{}^{200}Hg`$, $`{}_{}{}^{196}Pt`$, $`{}_{}{}^{198}Au`$, $`{}_{}{}^{192}Ir`$, $`{}_{}{}^{190}Os`$, $`{}_{}{}^{183}W`$, $`{}_{}{}^{182}Ta`$, $`{}_{}{}^{181}Hf`$, $`{}_{}{}^{177,176}Lu`$, $`{}_{}{}^{174}Yb`$, $`{}_{}{}^{170}Tm`$, $`{}_{}{}^{168}Er`$, $`{}_{}{}^{166}Ho`$, $`{}_{}{}^{164}Dy`$, $`{}_{}{}^{160}Tb`$, $`{}_{}{}^{158}Gd`$, $`{}_{}{}^{150}Sm`$, $`{}_{}{}^{146}Nd`$, $`{}_{}{}^{140}La`$, $`{}_{}{}^{139,138,137}Ba`$, $`{}_{}{}^{128}I`$, $`{}_{}{}^{125,124}Te`$, $`{}_{}{}^{114}Cd`$, $`{}_{}{}^{90}Zr`$, $`{}_{}{}^{80}Br`$. The experimental data were taken from Ref..

In Fig.2 the RS functions in $`{}_{}{}^{174}Yb`$ are shown. Here, the global parameter set is used for MLO2 model. It can be seen that MLO and GFL models are also usable in case of the $`E1+M1`$ emission.

## 4 Conclusions

The numerical studies led to the following conclusions on the description of the $`E1`$ and $`E1+M1`$ gamma-decay strength functions by the simple analytical models. The calculations by the MLO and GFL models are in more close agreement with experimental data than that ones within the EGLO and SLO methods at the energies $`ϵ_\gamma <4.5MeV`$. The overall comparison between calculations within the MLO, EGLO, SLO and GFL models and experimental data showed that the MLO and GFL approaches provide a rather reliable method of the $`\gamma `$ \- decay strength description in a relatively wide energy interval ranging from zero gamma-ray energy to values above GDR peak energy. The MLO and GFL methods are not time consuming. They can be applicable for calculations and predictions of the statistical contribution to the dipole strength functions as well as for extraction of the GDR parameters of heated nuclei with small errors with use of the $`\gamma `$-emission data.

This work is supported in part by the IAEA (Vienna) under contract 302-F4-UKR-11567.