2 Radiation-hydrodynamics - basic equations

  
2.1 Hydrodynamics and equation of state

For an ideal fluid characterized by the mass density $\rho$, the Cartesian components of the velocity vector $(v_1, v_2, v_3)^{\rm T}$, the specific energy density $\varepsilon = e + \frac{1}{2} v^2$ and the gas pressure p, the Eulerian, nonrelativistic equations of hydrodynamics in Cartesian coordinates read (sum over i implied):

   

$$\partial_t\rho + \partial_i(\rho v_i) = 0 ~,$$ (1)

$$\partial_t(\rho v_k) + \partial_i(\rho v_i v_k+ \delta_{ik}p) = - \rho\partial_k\Phi^{\rm Newt} + {Q_{\rm M}}_k ~,$$ (2)

$$\partial_t(\rho \varepsilon) + \partial_i(\{\rho \varepsilon +p \}v_i) = - \rho v_i \partial_i\Phi^{\rm Newt} + {Q_{\rm E}} + v_i{Q_{\rm M}}_i ~,$$ (3)

where $\Phi^{\rm Newt}$ denotes the Newtonian gravitational potential of the fluid, which can be determined by the Poisson equation $\partial_i\partial^i \Phi^{\rm Newt}=4\pi G~\rho$ (G is Newton's constant). $\vec{Q}_{\rm M}$ and $Q_{\rm E}$ are the neutrino source terms for momentum transfer and energy exchange, respectively, $\delta_{ik}$ is the Kronecker symbol, and $\partial_i:=\partial/\partial x^i$ is an abbreviation for the partial derivative with respect to the coordinate xⁱ. In order to describe the evolution of the chemical composition of the (electrically neutral) fluid, the hydrodynamical equations are supplemented by a conservation equation for the electron fraction ${Y_{\rm e}}$,

$$\partial_t (\rho {Y_{\rm e}}) + \partial_i(\rho {Y_{\rm e}}~ v_i) = Q_{\rm N} ~,$$ (4)

where the source term $Q_{\rm N}$ describes the change of the net electron number density (i.e. the density of electrons minus that of positrons) due to emission and absorption of electron-flavour neutrinos. Unless nuclear statistical equilibrium (NSE) holds, an equation like (4) has to be solved for the mass fraction X_k of each of the $N_{\rm nuc}$ individual (nuclear) species. In NSE, $X_k=X_k(\rho,T,{Y_{\rm e}})$ is determined by the Saha equations.

An equation of state is invoked in order to express the pressure as a function of the independent thermodynamical variables, i.e., $p=p(\rho,T,{Y_{\rm e}})$, if NSE holds, or $p=p(\rho,T,{Y_{\rm e}},{\{X_k\}}_{k=1\dots N_{\rm nuc}})$ otherwise (see Appendix B for the numerical handling of the equation of state).

In the following we will employ spherical coordinates and, unless otherwise stated, assume spherical symmetry.

2.2 Equations for the neutrino transport

  
2.2.1 General relativistic transport equation

Lindquist (1966) derived a covariant transfer equation and specialized it for particles of zero rest mass interacting in a spherically symmetric medium supplemented with the comoving frame metric (a is a Lagrangian coordinate) ${\rm d}s^2= -{\rm e}^{2\Phi(t,a)}c^2{\rm d}t^2+ {\rm e}^{2\Lambda(t,a)}{\rm d}a^2+ R(t,a)^2{\rm d}\Omega^2$.

The "Lindquist-equation'', which describes the evolution of the specific intensity ${\cal I}$ as measured in the comoving frame of reference, reads:

 

$$\frac{1}{c} {\cal D}_t ~{\cal I} + \mu {\cal D}_a~{\cal I} + \Gamma~\frac{1-\mu^2}{R} \frac{\partial}{\partial \mu}~~{\cal I}$$

$$\frac{\partial}{\partial \mu} \left[ \left(1-\mu^2\right)\left\{ ... ...R} -\frac{1}{c}{\cal D}_t\Lambda\Big) -{\cal D}_a\Phi \right\}~{\cal I} \right]$$ +

$$\frac{\partial}{\partial \epsilon} \left[ \epsilon~\Big( \left(1-... ...} +\mu^2\frac{1}{c}{\cal D}_t\Lambda +\mu {\cal D}_a\Phi \Big)~{\cal I} \right]$$ -

