Semiimplicit scheme for treating radiation under M1 closure in general relativistic conservative fluid dynamics codes
Abstract
A numerical scheme is described for including radiation in multidimensional generalrelativistic conservative fluid dynamics codes. In this method, a covariant form of the M1 closure scheme is used to close the radiation moments, and the radiative source terms are treated semiimplicitly in order to handle both optically thin and optically thick regimes. The scheme has been implemented in a conservative general relativistic radiation hydrodynamics code KORAL. The robustness of the code is demonstrated on a number of test problems, including radiative relativistic shock tubes, static radiation pressure supported atmosphere, shadows, beams of light in curved spacetime, and radiative Bondi accretion. The advantages of M1 closure relative to other approaches such as Eddington closure and fluxlimited diffusion are discussed, and its limitations are also highlighted.
keywords:
accretion, accretion disks, radiative transfer1 Introduction
Accretion disks are the power source behind many astrophysical systems. They span a wide range of central object mass, size and type, e.g., young stellar objects, neutron star binaries, black hole binaries, gammaray bursts, active galactic nuclei, to name a few, and exhibit a variety of different regimes of physics. In many systems, radiation is so intense that it strongly couples to the accreting gas and dramatically alters the flow structure and dynamics. To correctly infer the physics of such systems, we need models that properly take into account the interaction of gas and radiation.
Proper treatment of radiation is especially important for black holes (BHs) near the Eddington luminosity limit, . Transient black hole binaries (BHBs) approach nearEddington accretion rates near the peak of their outbursts (McClintock & Remillard, 2006; Remillard & McClintock, 2006; Done, Gierliński & Kubota, 2007), while some exceptional BHBs spend extended periods of time with (e.g., SS433, Margon et al. 1979; Margon 1984; GRS1915+105, Fender & Belloni 2004). Ultraluminous Xray sources have even larger luminosities and may conceivably be highly superEddington stellarmass BHs (Watarai et al., 2001), or intermediate mass BHs accreting at close to Eddington (Miller & Colbert, 2004). In either case, radiation must play an important role. Finally, luminous active galactic nuclei, especially those whose supermassive BHs (SMBHs) are growing rapidly in mass, may be perfect examples of systems with closely coupled gas and radiation (Collin et al., 2002).
Since most of the radiation from an accretion disk originates in the inner regions, general relativistic (GR) effects play an important role in determining the emergent radiation spectrum. Treating the interaction of fluid dynamics and radiation quantitatively in GR numerical codes is a daunting task that has been achieved only recently and only for optically thick flows, as we discuss below. However, the emergent spectrum is established at an optical depth of order unity, , which requires a proper treatment of the opticallythick disk (), the opticallythin corona () if one is present, and the optically thicktothin transition in between. In this paper, we present a GR numerical radiation hydrodynamics technique and a code which handles all three regimes of optical depth.
Depending on the mass accretion rate, a given accretion system can switch between different spectral states, with different radiation mechanisms dominating and with varying degrees of coupling between radiation and gas. At very low accretion rates, , e.g., Sgr A the SMBH at our Galactic Center (Narayan, Yi, & Mahadevan, 1995), the accretion flow advects most of the viscously released energy into the BH rather than radiating the energy (some energy probably also flows out in a wind). Such flows are optically thin, radiatively inefficient and geometrically thick, and are called advectiondominated accretion flows (ADAFs, see Narayan & Yi 1994, 1995; Abramowicz et al. 1995). Analytical and semianalytical models are reasonably successful in accounting for the main features in the spectra of ADAFs (e.g., Yuan, Quataert, & Narayan 2003). However, it was realized early on that numerical simulations are necessary to understand properly the physics of ADAFs.
The low optical depth and the extreme radiative inefficiency of an ADAF are a major advantage for simulations since they allow one to solve the fluid dynamics separately from the radiation field^{1}^{1}1However, Dibi et al. (2012) have shown that the accretion rate of Sgr A* is near the limit of the regime where radiative cooling may be important. A number of sophisticated GR magnetohydrodynamic (GRMHD) codes have been developed (e.g., De Villiers et al., 2003; Gammie et al., 2003; Anninos et al., 2005; Del Zanna et al., 2007), and ADAFlike models have been simulated using these (see Narayan et al. 2012; McKinney et al. 2012; Tchekhovskoy & McKinney 2012 for recent work). To compute observables, the output from such pure GRMHD codes are usually postprocessed with standalone radiation transfer schemes (e.g. Schnittman et al., 2006; Shcherbakov et al., 2010). Alternatively, a simple local cooling prescription is included directly into the code (e.g., Fragile & Meier, 2009; Dibi et al., 2012); this is not difficult since the gas is optically thin, so no radiative transfer is involved.
For higher accretion rates, , efficient cooling sets in, and the inner accretion disk collapses into an optically thick geometrically thin accretion disk (Narayan & Yi, 1995; Esin et al., 1997, 1998; MeyerHofmeister & Meyer, 2003; McClintock & Remillard, 2006; Done, Gierliński & Kubota, 2007). This state of accretion is the best understood of all states, and its thermal blackbodylike spectrum is welldescribed by the standard thin disk model with viscosity (Shakura & Sunyaev, 1973; Novikov & Thorne, 1973; Frank et al., 2002). However, despite its many successes, the disk model cannot describe important microphysical aspects of the flow, e.g., thermal and viscous stability, vertical structure of the disk, and the role of magnetic fields. Moreover, optically thick, geometrically thin disks are much more difficult to simulate numerically because of the strong coupling between gas and radiation.
Only a few timedependent GRMHD numerical simulations of thin disks have been computed so far, and all are based on an ad hoc (although physically motivated) local cooling function which artificially removes excess heat to keep the disk geometrically thin (Shafee et al. 2008; Penna et al. 2010; Noble et al. 2011; see also Reynolds & Fabian 2008). Although this approach has produced valuable results on the dynamics of the gas, it misses real physics involving propagation of photons along curved geodesics, radiation pressure support, radiative winds, etc.
Using a fluxlimited diffusion approximation, small patches of radiatively efficient thin accretion disks have been simulated using the local shearing box approximation (Turner et al., 2003; Krolik et al., 2007; Blaes et al., 2007, 2011; Hirose et al., 2009a, b), and even a few full disk simulations have been done with a nonrelativistic code using fluxlimited diffusion (Ohsuga et al., 2009; Ohsuga & Mineshige, 2011). However, to model a radiatively efficient disk selfconsistently it is necessary to handle the radiation field in both the optically thick (disk interior) and optically thin (corona) limits. This is difficult with the fluxlimited diffusion approximation which artificially enforces the flux to follow the local gradient of the radiative energy density. Fueled by advances in radiation algorithms, more advanced radiation moment closures are now becoming possible, which allow accurate treatment of both optically thick and thin photon fields in the “instant light” (nonrelativistic) approximation (Hayes & Norman, 2003; González et al., 2007; Jiang et al., 2012; Davis et al., 2012). However, as we discuss below, GR radiation transfer codes still continue to rely on fluxlimited diffusion or the Eddington approximation.
At superEddington accretion rates, , accretion flows again become radiatively inefficient. Here, the optical depth is so large that the photon diffusion time from the disk interior to the photosphere becomes longer than the accretion time. As a result, most of the photons are advected with the gas into the BH, leading to a “slim accretion disk” (Abramowicz et al., 1988). This important but poorly understood accretion state may be responsible for much of the SMBH mass growth in the Universe; for instance, it might explain the paradox of having BHs already at , when the Universe was less than Gyr old (Barth et al., 2003; Willott et al., 2005; Fan et al., 2006; Willott et al., 2010; Mortlock et al., 2011). SuperEddington accretion may also apply to ultraluminous Xray sources (e.g., Watarai et al., 2001; Kawashima et al., 2012). Clearly, to understand the accretion physics in these systems, radiation MHD models that selfconsistently couple gas, radiation and magnetic fields, are crucial. Such models are also important for measuring BH spins (McClintock et al., 2011; Straub et al., 2011), calculating the radiative efficiency of superEddington accretion disks (e.g., Sa̧dowski, 2009), and understanding largescale cosmological feedback by radiationdriven outflows from AGN (see Fabian 2012 for an observational review). Some efforts have already been made on attacking these important problems. In particular, using a nonrelativistic code and fluxlimited diffusion, superEddington accretion flows have been simulated and their spectra computed (Ohsuga et al., 2003; Ohsuga, 2006; Ohsuga & Mineshige, 2011; Kawashima et al., 2012).
Until a few years ago, all radiation hydrodynamic and MHD simulations of accretion disks were run with nonrelativistic or special relativistic codes. However, physics around BHs must be studied in GR because of the the strong gravity involved. Recently, progress has been made in implementing radiation into GR codes. Farris et al. (2008) developed a formalism for incorporating radiation under the Eddington approximation in a conservative MHD code. This method has been implemented in other codes (Zanotti et al., 2011; Fragile et al., 2012), but all these codes have problems because of the stiffness of radiative source terms at large optical depths. These stiff terms require implicit treatment, which is prohibitively complicated in curved spacetime. Roedig et al. (2012), extended the method by applying an implicitexplicit RungeKutta numerical scheme. However, the authors again relied on the Eddington approximation, which cannot handle optically thin flows accurately. A more advanced method has been recently described by Shibata & Sekiguchi (2012), who employ a truncated moment formalism, similar to our method, for neutrino transport in GR.
In the present paper, we describe a simple approach for including radiation transport in GR codes. The method makes use of a key simplification that is intrinsic to GR applications. Normally, in nonrelativistic codes, hydrodynamic or MHD signal speeds, which determine the time step via the Courant condition, are much slower than the speed of light. Correspondingly, the time step that one uses to evolve the fluid equations is much longer than the lightcrossing time across a cell. This is a problem when one wants to simulate radiation hydrodynamical systems, especially in optically thin regions, where one is faced with a large mismatch between the characteristic speeds of the fluid and the radiation. Either one must limit the time step to the lightcrossing time, which prohibitively increases the computational cost (as the fluid dynamics evolve on a much slower time scale), or one must handle all the radiation terms via an implicit method. The latter inevitably couples neighboring cells and makes the code very complicated. Moreover, it does not easily generalize to curved spacetime.
In contrast, a GR code is applied only in relativistic spacetimes near BHs and neutron stars. The time step in simulations is generally set by applying the Courant condition to the smallest grid cell, located at the innermost radius, where the flow is relativistic. Thus, the normal time step in a pure GR hydro or GRMHD code is already limited by the speed of light. Therefore, including radiation as an extra relativistic fluid is fairly easy. In particular, the advection terms in the radiation equations are no more difficult to compute than the corresponding terms in the fluid and magnetic field evolution equations. Nor does the time step need to be adjusted in any way. Because of this large simplification, the advective radiation operator can be treated via a standard explicit approach, just as one handles the corresponding hydro and MHD terms.
Of course, the interactions between radiation and gas via emission, absorption and scattering introduce their own time scales. These can sometimes be very short, requiring implicit handling. However, these interactions are local and can be handled via a local implicit scheme. This is a great simplification. It means that one can do implicit evolution independently in each grid cell, without coupling to other cells. Therefore, there are no spacetime curvature effects to contend with as in other more sophisticated multicell implicit schemes in GR.
In the work described here, we have implemented the above approach, closing the radiation moment equations using the M1 closure scheme (Levermore, 1984; Dubroca & Feugeas, 1999; González et al., 2007). M1 closure allows a limited treatment of anisotropic radiation fields and works well in both optically thick and thin regimes. We have implemented our method in a GR radiation hydrodynamics (GRRHD) code KORAL. The structure of the paper is as follows: In Section 2 we introduce the equations, in Section 3 we describe the numerical algorithm used by KORAL, in Section 4 we present a set of test problems which validate the scheme, and in Section 5 we discuss possible applications of the code.
2 Equations
2.1 Conservation laws
A pure hydrodynamic flow is described by the following conservation laws,
(1)  
(2) 
where is the gas density in the comoving fluid frame, is the gas fourvelocity as measured in the “lab frame”, and is the hydrodynamical stressenergy tensor in this frame,
(3) 
with and representing the internal energy and pressure of the gas in the comoving frame.
In the case of radiation hydrodynamics, it is convenient to introduce the radiation stressenergy tensor (e.g., Mihalas & Mihalas, 1984), and to replace the second equation above with the more general conservation law,
(4) 
The radiation stressenergy tensor in an orthonormal frame comprises various moments of the specific intensity , e.g., in the fluid frame it takes the following form,
(5) 
where the fluidframe quantities^{2}^{2}2Throughout the paper, “widehats” denote quantities in the fluid frame and “tildes” denote quantites in the ZAMO frame.
(6)  
(7)  
(8) 
are the radiation energy density, the radiation flux and the radiation pressure tensor, respectively, and is a unit vector in direction .
The fluid frame radiation stress tensor is related to the tensor defined in the locally flat nonrotating frame, or the zeroangular momentum frame (ZAMO) (Bardeen et al., 1972), by
(9) 
where is the Lorentz boost,
(10) 
, , and is the fourvelocity of the gas as measured by the locally nonrotating observer.
Quantities in the ZAMO frame (denoted with tildes) are related to those in the lab frame by tensors created from the components of the corresponding tetrads of the ZAMO, and , defined in Bardeen et al. (1972). The radiation tensor transforms as
(11)  
(12) 
while fourvectors transform as
(13)  
(14) 
The conservation law (4) may be rewritten with the help of the radiation fourforce density as
(15)  
where is given by (Mihalas & Mihalas, 1984),
(16) 
which takes a particularly simple form in the fluid frame,
(17) 
Here, is the integrated Planck function corresponding to the gas temperature , is the StefanBoltzmann constant, and denote the frequencydependent opacity and emissivity coefficients, respectively, while and are the frequency integrated absorption and total opacity coefficients, respectively. In Section 4, we occasionally refer to the scattering opacity , which is related to and by
(18) 
The fourforce may be transformed between frames as described above, e.g.,
(19) 
The rest mass conservation equation (1) and the energymomentum conservation equations (15) may be written in a coordinate basis in the following conservative form (Gammie et al., 2003),
(20)  
(21)  
(22) 
This formulation has a drawback for numerical computations: the terms involving Christoffel symbols on the right, when calculated at cell centers, will not balance the corresponding spatial derivatives on the left (approximated under a given reconstruction scheme). This is true even for particularly simple situations such as constant gas or radiation pressure, and can lead to catastrophic secular errors. To solve this issue we can either modify the values of the Christoffel symbols suitably (Appendix A of McKinney et al. 2012) or we can reformulate the equations so as to avoid the problem. We choose the second approach and make use of the following equations,
(23)  
(24)  
(25) 
where we assumed that the metric is static and moved its determinant out of the derivatives on the left. In this formulation, the two terms that are expected to cancel each other both appear as source terms and therefore are calculated at the same location^{3}^{3}3In principle, it is sufficient to move the determinant out of the spatial derivatives in and components of equations (21) and (22). .
2.2 Closure scheme
To close the above set of equations we need a prescription to compute the second moments of the angular radiation intensity distribution. Specifically, we need a prescription to write down the full radiation stress tensor knowing only the radiative energy density and the fluxes .
The simplest approach, which corresponds to assuming a nearly isotropic radiation field, is the Eddington approximation, which in the fluid frame gives
(26) 
However, the assumption of isotropic specific intensity is good only in optically thick media. In many astrophysical applications we are interested in radiation that escapes from the photosphere to infinity, for which we need a better closure scheme.
Following Levermore (1984), we assume that the radiation tensor is isotropic and satisfies the Eddington closure, not in the fluid frame, but in the orthonormal “rest frame” of the radiation. The latter is defined as the frame in which the radiative flux vanishes. Thus, in this frame, we assume that , , and all other components of are zero. This leads to the M1 closure scheme.
In the radiation rest frame, the radiation stress tensor can be written in a compact form as
(27) 
where and is the contravariant metric tensor, which in this frame is given by the flat space Minkowski metric. Since equation (27) is in a covariant form, it is also valid in the lab frame (as well as all other frames), with being the fourvelocity of the radiation rest frame as measured by an observer in the lab frame. Note that, regardless of which frame one works in, the quantity should be interpreted as the radiation energy density as measured in the radiation rest frame.
Each time step in the numerical integration in any particular cell gives an update to the “time” row of the radiation tensor, , for that cell in the lab frame. Thus, we obtain numerical values of these four particular components of the tensor. According to equation (27), the full tensor is a function of 5 numbers, , , though only four of these are independent since the norm of the fourvelocity is equal to . Hence, we can use the four given tensor elements to solve for the four unknowns and thereby compute the full matrix. Below we give an algorithm for doing this analytically.
Consider the quantity which can be expressed as (equation 27),
(28) 
If we contract this with and use the following results,
(29)  
(30)  
(31) 
we obtain
(32) 
The lefthand side is computable from the four given tensor elements and the righthand side involves two of the unknowns: , . We also have the following expression for ,
(33) 
which again involves the same two unknowns. Thus, we can solve equations (32) and (33) to obtain and (it reduces to a quadratic equation). It is then straightforward to calculate the remaining from the other time components of equation (27) and to calculate the entire radiation stress tensor.
For flat spacetime, the above formulation reduces to the standard formulae (Levermore, 1984; Dubroca & Feugeas, 1999; González et al., 2007). For instance, the radiation pressure tensor in the fluid frame has the form,
(34) 
where is the reduced radiative flux and is the Eddington factor given by (Levermore, 1984),
(35) 
In the extreme “optically thick limit”, , and we find , and , which corresponds to the correct answer, viz., the Eddington approximation,
(36) 
In the opposite extreme “optically thin limit”, , i.e., a unidirectional radiation field directed along the xaxis, we have , and , which gives
(37) 
This corresponds to an intensity distribution in the form of a Dirac function parallel to the flux vector, which is again the correct answer. The M1 closure scheme thus handles both optical depth extremes well and it is found to be fairly good at intermediate optical depths as well.
As explained, the M1 closure scheme assumes that radiation is isotropic in the radiation “rest frame”. The stress tensor in an arbitrary frame is obtained by applying a Lorentz boost to the isotropic “rest frame” tensor. As a result, only one direction, the direction of the boost, is distinguished. In other words, the specific intensity is always symmetric with respect to the mean flux. The M1 closure scheme is thus expected to be only approximate when multiple sources of light are involved (see Section 4.5). In problems involving accretion disks, which are the primary area of interest of the present authors, highly anisotropic configurations with multiple beams are not very common, and the M1 scheme is probably adequate. In any case, M1 closure will provide a significantly superior treatment of radiation in the optically thin regions near and above the disk photosphere, compared to the Eddington approximation or fluxlimited diffusion.
3 The Koral code
The scheme described in this paper has been implemented into a GRRHD code KORAL which solves equations (23)–(25) in an arbitrary metric. The code uses a finite difference scheme with either linear slopelimited reconstruction (Kurganov & Tadmor, 2000) or fifthorder nonlinear monotonizing filter reconstruction (MP5, see Suresh & Huynh, 1999; Del Zanna et al., 2007). The fluxes at the cell faces are calculated using the LaxFriedrichs scheme. The source terms are applied at the cell centers and the time stepping is performed using the optimal RungeKutta method of third order (Shu & Osher, 1984). The vector of conserved quantities is (Section 3.5)
(38) 
while the primitive quantities are,
(39) 
Conversion from conserved to primitive quantities is described in Section 3.4, while the algorithm itself is described in the next one.
3.1 The algorithm
During each substep of the RungeKutta time integration, the code carries out the following steps in the given order:

