{'fraction_non_alphanumeric': 0.04386555160370406, 'fraction_numerical': 0.015027716363869396, 'mean_word_length': 4.249264149356656, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Relativistic Radiation Magnetohydrodynamics in Dynamical Spacetimes: Numerical Methods and Tests

Many systems of current interest in relativistic astrophysics require a knowledge of radiative transfer in a magnetized gas flowing in a strongly-curved, dynamical spacetime. Such systems include coalescing compact binaries containing neutron stars or white dwarfs, disks around merging black holes, core collapse supernovae, collapsars, and gamma-ray burst sources. To model these phenomena, all of which involve general relativity, radiation (photon and/or neutrino), and magnetohydrodynamics, we have developed a general relativistic code capable of evolving MHD fluids and radiation in dynamical spacetimes. Our code solves the coupled Einstein-Maxwell-MHD-Radiation system of equations both in axisymmetry and in full 3+1 dimensions. We evolve the metric by integrating the BSSN equations, and use a conservative, high-resolution shock-capturing scheme to evolve both the MHD and radiation moment equations. In this paper, we implement our scheme for optically thick gases and grey-body opacities. Our code gives accurate results in a suite of tests involving radiating shocks and nonlinear waves propagating in Minkowski spacetime. In addition, to test our code's ability to evolve the relativistic radiation-MHD equations in strong-field dynamical spacetimes, we study"thermal Oppenheimer-Snyder collapse"to a black hole, and find good agreement between analytic and numerical solutions.


I. INTRODUCTION
Many relativistic systems of current astrophysical interest are characterized by the dynamical coupling of strong-field gravitation, high magnetic fields and intense radiation (where the latter may be photons or neutrinos). Quasars, active galactic nuclei (AGNs), galactic "microquasars", core-collapse supernovae, collapsars, gammaray burst sources (GRBs), merging neutron star binaries (NSNSs), merging black hole-neutron star binaries (BHNSs), and merging neutron star-white dwarf binaries (NSWDs) are all examples of such systems. Developing robust computational methods that can treat simultaneously the different dynamical phenomena that govern these systems is necessary in order to simulate their physical behavior reliably and identify their observational signatures.
Many of the systems listed above involve compact objects, such as black holes and neutron stars. Hence general relativity is required to describe their dynamical evolution accurately. Both observations and theory strongly suggest that magnetic fields play an important role in many of these systems. For example, magnetic fields are crucial in launching jets from black holes in AGNs and GRBs (see, e.g., [1,2]), driving accretion onto black holes in disks (see, e.g., [3,4]), and inducing 'delayed' collapse in hypermassive neutron stars that may form following NSNS mergers [5,6,7,8]. Radiation, apart from its role as an observational tracer and diagnostic probe, also can play an important dynamical role in many relativistic * Also at Department of Astronomy & NCSA, University of Illinois at Urbana-Champaign, Urbana, IL 61801 systems. For example, the role of neutrino transport may be essential to understanding core-collapse supernovae (see, e.g. [9,10,11]). As a second example, consider that the interior pressure of supermassive stars and massive Population III stars is dominated by thermal radiation pressure. These objects may collapse in the early universe to form the seeds of the supermassive black holes that reside in the centers of many, and perhaps most, galaxies [12,13]. Radiation thus plays a crucial role in determining the onset and dynamics of the collapse of these stars and the masses and spins of the black holes that are formed [14,15,16]. As a final example, accretion onto compact objects leading to outgoing radiation near and above the Eddington value is controlled by the competition between inward gravitational forces and outward radiation pressure forces. All of these systems need to be handled in a computational scheme designed to probe these physical phenomena self-consistently, simultaneously accounting for radiation, magnetic fields and relativistic gravitation.
We have developed previously a robust numerical scheme in 3+1 dimensions that simultaneously evolves the Einstein equations of general relativity for the gravitational field (metric), the equations of relativistic magnetohydrodynamics (MHD) for the matter, and Maxwell's equations for a magnetic field [17]. Our approach is based on the BSSN (Baumgarte-Shapiro-Shibata-Nakamura) formalism to treat the gravitational field [18,19], a high-resolution, shock-capturing (HRSC) scheme to handle the fluid and a constrained-transport scheme to treat magnetic induction [20]. Our resulting GRMHD code has been subjected to a rigorous suite of numerical tests to check and calibrate its validity [17]. We have applied our code to explore a number of dynamical scenarios, including the collapse of magnetized, dif-ferentially rotating, hypermassive neutron stars to black holes [7,8], the collapse of rotating stellar cores to neutron stars [21], the collapse of rotating, supermassive stars and massive Pop III stars to black holes [22], and the merger of binary black holes [23] and binary black hole-neutron stars [24]. The purpose of this paper is to present a generalization of our current GRMHD scheme that accounts for the presence of radiation (photon or neutrino).
Our approach for handling the radiation follows in the long tradition of formalisms designed to treat radiation transport in the framework of general relativity. However, we have developed a new version specifically designed to fit neatly onto our existing 3 + 1 GRMHD scheme. The general relativistic radiative transfer equation has been derived in full detail by Lindquist in 1966 [25]. His treatment has been followed by numerous adaptations and implementations in various approximations. For example, Thorne has derived a set of radiation moment equations to arbitrary order by the technique of projected symmetric trace-free (PSTF) tensors [26]. So far, most GR radiation hydrodynamics calculations (e.g. [27,28,29,30,31,32,33,34]), including those based on the PSTF scheme (e.g. [35,36,37,38,39,40]), have been implemented in spherical symmetry only. Once spherical symmetry is broken, most radiation schemes become quite difficult to implement, given the large number of phase space degrees of freedom that need to be tracked for the radiation field.
In this paper, we formulate the radiation transport equations in the framework of our 3+1 GRMHD scheme, which operates without any restrictions regarding the spatial symmetry of the system. However, our implementation focuses on the optically thick limit for the radiation field, which simplifies the analysis by allowing us to assume that the radiation field in the comoving frame of the fluid is nearly isotropic. Our emphasis is geared to treating systems in which the radiation has a strong dynamical influence on the matter flow and, in some cases, on the spacetime geometry itself. We are less concerned in this initial treatment with the radiation that escapes from the matter surface, or with the radiation spectrum measured by a distant observer. It is in the interiors of collapsing stars, neutron stars in merging compact binaries, and dense accretion disks orbiting black holes where the dynamical influence of the radiation field is likely to play its most significant role. In these interior regions the optically thick assumption should be quite reliable in many cases. In the implementation presented here we also adopt a grey-body opacity law, which, though simple, suffices to illustrate our method. However, the formalism makes no assumptions regarding the spatial symmetry of the matter source, radiation field, or spacetime.
We present two sets of tests to check our new radiation GRMHD code. The first set of tests involves radiation shocks and nonlinear waves propagating in a fixed Minkowski spacetime. The second set of tests is the "thermal Oppenheimer-Snyder collapse" problem origi-nally proposed and solved by Shapiro [29,30], wherein radiation propagates in a spherical spacetime that, though simple, is highly dynamical and characterized by a strong gravitational field (i.e. one in which a black hole forms). In both sets of tests, we compare our numerical results with analytic solutions and perform convergence tests.
The structure of the paper is as follows: In Sec. II, we formulate the system of coupled Einstein-Maxwell-MHDradiation equations in 3+1 form, with the Maxwell, MHD and radiation equations written in conservative form. In Sec. III, we describe techniques for evolving this system of equations. In Sec. IV, we present the new code tests and their results. Finally, we summarize our results in Sec. V.

A. Evolution of gravitational fields
We write the spacetime metric in the standard 3+1 form: where α, β i , and γ ij are the lapse, shift, and spatial metric, respectively. The extrinsic curvature K ij is defined by where L β is the Lie derivative with respect to β i . The evolution of γ ij and K ij is governed by the Einstein equation G µν = 8πT µν , where G µν is the Einstein tensor and T µν is the stress-energy tensor. We evolve γ ij and K ij using the BSSN formulation [18,19]. The fundamental variables for BSSN evolution are The Einstein equation G µν = 8πT µν gives rise to the evolution equations and constraint equations for these fields, which are summarized in [19]. In this paper, we use the same field evolution equations as Eqs. (11)- (15) of [41]: +αÃ ijÃ ij + 4πα(ρ + S), and where D denotes covariant derivative operator associated with γ ij , and T F denotes the trace-free part of a tensor. The constraint equations, expressed in terms of the BSSN variables, are whereD denotes covariant derivative operator associated withγ ij . The matter-energy source terms are given by ρ = n α n β T αβ , Here n α = (α −1 , −α −1 β i ) is the time-like unit vector normal to the t = constant time slices. In this paper T αβ contains three components: where T αβ (hydro) , T αβ (em) and R αβ are the stress-energy tensor for the hydrodynamic matter field, (large scale) electrodynamic field and the radiation field, respectively. Hence all components here contribute to the BSSN source terms in Eq. (15).
In order to evolve the 3+1 Einstein equations forward in time, one must choose lapse α and shift β i functions, which specify how the spacetime is foliated. The lapse and shift must be chosen in such a way that the total system of evolution equations is stable. In the past few years, we have experimented with several gauge conditions. We find that, in general, the most useful gauge choices are the hyperbolic driver conditions [42,43], and the puncture gauge conditions (see e.g. [44,45]). In this paper, we use the hyperbolic driver conditions as in [43] when evolving a dynamical spacetime: where a 1 , a 2 , a 3 , b 1 , and b 2 are freely specifiable constants.
B. Evolution of radiation fields

Radiation fields
The equations governing the dynamics of the radiation can be expressed as where R αβ is the radiation stress-energy tensor, and G α is the radiation four-force density which describes the interaction of the matter with the radiation [30,46]. The radiation stress-energy tensor R αβ is defined as where ν is the frequency, I ν = I(x α ; N i , ν) is the specific intensity of radiation at x α moving in the direction N α ≡ p α /hν, p α is the photon 4-momentum, h is the Planck constant, and dΩ is the differential solid angle. Here ν, I ν and dΩ are all measured in the local Lorentz frame of a fiducial observer with 4-velocity u α (fid) , i.e. hν = −p α u α (fid) . The integral is evaluated over all frequency and solid angles.
We now choose our fiducial observer to be comoving with the fluid. In the comoving frame of the fluid the radiation stress-energy tensor R αβ takes the form is the comoving radiation energy density, is the comoving radiation flux, and is the comoving radiation stress tensor. We are interested in the optically thick regime, in which the radiation is very nearly isotropic in the comoving frame of the fluid. In the limit of strict isotropy, independent of the propagation direction Nî, the intensity is I ν = I(x α ; ν). Using this fact and the expression of N α in the comoving frame one can show that Fî = 0 and Pî = 1 3 δîE ≡ δîP, where P is the radiation pressure, θ is the polar angle measured from theẑ-axis, and ϕ is the azimuthal angle (i.e. tan ϕ = Nŷ/Nx) . Henceforth, we include the effect of a small anisotropy by allowing a small non-zero radiation flux Fî, but we retain the closure relation P = E/3. That is, we adopt an Eddington factor equal to 1/3.
The radiation stress-energy tensor R αβ can be written in covariant form as where u α is the fluid 4-velocity. This expression reduces to the same form as Eq. (21) in the comoving frame. Here we have introduced the projection tensor, h αβ , defined as and the radiation flux four-vector defined as Note that with this definition, the flux satisfies Following [30], the radiation four-force density is given by where χ ν = χ a ν + χ s ν is the total opacity (the superscript a and s denote the absorption and scattering opacities respectively) and η ν = η a ν + η s ν is the total emissivity. By assuming isotropic and coherent scattering, and that the thermal emissivity η a ν and absorption coefficient χ a ν are related by Kirchhoff's law η a ν = χ a ν B ν , we can write, in the fluid comoving frame, where B ν is the intensity in thermal equilibrium (e.g. the Planck function for photons, the analogous Fermi-Dirac function for neutrinos, etc.). We further assume a grey-body form for all opacities, χ ν = κρ 0 , where κ is a frequency independent opacity, and ρ 0 is the rest mass-energy density. Then we may write [30] It is straightforward to express G α in covariant form as Note that the frequency integrated equilibrium intensity B(T ) can be written as where T is the temperature of the fluid, and a R is a constant depending on the type of radiation: for thermal photons it equals the usual radiation constant a; for each flavor of nondegenerate thermal neutrino or antineutrino (chemical potentials = 0) it is (7/16)a; lumping the contributions of all neutrinos and antineutrinos together, it is (7N ν /8)a, where N ν is the number of neutrino flavors which contribute to thermal radiation. We emphasize here that our method allows for situations in which the gas may be out of thermal equilibrium with the radiation (E = 4πB). Our formalism is equivalent to keeping the first two radiative moment equations, and using an Eddington factor to close the set. Our choice of P = 1/3E serves as the necessary closure relation for these equations. We demonstrate in Appendix A 1 that our formalism, while more general, reduces to the relativistic diffusion approximation in a simplifying limit. In this limit, the diffusion approximation transforms the radiation moment evolution equations from a hyperbolic to a parabolic (i.e. diffusion) form, which does not have the same causal structure as the original system of equations. Moreover, the parabolic form is not suitable for implementing the conservative HRSC scheme used to integrate the combined MHD-radiation equations (see Sec. III). In any case, we do not adopt the diffusion approximation here, but treat the original set without simplification.

Radiation evolution
We can decompose the radiation evolution equations given by (19) in a manner analogous to the way we decompose the MHD evolution equations (e.g., see Sec. IIC in [17] and Sec. II D below). The resulting equations are therefore cast in conservative form, as are the MHD evolution equations. Taking the scalar product of Eq. (19) with n α on both sides gives the energy equation where the radiation energy density variableτ is defined asτ (36) and the source terms is Here γ = e 12φ denotes the determinant of the spatial metric γ ij . The spatial components of Eq. (19) give the momentum equation, where the radiation momentum density variable is defined asS

C. Evolution of large-scale electromagnetic fields
The evolution equation for the electromagnetic field in a perfectly conducting MHD fluid (F µν u ν = 0) can be obtained in conservative form by taking the dual of Maxwell's equation F [µν,λ] = 0. One finds where F αβ is the Faraday tensor, and * F αβ = ǫ αβµν F µν /2 is its dual. Using the fact that the magnetic field as measured by a normal observer n α is given by B i = n µ * F µi , the time component of Eq. (40) gives the no-monopole constraint ∂ jB j = 0, whereB j = √ γ B j . The spatial components of Eq. (40) give the magnetic induction equation, which can be written as where v i ≡ u i /u 0 .

D. Evolution of the MHD field
In the MHD limit, T µν (em) can be expressed as is the magnetic field measured by an observer comoving with the fluid. The stress-energy tensor associated with the perfect fluid can be expressed as where ρ 0 is the (baryon) rest-mass density, P is matter pressure, h = 1 + ǫ + P/ρ 0 is the specific enthalpy, and ǫ is the specific internal energy density of the matter. For brevity, we denote Thus, we see that the conservation of the total stressenergy tensor can be written as This can be combined with (19) to give Additionally, we have the continuity equation expressing baryon number conservation, Rewriting Eqs. (47) and (48) in conservative form gives (cf. Section IIC in [17]) where the MHD evolution variables are and the source term s is Note that these evolution variables are very similar to those in [17]. The only difference is that there are new radiative source terms G i and G 0 in the momentum and energy equations (50) and (51), respectively.
To complete the system of equations, it remains only to specify the equation of state (EOS) of the fluid. In this paper, we adopt a Γ-law EOS, where Γ is the adiabatic gas constant. We choose a Γ-law EOS because it simplifies some of the calculations, it is applicable to many cases of interest, and it is a standard choice for demonstrating new computational techniques in the numerical relativity literature. Also, the analytic solutions we are going to use as code tests also use this EOS. Nevertheless, all evolution equations derived in this section apply for any equation of state, and generalization to a more realistic EOS is straightforward. In fact, our code is currently capable of handling the general class of EOSs of the form P = P (ρ 0 , ǫ).
The fluid temperature T is required in the radiation force density term G µ (Eq. (34)). In this paper, we compute it by using the ideal gas law P = nk B T = ρ 0 k B T /m, where n is the baryon number density, m = ρ 0 /n is the mean mass of the baryons in the fluid, and k B is Boltzmann's constant. Hereafter we set k B = 1.
We point out that the fluid flow is nonadiabatic in general. In particular, there is energy exchange between the matter and radiation fields. Also, shocks may be present in some applications.

E. Summary of equations
To reiterate, the system of coupled Einstein-radiation-Maxwell-MHD equations we consider are the BSSN equations (8)- (12), the radiation transport equations (35) and (38), the magnetic induction equation (41), and the MHD equations (49)- (51). In Appendix A, we demonstrate that our equations reduce to the more familiar Newtonian form in the weak-field, slow-velocity limit. The evolution variables are φ,γ ij , K,Ã ij ,Γ i ,τ ,S i ,B i , ρ * ,S i andτ . These variables are not completely independent: the BSSN variables φ,γ ij , K,Ã ij , andΓ i have to satisfy the Hamiltonian constraint (13) and the momentum constraint (14); the magnetic field variablesB i have to satisfy the no-monopole constraint ∂ iB i = 0. The total stress-energy tensor T µν is given by The BSSN matter-energy source terms [Eq. (15)] can be expressed as

III. IMPLEMENTATION
We use a cell-centered Cartesian grid in our threedimensional simulations. Sometimes, symmetries can be invoked to reduce the integration domain. For octant symmetric systems, we evolve only the upper octant; for equatorially symmetric systems, we evolve only the upper half-plane. For axisymmetric systems, we evolve only the x-z plane (a 2+1 dimensional problem). In axisymmetric evolutions, we adopt the Cartoon method [47] for evolving the BSSN equations, and use a cylindrical grid for evolving the induction, MHD, and radiation equations [48].
Our code uses the Cactus parallelization framework [49], with the time-stepping algorithm based on the MoL, or method of lines, thorn. In the metric evolution (BSSN sector), spatial derivatives can be calculated using second-order or fourth-order finite differencing schemes. The Cactus MoL thorn allows us to switch to a higher-order time-stepping scheme easily. Higher order schemes are very useful for evolving spacetimes containing black holes using the moving puncture techniques (see, e.g., [44,45]). However, we do not treat puncture black holes here and we are currently using a HRSC scheme which is at most second-order accurate to evolve the Maxwell, radiation, and MHD equations. Hence we use second-order finite differencing scheme in the BSSN sector and (second-order) iterated Crank-Nicholson timestepping in our calculations.
Our technique for metric evolution is described in our earlier papers [23,41,50], so we focus here on our MHD, induction and radiation algorithms. The goal of this part of the numerical evolution is to determine the fundamental "primitive" variables P ≡ (ρ 0 , P, v i , B i , E, F i ) at future times, given initial values of P. The evolution equations (35), (38), (41), (49), (50) and (51) are written in conservative form: where the evolution variables U(P) ≡ (ρ * ,τ ,S i ,B i ,τ ,S i ), the fluxes F(P) and the sources S(P) are not explicit functions of derivatives of the primitive variables, although they are explicit functions of the metric and its derivatives. As mentioned above, we evolve Eq. (61) using the iterated Crank-Nicholson scheme. This scheme is second order in time and will be stable if ∆t < min(∆x i )/c max , where in our case c max is the speed of light. For each Crank-Nicholson substep, we first update the gravitational field variables (the BSSN variables). We then update the electromagnetic fields B i by integrating the induction equation. Next, the MHD variables (ρ ⋆ ,τ , andS i ) are updated. Then we update the radiation variables (τ , andS i ). Finally, we use these updated values to recover the primitive variables on the new timestep. Below, we briefly summarize some of the important techniques we utilize during the evolution.

A. Reconstruction step
We implement an approximate Riemann solver to handle the advection in Eq. (61). For simplicity, we consider the one-dimensional case here. The generalization to multi-dimension is straightforward. The first step in calculating this flux is to compute P L = P i+1/2−ǫ and P R = P i+1/2+ǫ , i.e. the primitive variables to the left and right of the grid cell interface. As in [17], we use the Monotonized central (MC) scheme [51] to compute the primitive variables at the cell interface. This scheme is second-order accurate at most points when the data are smooth, but becomes first-order accurate across a discontinuity (e.g. shock). (See [17] for other reconstruction methods in our MHD code.)

B. Riemann solver step
Next, we take the reconstructed data as initial data for a piecewise constant Riemann problem, with P = P L on the left of the interface, and P = P R on the right of the interface. The net flux at the cell interface is given by the solution to this Riemann problem.
We use the HLL (Harten, Lax, and van Leer) approximate Riemann solver [52]. Our implementation has been described in [17]. To summarize, HLL fluxes are given by Here where c + is the maximum right-going wave speed and c − is the maximum left-going wave speed. We obtain c ± by solving the dispersion relation for waves with wave vectors of the form The wave speed is simply the phase speed ω/k 1 . We find the dispersion relation in the comoving frame of the fluid (denoted by subscript cm), and hence ω cm /k cm , as described in Appendix B. To obtain ω/k 1 in the grid frame, we use the dispersion relation (B4) and substitute the ω cm = −k µ u µ , and k 2 cm = K µ K µ , where K µ = (g µν + u µ u ν )k ν . Wave speeds in the y-direction and z-direction are found analogously.

C. Recovery of primitive variables
Having computed U at the new timestep, we must use these values to recover P, the primitive variables on the new time level. We can recover the hydrodynamics primitive variables ρ 0 , P, v i from the MHD evolution variables ρ * ,τ ,S i numerically, as described in Section III C in [17]. Once the fluid velocity v i is found, the radiation primitive variables E and F i can be computed from the radiation evolution variablesτ andS i analytically using Eqs. (36) and (39). We solve the following set of two coupled linear equations to recover E and F 0 : The first one is just Eq. (36), while the second one is obtained by using u µ F µ = 0 to eliminate F i in Eq. (39). After solving for E and F 0 , we compute F i by which is derived by raising the index of F i in Eq. (39).

D. Constrained Transport
The Maxwell equation demands that the magnetic fieldsB i satisfy the no-monopole constraint ∂ iB i = 0. Unphysical behavior may arise if this constraint is violated. Thus, "constrained transport schemes" have been designed to evolve the induction equation while maintaining ∂ iB i = 0 to roundoff precision [53]. We use the flux-interpolated constrained transport (flux-CT) scheme introduced by Tóth [20] and used by Gammie et al [54]. This scheme involves replacing the induction equation flux computed at each point with linear combinations of the fluxes computed at that point and neighboring points. The combination assures both that second-order accuracy is maintained, and a particular finite-difference representation of ∂ iB i = 0 is enforced to machine precision.
E. Low-density regions and boundary conditions

Low-density regions
If vacuum exists anywhere in our computational domain, the MHD approximation will not apply in this region, and we will have to solve the vacuum Maxwell equations there (see e.g. [55]). In addition, the optically thick assumption on the radiation field in Sec. II B also breaks down in sufficiently low density regions. In many astrophysical scenarios, however, a sufficiently dense, ionized plasma will exist outside the stars or disks, whereby MHD will remain valid in its force-free limit. A similar situation may arise for the radiation field, where the ambient gas in our computational domain may be sufficiently dense to maintain an optical depth above unity. However, in some applications we may need to take into account the transition from optically thick to optically thin limits in the low-density regions, depending on the magnitude of the opacity. A precise treatment of this problem requires solving the full Boltzmann radiative transfer equation (see e.g. [30,56,57,58]; See also [59,60] for approximation schemes.) In this paper, however, we avoid this issue. For the code tests that do not have low density regions (Sec. IV A), no special treatment is required. We do, however, present a test involving Oppenheimer-Snyder collapse (Sec. IV B) where there is a vacuum outside the star. As in many hydrodynamics simulations in astrophysics, we impose a low-density "atmosphere" outside the star to facilitate the integration of hydrodynamics equations. It turns out that our atmosphere scheme suffices to mimic the correct ("zero temperature") radiation boundary conditions that we wish to impose at the surface of the star (see Sec. IV B).
In the low-density regions near the surface of the star, we sometimes encounter problems when recovering the primitive variables; in particular, the equations U = U(P) occasionally have no physical solution. Usually, unphysical U are those values corresponding to negative pressure. As in [17], we apply a fix at these points, first suggested by Font et al [61]. In the system of equations (52)-(54) to be solved, we replace Eq. (54) with the adiabatic relation P = κρ Γ 0 , where κ is set equal to its initial value. This substitution guarantees a positive pressure. Typically, these low density regions have little influence on the dynamical evolution of the system, which is the principal target of our current investigations.

Boundary conditions
For the code test in Sec. IV A, we evolve onedimensional, radiation-hydrodynamics equations in a fixed, Minkowski spacetime. We impose the "copy" boundary condition on all the evolution variables, i.e. variables at the boundaries are copied from the closest grid point.
For the Oppenheimer-Snyder code test in Sec. IV B, we evolve the system of coupled Einstein-radiationhydrodynamics equations. In this case, we employ Sommerfeld outgoing wave boundary conditions for all BSSN and gauge variables: where ∆T is the timestep and ∆r = αe −2φ ∆T . For the radiation hydrodynamics, we impose the outflow boundary condition on the primitive variables ρ, P , v i , E, and F i (i.e., the variables are copied along the grid directions with the condition that the velocities be positive or zero in the outer grid zones). We note that the radiation in this test is initially confined inside the star, but escapes from the stellar surface during the evolution. In the end of the simulation, the total emitted radiation remains small and the dynamics of the system is insensitive to the boundary condition employed.

IV. CODE TESTS
Our GRMHD code has previously been thoroughly tested by maintaining stable rotating stars in stationary equilibrium, by reproducing Oppenheimer-Snyder collapse to a black hole, and by reproducing analytic solutions involving MHD shocks, nonlinear MHD wave propagation, magnetized Bondi accretion, and MHD waves induced by linear gravitational waves [17]. It has also been compared with the GRMHD code of Shibata and Sekiguchi [62] by performing simulations of the evolution of magnetized hypermassive neutron stars [7,8], and of magnetorotational collapse of stellar cores [21]. We obtain good agreement between these two independent codes. Our code has also been used to study the evolution of BHBH and BHNS binaries [24], and the evolution of relativistic hydrodynamic matter in the presence of puncture black holes [63]. Here we restrict our attention to testing the new radiation-hydrodynamics sector, setting large-scale magnetic fields to zero. We also choose the grey body absorption opacity κ a to be a constant and set the scattering opacity κ s to zero.

A. Minkowski radiation-hydrodynamics tests
We present here a series of tests of nonlinear radiationhydrodynamic waves in Minkowski spacetime with planar symmetry. These tests are summarized in Table I. For our initial data, we generate semi-analytic, stationary configurations using the method outlined in Appendix C.
To test the ability of our code to handle shocks and waves moving across the grid, we boost the stationary solutions derived in Appendix C for tests 2-4. In each case, our computational domain is x ∈ (−20, 20). We choose the opacities in each case to ensure that the grid boundaries at x = ±20 reside in the asymptotic region where all  hydrodynamic and radiation quantities approach their asymptotic values, and that the total optical depth across the grid is τ ∼ 10 (see Appendix C). We evolve the system with a timestep ∆t = ∆x for test 1 and test 2, and ∆t = 0.1∆x for test 3 and test 4. We use resolutions ranging from ∆x = 0.0125 to ∆x = 0.1 in order to perform convergence tests. To demonstrate convergence, we consider a grid function g with error δg = g − g exact . We calculate the L1 norm of δg (the "average" of δg) by summing over every grid point i: where N ∝ 1/∆x is the number of grid points. We find that for our continuous configurations (tests 3 and 4), we achieve second order convergence. Because our shock capturing scheme becomes 1st order when discontinuities are present, we achieve the expected first order convergence for our discontinuous configurations (tests 1 and 2). The initial configurations listed in Table I are chosen to test our code in a variety of regimes, including gaspressure dominated, radiation-pressure dominated, Newtonian, relativistic, continuous, and discontinuous matter and radiation profiles. In each test, the equation of state of the gas is given by a Γ-law EOS. We choose Γ = 5/3 for each test except our highly-relativistic case, in which we choose Γ = 2. This latter choice is adopted because the sound speed (c s = ΓP/ρ 0 h) for a Γ-law EOS is limited by c s < √ Γ − 1. Highly-relativistic sound speeds (c s → 1) can only be achieved for Γ ≥ 2. Below, we provide a brief description of each test.
1. Nonrelativistic strong shock. For this test, we set up a strong, gas-pressure dominated, Newtonian (u x max = 0.015 ≪ 1) shock propagating into a cold gas. We have chosen to simulate this scenario because it can be compared to the analytic solution for a subcritical radiating shock first derived by Zel'dovich and Raizer [64] and summarized in [46]. We find very good agreement with this analytic result (see Fig. 1). We note that the radiative shock junction conditions (see Appendix C) require that R 0x and R 00 be continuous at the shock front, even though E and F x are, in general, discontinuous at the shock. In the Newtonian limit, however, the continuity of R 0x and R 00 is equivalent to the continuity of E and F x .
2. Mildly-relativistic strong shock. In this test, we set up a mildly-relativistic (u x max = 0.25), gas-pressure dominated shock. In this case, we see that E and F x no longer appear continuous. We boost this shock so that the shock speed is µ = 0.1. We find that the discontinuity is able to retain its shape very well as the shock travels and matches very well with the analytic solution (see Fig. 2).
3. Highly-relativistic wave. In this test, we simulate a highly-relativistic (u x max = 10), gas-pressure dominated configuration in which all quantities are continuous, but asymptote to different values on either side of the computational domain. We boost this configuration so that it travels across the grid with velocity µ = 0.8. The numerical results agree very well with semi-analytic solution (see Fig. 3). Figure 4 shows the L1 norms of the errors in F x , E, v x , P and ρ 0 at t = t final = 20. We find that all errors converge to zero at second order in ∆x.

Radiation-pressure dominated, mildly-relativistic
wave. In this test, we study the performance of our code in the radiation-pressure dominated (P ≪ P), mildly-relativistic (u x max = 0.69) regime. We boost this configuration so that it travels across the grid with velocity µ = 0.1. The numerical results again agree with the semi-analytic solution (see Fig. 5).

B. Dynamical Spacetime test: Thermal Oppenheimer-Snyder Collapse
The collapse from rest of a homogeneous dust ball (P = 0) in general relativity can be described by the analytic Oppenheimer-Snyder solution [65]. The collapse results in the formation of a Schwarzschild black hole. The evolution of thermal radiation within the dust ball has been considered by Shapiro in [29,30]. In both papers, the radiation is assumed to be a small perturbation, so that the dynamics are unaffected by the presence of radiation, and the matter and metric profiles can still be described by the Oppenheimer-Snyder solution. The first paper employs the relativistic thermal diffusion approximation for the radiation, and derives analytic solutions for both the Newtonian and general relativistic cases. In the second paper this approximation is removed and replaced by solving the exact radiative transfer (Boltzmann) equation for the intensity, coupled to the radiation moment equations for the radiation flux and energy density. It is found that the results obtained by solving the Boltzmann transport equation agree very well with the analytic solutions assuming diffusion approximations, provided the optical depth of the star is sufficiently large (≫ 1). Here we perform a numerical simulation of "thermal Oppenheimer-Snyder collapse" using our radiation GRMHD code, and compare our results to the analytic solution in the diffusion approximation limit given in [29,30]. For convenience, we summarize the analytic solution in Appendix D.
For our initial data, the areal radius of the star is set to be R i = 3M , where M is the ADM mass of the star. We choose the initial profiles for all hydrodynamic and radiation quantities to be homogeneous throughout the star, in compliance with the analytic solution in Appendix D. The analytic solution assumes that (1) the matter and radiation pressure is small enough to be dynamically unimportant (i.e. P/ρ 0 ≪ M/R and P/ρ 0 = E/3ρ 0 ≪ M/R), (2) radiation pressure dominates over gas pressure (P ≫ P ), (3) gas and radiation are in local thermal equilibrium (LTE), i.e. E = 4πB = a R T 4 , and (4) the star is optically thick. To satisfy these conditions, we choose the following initial data: ρ 0 = M/( 4 3 πR 3 i ), P = 10 −4 ρ 0 , E = 10 −3 ρ 0 , v i = 0, F i = 0. LTE is achieved in the initial data by fixing the constant a R m 4 = m 4 E/T 4 = E(ρ 0 /P ) 4 = 10 13 M/( 4 3 πR 3 i ). We note that in our formalism the system is allowed to deviate from the LTE during the evolution. However, we find that the system remains close to the LTE during the entire evolution and the numerical data agree well with the analytic solution (see below). We choose κ a so that τ a = κ a ρR = 50 initially. This guarantees that the star is optically thick initially. As the collapse proceeds, the optical depth increases as τ ∝ 1/R 2 , so the star remains optically thick. We construct the initial data for the spatial metric by transforming the analytic Oppenheimer-Snyder metric from Friedmann to isotropic coordinates, following the procedure described in [66]. We use the analytic solution only at t=0. The metric at later times is evolved, together with hydrodynamics and radiation. The lapse and shift are determined by the hyperbolic driver conditions (Eqs. (17) and (18)). These are gauge conditions that have been widely used in stellar collapse calculations using the BSSN scheme. We choose a 1 = 0.75, b 1 = 0.15, a 2 = b 2 = 2M −1 , a 3 = 1. A smaller b 1 prevents "blowing out" of the coordinate system, a well-known effect [41,67] which can spoil grid resolution in the center of the collapsing object. We perform our numerical simulation in axisymmetry, with 200 2 , 400 2 , 800 2 and 1600 2 grid points. We choose ∆t = 0.1∆x in these simulations. The outer boundary is placed at 4M in isotropic coordinates (R out = 5.06M in areal radius). Note that we do not impose any special boundary condition at the stellar surface, in contrast to the zero temperature boundary condition (E = 0) used in the derivation of the analytic solution [29]. The low density region outside the star mimics this surface boundary condition, as the atmosphere is made to be much colder than the interior of the star, and hence the thermal emission and build-up of radiation energy density in the atmosphere is negligible.
The analytic solution given in Appendix D is expressed in Friedmann coordinates (i.e. Gaussian normal coordinates comoving with the fluid), which is equivalent to using the gauge conditions α = 1 and β i = 0 (geodesic slicing and zero shift), which are different from the gauge conditions we adopt in our numerical simulations. In order to compare our numerical result to the analytic solution, we perform a mapping between these two different gauges. This is achieved first by following a set of Lagrangian fluid elements inside the star, and calculating the proper time and position of these elements by integrating the equations Next, we use the metric and the positions of the fluid elements to compute their areal radii r s . Finally, knowing the proper times τ and areal radii r s of the fluid elements, we use Eqs. (D2) and (D3) and r s = a(τ ) sin χ to compute their Friedmann coordinates (τ, χ). The mapping between these two gauges is thus established. The pair (τ j , χ j ) for each element j uniquely specifies the fluid and radiation parameters. Figures 6 and 7 show the profiles of ρ 0 , P , E and F at different times during the collapse. Note that while the density remains spatially constant during the collapse in comoving Friedmann coordinates, this is not true in our gauge. We see that the numerical results agree very well with the analytical solution, even after the apparent horizon appears at t = 6.78M and all of the stellar material is inside the horizon.
We next perform a convergence test for the radiation quantities. We find that our numerical data converge to a solution slightly different from the analytic solution. This can be explained by the fact that the analytic solution is strictly valid only in the perturbative limits that P/ρ 0 → 0 and P/P → 0. While we set up our initial data to approximate these limits, the slight deviation from the analytic solution is still detectable with our resolutions. In the absence of radiation (E = 0, F i = 0), we have checked that the deviation in ρ 0 between numerical and analytic values is reduced by a factor of 10 if we reduce the ratio P/ρ 0 by a factor of 10. In the presence of radiation, however, decreasing the ratios P/P and P/ρ 0 arbitrarily small makes the numerical simulations quite challenging, as accurate evolution for the radiation quantities E and F i requires accurate evolution of the temperature T ∝ P/ρ 0 , which in turn requires accurate determination of the pressure P . However, accurate computation of P in the limit P/ρ 0 → 0 is difficult. Since the evolution variables are dominated by the rest mass density ρ 0 , in order to recover the tiny P accurately from them, the numerical evolution has to be very accurate. This requires very high resolution. Thus, we perform a convergence test in which we compare numerical solutions with small but finite P/P and P/ρ 0 for different resolutions, rather than comparing with the analytic solution. Figure 8 shows the result of the convergence test for E and F , with the differences scaled for second order convergence. We follow a Lagrangian fluid element halfway between the center and the surface of the star (in terms of areal radius), determine the radiation parameters at Resolutions with 200 2 , 400 2 , 800 2 and 1600 2 grids are used here. The differences between lower and the highest solution (1600 2 ) are rescaled to demonstrate second-order convergence.
the position of the element versus time, and subtract the numerical results from different resolutions. Since we use a HRSC scheme that is second-order accurate except at positions where discontinuities appear, e.g. at the surface of the star where the density falls abruptly, we expect that the order of convergence depends on how much the physical quantity is affected by the discontinuity at the stellar surface. In principle, the first-order error will propagate everywhere inside the star, but its effect may be small (depending on the physical quantity under consideration) and may only be detectable at very high resolution. Fig. 8 shows that E converges at second order, whereas F converges at less than second order but better than first order. This shows that F is more susceptible to propagation of the first order behavior at the surface of the star, which can be anticipated by looking at the shape of the profiles in Fig. 7: E drops abruptly before reaching the surface, while F increases monotonically up to the surface. We see that convergence of F deviates further from second order as the resolution is increased. This is consistent with the presence of a first order term with small coefficient due to the discontinuity at the stellar surface: Here F (∆, t) is the value of F at time t evolved with a grid size ∆, F exact (t) is the exact solution, and c 1 (t) and c 2 (t) are resolution-independent functions. The firstorder term c 1 (t)∆ results from the discontinuity. We expect that c 1 (t) ≪ c 2 (t) since we look at a point far away from the discontinuity. With a lower resolution, and hence a larger grid size ∆, the first-order term is not as significant relative to the second order term because of the small coefficient. Once we decrease ∆ by increasing the resolution, the second-order term diminishes as ∆ 2 , while the first order term shrinks as ∆ only, making the first order term more conspicuous. Similar behavior for E is not evident since E drops to a very low value at the surface.

V. CONCLUSIONS
We have developed a code which can evolve the coupled Einstein-Maxwell-MHD-Radiation equations in 3 + 1 dimensions. In the implementation presented here this code is able to model the behavior of magnetized, perfectly conducting, radiating fluids in dynamical spacetimes in which the characteristic length scales of the system are longer than the mean free path of the radiation and the opacity has a grey-body form. Our formalism allows us to evolve the radiation fields using a HRSC scheme which is analogous to the method we use for the hydrodynamic fields. In this paper, we have tested the shock-capturing capabilities of our code by simulating both continuous and discontinuous onedimensional radiating hydrodynamic waves. We have been successful in evolving highly relativistic, radiationpressure dominated, gas-pressure dominated, and Newtonian waves. We have treated both stationary waves and boosted waves that propagate across our computational grid in Minkowski spacetime. We have also confirmed our ability to accurately capture the behavior of radiation in a strong-field dynamical spacetime by simulating a thermal Oppenheimer-Snyder collapse. Our numerical results agree well with the analytic solutions. We perform convergence tests on our test problems and find the expected order of convergence in all cases.
We plan to use our radiation GRMHD code to study many interesting systems, and revisit some of the problems we have considered before, such as core-collapse supernovae, accretion onto a black hole, merging NSNSs and BHNSs, etc. By taking into account the effect of radiation, we hope to gain more insights and provide some answers to questions that are relevant to observations. For example, we can calculate the radiation luminosity, and study the radiation feedback to the dynamics of the systems.

VI. ACKNOWLEDGMENTS
It is a pleasure to thank B. Stephens and Z. Etienne for useful suggestions and discussions. Some of these calculations were performed at the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign (UIUC). This paper was supported in part by NSF Grants PHY02-05155, PHY03-45151 and PHY06-50377 as well as NASA Grants NNG04GK54G and NNX07AG96G at UIUC. APPENDIX A: LIMITS

Diffusion Approximation
Using the form of R αβ given by Eq. (26) in the radiation moment equation R αβ ;β = −G α gives where P = E/3. In our formalism, we have already assumed near isotropy, which implies that F α /E ≪ 1. We can make the further approximation that the radiation flux F α can be neglected in the left hand side of Eq. (A1) to get Operating on both sides of the above equation with the projection tensor h γ α , using Eq. (33) for G α and the fact that F α u α = 0, we get where κ ≡ κ a + κ s is the total opacity, a α ≡ u α ;β u β is the 4-acceleration, and where we have used the fact that h γ α u α = 0 and P = E/3. Thus, we arrive at the expression for the radiation flux in the diffusion approximation, This is the relativistic diffusion equation relating the radiation flux F α to the local radiation energy density E.
If we further assume LTE, then E = a R T 4 , in which case, where This familiar result is in agreement with Eq. (3.2) of [29], and with Eq. (2.5.28) of [68]. For a further discussion of the simplification in Eq. (A1) leading to the diffusion approximation, see [30].

Newtonian Limit
It is instructive to consider the weak-field, slowvelocity (Newtonian) limit of the GR radiation hydrodynamic equations, and show that our equations reduce to the familiar expressions of Newtonian radiation hydrodynamics. For simplicity, we set all large-scale electromagnetic fields to zero (i.e. T µν (em) = 0 = B i ).

a. Continuity equation
From Eq. (48), we have (ρ 0 u ν ) ;ν = 0 (A7) In Newtonian limit, the covariant derivative reduces to partial derivative (in Cartesian coordinates), and the 4vector u α reduces to u α ≈ (1, v i ). Hence the continuity equation reduces to the familiar expression: In Newtonian limit, the metric can be approximated by where Φ ≪ 1 is the Newtonian gravitational potential. Keeping only the lowest order terms, we obtain Using Eq. (57) and assuming ρ 0 ≫ P , ρ 0 ≫ E, and Combining Eq. (A12) with the continuity equation (A8) yields which is the familiar Newtonian Euler equation, allowing for gas plus radiation pressure, (P + P).

c. Energy equation
The energy equation in the Newtonian limit is derived by contracting u µ with the equation ∇ ν T µν = 0. Using Eq. (46) and the continuity equation ∇ µ (ρ 0 u µ ) = 0, we find, after some algebra, that where a µ = u ν ∇ ν u µ is the 4-acceleration. We note that u µ a µ = 0 and F µ u µ = 0, which implies that so we may neglect F µ a µ in the Newtonian limit. Also, we may neglect this term as well.
Hence in the Newtonian limit, we obtain This can be identified as the Newtonian energy equation for the coupled fluid.
Equations (A14) and (A17) can be expressed in more familiar forms by introducing the total (Lagrangian) time derivative comoving with the fluid: The continuity equation ∇ µ (ρ 0 u µ ) = 0 gives dρ 0 /dτ = −ρ 0 ∇ µ u µ . Hence, Combining Eqs. (A14), (A18) and (A19) yields Dividing both sides of Eq. (A20) by ρ 0 and writing E tot = ρ 0 ǫ + E and P tot = P + P, we get d dτ (A21) In Newtonian limit, d/dτ → d/dt, and we may neglect F µ a µ and ∂ t F 0 as explained above, and Eq. (A21) reduces to d dt which is the familiar first-law of thermodynamics in the case where entropy is generated by radiation.

d. Radiation equations
From Eqs. (38) and (35), we get In Newtonian limit,S Inserting these into (A23) and (A24), and dropping higher order terms, we obtain In order to compute HLL fluxes, we must compute the maximum left-going wave speed c − and maximum right-going wave speed c + on both sides of the interface. We estimate the wave speed by computing the dispersion relation due to a small perturbation on a magnetized, radiating plasma of uniform ρ 0 , P , B i , E and F i . In the comoving frame, u i = 0. We further choose our coordinate system so that the spacetime is locally Minkowski (i.e. g µν = η µν ). To compute the dispersion relation, we consider a perturbation of the form ρ 0 =ρ 0 + δρ 0 e i(kcm·x−ωcmt) , P =P + δP e i(kcm·x−ωcmt) , where bar denotes the unperturbed quantity, and the subscript "cm" refers to comoving frame values. Substituting these 12 expressions into our 12 radiation-MHD equations (35), (38), (49), (50) and (51) and keeping terms linear in the perturbation, we obtain a matrix equation of the form M X = 0. Here M is a 12 × 12 matrix, and X = (δρ 0 δP δu i δB i δE δF i ) t , where superscript t denotes the transpose. For simplicity, we drop the radiative source terms G µ in deriving the dispersion relation (the curvature source terms vanish in Minkowski spacetime). The dispersion relation is obtained by setting det(M ) = 0, which leads, after some algebra, to the following equation: Here . The solution ω cm = 0 corresponds to pure, stationary density perturbations, ω 2 cm = k 2 cm /3 is related to the propagation speed of the radiation flux for the nearly isotropic radiative diffusion, ω 2 cm = (k cm ·v A ) 2 corresponds to the Alfvén waves, and Q(ω cm , k cm ) = 0 corresponds to magnetosonic waves. As in [54], we replace the dispersion relation Q(ω cm , k cm ) = 0 by ω 2 cm − c 2 m k 2 cm = 0 as it is more convenient when calculating the characteristic speed in the grid frame. As pointed out in [54], this modified dispersion relation overestimates the maximum wave speed by a factor of ≤ 2 in the comoving frame.
Since the HLL scheme only requires the information on the maximum and minimum characteristic speeds, we use the following dispersion relation to estimate the characteristic speeds: .

(B4)
APPENDIX C: ANALYTIC SOLUTIONS FOR RADIATION TESTS IN MINKOWSKI SPACETIME

One-Dimensional waves in Minkowski Spacetime
We begin by assuming a Γ-law EOS and write Γ = 1 + 1/n, where n is the polytropic index. Hence the EOS (56) becomes ρ 0 ǫ = nP .
We consider a stationary, infinite fluid in Minkowski spacetime with planar symmetry, hence we drop all time derivatives, all y and z components, and all y and z derivatives. It follows from ∇ µ (ρ 0 u µ ) = 0, ∇ ν T µν = 0 and ∇ ν R µν = 0 that It is convenient to define The system of ODEs in (C1)-(C5) becomes ∂ x U(P) = S(P) .
The first three equations [Eqs. (C1)-(C3)] are readily integrated, giving three "conserved quantities" U 1 , U 2 , and U 3 . In the presence of a shock, across which P is discontinuous, these "conserved quantities" give the Rankine-Hugoniot conditions, relating P 's on both sides of the shock. The remaining two ODEs [Eqs. (C4) and (C5)] can be integrated numerically using, for example, a fourth-order Runge-Kutta integrator. During the integration, we need to compute P from U , which we outline as follows. For simplicity, we consider the case without a large-scale electromagnetic field (T αβ (em) = 0 = B i ). Using Eq. (57) for T µν , Eq. (26) for R µν , and combining Eqs. (C1)-(C5), we obtain Eliminating ρ 0 and P from Eqs. (C6)-(C8), we get an expression for u x : Hence u x can be determined by solving the above algebraic equation. Having obtained u x , the quantities ρ 0 and P are computed from The values of E and F x are obtained by solving the linear Eqs. (C9) and (C10). The result is where We note that Eq. (C11) will, in general, have more than one real root. Indeed, this must be the case in order for shocks to exist. In the absence of a shock, the appropriate root is chosen by continuity.
Our one-dimensional tests can be divided into two groups: fully continuous configurations, and discontinuous configurations (i.e. shocks are present). In either case, we begin by specifying boundary conditions on the asymptotic left side (x = −∞). In practice, we set up our computational domain with x ∈ [−L, L]. We specify ρ 0 , P , and u x at x = −L, denoting them as ρ 0L , P L and u x L . We also impose that the radiation be in LTE with the gas at x = −L [E L = a R T 4 L = a R m 4 (P L /ρ 0L ) 4 ] and set F x L = f L , where f L is a parameter chosen so that F x L /E L ≪ 1. Thus, we have specified 5 boundary conditions for our 5 ODE's. The values of U 1 , U 2 and U 3 , which are independent of x, are determined. If all quantities in the configuration are continuous throughout the computational domain, these boundary conditions at x = −L are sufficient for us to integrate Eqs. (C9) and (C10) from x = −L to x = L.
If however, we wish to determine a configuration that contains a discontinuity, which we set at x = 0, we use a shooting method, described as follows. In addition to the boundary conditions set at x = −L, we also demand that the radiation be in LTE with the gas, and that F x = f R at x = L. The parameter f R is chosen so that F x /E ≪ 1 at x = L. The constants U 1 , U 2 and U 3 are fixed, giving the values of P at x = −L. Denote the values of ρ 0 , P , u x and E at x = L by ρ 0R , P R , and u x R , and E R , respectively. The LTE condition at x = L gives E R = a R m 4 (P R /ρ 0R ) 4 . From the definition of U 1 , U 2 and U 3 , it is straightforward to show that where u 0 R = 1 + (u x R ) 2 . Substituting Eq. (C19) into Eq. (C20) gives an algebraic equation for u x R , which can be solved numerically. Obviously, u x R = u x L is a solution, but we look for another solution in order to obtain a configuration containing a shock. Having determined u x R , the other quantities are computed as follows: where P R /ρ 0R in above equations are computed from Eq. (C19). To generate the 1D shock configuration, we first specify ρ L , P L , u x L , and consider f L and f R as free parameters. The other quantities at x = ±L are fixed by the LTE condition at x = ±L and Eqs. (C19)-(C23). We then integrate Eqs. (C4) and (C5) from both x = ±L to x = 0. Since U 1 , U 2 and U 3 are constants, the Rankine-Hugoniot junction conditions are automatically satisfied at the shock front (x = 0). Hence, we only have to impose the junction conditions for the radiation variables, which are the continuity of R 0x and R xx . In the Newtonian limit, these conditions reduce to the continuity of E and F x at the shock front, but this is not the case in general. In any case, we need two junction conditions at the shock front, so a well-posed shooting problem requires varying two boundary condition parameters until the solution satisfies the two junction conditions at the shock front. We use f L and f R as such two parameters.
For a pure hydrodynamic (or MHD) shock, the profiles of P are constants on each side of the shock front. This is not the case for a radiating hydrodynamic shock, where P vary with x and approach constants only in the asymptotic regions (x → ±∞). This variation results from the radiative source terms G 0 and G x , given by Eq. (33), which vanishes when the radiation and fluid are in strict LTE and the radiation flux vanishes (i.e. in the asymptotic regions). The length scale over which the parameters vary between the two asymptotic regions is a few optical depths. We need to choose κ to ensure that x = ±L are in the asymptotic regions. In practice, one can estimate the optical depth τ = L −L ρ 0 κdx ∼ (ρ 0L + ρ 0R )κL and choose κ so that τ ≫ 1. For our tests, we choose τ ∼ 10. Choosing a larger κ does not change the profile (plotted against the rescaled coordinate κx) significantly.

Special Analytic Case
In Newtonian limit, the solutions of (C6)-(C10) can be written in analytic form under special conditions, as stated in [46,64]. Here we briefly summarize this solution.
In Newtonian limit, Eqs. (C6)-(C10) become while by dropping time derivatives in Eqs. (A29)-(A30), we get where v = v x and τ is the optical depth, given by dτ = ρ 0 κdx. As in [46,64], we consider a strong shock propagating into cold gas, so that the pressure and internal energy of the unshocked gas (x < 0) can be neglected. We also assume that the gas is optically thick so that we may use the diffusion approximation, and that it is sufficient to account for the radiation flux, while neglecting radiation energy density and radiation pressure. Under these assumptions, the above equations can be rewritten as: Combining Eqs. (C28) and (C27) gives Solving these coupled equations using methods outlined in [46,64] one arrives at: Preshock Medium (x < 0): Postshock Medium (x > 0): Here T R is the asymptotic temperature in the postshock region and τ is measured from the shock front (i.e. τ (x) = x 0 κρ(x ′ )dx ′ ).

APPENDIX D: THERMAL OPPENHEIMER-SNYDER SOLUTION
Here we summarize the analytic solutions derived in [29] for a strong field (black hole) dynamical scenario, used to compare with our numerical results in Section IV B.
The standard Oppenheimer-Snyder (OS) collapse solution for a homogeneous dust ball was first derived in [65]. In [29], the collapsing sphere is subjected to radiation and gas pressure perturbations which are assumed to be sufficiently small that the spacetime metric and density evolution are well-approximated by the OS solution. For this to be true, we require P/ρ 0 ≪ M/R and P/ρ 0 = E/3ρ 0 ≪ M/R for a star with mass M and radius R. We are interested in the evolution of the radiation quantities E and F i inside the star.
Here η is the conformal time and a m is a constant which is related to the initial areal radius R i of the star. (The subscript i denotes initial values.) The radius is given by The exterior Schwarzschild line element is dr 2 s + r 2 s dΩ 2 (D5) where r s is the Schwarzschild (areal) radius. Matching the two metrics at r s = R gives Denote χ 0 as the Lagrangian radial coordinate χ at the stellar surface. It follows from Eqs. (D1), (D5), (D2) and (D4) that In Friedmann comoving coordinates the density ρ 0 is always homogeneous and given by where Q = a a m = 1 2 (1 + cos η) .
An analytic solution for the E and F is derived in [29,30] by assuming (1) diffusion approximation, (2) that the radiation and fluid are in LTE (E = a R T 4 ), and (3) that radiation pressure is much greater than gas pressure (P ≫ P ). We summarize the solution below.
Define the radiation energy density and flux "corrected" for adiabatic contraction as E c = Q 4 E and F c = Q 4 F respectively. Define a time parameterτ as τ = 1 4(τ a + τ s ) (D10) where τ a and τ s are the initial absorption and scattering optical depths respectively, related to the absorption and scattering opacities κ a and κ s by τ a = κ a ρ 0i R i and τ s = κ s ρ 0i R i . (Note that the analytic solution assumes κ a and κ s to be constant throughout the collapse.) Define also a normalized Lagrangian radius z = χ/χ 0 , such that 0 ≤ z ≤ 1 within the star. Then the interior corrected energy density is given by E c (τ , z) = 2E i sin χ 0 sin(χ 0 z) e χ 2 0τ ∞ n=1 (−1) n+1 e −n 2 π 2τ sin(nπz) nπ n 2 π 2 − χ 2 0 (D11) where E i is the initial value of E, assumed to be constant throughout the star. Note that in the derivation of Eq. (D11) the "zero temperature" boundary condition (E = 0 at the stellar surface) has been used.
The diffusion approximation gives the expression for the corrected radiation flux F c in terms of the gradient of E c : Finally, an analytic expression for the ideal gas pressure P = ρ 0 T /m is obtained from the LTE condition E = a R T 4 = a R m 4 P ρ0 4 . Hence where ρ 0 and E are given by the analytic solutions above.
A comparison between the analytic solution for thermal OS collapse and the exact solution of the Boltzmann equation of radiative transfer for the same problem is given in [30].