$$\Big( \left(3-\mu^2\right)\frac{U}{R}+ \left(1+\mu^2\right)\frac{1}{c}{\cal D}_t\Lambda+ 2\mu {\cal D}_a\Phi \Big)~{\cal I}=C ~,$$ + (5)

where we use the classical abbreviations $\Gamma(t,a):={\cal D}_a R$ and $U(t,a):=c^{-1}{\cal D}_t R$ with ${\cal D}_a:={\rm e}^{-\Lambda} \partial_a$ and ${\cal D}_t:={\rm e}^{-\Phi}\partial_t$. The latter definition has been used in Eq. (5) also in order to emphasize that the time derivative has to be taken at fixed Lagrangian coordinate a.

The functional dependences of the metric functions $\Phi(t,a), \Lambda(t,a)$, R(t,a), the specific intensity ${\cal I}(t,a,\epsilon,\mu)$, and the collision integral $C(t,a,\epsilon,\mu)$ were suppressed for brevity. Momentum space is described by the coordinates $\epsilon$ and $\mu$, which are the energy and the cosine of the angle of propagation of the neutrino with respect to the radial direction, both measured in the locally comoving frame of reference. Note that the opacity $\chi$ and the emissivity $\eta$, and thus the collision integral $C=\eta - \chi{\cal I}$ in general depend also explicitly on momentum-space integrals of ${\cal I}$, which makes the transfer equation an integro-partial differential equation. Examples of the actual computation of the collision integral for a number of interaction processes of neutrinos with matter can be found in Appendix A.

  
2.2.2 $\mathsf{{\cal O}(v/c)}$ transport equations

In general, the metric functions $\Phi(t,a), \Lambda(t,a)$ and R(t,a) have to be computed numerically from the Einstein field equations. When working to order ${\cal O}(\beta:=v/c)$ and in a flat spacetime (usually called the "Newtonian approximation''), it is however possible to express these functions analytically in terms of only the velocity field and its first time derivative (the fluid acceleration). Details of the derivation can be found in Castor (1972). Alternatively one can simply reduce the special relativistic transfer equation (Mihalas 1980) to order ${\cal O}(v/c)$.

This transfer equation, together with its angular moment equations of zeroth and first order reads (e.g., Mihalas & Mihalas 1984, see also Lowrie et al. 2001):

 

$$\Big(\frac{1}{c}\frac{\partial}{\partial t} \beta \frac{\partial}{\partial r}\Big)~{\cal I} + \mu\frac{\partial}{\partial r}~{\cal I} + \frac{1-\mu^2}{r}\frac{\partial}{\partial \mu}~{\cal I}$$ +

$$\frac{\partial}{\partial \mu}\left[(1-\mu^2) \left\{\mu\Big(\frac... ...r}\Big)-\frac{1}{c}\frac{\partial \beta}{\partial t} \right\}~ {\cal I} \right]$$ +

$$\frac{\partial}{\partial \epsilon}\left[\epsilon~\Big((1-\mu^2) \... ...ial r}+\mu \frac{1}{c}\frac{\partial \beta}{\partial t} \Big) ~{\cal I} \right]$$ -

$$\Big((3-\mu^2)\frac{\beta}{r}+(1+\mu^2) \frac{\partial \beta}{\partial r} + \mu \frac{2}{c}\frac{\partial \beta}{\partial t} \Big)~{\cal I} = C ~,$$ + (6)

 

$$\Big(\frac{1}{c}\frac{\partial}{\partial t} \beta\frac{\partial}{{\partial r}}\Big) J + \frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2H\right)$$ +

$$\frac{\partial}{\partial \epsilon}\left[\epsilon\left(\frac{\beta... ...eta}{\partial r}K+\frac{1}{c}\frac{\partial \beta}{\partial t}H \right) \right]$$ -

$$\frac{\beta}{r}(3J-K)+\frac{\partial \beta}{\partial r}(J+K)+\frac{2}{c}\frac{\partial \beta}{\partial t}H=C^{(0)} ~,$$ + (7)

 

$$\Big(\frac{1}{c}\frac{\partial}{\partial t} \beta\frac{\partial}{{\partial r}}\Big) H + \frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2K\right) + \frac{K-J}{r}$$ +

$$\frac{\partial}{\partial \epsilon}\left[\epsilon~\left(\frac{\bet... ...eta}{\partial r}L+\frac{1}{c}\frac{\partial \beta}{\partial t}K \right) \right]$$ -

$$2\left(\frac{\partial \beta}{\partial r}+\frac{\beta}{r}\right) H +\frac{1}{c}\frac{\partial \beta}{\partial t}(J+K)=C^{(1)} ~,$$ + (8)