The vector of conserved quantities in each cell is used to calculate the primitive quantities at the cell center (see Section 3.4).

Ghost cells at the boundaries of the computational domain are assigned primitives appropriate to the boundary conditions of the particular problem of interest.

For each cell, the maximal characteristic left and rightgoing wave speeds are calculated, following the algorithm described in Section 3.2.

For each dimension, primitives are interpolated using the chosen reconstruction scheme (linear slopelimiter or nonlinear monotonizing filter) to obtain their left and rightbiased values at cell faces: and .

From and , left and rightbiased fluxes and are calculated at cell faces.

The flux at a given cell face is calculated using the LaxFriedrichs formula,
(40) where is the maximal absolute value of the characteristic speeds at the centers of the two cells on either side of the face, and and are the conserved quantities calculated at the cell face based on and .

The advective time derivative is calculated using an unsplit scheme,
(41) where denotes cell size in the direction . Note that all primitives, including the radiation density and radiation flux , are treated identically as far as the advective term is concerned.

The advective and geometrical operators are used to update the conserved quantities according to
(42) That is, all these terms are treated in an explicit fashion.

The updated vectors of conserved quantities are used to calculate the corresponding updated primitive quantities at cell centers.

Finally, the remaining terms, viz., those involving the fourforce density , are handled implicitly using the method described in Section 3.3. This results in a final update of the vector of conserved quantities at each cell center.
3.2 Characteristic wavespeeds
The LaxFriedrichs scheme requires knowledge of the maximal characteristic wave speeds of the system ( in equation 40). The hydrodynamical and the radiative components of equations (20)–(22) are coupled only through the radiative source term . Therefore, for the purpose of calculating the fluxes at cell faces and the advective time derivative, we are allowed to separate the hydrodynamical and radiative wave speeds. We can calculate each separately, and use it for evaluating its corresponding flux. Such an approach avoids excessive artificial numerical viscosity which appears when the characteristic wavespeeds are not separated. The radiative wavespeed, as described below, never drops below . If such a high value was used in equation (40) for the hydrodynamical subsystem, it would result in strong artificial diffusion of the gas.
The hydrodynamical characteristic velocity in the fluid frame is the local speed of sound,
(43) 
To get the left and rightgoing wave speeds in the lab frame we transform this velocity in the standard way (e.g., Gammie et al. 2003).
The evolution of the radiation field is described in the fluid frame by the following set of equations,
(44)  
which is the nonrelativistic limit of equations (22) in flat spacetime. Given a closure scheme, it is possible to calculate the Jacobi matrices,
(45) 
whose eigenvalues give the wave speeds of interest in a given direction . The maximal left and right going speeds can then be transformed to the lab frame following the same method as for fluid velocities.
For the M1closure scheme, the equation for the eigenvalues is quartic and may be solved efficiently using standard algorithms. Another approach, which will likely improve code performance, is to precalculate tables of radiative wave speeds as a function of the reduced flux vector components (González et al., 2007). At present, KORAL uses the analytical approach.
In the limit of large optical depths, the radiative energy density, when decoupled from gas (e.g., for but ), has a diffusion coefficient given by (see Section 4.4)
(46) 
In this limit the distribution of radiative energy density should remain stationary (). On the other hand, the maximal eigenvalue of the Jacobi matrices (equation 45) is (González et al., 2007). If such large wave speeds are incorporated into a numerical scheme they will result in large, unphysical, numerical diffusion. To limit this effect, we modify the radiative wave speeds in the fluid frame according to
(47)  
where and are the maximal right and leftgoing radiative wave speeds in the fluid frame in the direction , and is the total optical depth of a given cell in that direction.
The smaller the characteristic wave speed in equation (41), the weaker the numerical diffusion. Thus, one may be tempted to set the wave speed to zero. However, the numerical scheme will then no longer satisfy the total variation diminishing (TVD) condition and the algorithm will be unstable. Our choice of the wave speed limiter (equation 47) is motivated by the fact that, for a diffusion equation of the form , the maximum allowed time step for an explicit numerical solver is
(48) 
This expression, combined with equation (46), gives
(49) 
which is the limiter introduced in equation (47).
3.3 Implicit treatment of radiative source terms
It is wellknown that, under some circumstances, e.g. for large optical depths, the radiative source terms in equations (21) and (22) become stiff, making explicit integration practically impossible (e.g., Zanotti et al., 2011). We then need to treat these terms via implicit time integration. In principle, we could try to solve the whole system of partial differential equations (23)–(25) implicitly, as done for instance in nonrelativistic or special relativistic radiation hydrodynamics codes (e.g., Krumholz et al., 2007; Jiang et al., 2012). However, this approach is very difficult in GR, where the curvature of spacetime makes the problem highly complicated and it is nontrivial to ensure that an implicit code is conservative.
As already explained, our approach is to split the advective derivative operator from the radiative source terms operator. The former is applied explicitly in the usual way while the latter is handled implicitly. This approach is possible because the time step is already limited by the speed of light just from the fluid dynamics, so radiation advection is also guaranteed to be stable in an explicit scheme.^{4}^{4}4 Note that if we were to apply the code to nonrelativistic problems, or situations involving weak gravity, there would be a large disparity in the advective time scales of fluid and radiation, and the code would no longer be efficient. The main advantage is that the radiation source terms are local, so they can be treated semiimplicitly and pointwise in the fluid frame, without having to deal with curvature effects. In our experience, this approach is both simple and robust.
The radiative source term operator describes the action of the radiative fourfource on the energy and momentum density of the gas and the radiation. The corresponding equations are
(50)  
(51) 
In an explicit scheme, updates would be calculated very simply as
(52)  
(53) 
where the subscripts and denote values at the beginning and end of a time step of length , respectively. This approach, though simple, is numerically very unstable whenever the terms in the force vector are large. Our scheme avoids the instability by computing the updates implicitly via
(54)  
(55) 
i.e., using quantities at time rather than to compute the force vector on the righthand side. It is wellknown that this simple change has a profound effect on stability.
KORAL solves equations (54) and (55) numerically. Because of the symmetry of the problem, specifically, the righthand sides of eqs. (54) and (55) differing only by sign, the system of equations may be reduced to four nonlinear equations, e.g., for . We use the standard Newton method to solve these equations and estimate the Jacobian matrix numerically. During each iteration we use the current guess of to calculate the corresponding values of ,
(56) 
We then convert the current vector of conserved quantities to primitives (see Section 3.4) and calculate the radiative four force in the fluid frame, . This is then boosted to the laboratory frame to obtain .
The method described above is numerical and can occasionally fail. In Appendix A we describe an alternate, fully analytical, but approximate, method for applying the radiative source operator. KORAL uses that algorithm as a failsafe backup whenever the numerical scheme descibed in this section fails.
3.4 Conversion between conserved and primitive quantities
In conservative GR numerical codes it is necessary to convert conserved quantities to primitives at least once per subtime step (twice in the case of our algorithm). In our problem, the hydrodynamical and radiative variables decouple, so the conversion may be done separately for each.
The conversion of hydrodynamical quantities is performed in the usual manner (e.g., Noble et al., 2006; Del Zanna et al., 2007). We take the explicit forms of the conserved quantities (equation 3), use the equation of state , the fourvelocity normalization , and combine all the terms into a nonlinear equation for the internal energy density . We solve this equation numerically by the Newton method using the value of from the end of the previous time step as the initial guess.
The radiative conserved variables () may be easily converted to the radiative primitives (, ). First, one has to calculate all the components of the radiation stress tensor in the lab frame following the algorithm described in Sect. 2.2. This tensor is then boosted to the fluid frame (Sect. 2.1),
(57) 
The fluidframe radiative energy density and fluxes are then given by,
(58)  
(59) 
3.5 Implementation notes

