Radiation hydrodynamics in spherical symmetry

Many spherically symmetric simulations of compact objects have been approached in comoving orthogonal coordinates Misner_Sharp_64,May_White_66. Finite difference schemes of varying complexity were designed in May_White_67,VanRiper_79,Bruenn_85,Rezzolla_Miller_94,Swesty_95,Liebendoerfer_Rosswog_Thielemann_02, culminating in an approximate Riemann solver Yamada_97. The left-hand side of the Einstein field equation, the Einstein tensor, is based on the metric

$$ds^{2}=-\alpha ^{2}dt^{2}+\left( \frac{r'}{\Gamma }\right) ^... ...eft( d\vartheta ^{2}+\sin ^{2}\vartheta d\varphi ^{2}\right) ,$$ (1)

where is the areal radius and is a label corresponding to an enclosed rest mass (the prime denotes a derivative with respect to : ). The proper time lapse of a comoving observer is related to the coordinate time by the lapse function . We have made the substitution , based on a function , for the space-space component of the metric. The angles and describe a two-sphere. We use natural units such that the velocity of light, , and the gravitational constant, , become .

The right-hand side of the Einstein equations is given by the fluid- and radiation stress-energy tensor, . In a comoving orthonormal basis, it has the components Lindquist_66

$$T^{tt} = \rho \left( 1+e+J\right)$$

$$T^{ta}=T^{at} = H$$

$$T^{aa} = p+\rho K$$

$$T^{\vartheta \vartheta }=T^{\varphi \varphi } = p+\frac{1}{2}\rho \left( J-K\right) .$$ (2)

The total energy is expressed in terms of the rest mass density, , the specific internal fluid energy, , and the specific radiation energy, . The isotropic fluid pressure is denoted by , and the radiation stress is composed from the zeroth () and second () angular moments of the radiation intensity. Radial net energy transport is accounted for by the nondiagonal component of the stress-energy tensor, the first angular moment () of the specific radiation intensity.

We define a velocity , equivalent to the component of the fluid four-velocity as observed from a frame at constant areal radius May_White_67, and identify the total energy enclosed in a sphere with the gravitational mass, . In the special relativistic limit, then becomes the Lorentz factor corresponding to the boost between inertial and comoving observers. As in nonrelativistic hydrodynamics we can define a specific volume, , specific energy, , and specific radial momentum, , by

$$\frac{1}{D} = \frac{\Gamma }{\rho }$$ (3)

$$\tau = \Gamma \left( e+J\right) +\frac{2}{\Gamma +1}\left( \frac{1}{2}u^{2}-\frac{m}{r}\right) +uH$$ (4)

$$S = u\left( 1+e+J\right) +\Gamma H.$$ (5)

It has been shown in Liebendoerfer_Mezzacappa_Thielemann_01 that these definitions lead to conservation equations (6)-(8) that are analogous to the continuity equation, the conservation of total energy, and the conservation of radial momentum :

$$\frac{\partial }{\partial t}\left[ \frac{1}{D}\right] = \frac{\partial }{\partial a}\left[ 4\pi r^{2}\alpha u\right]$$ (6)

$$\frac{\partial \tau }{\partial t} = -\frac{\partial }{\partial a}\left[ 4\pi r^{2}\alpha \left( up+u\rho K+\Gamma \rho H\right) \right]$$ (7)

$$\frac{\partial S}{\partial t} = -\frac{\partial }{\partial a}\left[ 4\pi r^{2}\alpha \left( \Gamma p+\Gamma \rho K+u\rho H\right) \right]$$

$$- \frac{\alpha }{r}\left[ \left( 1+e+\frac{3p}{\rho }+J+3K\right) \frac{m}{r}-\left( 1-\frac{2m}{r}\right) \left( J-3K\right) \right.$$

$$+ \left. 8\pi r^{2}\left( \left( 1+e+J\right) \left( p+\rho K\right) -\rho H^{2}\right) -2\left( \frac{p}{\rho }+K\right) \right]$$ (8)

$$\frac{\partial V}{\partial a} = \frac{1}{D}$$ (9)

$$\frac{\partial m}{\partial a} = 1+\tau$$ (10)

$$\frac{\partial }{\partial t}\left[ \frac{1}{4\pi r^{2}\rho }H\right] = -\left( 1+e+J\right) \frac{\partial \alpha }{\partial a}-\frac{1}... ... \alpha \left( p+\rho K\right) \right] +\frac{\alpha }{3VD}\left( J-3K\right) .$$ (11)

The change of the specific volume in Eq. (6) is given by the balance in the displacement of the zone boundaries. The rate of change of total energy in Eq. (7) is determined by the surface luminosity, , and the work on the surface of the mass shell against the pressure, . Of leading order in the momentum equation (8) are the pressure gradient and the gravitational force, . The constraints (9) and (10) are most easily understood in the Newtonian limit (the enclosed volume is defined by ), where the first becomes the definition of the rest mass density and the second the Poisson equation for the gravitational potential. The time derivative in equation (11) is very small; therefore, this equation essentially acts as a constraint on the lapse function, . This equation derives from the space component of the four-divergence of the stress-energy tensor. In addition to the evolution of the total energy, we also need an equation for the evolution of the internal energy that we may derive from the time component of the four-divergence of the stress-energy tensor:

$$\frac{\partial }{\partial t}\left[ e+J\right] =-\frac{1}{\al... ...\frac{1}{\rho }\right) -\frac{\alpha u}{r}\left( J-3K\right) .$$ (12)

Next, we detail the description of the radiation field. We identify the energy flux with a particle flux that is determined by a Boltzmann transport equation. The transport equation is split into a left-hand side and a right-hand side. The left-hand side is the directional derivative of the particle distribution function along trajectories of free particle propagation. This derivative is equated to the changes in the distribution function due to collisions, which are described by the right hand side of the equation. Once a 1+1 decomposition of space-time Arnowitt_Deser_Misner_62,Smarr_York_78 and a basis in the momentum phase space for the particle four-momentum have been chosen, the directional derivative along the phase flow can be expressed in terms of partial derivatives of the distribution function with respect to the space-time coordinates and momenta Lindquist_66,Mezzacappa_Matzner_89. We measure the particle four-momentum in a comoving orthonormal frame, with components

$$p^{a}=p\cos \theta ,\; p^{\vartheta }=p\sin \theta \cos \phi ,\; p^{\varphi }=p\sin \theta \sin \phi .$$ (13)

In spherical symmetry, the particle energy, , measured in a comoving frame, and the cosine of the angle between the particle momentum and the radial direction, , completely describe the particle phase space. The neutrinos are assumed to have no mass. In spherical symmetry, the distribution function does not depend on the three-momentum azimuth angle . Thus, the specific particle distribution function depends on four arguments and describes the number of particles at a given time, , in the phase space volume by

$$dN=F(t,a,\mu ,E)E^{2}dEd\mu da.$$ (14)

With the metric of Eq. (1), the Boltzmann equation reads Yamada_Janka_Suzuki_99,Liebendoerfer_Mezzacappa_Thielemann_01,

$$C_{t}+D_{a}+D_{\mu }+D_{E}+O_{\mu }+O_{E}=C_{c},$$ (15)

with

$$C_{t} = \frac{\partial F}{\alpha \partial t}$$ (16)

$$D_{a} = \frac{\mu }{\alpha }\frac{\partial }{\partial a}\left[ 4\pi r^{2}\alpha \rho F\right]$$ (17)

$$D_{\mu } = \Gamma \left( \frac{1}{r}-\frac{1}{\alpha }\frac{\partial \alpha ... ...\right) \frac{\partial }{\partial \mu }\left[ \left( 1-\mu ^{2}\right) F\right]$$ (18)

$$D_{E} = -\mu \Gamma \frac{1}{\alpha }\frac{\partial \alpha }{\partial r}\frac{1}{E^{2}}\frac{\partial }{\partial E}\left[ E^{3}F\right]$$ (19)

$$O_{E} = \left( \mu ^{2}\left( \frac{\partial \ln \rho }{\alpha \partial t... ...ac{u}{r}\right) \frac{1}{E^{2}}\frac{\partial }{\partial E}\left[ E^{3}F\right]$$ (20)

$$O_{\mu } = \left( \frac{\partial \ln \rho }{\alpha \partial t}+\frac{3u}{r}\right) \frac{\partial }{\partial \mu }\left[ \mu \left( 1-\mu ^{2}\right) F\right]$$ (21)

$$C_{c} = \frac{j}{\rho }-\chi F.$$ (22)

The source on the right-hand side, , is the collision term that describes changes in the particle distribution function due to local interactions with matter. It is represented here by an emissivity and an opacity . All other terms stem from the partial derivatives of the distribution function with respect to the phase-space coordinates in the directional derivative along the phase flow. They can all be physically interpreted. The first term on the left hand side of the equation, , is the temporal change of the particle distribution function. The second term, , counts the particles that are propagating into or out of an infinitesimal mass shell. The third term, , accounts for the change in the neutrino distribution function in an angle interval owing to the propagation of the neutrinos along geodesics with changing local angle cosine . The curved particle trajectories in general relativity are accounted for by the term proportional to the gradient of the gravitational potential, ,

$$\frac{1}{\alpha }\frac{\partial \alpha }{\partial r}=\frac{\partial \Phi }{\partial r}.$$

The fourth term, , expresses the redshift or blueshift of the particle energy that applies when the particles have a velocity component in the radial direction () and, therefore, change their position in the gravitational well. The fifth and sixth term, and , account for the Doppler shift and the angular aberration between adjacent comoving observers.

The integration of the Boltzmann equation over momentum space, spanned by the particle direction cosine and energy, gives the local conservation laws for particle number and energy. We define and to represent the zeroth and first moments of the distribution function:

$$J^{N} = \int ^{1}_{-1}\int ^{\infty }_{0}FE^{2}dEd\mu ,$$

$$H^{N} = \int ^{1}_{-1}\int ^{\infty }_{0}FE^{2}dE\mu d\mu .$$ (23)

Integration of Eq. (15) over and with as the measure of integration gives the following evolution equation for :

$$\frac{\partial J^{N}}{\partial t}+\frac{\partial }{\partial ... ... \frac{j}{\rho }E^{2}dEd\mu +\alpha \int \chi FE^{2}dEd\mu =0.$$ (24)

The derivatives with respect to the momentum phase space in Eq. (15) do not contribute because vanishes at and is zero for and . Eq. (24) is a continuity equation analogous to Eq. (6), extended by source and sink terms for the radiation particles. One more integration over the rest mass from the center of the star to its surface gives the evolution equation of the total particle number.

Slightly less straightforward is the derivation of total radiation energy conservation. We define the energy moments

$$J = \int FE^{3}dEd\mu$$

$$H = \int FE^{3}dE\mu d\mu$$

$$K = \int FE^{3}dE\mu ^{2}d\mu$$

$$Q = \int FE^{3}dE\mu ^{3}d\mu$$ (25)

and evaluate the evolution of the radiation energy as measured by an observer at infinity,

$$\frac{\partial }{\partial t}\int \left( \Gamma +u\mu \right) FE^{3}dEd\mu .$$ (26)

To this purpose, we integrate Eq. (15) again over phase space, but this time with measure of integration . After performing some integrations by parts to account for the time and space dependence of and , this leads to the concise result Liebendoerfer_Mezzacappa_Thielemann_01

$$0 = \frac{\partial }{\partial t}\left( \Gamma J+uH\right) +\frac{\par... ...+\Gamma H\right) \right] +4\pi r\alpha \rho \left( 1+e+\frac{p}{\rho }\right) H$$

$$- \alpha \Gamma \int \left( \frac{j}{\rho }-\chi \right) E^{3}dEd\mu +\alpha u\int \chi FE^{3}dE\mu d\mu .$$ (27)

Note that the conserved quantity is the radiation energy density in the frame of an observer at infinity. It is expressed in terms of the momentum moments and in the comoving frame. The second term describes the surface work by the radiation pressure, , and the energy loss or gain due to the luminosity at the boundary. The third term contains a gravitational term coupling the matter enthalpy with the luminosity that we neglected in previous work Liebendoerfer_Mezzacappa_Thielemann_01. The source terms in Eq. (27) describe the energy exchange with matter by particle emission, absorption, and radiation stress. The omitted terms from neutrino scattering enter the equation in a similar form. Beforehand, we found that Eq. (7) describes the evolution of the total energy. Then, from the Boltzmann equation, we derived Eq. (27) for the evolution of the radiation energy. Thus, we will find a consistent equation for the evolution of the hydrodynamics part by the subtraction of Eq. (27) from Eq. (7). The result is Eq. () in section where we discuss the implementation of the hydrodynamics part. The same procedure applied to other conserved quantities leads the full set of consistent hydrodynamics equations. By taking moments of the Boltzmann equation with the measures of integration , , and , respectively, we derive equations for the evolution of the radiation momentum, Eq. (28), a radiative contribution to the lapse function, Eq. (29), and the radiation energy in the comoving frame, Eq. (30):

$$\frac{\partial }{\partial t}\left[ uJ+\Gamma H\right] = -\frac{\partial }{\partial a}\left[ 4\pi r^{2}\alpha \rho \left( \Gamma K+uH\right) \right]$$

$$- \frac{\alpha }{r}\left[ \left( J+3K\right) \frac{m}{r}-\left( 1-\frac{2m}{r}\right) \left( J-3K\right) -2K\right.$$

$$+ \left. 4\pi r^{2}\left( J\left( p+\rho K\right) +\left( 1+e+J\right) \rho K-2\rho H^{2}\right) \right]$$

$$- \alpha \Gamma \int \chi FE^{3}dE\mu d\mu +\alpha u\int \left( \frac{j}{\rho }-\chi F\right) E^{3}dEd\mu$$ (28)

$$\frac{\partial }{\partial t}\left[ \frac{1}{4\pi r^{2}\rho }H\right] = -J\frac{\partial \alpha }{\partial a}-\frac{1}{\rho }\frac{\partial }{\partial a}\left[ \alpha \rho K\right] +\frac{\alpha }{3VD}\left( J-3K\right)$$

$$- \frac{1}{4\pi r^{2}\rho }\int \chi FE^{3}dE\mu d\mu$$ (29)

$$\frac{\partial J}{\partial t} = -\frac{1}{\alpha }\left[ 4\pi r^{2}\alpha ^{2}\rho H\right] -\rho... ...{\partial t}\left( \frac{1}{\rho }\right) -\frac{\alpha u}{r}\left( J-3K\right)$$

$$+ \alpha \int \left( \frac{j}{\rho }-\chi F\right) E^{3}dEd\mu .$$ (30)

The subtraction of Eq. (28) from Eq. (8) leads to a hydrodynamics equation for the evolution of the momentum, Eq. (). The subtraction of Eq. (29) from Eq. (11) leads to a hydrodynamics equation for the update of the lapse function, Eq. (). Finally, the substraction of Eq. (30) from Eq. (12) leads to a hydrodynamics equation for the evolution of the internal energy, Eq. (). We will pay attention to preserve this consistency also in our finite difference representation of the equations of radiation hydrodynamics. However, before we proceed with the technical details in section , we provide in the next subsection an overview of two exemplary simulation runs to complete the physical context and to illustrate the numerical challenges we face.