where, in spherical symmetry, the angular moments of the specific intensity are given by

$$\{J,H,K,L,\dots\}(t,r,\epsilon):= \frac{1}{2}\int\limits_{-1}... ...rm d}\mu~\mu^{\{0,1,2,3,\dots\}} {\cal I}(t,r,\epsilon,\mu) ~.$$ (9)

Note that in Eqs. (6)-(8) all physical quantities, in particular also the collision integral and its angular moments $C^{(k)}(t,r,\epsilon):= \frac{1}{2}\int_{-1}^{+1}{\rm d}\mu~\mu^{k}C(t,r,\epsilon,\mu)$ are measured in the comoving frame, but the choice of coordinates (r,t) is Eulerian. The simple replacement $\partial/\partial t+v\partial/\partial r\to {\rm D}/{\rm D}t$ yields the conversion to Lagrangian coordinates.

For reference we also write down the transformations (correct to ${\cal O}(\beta)$) which allow one to relate the frequency-integrated moments in the comoving ("Lagrangian'') and in the inertial ("Eulerian'') frame of reference (indicated by the superscript "Eul'').

 

$$J^{\rm Eul} J + 2\beta~ H ~,$$ =

$$H^{\rm Eul} H + \beta~ (J + K ) ~,$$ = (10)

$$K^{\rm Eul} K + 2\beta~ H ~.$$ =

Equation (10) can easily be deduced from a Lorentz-transformation of the radiation stress-energy tensor (e.g., Mihalas & Mihalas 1984). In principle also relations for monochromatic moments can be derived by transforming the specific intensity ${\cal I}$, the angle cosine $\mu$ and the energy $\epsilon$, which, however, leads to more complicated expressions.

The system of Eqs. (6)-(8) is coupled to the evolution equations of the fluid (Eqs. (1)-(4)) in spherical coordinates and symmetry) by virtue of the definitions of the source terms

   

$$Q_{\rm E} -4\pi\int_0^\infty{\rm d}\epsilon~C^{(0)}(\epsilon)~,$$ = (11)

$$Q_{\rm M} -\frac{4\pi}{c}\int_0^\infty{\rm d}\epsilon~C^{(1)}(\epsilon)~,$$ = (12)

$$Q_{\rm N} -4\pi~ m_{\rm B}\int_0^\infty{\rm d}\epsilon~\frak{C}^{(0)}(\epsilon)~,$$ = (13)

where $m_{\rm B}$ denotes the baryonic mass, and $\frak{C}^{(0)}(\epsilon):=\epsilon^{-1}C^{(0)}(\epsilon)$.

Recently, Lowrie et al. (2001) emphasized the fundamental significance of a term

$$\beta\mu\cdot\frac{1}{c}\frac{\partial}{{\partial t}}~{\cal I}$$ (14)

on the left hand side of the special relativistic transport equation. This term has traditionally been dropped in deriving the ${\cal O}(v/c)$-approximation (Eq. (6)) by assuming the time variation of all dependent variables (like, e.g., ${\cal I}$) to be on a fluid time scale (given by l/v, where l is the characteristic length scale of the system and v a typical fluid velocity). In this case the term given by Eq. (14) is of order v²/c² and can be omitted (Mihalas & Mihalas 1984). If, on the other hand, appreciable temporal variations of, e.g., the specific intensity occur on radiation time scales (given by l/c), Eq. (6) is no longer correct to ${\cal O}(v/c)$ (Lowrie et al. 2001).

In core-collapse supernova simulations carried out so far, the dynamics of the stellar fluid presumably was not affected by neglecting the term in Eq. (6). However, our tests with Eq. (6), including the additional time derivative of Eq. (14) and the corresponding changes in the moment equations (Eqs. (7), (8)), have shown that the neutrino signal computed in a supernova simulation is indeed altered compared to the traditional treatment. We will therefore take the term of Eq. (14) into account in future simulations.

For calculating nonrelativistic problems Mihalas & Mihalas (1984) suggested a form of the radiation momentum equation (Eq. (8)), in which all velocity-dependent terms except for the $\beta\partial/\partial r$-term in the first line of Eq. (8) are dropped. When the velocities become sizeable, however, it may be advisable to solve the momentum equation in its general form (Eq. (8)). Doing so, we indeed found that the terms omitted by Mihalas & Mihalas (1984) can have an effect on the solution of the neutrino transport in supernovae, in particular on the neutrino energy spectrum.