The mass conservation equation (23) and the gas internal energy conservation law, i.e., the component of equation (24), are aggregated to give (Gammie et al., 2003)
(60) which replaces the component of equation (24). Then, becomes the relevant conserved quantity, which reduces in the nonrelativistic limit to .

In cold relativistic flows, where , the numerical accuracy is not suficient to evolve the internal energy reliably. As a result, negative internal energy densities may be occasionally found. Whenever this happens, we recalculate the gas properties by evolving the gas entropy per unit volume,
(61) That is, we compute the entropy at the end of the last successful time step and evolve it according to .
4 Test problems
4.1 Relativistic shock tube
To validate the hydrodynamics part of our code, we tested it using the relativistic shock tube problem introduced by Hawley et al. (1984). The modeled system consists of two states (left and right), separated initially by a membrane. The gas on the left is dense and hot, while that on the right is rarefied and cool (Table 1). At time the membrane is removed and the hot gas pushes the cool gas to the right causing a shock to travel to the right. Meanwhile, a rarefaction wave moves at the local sound speed back to the left. An analytical solution for this shock tube problem may be obtained by solving the appropriate jump conditions.
In Figure 1 we show our numerical solutions at time , which may be compared with the exact solution (black dashed line) obtained using an exact Riemann solver (Giacomazzo & Rezzolla, 2006). The upper panel shows profiles of density obtained using three different reconstruction schemes: linear slope limiter with , linear slope limiter with ,^{5}^{5}5The parameter in the generalized minmod limiter determines the diffusivity of the scheme (Kurganov & Tadmor, 2000); corresponds to the MC (monotonized central) scheme and corresponds to the MINMOD scheme. and fifthorder monotonizing filter, MP5. The rarefaction wave and the plateau are resolved equally well in all three schemes. However, there are differences in the post shock region. The linear scheme (red lines) is most diffusive and smooths the edges of the postshock region, while the MP5 scheme (blue) somewhat underestimates the density here. The lower panel shows the velocity. This is reproduced well in all three schemes, though the MP5 scheme produces lowlevel oscillations near the edge of the rarefaction wave.
Left state:  Right state:  
5/3  10.0  13.33  0.0  1.0  0.0 
This test shows that all the reconstruction schemes currently implemented in KORAL work reasonably well on relativistic hydrodynamic shocks. The MP5 scheme is most accurate in smooth regions, but it is also prone in some circumstances to give unphysical oscillations. The linear scheme is most diffusive, and at the same time is also the most stable. In the following tests, if not stated otherwise, we use the linear reconstruction scheme with .
4.2 Polish doughnuts
To test the ability of the code to handle multidimensional hydrodynamics in curved spacetime, we set up analytical equilibrium torii (Polish doughnuts, Abramowicz et al., 1978) in the Schwarzschild metric as initial conditions and let the code evolve to a numerical equilibrium configuration. For the analytical model, we assume a constant specific angular momentum, constant. From the condition , it follows that
(62) 
We choose the specific internal energy at the inner edge of the torus, , which determines the radius of the inner edge of the torus, and we then calculate the fluid enthalpy, (e.g., Hawley et al., 1984),
(63) 
Using an equation of state (where the constant determines the entropy of the torus gas), we obtain
(64)  
(65) 
We set the initial velocity to , , and choose .
Table 2 gives parameter values corresponding to three models that we ran to test the code. Models 1 and 2 have the same value of the specific angular momentum, . Model 1 has , corresponding to a torus inner radius , while Model 2 has , corresponding to . The specific angular momentum for the third model, , lies in between the Keplerian values of at the marginally stable and marginally bound orbits (). Therefore, this torus has a cusp (selfcrossing equipotential surface) near the marginally bound orbit (). All the equipotential surfaces outside the body of the torus (defined by the critical surface producing the cusp) are open and continue into the BH. Therefore, the assumption of zero poloidal velocity cannot be applied to this region.
Model  

1  4.5  1.0  0.03 
2  4.5  0.98  0.03 
3  3.77  1.0  0.03 
The above three torus models were simulated on a 50x50 grid in BoyerLindquist () coordinates, linearly spanning the range and . At the spin axis and the equatorial plane, reflection boundary conditions were assumed. The analytical solution was imposed as boundary conditions in the ghost cells outside . At the inner edge of the grid , outflow boundary conditions were implemented. Linear reconstruction with was used.
Figure 2 compares the relaxed, stationary solutions obtained by running the code with the corresponding analytical Polish doughnut solution. Colors and orange solid lines show the numerical solution (density), while dashed green lines show the analytical model. Vectors indicate the poloidal component of velocity.
The top panel in Figure 2 corresponds to Model 1, where the analytical solution extends from to the outer boundary. The analytical and numerical contours agree very well. The only visible discrepancy is in the lowdensity region near the torus inner edge, where there is some numerical dissipation.
The middle panel in Figure 2 corresponds to Model 2, where the analytical torus is entirely confined within the grid: , . Again, the numerical solution agrees very well with the analytical solution except for the innermost region of the torus where the numerical solution is slightly stretched inward.
The bottom panel in Figure 2 shows Model 3. The analytical contours plotted correspond to open equipressure surfaces which intersect the BH horizon. Nonzero poloidal velocities are thus expected in this model and the analytical solution is not fully consistent. Nevertheless, there is a close similarity between the numerical and analytical solutions, which means that the code is able to represent even this torus quite well.
4.3 Radiative shock tubes
The previous two tests did not involve radiation. For our first test with radiation, we set up a number of radiative shock tube problems as described in Farris et al. (2008) and Roedig et al. (2012). These shock tubes are similar to the pure hydrodynamic shock tube described in Section 4.1, i.e., the system begins with gas in two different states (left and right), separated by a membrane. The membrane is removed at and the system is allowed to evolve. The difference here is that the evolution is described by the full set of equations (23)–(25). In addition, the left and rightstates of all the tests except test No. 5 are set up in such a way that the shock asymptotically becomes stationary (see Appendix C of Farris et al., 2008).
Table 3 lists the parameters describing the initial states of seven test problems that we have simulated. The scattering opacity in all the tests is set to zero, so (equation 18). The value of the radiative constant in code units is given in the table. All the tests were solved on a grid of uniformly spaced points. For consistency with previous work by other authors, the Eddington approximation was used (equation 26).
Test  Left state:  Right state:  

No.  
1  5/3  0.4  1.0  0.015  2.4  
2  5/3  0.2  1.0  0.25  3.11  
3a  2  0.3  1.0  10.0  8.0  
3b  2  25.0  1.0  10.0  8.0  
4a  5/3  0.08  1.0  0.69  3.65  
4b  5/3  0.7  1.0  0.69  3.65  
5  2  1000.0  1.0  1.25  1.0 
Figure 3 shows the numerical (solid) and analytical (dashed lines) solutions for radiative shock tube problem No. 1, which corresponds to a nonrelativistic strong shock. The panels show (top to bottom) density, proper velocity, fluidframe radiative energy density and fluidframe flux, and may be directly compared (but for the flux which the other authors plot in the lab frame) to the corresponding figures and analytical solutions of (Farris et al., 2008; Zanotti et al., 2011; Fragile et al., 2012). The agreement is good.
Figure 4 shows results for radiative shock tube test No. 2, which corresponds to a mildly relativistic strong shock. Again, the agreement between the numerical and semianalytical (Farris et al., 2008) profiles is good, except for a slight smoothing of the numerical profiles at the position of the shock (see the bottom panel).
Figure 5 shows results corresponding to radiative shock tube tests No. 3a and 3b. These are strongly relativistic shocks with upstream . Test No. 3a corresponds to shock tube test 3 of Farris et al. (2008), while test 3b is the optically thick version of the same test which was proposed and solved by Roedig et al. (2012). These two tests verify that the code is able to resolve a highly relativistic wave in two very different optical depth limits. In both cases, the numerical solution reaches a steady state and closely follows the corresponding semianalytical solution. The only noticeable difference is that in the highopacity case the discontinuity in the numerical solution is less steep than in the exact analytical profile. This discrepancy measures the amount of numerical dissipation introduced by the algorithm.
Figure 6 shows results for radiative shock tube tests No. 4a and 4b. These tests correspond to radiation pressure dominated mildly relativistic waves. Test 4b is the optically thick version of test 4a that was proposed by Roedig et al. (2012). In both tests, the numerical solution reaches a stationary state and agrees well with the semianalytical solution. Note that the values of the opacity coefficient in tests 3b and 4b are the maximum values that the numerical scheme of Roedig et al. (2012) could handle. The algorithm implemented in KORAL has no such a limitation. We could increase to much larger values and the scheme would remain stable.
Figure 7 corresponds to radiative shock tube test No. 5. This is the only test that does not asymptote to a stationary solution. This test was proposed and solved by Roedig et al. (2012) and represents an optically thick flow with mildly relativistic velocities. The left and rightstates are identical except that they have different velocities. As a result, two shock waves propagate in opposite directions. This test does not have an analytical solution. However, by comparing our numerical solution with that presented in Roedig et al. (2012), we confirm that our scheme performs satisfactorily.
4.4 Radiative pulse
We now test the ability of our scheme to handle the evolution of a radiation pulse in the optically thick and optically thin limits. We start with the optically thin case. We set up a Gaussian distribution of radiative energy density at the center of a 3D Cartesian coordinate system. The pulse radiative temperature is set according to,
(66) 
with , and the radiative constant . We assume zero absorption opacity () and nonzero but negligible scattering opacity (). The background fluid field has constant density and temperature . We solve the problem in three dimensions on a coarse Cartesian grid of cells using the linear reconstruction with .
The initial pulse in radiative energy density is expected to spread isotropically with the speed of light (optically thin medium) and to decrease inversely proportionally to the square of radius (energy conservation). Such behavior is visible in Fig. 8 showing the radiative energy distribution in the plane (top panels) and its crosssection along (bottom panels). The orange circles in the top set of panels show the expected size of the pulse. It is clear that the propagation speed of the pulse is consistent. This problem was solved on a relatively coarse Cartesian grid. This results in deviations from the perfectly spherical shape — the radiative energy density of the pulse along the axes is higher than along the diagonals. This effect is reduced by choosing larger resolution or a more suitable grid (e.g., spherical). The bottom set of panels in Fig. 8 shows the profiles of the energy density along the axis. The black dotted lines show the expected rate of energy decrease with increasing distance from the center. The numerical solution perfectly follows this trend.
To test the optically thick limit we choose to set up a similar pulse but this time planar instead of a pointlike, i.e., according to,
(67) 
This time we set the scattering opacity to and solve the problem as onedimensional on grid points distributed uniformly between and with periodic boundary conditions in and . The total optical depths per cell and per pulse are therefore and , respectively.
In the optically thick limit the evolution of such a system is described by a diffusion equation,
(68) 
which can be derived from the nonrelativistic limit of Eq. 25 assuming the time derivative of the component of the flux vanishes. An initially Gaussian pulse of radiative internal energy will therefore diffuse according to the value of the diffusion coefficient .
In Fig. 9 we plot profiles of the radiative energy at various moments (solid lines) and compare them to the exact solution of Eq. 68 (dotted lines). The numerical solution diffuses slightly faster due to the additional numerical dissipation introduced by the scheme^{6}^{6}6This artificial numerical diffusion may be further reduced by stronger damping of the characteristic radiative wavespeed in the optically thick limit (Eq. 47). However, it might cause problems in the intermediate regime.. At later times this difference becomes insignificant.
4.5 Shadow test
Here we test the ability of the M1 closure scheme, as incorporated in KORAL, to resolve shadows. We set up a blob of dense, optically thick gas in flat spacetime, surrounded by an optically thin medium, and we illuminate this system.
We start with a single source of light imposed on the left boundary. We solve the problem in two dimensions on a 100x50 grid, with the density distribution set to be
(69) 
where , and . The gas temperature is adjusted so as to give constant pressure throughout the domain,
(70) 
The initial radiative energy density is set to the local thermal equilibrium value, and the initial velocities and radiative fluxes are zero. We apply periodic boundary conditions at the top and bottom and outflow boundary conditions at the right border of the domain. At the left border we have the external source of light, which we specify with , , . All other quantities are set to match the ambient gas. We evolve the system with both the Eddington approximation and M1 closure, assuming .
Figure 10 shows the results at . By this time, the initial radiation wave has passed through the domain and the system has reached a stationary state. The upper panel shows the solution we obtain with the Eddington approximation. Because this closure treats all directions equally, radiation readily diffuses into the region behind the blob. As a result, there is no shadow behind the optically thick blob. The lower panel in Figure 10 shows results with the the M1closure. In contrast to the case of Eddington closure, here the direction is distinguished because the incoming radiation at the left boundary moves in this direction. The M1 closure is designed to keep flux moving parallel to itself in optically thin regions for . As a result, a strong shadow develops behind the optically thick blob. This test illustrates one of the key differences between the Eddington and M1 closure schemes.
It is appropriate to mention that the excellent performance of the M1 closure scheme in this shadow test problem is partly because the setup is particularly favorable. First, we have only a single source of radiation. Second, the shadow here is aligned well with the grid, which helps to minimize diffusion. For other grid orientations, the numerical results would exhibit more diffusion, e.g., see Section 4.7.
The M1 closure assumes that the specific intensity is symmetric with respect to the radiative flux, i.e., only one direction is distinguished. It means that this approach is supposed to be less efficient when multiple sources of light are involved. To test its performance for such a setup we implemented a twobeam test problem similar to the one described in Jiang et al. (2012)^{7}^{7}7The key principle behind the nonrelativistic algorithm described in Jiang et al. (2012) is the use of a “Variable Eddington Tensor” (VET). The VET is used to close the radiative equations, relating radiative pressure to the local radiative energy density. The VET is computed through a separate radiative transfer solver, which calculates (at each time step) the timeindependent radiation field (using the fixed fluid background of the previous timestep as its input). The authors solve the radiative transfer equations along a discrete set of rays, and so their scheme accurately captures all shadows that can be resolved by these ray angles. The radiation pressure obtained from the VET is then used to evolve the MHD fluid equations.. We set up exactly the same initial conditions for gas and radiation as in the two previous tests. This time, however, we set up a reflection symmetry at the lower boundary () and we impose an inclined (, ) beam on the upper boundary and on the part () of the left boundary. As a result, the domain is effectivelly enlighted by two selfcrossing beams of light. We plot the result of a numerical simulation in Fig. 11 where the region of negative coordinates was plotted by reflecting the positive data. In the region near the left top and bottom corners, where the beams do not overlap, the direction of the flux follows the imposed boundary condition. In the region of the overlap the radiative energy density increases twice () while the flux becomes equivalent to the superposition of the beamintrinsic fluxes, i.e., it is purely horizontal and its component equals . The clump of optically thick gas is, therefore, effectively illuminated by a purely horizontal beam. One could expect it creates a parallel shadow similar to the one obtained in the single beam problem. This is, however, not the case. There is an important difference between the beam we imposed on the left boundary in the single beam problem and the one which develops in the overlapping region. The former had what implied that the specific intensity was almost a function in the direction of the flux. The latter, however, has which means that the implied distribution of the specific intensity is only an elongated ellipsoid pointing in the direction of the flux, i.e., photons moving in other directions than the direction of the flux are allowed. This has an effect on the shadow produced behind the clump. Instead of sharp parallel edges we get regions of the partial shadow (penumbra) resulting from these perpendicular photons allowed by the closure when . The region of the total shadow (umbra) is therefore limited by the edges of the penumbra and follows the expected shape (compare Fig. 11 in Jiang et al. (2012)) to a good accuracy. A significant difference between the exact solution and our numerical one arises in the region where the penumbrae overlap. One could expect a uniform, triangular region of no shadow. The M1 closure, however, produces an extra narrow horizontal shadow along the axis.
This test shows limits of the M1 closure approach but at the same time stresses the fact that, in principle, it does not limit specific intensity to one particular direction (assuming only its symmetry with respect to the flux). It performs much better than the Eddington approximation but in the case of multiple sources of light it must be used with caution.
4.6 Static atmosphere
An important aspect of radiation in accretion disks is momentum transfer between radiation and gas. The Eddington luminosity limit, for instance, arises from this interaction. To validate the treatment of gasradiation momentum exchange in our method, we study a static atmosphere which is in equilibrium under the combined action of gravity, gas pressure gradient and radiation force. We consider a polytropic atmosphere on a spherical object. We take the optically thin limit and assume that gasradiation interactions occur only through a scattering coefficient, i.e., , (equation 18).
An analytical solution is easily derived for this model problem. Because we assume a polytropic equation of state and set , there is no energy equation, and the radial component of the momentum equation ( component of equation 24) is all that matters. In the nonrelativistic limit (), assuming stationarity () and zero velocity (), this equation takes the form
(71) 
where
(72) 
Here is the radiative flux imposed as a boundary condition at the bottom of the atmosphere, , and gives the ratio of the radiative to gravitational (or geometrical) forces; corresponds to the Eddington limit, where the luminosity is and the radiative flux is . Since radiative energy must be conserved, in the stationary state the flux must satisfy (nonrelativistic limit).
Equation (71) may be solved with the help of the polytropic equation of state to give,
(73) 
where
(74) 
and is the assumed density at . The entropy constant is calculated at the bottom of the atmosphere from the assumed gas temperature .
We set up a linear grid of points between and gravitational radii and we solved the problem using MP5 reconstruction scheme and the standard M1 closure. We scaled all quantities to physical units assuming and . At the innermost radius we set (optically thin atmosphere) and . All the velocities were initially zero and the radiative energy density . Initial values of the gas density and temperature in the domain and in the ghost cells were assigned based on the analytical solution. We ran four models corresponding to four luminosities: , , and . Each model was run up to a time , which is sufficient to reach relaxed steady state for these optically thin atmospheres.
Figure 12 shows the results. The top panel shows the density profiles corresponding to the four models. Solid lines are the analytical solutions and filled circles correspond to the numerical solutions. The agreement is very good. The higher the luminosity, the flatter is the density profile, indicating the effect of the outward force due to radiation. For the particular case of the Eddington luminosity, the density is perfectly constant, reflecting the fact that the gravitational force is exactly balanced by radiation and no pressure gradient in required. We see that the relaxed numerical solution is indistinguishable from the analytical solution, confirming that KORAL properly handles gasradiation momentum exchange. The plot of residuals at the bottom of the panel indicates that fractional deviations in the density are below . Even this small discrepancy is at least in part because we are comparing a numerical solution from a GR code with an analytical solution derived under Newtonian physics.
The middle panel in Figure 12 shows our results for the radial radiative flux. Once again, the models behave very well and the agreement with the analytical solution is excellent. Finally, the bottom panel shows the residual radial velocities (). These are of the order of (they should be zero), and appear to be mostly driven by slight inconsistencies near the boundaries (possibly again because of a slight tension between GR and Newtonian physics).
4.7 Beam of light near BH
To test the performance of the code for radiation in strong gravitational field, we study propagation of a beam of light in the Schwarzschild metric. We consider three models, in each of which a beam of light is emitted in the azimuthal direction at a different radius. We decouple gas and radiation by neglecting absorptions and scatterings (). We run the models on a twodimensional grid with 30 points distributed uniformly in between and (see Table 4 for values) and 60 points distributed uniformly in azimuthal angle between and . Initially, we assign negligibly small values for all primitive quantities, including the radiation energy density and flux. We use outflow boundary conditions on all borders except the region covered by the beam at the equatorial plane (see the range of in Table 4), where we set the radiation temperature to and the flux to . Here is the initial gas and radiation temperature of the ambient medium. We use linear reconstruction with .
Model  

1  2.5  3.5  
2  5.3  7.5  
3  14.0  20.5 
The panels in Figure 13 show the results for the three models and geodesics of photons at beam boundaries. Consider the right panel, which corresponds to Model 3 (Table 4) with the beam centered at . At such a large radius we do not expect significant bending of photon geodesics and this is indeed the case — the beam is only slightly bent towards the BH. We also expect the beam to be tighly confined, i.e., it should propagate with a nearly constant width, as indicated by the two solid green lines, the true geodesics of photons at the boundaries of the beam. However, the numerical solution shows significant artificial broadening. This is caused by numerical diffusion which is significant whenever the radiative flux vector is not aligned with the grid geometry, i.e., the beam is tilted with respect to the grid axes.
The middle panel in Figure 13 shows Model 2, where the beam is centered at the marginally stable orbit: . At this radius, photon geodesics are significantly deviated by gravity, resulting in strong curvature in the beam. The numerical beam follows the correct trajectory. Moreover, numerical diffusion is lower in this case because the curvature brings the beam into closer alignment with the grid. There is in addition some real beam divergence because the geodesics at different radii within the beam have different curvatures (see the solid green lines), but this effect is not very significant.
Finally, the left panel in Figure 13 shows Model 1, where the center of the beam is exactly at the photon orbit: . An azimuthally oriented ray at this radius is expected to orbit around the BH at a constant . This is seen clearly in the numerical solution. Moreover, since the photon geodesic follows the grid, there is practically no numerical diffusion. Indeed, there is less diffusion than there should be (compare the numerical beam with the two solid green lines). The beam should have some divergence as it propagates around the BH because photons emitted inside curve inward and will ultimately fall into the BH, while those emitted outside curve outward and will move towards infinity. The simulated beam does not reproduce this physical broadening very well.
4.8 Radiative spherical accretion
Our last test problem considers radiative spherical accretion onto a nonrotating BH. This problem has been studied in the past by Vitello (1984) and Nobili et al. (1991) and more recently by Roedig et al. (2012) and Fragile et al. (2012). We follow Fragile et al. (2012) in the setup of our simulations to facilitate comparison with their results. As in their work, we consider Thomson scattering and thermal bremsstrahlung, which give the following opacity coefficients,
(75)  
(76) 
where is in and is the mass of the proton. Our numerical grid spans from to and is resolved by 512 grid points spaced logarithmically following where the auxiliary variable is spaced linearly between values corresponding to and . We assume a BH mass of . For the initial state, we choose the mass accretion rate