A&A 396, 361-392 (2002)
DOI: 10.1051/0004-6361:20021398

Radiation hydrodynamics with neutrinos

Variable Eddington factor method for core-collapse supernova simulations

M. Rampp - H.-T. Janka

Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85741 Garching, Germany

Received 8 March 2002 / Accepted 19 September 2002

Abstract
Neutrino transport and neutrino interactions in dense matter play a crucial role in stellar core collapse, supernova explosions and neutron star formation. Here we present a detailed description of a new numerical code for treating the time and energy dependent neutrino transport in hydrodynamical simulations of such events. The code is based on a variable Eddington factor method to deal with the integro-differential character of the Boltzmann equation. The moments of the neutrino distribution function and the energy and lepton number exchange with the stellar medium are determined by iteratively solving the zeroth and first order moment equations in combination with a model Boltzmann equation. The latter is discretized on a grid of tangent rays. The integration of the transport equations and the neutrino source terms is performed in a time-implicit way. In the present version of the program, the transport part is coupled to an explicit hydrodynamics code which follows the evolution of the stellar plasma by a finite-volume method with piecewise parabolic interpolation, using a Riemann solver to calculate the hydrodynamic states. The neutrino source terms are implemented in an operator-split step. Neutrino transport and hydrodynamics can be calculated with different spatial grids and different time steps. The structure of the described code is modular and offers a high degree of flexibility for an application to relativistic and multi-dimensional problems at different levels of refinement and accuracy. We critically evaluate results for a number of test cases, including neutrino transport in rapidly moving stellar media and approximate relativistic core collapse, and suggest a path for generalizing the code to be used in multi-dimensional simulations of convection in neutron stars and supernovae.

Key words: stars: supernovae: general - elementary particles - hydrodynamics - neutrinos

1 Introduction

When the core of a massive star collapses to a neutron star, a huge amount of gravitational binding energy is released mainly in the form of neutrinos which are abundantly produced by particle reactions in the hot and dense plasma. In fact, the emission of neutrinos determines the sequence of dynamical events which precede the death of the star in a supernova explosion. Electron neutrinos from electron captures on protons and nuclei reduce the electron fraction and the pressure and thus accelerate the contraction of the stellar iron core to a catastrophic implosion. The position of the formation of the supernova shock at the moment of core bounce depends on the degree of deleptonization during the collapse phase. The outward propagation of the prompt shock is severely damped by energy losses due to the photodisintegration of iron nuclei, and finally stopped only a few milliseconds later when additional energy is lost in a luminous outburst of electron neutrinos. These neutrinos are created in the shock-heated matter and leave the star as soon as their diffusion is faster than the expansion of the shock. "Prompt'' supernova explosions generated by a direct propagation of the hydrodynamical shock have been obtained in simulations only when the stellar iron core is very small (Baron & Cooperstein 1990) and/or the equation of state of neutron-rich nuclear matter is extraordinarily soft (Baron et al. 1985; Baron et al. 1987).

At a later stage (roughly a hundred milliseconds after the shock formation) the situation has changed. The temperature behind the stalled shock has dropped such that increasingly energetic neutrinos diffusing out from deeper layers start to transfer energy to the stellar gas around the nascent neutron star. If this energy deposition is sufficiently strong, the stalled shock can be revived and a "delayed'' explosion can be triggered (Wilson 1985; Bethe & Wilson 1985; the idea that neutrinos provide the energy of the supernova explosion was originally brought up by Colgate & White 1966). A small fraction of less than one per cent of the energy released in neutrinos can account for the kinetic energy of a typical type II supernova. This explanation is currently favored for the explosion of massive stars in the mass range between about 10 $~M_{\odot}$ and roughly 25 $~M_{\odot}$. It is supported by numerical simulations and analytic considerations. A finally convincing hydrodynamical simulation, however, is still missing (a summary of our current understanding of the explosion mechanism can be found in Janka 2001). The measurement of a neutrino signal from a Galactic supernova would offer the most direct observational test for our theoretical perception of the onset of the explosion. Due to the central role of neutrinos in supernovae, the neutrino transport deserves particular attention in numerical models.

1.1 Brief review of the status of supernova models

Current hydrodynamic models of neutrino-driven supernova explosions leave an ambiguous impression and have caused confusion about the status of the field outside the small community of supernova modelers. Simulations by different groups seem to be contradictory because some models show successful explosions by the neutrino-driven mechanism whereas others have found the mechanism to fail.

Wilson and collaborators have performed successful simulations for more than 15 years now (e.g., Wilson 1985; Wilson et al. 1986; Wilson & Mayle 1988; Wilson & Mayle 1993; Mayle et al. 1993; Totani et al. 1998). Their models were calculated in spherical symmetry, but shock revival was obtained by making the assumption that the neutrino luminosity is boosted by neutron-finger instabilities in the nascent neutron star (Wilson & Mayle 1988, 1993). In fact, in their calculations explosion energies comparable to typically observed values (of the order of 10⁵¹ erg) require pions to be abundant in the nuclear matter already at moderately high densities (Mayle et al. 1993). The corresponding equation of state (EoS) with pions leads to higher core temperatures and the emission of more energetic neutrinos. Both the convectively boosted luminosities and the harder spectra enhance the neutrino heating behind the stalled shock. The relevance of neutron-finger instabilities for efficient energy transport on large scales, however, has not been demonstrated by direct simulations. The existence of neutron-finger instabilities (in the sense of the definition introduced by Wilson & Mayle 1993) is indeed questioned by other investigations (Bruenn & Dineva 1996) and might be a consequence of specific properties of the high-density EoS or the treatment of the neutrino transport by Wilson and collaborators. Also the implementation of a pionic component in their EoS is at least controversial and not in agreement with other, more conventional descriptions of nuclear matter (Pethick & Pandharipande, personal communication).

Rather than finding neutron-finger instabilities, two-dimensional hydrodynamic simulations have shown that regions inside the nascent neutron star can exist where Ledoux or quasi-Ledoux convection (as defined by Wilson & Mayle 1993) develops on a time scale of tens of milliseconds after core bounce (Keil et al. 1996; Keil 1997; also Janka et al. 2001). The models were constrained to a simulation of the neutrino cooling of the nascent neutron star, but it must be expected that the convective enhancement of the neutrino luminosity, which becomes appreciable after about 200 milliseconds (see also Janka et al. 2001), can have important consequences for the revival of the stalled supernova shock (Burrows 1987). Ledoux unstable regions were also detected in spherical cooling models of newly formed neutron stars (Burrows 1987; Burrows & Lattimer 1988; Pons et al. 1999; Miralles et al. 2000). Nevertheless, their significance is controversial, because Bruenn et al. (1995) in spherically symmetric models and Mezzacappa et al. (1998a) in two-dimensional simulations observed convection inside the neutrinosphere only as a transient phenomenon during a relatively short period after core bounce.

The question is therefore undecided whether convection plays an important role in newly formed neutron stars and if so, what its implications for the supernova explosion are. Unfortunately, the various calculations were performed with different treatments of the neutrino transport, one-dimensional vs. two-dimensional hydrodynamics, general relativistic gravitational potential vs. Newtonian gravity, or even an inner boundary condition to replace the central, dense part of the neutron star core (Mezzacappa et al. 1998a). It is this inner region, however, where Keil et al. (1996) found convection to develop roughly 100 milliseconds after bounce. Convection at larger radii and thus closer below the neutrinosphere had indeed died out within only 20-30 milliseconds after shock formation in agreement with the results of Bruenn et al. (1995). Future studies with a better and more consistent handling of the different aspects of the physics are definitely needed to clarify the situation.

A second hydrodynamically unstable region develops exterior to the nascent neutron star in the neutrino-heated layer behind the stalled supernova shock. Convective overturn in this region is helpful and can lead to explosions in cases which otherwise fail (Herant et al. 1992; Herant et al. 1994; Janka & Müller 1995; Janka & Müller 1996; Burrows et al. 1995). In two-dimensional supernova models computed recently by Fryer (1999); Fryer & Heger (2000), and Fryer et al. (2002) these hydrodynamical instabilities in the postshock region are crucial for the success of the neutrino-driven mechanism, because they help transporting hot gas from the neutrino-heating region directly to the shock, while downflows simultaneously carry cold, accreted matter to the layer of strongest neutrino heating where a part of this gas readily absorbs more energy from the neutrinos. The existence of this multi-dimensional phenomenon seems to be generic for the situation which builds up in the stellar core some time after shock formation. It is therefore no matter of dispute.

All simulations showing explosions as a consequence of post-shock convection, however, have so far been performed with a strongly simplified treatment of the crucial neutrino physics, e.g., with grey, flux-limited diffusion which is matched to a "light bulb'' description at some "low'' value of the optical depth (Herant et al. 1994; Burrows et al. 1995). Alternatively, an inner boundary near the neutrinosphere had been used where the spectra and luminosities of neutrinos were prescribed to parameterize our ignorance and the potential uncertainties of the exact properties of the neutrino emission from the newly formed neutron star (Janka & Müller 1995, 1996). It is therefore not clear whether the instabilities and their associated effects in fully self-consistent and more accurate simulations will be strong enough to cause successful explosions. Doubts in that respect were raised by Mezzacappa et al. (1998b), whose two-dimensional models showed convective overturn in the neutrino-heating region but still no explosion. Mezzacappa et al. (1998b) combined two-dimensional hydrodynamics with neutrino transport results obtained by multi-group flux-limited diffusion in spherical symmetry. Although not self-consistent, their approach is nevertheless an improvement compared to previous treatments of the neutrino physics by other groups.

All computed models bear some pieces of truth. Essentially they demonstrate the remarkable sensitivity of the supernova dynamics to the different physical aspects of the problem, in particular the treatment of neutrino transport and neutrino-matter interactions, the properties of the nuclear EoS, multi-dimensional hydrodynamical processes, and general relativity. Considering the huge energy reservoir carried away by neutrinos, the neutrino-driven mechanism appears rather inefficient (the often quoted value of 1% efficiency, however, misjudges the true situation, because neutrinos can transfer between 5% and 10% of their energy to the stellar gas during the critical period of shock revival; see Janka 2001). Nevertheless, neutrino-driven explosions are "marginal'' in the sense that the energy of a standard supernova explosion is of the same order as the gravitational binding energy of the ejected progenitor mass. The final success of the supernova shock is the result of different physical processes which compete against each other. On the one hand, neutrino heating tries to drive the shock expansion, on the other hand energy losses, e.g., by neutrinos that are reemitted from the inward flow of neutrino-heated matter which enters the neutrino-cooling zone below the heating layer, extract energy and thus damp the shock revival. It is therefore not astonishing that different approximations or descriptions for one or more physical components of the problem can decide between an explosion or failure of a simulation.

None of the current supernova models includes all relevant aspects to a satisfactory level of accuracy, but all of these models are deficient in one or more respects. Wilson and collaborators obtain explosions, but their input physics is unique and cannot be considered as generally accepted. Two-dimensional (Herant et al. 1994; Burrows et al. 1995; 1999; Fryer & Heger 2000; Fryer et al. 2002) and three-dimensional (Fryer, personal communication) models show explosions due to strong postshock convection, but the neutrino transport and neutrino-matter interactions are handled at a level of accuracy which falls back behind the most elaborate treatments that have been applied in spherical symmetry. The sensitivity of the outcome of numerical calculations to details of the neutrino transport demands improvements. Self-consistent multi-dimensional simulations with a sophisticated and quantitatively reliable treatment of the neutrino physics have yet to be performed.

1.2 Boltzmann solvers for neutrino transport

Most recently, spherically symmetric Newtonian (Rampp & Janka 2000; Mezzacappa et al. 2001), and general relativistic (Liebendörfer et al. 2001b) hydrodynamical simulations of stellar core-collapse and post-bounce evolution including a Boltzmann solver for the neutrino transport have become possible. Although the models do not yield explosions, they must be considered as a major achievement for the modeling of supernovae. Before these calculations only the collapse phase of the stellar core had been investigated with solving the Boltzmann equation (Mezzacappa & Bruenn 1993c). Boltzmann transport has now superseded multi-group flux-limited diffusion (Arnett 1977; Bowers & Wilson 1982; Bruenn 1985; Myra et al. 1987; Baron et al. 1989; Bruenn 1989a,b, 1993; Cooperstein & Baron 1992) as the most elaborate treatment of neutrinos in supernova models.

The differences between both methods in dynamical calculations have still to be figured out in detail, but the possibility of solving the Boltzmann equation removes imponderabilities and inaccuracies associated with the use of flux-limiters in particular in the region of semi-transparency, where neutrinos decouple from the stellar background and also deposit the energy for an explosion (Janka 1991; Janka 1992; Messer et al. 1998; Yamada et al. 1999). For the first time the transport can now be handled at a level of sophistication where the technical treatment is more accurate than our standard description of neutrino-matter interactions, which includes various approximations and simplifications (for an overview of the status of the handling of neutrino-nucleon interactions, see Horowitz 2002 and the references therein).

In this paper we describe our new numerical code for solving the energy and time dependent Boltzmann transport equation for neutrinos coupled to the hydrodynamics of the stellar medium. First results from supernova calculations with this code were published before (Rampp & Janka 2000), but a detailed documentation of the method will be given here. It is based on a variable Eddington factor technique where the coupled set of Boltzmann equation and neutrino energy and momentum equations is iterated to convergence. Variable Eddington moments of the neutrino phase space distribution are used for closing the moment equations, and the integro-differential character of the Boltzmann equation is tamed by using the zeroth and first order angular moments (neutrino density and flux) in the source terms on the right hand side of the Boltzmann equation. This numerical approach is fundamentally different from the so-called $S_{\rm N}$ methods (Carlson 1967; Yueh & Buchler 1977; Mezzacappa & Bruenn 1993a; Yamada et al. 1999) which employ a direct discretization of the Boltzmann equation in all variables including the dependence on the angular direction of the radiation propagation.

Some basic characteristics of our code are similar to elements described by Burrows et al. (2000). Different from the latter paper we shall discuss the details of the numerical implementation of the transport scheme and its coupling to a hydrodynamics code (Rampp 2000), which in our case is the PROMETHEUS code with the potential of performing multi-dimensional simulations. The variable Eddington factor technique was our method of choice for the neutrino transport because of its modularity and flexibility, which offer significant advantages for a generalization to multi-dimensional problems. We shall suggest and motivate corresponding approximations which we consider as reasonable in the supernova case, at least as a first step towards multi-dimensional hydrodynamics with a Boltzmann treatment of the neutrino transport.

This paper is arranged as follows: in Sect. 2 the equations of radiation hydrodynamics are introduced. Section 3 provides a general overview of the numerical methods used to solve these equations and contains details of their practical implementation. Results for a number of idealized test problems as well as applications of the new method to the supernova problem are presented in Sect. 4. A summary will be given in Sect. 5. In the Appendix the numerical implementation of neutrino opacities and the equation of state used for core-collapse and supernova simulations is described.

  
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.

  
3 Numerical implementation

3.1 Overview of the basic numerical strategy

For solving the neutrino radiation hydrodynamics problem we employ the well-known "operator splitting''-approach, a popular method to make large problems computationally tractable by considering them as a combination of smaller subproblems. The coupled set of equations of radiation hydrodynamics (Eqs. (1)-(4), (6)-(8)) is split into a "hydrodynamics part'' and a "neutrino part'', which are solved independently in subsequent ("fractional'') steps. In the "hydrodynamics part'' (Sect. 3.2) a solution of the equations of hydrodynamics including the effects of gravity is computed without taking into account neutrino-matter interactions.

In what we call the "neutrino part'' (Sects. 3.3-3.4) the coupled system of equations describing the neutrino transport and the changes of the specific internal energy e and the electron fraction ${Y_{\rm e}}$ of the stellar fluid due to the emission and absorption of neutrinos are solved. The latter changes are given by the equations

  

$$\frac{\delta}{\delta t} e= -\frac{4\pi}{\rho}\int_{0}^{\infty} {\... ...\sum_{\nu\in\{{\nu_{\rm e}},{\bar\nu_{\rm e}},\dots\}} C^{(0)}_\nu(\epsilon) ~,$$ (15)

$$\frac{\delta}{\delta t} {Y_{\rm e}}= -\frac{4\pi m_{\rm B}}{\rho}... ...}_{\nu_{\rm e}}(\epsilon)- \frak{C}^{(0)}_{\bar\nu_{\rm e}}(\epsilon)\right) ~.$$ (16)

Here, the notation $\delta/\delta t$ indicates that the full problem as given, e.g., by Eq. (4), is split into the equations $\partial_t(\rho {Y_{\rm e}}) + r^{-2} \partial_r~(r^2\rho {Y_{\rm e}}~v) = 0$ and $\delta_t~{Y_{\rm e}}=Q_{\rm N}/\rho$.

Within the neutrino part, we calculate in a first step the transport of energy and electron-type lepton number by neutrinos and the corresponding exchange with the stellar fluid for electron neutrinos ( ${\nu_{\rm e}}$) and antineutrinos ( ${\bar\nu_{\rm e}}$) only, which also means that the sum in Eq. (15) is restricted to $\nu \in \{{\nu_{\rm e}}, {\bar\nu_{\rm e}}\}$. The $\mu$ and $\tau$ neutrinos and their antiparticles, all represented by a single species ("${\nu_x}$'') are handled together with the remaining terms of the sum in Eq. (15) in a second step.

For the neutrino transport we use the so-called "variable Eddington factor'' method, which determines the neutrino phase-space distribution by iteratively solving the radiation moment equations for neutrino energy and momentum coupled to the Boltzmann transport equation. Closure of the set of moment equations is achieved by a variable Eddington factor calculated from a formal solution of the Boltzmann equation, and the integro-differential character of the latter is tamed by making use of the integral moments of the neutrino distribution as obtained from the moment equations. The method is similar to the one described by Burrows et al. (2000).

  
3.2 The hydrodynamics code

For the integration of the equations of hydrodynamics we employ the Newtonian finite-volume hydrodynamics code PROMETHEUS developed by Fryxell et al. (1989). PROMETHEUS is a direct Eulerian, time-explicit implementation of the Piecewise Parabolic Method (PPM) of Colella & Woodward (1984), which is a second-order Godunov scheme employing a Riemann solver. PROMETHEUS is particularly well suited for following discontinuities in the fluid flow like shocks or boundaries between layers of different chemical composition. It is capable of tackling multi-dimensional problems with high computational efficiency and numerical accuracy.

The version of PROMETHEUS used in our radiation hydrodynamics code was provided to us by Keil (1997). It offers an optional generalization of the Newtonian gravitational potential to include general relativistic corrections. The hydrodynamics code was considerably updated and augmented by K. Kifonidis who added a simplified version of the "Consistent Multifluid Advection (CMA)'' method (Plewa & Müller 1999) which ensures an accurate advection of the individual chemical components of the fluid. In order to avoid spurious oscillations arising in multidimensional simulations, when strong shocks are aligned with one of the coordinate lines (the so-called "odd-even decoupling'' phenomenon), a hybrid scheme was implemented (Kifonidis 2000; Plewa & Müller 2001) which, in the vicinity of such shocks, selectively switches from the original PPM method to the more diffusive HLLE solver of Einfeldt (1988).

  
3.3 Multi-group formulation of the $\mathsf{{\cal O}(v/c)}$ moment equations

In what follows we describe the numerical implementation of the "Newtonian'' limit of the transport equations including ${\cal O}(v/c)$-terms (cf. Sect. 2.2.2). The general relativistic version of the code is based on almost exactly the same considerations.

  
3.3.1 General considerations

In order to construct a conservative numerical scheme for the system of moment equations, we integrate the moment equations (Eqs. (7), (8)) over a zone of the numerical $(r,\epsilon)$-mesh.

Spatial discretization:

After having performed the integral over the volumes $\Delta V$ of the radial mesh zones, the moment equations can be rewritten as evolution equations for the volume-integrated moments (e.g.,  $\int_{\Delta V}J~{\rm d}V$). In case of a moving radial grid, where the volume $\Delta V$ of a grid cell is time dependent, one has to apply the "moving grid transport theorem'' (see e.g. Müller 1991) in order to interchange the operators D_t and $\int_{\Delta V}{\rm d}V$. For the special case of a Lagrangian grid, where the grid moves with the velocity of the stellar fluid, the latter reads:

$$\int_{\Delta V}{\rm d}V~(D_t~\Xi)= D_t\int_{\Delta V} {\rm d}V~ \Xi - \int_{\Delta V}{\rm d}V~(\Xi~ {\rm div}~ \vec{v}) ~.$$ (17)

Here and in the following the symbol $\Xi$ represents one of the moments J or H (and also ${\cal J}:=\epsilon^{-1}J$ or ${\cal H}:=\epsilon^{-1}H$).

The computational domain $r_{\rm min} \le r \le r_{\rm max}$ is covered by N_r radial zones. As the zone center $r_{i+{\frac{1}{2}}}$ we define the volume-weighted mean of the interface coordinates r_i and r_i+1:

$$r_{i+{\frac{1}{2}}}:=\left(\frac{1}{2}~\left(r_i^3+r_{i+1}^3\right)\right)^{1/3},\quad i=0,\dots,N_r-1 ~.$$ (18)

Scalar quantities like $J_{i+{\frac{1}{2}}}$, $K_{i+{\frac{1}{2}}}$ are defined on zone centers, whereas vectors like H_i and L_ilive on the interfaces of radial zones. "Miscentered'' values like, e.g,. J_i or $H_{i+{\frac{1}{2}}}$, are in general computed by linear interpolation, using the values of the nearest neighbouring coordinates ( $J_{i-{\frac{1}{2}}}$ and $J_{i+{\frac{1}{2}}}$, or H_i and H_i+1, respectively).

Spectral discretization:

The spectral range $0\le \epsilon \le \epsilon_{\rm max}$ is covered by a discrete energy grid consisting of $N_\epsilon$ energy "bins'', where the centers of these zones are given in terms of the interface values as

$$\epsilon_{j+{\frac{1}{2}}}:=\frac{1}{2}~(\epsilon_{j}+\epsilon_{j+1}),\quad j=0,\dots,N_\epsilon-1 ~,$$ (19)

and the width of an energy bin is denoted by ${\Delta\epsilon}_{j+{\frac{1}{2}}} :=\epsilon_{j+1}-\epsilon_{j}$. The radiation moments are interpreted as average values within energy bins:

$$\Xi_{i,j+{\frac{1}{2}}}:= \frac{1}{{\Delta \epsilon}_{j+{\fr... ...\epsilon_j}^{\epsilon_{j+1}}{\rm d}\epsilon~\Xi_i(\epsilon) ~.$$ (20)

Time differencing:

We employ a time-implicit finite differencing scheme for the moment equations (Eqs. (7), (8)) in combination with the equations which describe the exchange of internal energy and electron fraction with the stellar medium (Eqs. (15), (16)). Doing so we avoid the restrictive Courant Friedrichs Lewy (CFL) condition which explicit numerical schemes have to obey for stability reasons: for neutrino transport in the optically thin regime the CFL condition would require $c\Delta t < \Delta r$, where $\Delta r$ is the size of the smallest zone and $\Delta t$ the numerical time step. Even more importantly, the time-implicit discretization allows one to efficiently cope with the different time scales of the problem: the "stiff'' source terms introduce a characteristic time scale proportional to the mean free path $\lambda$ of the neutrinos, $c\Delta t = \lambda=\Delta r/\Delta\tau$, which can be orders of magnitude smaller than the CFL time step in the opaque interior of a protoneutron star, where the optical depth $\Delta\tau$ of individual grid cells becomes very large.

  
3.3.2 Finite differencing of the moment equations

Applying the procedures described in Sect. 3.3.1 to the moment equations (Eqs. (7), (8)) we obtain for each neutrino species $\nu\in\{{\nu_{\rm e}},{\bar\nu_{\rm e}},\nu_{\mu},\dots\}$ (the corresponding index is suppressed in the notation below) the finite differenced moment equations. On an Eulerian grid they read:

 

$$\frac{J_{i+{\frac{1}{2}},j+{\frac{1}{2}}}^{n+1}-J_{i+{\frac{1}{2}},j+{\frac{1}{2}}}^{n}} {ct^{n+1}-ct^{n}} \frac{4\pi}{\Delta V_{i+{\frac{1}{2}}}} \left(r_{i+1}^2H_{i+1,j+{\frac{1}{2}}}^{n+1}-r_{i}^2H_{i,j+{\frac{1}{2}}}^{n+1}\right)$$ +

$$\frac{4\pi}{\Delta V_{i+{\frac{1}{2}}}} \left(r_{i+1}^2\beta_{i+1... ...1}{2}}}- r_{i}^2\beta_{i}J_{\iota(i),j+{\frac{1}{2}}}^{ n+{\frac{1}{2}}}\right)$$ +

$$\frac{1}{{\Delta \epsilon}_{j+{\frac{1}{2}}}} \left\{ \epsilon_{j... ..._{i+{\frac{1}{2}},j+1} -\epsilon_j {[u\;J^{n+1}]}_{ i+{\frac{1}{2}},j} \right\}$$ -

$$\left[\frac{\beta}{r}\right]_{i+{\frac{1}{2}}}{[(1-f_K)\; J^{n+1}]}_{i+{\frac{1}{2}},j+{\frac{1}{2}}}$$ +

$$\left[\frac{\partial \beta}{\partial r}\right]_{ i+{\frac{1}{2}}}... ...{ i+{\frac{1}{2}},j+{\frac{1}{2}}}\; J_{ i+{\frac{1}{2}},j+{\frac{1}{2}}}^{n+1}$$ +

$$\frac{2}{c}\left[\frac{\partial \beta}{\partial t}\right]_{i+{\fr... ...{1}{2}},j+{\frac{1}{2}}} = {C^{(0)}}^{n+1}_{i+{\frac{1}{2}},j+{\frac{1}{2}}} ~,$$ + (21)

 

$$\frac{H_{i,j+{\frac{1}{2}}}^{n+1}-H_{i,j+{\frac{1}{2}}}^{n}} {ct^{n+1}-ct^{n}} \frac{4\pi}{\Delta V_{i}} (r_{i+{\frac{1}{2}}}^2{[f_K J^{n+1}]}_{... ...}{2}}}- r_{i-{\frac{1}{2}}}^2{[f_K J^{n+1}]}_{i-{\frac{1}{2}},j+{\frac{1}{2}}})$$ +

$$\frac{2\pi}{\Delta V_{i}}(r_{i+{\frac{1}{2}}}^2-r_{i-{\frac{1}{2}}}^2 )~ {[(f_K-1)~J^{n+1}]}_{i,j+{\frac{1}{2}}}$$ +

$$\frac{4\pi}{\Delta V_{i}}\left( r_{i+{\frac{1}{2}}}^2\beta_{i+{\f... ...rac{1}{2}}}H_{\iota(i-{\frac{1}{2}}),j+{\frac{1}{2}}}^{ n+{\frac{1}{2}}}\right)$$ +

$$\frac{1}{{\Delta \epsilon}_{j+{\frac{1}{2}}}} \left\{ \epsilon_{j+1} {[w\;J^{n+1}]}_{i,j+1} -\epsilon_j {[w\;J^{n+1}]}_{ i,j} \right\}$$ -

$$\left[\frac{\partial \beta}{\partial r}\right]_{i}\;H_{ i,j+{\fra... ...{f_K})\; J^{n+1}]}_{ i,j+{\frac{1}{2}}} ={C^{(1)}}^{n+1}_{i,j+{\frac{1}{2}}} ~.$$ + (22)

Note that by virtue of introducing the "Eddington factors'' (see the following section for their actual computation) f_H:=H/J, f_K:=K/J and f_L:=L/J, the system of moment equations has been closed. The "flux'' in energy space across the boundaries between energy zones appearing in the third and fifth line of Eq. (21) and Eq. (22), respectively, is given in terms of the interface values $J_{i+{\frac{1}{2}},j}$, which are defined as the geometrical mean of $J_{i+{\frac{1}{2}},j-{\frac{1}{2}}}$ and $J_{i+{\frac{1}{2}},j+{\frac{1}{2}}}$, and the advection velocities

 

$$u_{i+{\frac{1}{2}},j}:= \left[\frac{\beta}{r}\right]_{i+{\frac{1}... ...rtial \beta}{\partial t}\right]_{i+{\frac{1}{2}}} {[f_H]}_{i+{\frac{1}{2}},j}~,$$ (23)

 

$$w_{i,j}:= \left[\frac{\beta}{r}\right]_{i}{[f_H-f_L]}_{i,j} + \le... ... \frac{1}{c}\left[\frac{\partial \beta}{\partial t}\right]_{i} {[f_K]}_{i,j} ~.$$ (24)

The nonlinear interpolation of J in energy space turned out to be crucial for obtaining a numerically stable algorithm. The Eddington factors f_H, f_K and f_L are interpolated linearly in energy space in order to avoid problems in case the odd moments (H or L) change their signs.

Since the radiation moments are defined on a staggered mesh, the finite-difference equations are second-order accurate in space. First-order accuracy in time is achieved by employing fully implicit time-differences ("backward Euler''). Only the advection terms in the second and forth line of Eq. (21) and Eq. (22), respectively, are treated differently (see Sect. 4.3.1 for a motivation):

$$\Xi^{n+{\frac{1}{2}}}:= (1-\zeta)~\Xi^{n}+\zeta~\Xi^{n+1} ~.$$ (25)

For stability reasons, the interpolation parameter $\zeta=0.51$ was chosen slightly larger than 0.5, which makes the treatment of the advection terms "almost'' second-order accurate in time (see, e.g. Dorfi 1998). The interface values $J_{\iota}$ and $H_{\iota}$ are taken as the corresponding upwind quantities with

$$\iota(i):= \left\{ \begin{array}{ll} i-\frac{1}{2}&\quad {\rm... ...[2mm] i+\frac{1}{2}&\quad {\rm else}~. \\ \end{array}\right.$$ (26)

In order to damp spurious oscillations behind sharp neutrino pulses travelling through nearly transparent zones, an additional diffusive term is added to the left-hand side of the first moment equation (Eq. (22)). Its finite-differenced form is (Blinnikov et al. 1998):

 

$${\cal D}_{\rm visc}=\frac{{\cal A}^n_{i,j+{\frac{1}{2}}}}{2~r_i^2... ...1}{2}}}r_{i}^2-H^{n+1}_{i-1,j+{\frac{1}{2}}}r_{i-1}^2}{r_{i}-r_{i-1}} \Bigg) ~.$$ (27)

In effect this term resembles the diffusivity of the donor cell scheme and hence allows one to selectively switch from the second order radial discretization to a first order scheme. Following Blinnikov et al. (1998) the (energy dependent) "artificial radiative viscosity'' ${\cal A}$ is set to unity in transparent zones and goes to zero for optical depth $\Delta \tau \gtrsim 1$. By construction the diffusive term vanishes when the luminosity H r² does not depend on radius and hence is only active during transient phases when sharp neutino pulses propagate across the grid.

Using a similar diffusive term in the zeroth moment equation (Eq. (21)) allowed us also to overcome stability problems with the numerical handling of the velocity-dependent terms in the general form of the radiation momentum equation, Eq. (8).

  
3.3.3 Boundary conditions

For the solution of the moment equations (Eqs. (7), (8)), boundary conditions have to be specified at $r=r_{\rm min}$, $r=r_{\rm max}$, $\epsilon=0$ and $\epsilon=\epsilon_{\rm max}$. In the radial domain the values of quantities at the boundaries are iteratively obtained from the solution of the Boltzmann equation (see Sect. 3.4.3 for the boundary conditions employed there), which has the advantage that in the moment equations no assumptions have to be made about the angular distribution of the specific intensity at the boundaries. Specifically, at the inner boundary $H_{0,j+{\frac{1}{2}}}$ is given by $H(t,r=r_{\rm min},\epsilon)$ as computed from the Boltzmann equation. To handle the outer boundary, we apply Eq. (22) on the "half shell'' between $r_{N_r-{\frac{1}{2}}}$ and r_Nr (Mihalas & Mihalas 1984). Where necessary, $J_{N_r+{\frac{1}{2}},j+{\frac{1}{2}}}$ is replaced by $H_{N_r,j+{\frac{1}{2}}}/[{f_H}]_{N_r,j+{\frac{1}{2}}}$ in Eqs. ((21), (22)).

At $\epsilon=0$ the flux in energy space vanishes and hence we simply set $\epsilon_0 u_{i+{\frac{1}{2}},0}J_{ i+{\frac{1}{2}},0}^{n+1}=0$ in Eq. (21). At the upper boundary of the spectrum the flux in energy space, $\epsilon_{N_\epsilon} u_{i+{\frac{1}{2}},N_\epsilon}J_{i+{\frac{1}{2}},N_\epsilon}^{n+1}$ is computed by a (geometrical) extrapolation of the moments $J_{i+{\frac{1}{2}},N_\epsilon-\frac{3}{2}}$ and $J_{i+{\frac{1}{2}},N_\epsilon-\frac{1}{2}}$ and by a linear extrapolation of the advection velocity using $u_{i+{\frac{1}{2}},N_\epsilon-\frac{3}{2}}$ and $u_{i+{\frac{1}{2}},N_\epsilon-\frac{1}{2}}$. Analogous expressions are used for Eq. (22).

3.3.4 Finite differencing of the source term

The finite differenced versions of the operator-splitted source terms in the energy and lepton number equations (Eqs. (15), (16)) of the stellar fluid read:

  

$$\frac{e_{i+{\frac{1}{2}}}^{n+1}-e_{i+{\frac{1}{2}}}^{n}} {t^{n+1}... ...\nu_{\rm e}},\dots\}}~ {C^{(0)}}^{n+1}_{i+{\frac{1}{2}},j+{\frac{1}{2}},\nu} ~,$$ (28)

$$\frac{{{Y_{\rm e}}}_{i+{\frac{1}{2}}}^{n+1}-{{Y_{\rm e}}}_{i+{\fr... ...k{C}^{(0)}}^{n+1}_{i+{\frac{1}{2}},j+{\frac{1}{2}},{\bar\nu_{\rm e}}} \right)~.$$ (29)

  
3.3.5 Neutrino number transport

In the finite-differenced version of the (monochromatic) neutrino energy equation (Eq. (21)) the derivative with respect to energy, $\partial/\partial\epsilon$, has been written in a conservative form. When a summation over all energy bins is performed in Eq. (21), the terms containing fluxes across the boundaries of the energy zones telescope and an appropriate finite differenced version of the total (i.e. spectrally integrated) neutrino energy equation is recovered. This essential property does, however, not hold automatically also for the neutrino number density ${\cal N}:= 4\pi/c ~ \int_0^\infty{\rm d}\epsilon~ {\cal J}(\epsilon)$, when the naive definition ${\cal J}_{i+{\frac{1}{2}},j+{\frac{1}{2}},\nu}:=J_{i+{\frac{1}{2}},j+{\frac{1}{2}},\nu}/\epsilon_{j+{\frac{1}{2}}}$is adopted. By inserting the latter expression into Eq. (21) and summing over all energy bins, it can easily be verified that the corresponding fluxes across the boundaries of the energy bins do not cancel anymore due to the finite cell size $\Delta\epsilon_{j+{\frac{1}{2}}}$.

In order to avoid this problem, we instead derive the moment equations for ${\cal J}$ and ${\cal H}$ (by multiplying Eqs. (7), (8) with $\epsilon^{-1}$) and recast them into a conservative form:

 

$$\Bigg(\frac{1}{c}\frac{\partial}{{\partial t}} \beta\frac{\partial}{{\partial r}}\Bigg)~ {\cal J} + \frac{1}{r^2... ...cal J} +\frac{1}{c}\frac{\partial \beta}{\partial t}{\cal H} =\frak{C}^{(0)} ~,$$ + (30)

 

$$\Bigg(\frac{1}{c}\frac{\partial}{{\partial t}} \beta\frac{\partial}{{\partial r}}\Bigg)~ {\cal H} + \frac{1}{r^2... ...\right) +\frac{1}{c}\frac{\partial \beta}{\partial t}{\cal J} =\frak{C}^{(1)} .$$ + (31)

Equations (30), (31) are then discretized as a separate set of equations. Consequently, the finite-difference representation of the monochromatic number density ${\cal J}_{i+{\frac{1}{2}},j+{\frac{1}{2}}}$ and the monochromatic number flux ${\cal H}_{i,j+{\frac{1}{2}}}$ are now considered as new variables, which are determined by solving an extra set of two moment equations. Note that treating ${\cal J}_{i+{\frac{1}{2}},j+{\frac{1}{2}}}$ and ${\cal H}_{i,j+{\frac{1}{2}}}$ as additional variables is not necessarily in conflict with the relations ${\cal J}(\epsilon)\equiv J(\epsilon)/\epsilon$ and ${\cal H}(\epsilon)\equiv H(\epsilon)/\epsilon$ (${\cal K}$ and ${\cal L}$are defined likewise), which constitute dependences between the quantities. Since $J_{i+{\frac{1}{2}},j+{\frac{1}{2}}}$ and ${\cal J}_{i+{\frac{1}{2}},j+{\frac{1}{2}}}$, $H_{i,j+{\frac{1}{2}}}$ and ${\cal H}_{i,j+{\frac{1}{2}}}$ are defined as average values within the energy interval $[\epsilon_j,\epsilon_{j+1}]$ of different energy moments of the same function (rather than values of a function at point $\epsilon_{j+{\frac{1}{2}}}$; cf. Eq. (20)) they can be considered as variables that are independent of each other, provided that the inequalities $\epsilon_j \le \langle\epsilon\rangle_{i+{\frac{1}{2}},j+{\frac{1}{2}}} \le \epsilon_{j+1}$ with $\langle\epsilon\rangle_{i+{\frac{1}{2}},j+{\frac{1}{2}}} = J_{i+{\frac{1}{2}},j+{\frac{1}{2}}}/{\cal J}_{i+{\frac{1}{2}},j+{\frac{1}{2}}}$ or $\langle\epsilon\rangle_{i,j+{\frac{1}{2}}} = \vert H_{i,j+{\frac{1}{2}}}\vert/\vert{\cal H}_{i,j+{\frac{1}{2}}}\vert$are fulfilled (we thank our referee, A. Mezzacappa, for stressing the importance of this point). In practice, we find that this is always well fulfilled at high optical depths and low velocities. The constraint is violated only in the steeply falling, high-energy tail of the neutrino spectrum. Systematically inspecting results from simulations we find that a part of the spectrum is affected where $J_{i+{\frac{1}{2}},j+{\frac{1}{2}}}\Delta\epsilon_{i+{\frac{1}{2}},j+{\frac{1}{2}}}$ is less than one percent of the maximum value, which means that the energy integral over this spectral tail accounts for a negligible fraction of the total neutrino (energy) density.

3.3.6 Numerical solution

With the Eddington factors f_H, f_K, f_L and flow field $\beta$ being given, the nonlinear system of equations (Eqs. (21), (22), (30), (31), (28), (29)) is solved for the variables $\{J_{i+{\frac{1}{2}},j+{\frac{1}{2}}}, H_{i,j+{\frac{1}{2}}}, {\cal J}_{i+{\fr... ...1}{2}},j+{\frac{1}{2}}}, e_{i+{\frac{1}{2}}}, {{Y_{\rm e}}}_{i+{\frac{1}{2}}}\}$ by a Newton-Raphson iteration (e.g. Press et al. 1992), using the aforementioned boundary conditions.

This requires the inversion of a block-pentadiagonal system with 2 N_r+1 rows of blocks. The dimension of an individual block of the Jacobian is $(4N_\epsilon +2)$ in case of the transport of electron neutrinos and antineutrinos (the number of energy bins $N_\epsilon$ is multiplied by a factor of two because ${\nu_{\rm e}}$ and ${\bar\nu_{\rm e}}$ are treated as separate particles, the other factor of two comes from the additional transport of neutrino number). In contrast, muon and tau neutrinos and their antiparticles are considered here as one neutrino species because their interactions with supernova matter are roughly the same. In the absence of the corresponding charged leptons, muon and tau neutrinos are only produced (or annihilated) as $\nu\bar\nu$ pairs and hence do not transport lepton number. Therefore the blocks have a dimension of $(N_\epsilon +1)$ for these flavours.

For setting up the Jacobian all partial derivatives with respect to the radiation moments can be calculated analytically, whereas the derivatives with respect to electron fraction and specific internal energy are approximated by finite differences.

  
3.4 Formal solution of the Boltzmann equation

In order to provide the closure relations (the "variable Eddington factors'') for the truncated system of moment equations, we have to solve the Boltzmann equation. Since the emissivity $\eta$ and the opacity $\chi$ depend in general also on angular moments of ${\cal I}$, this is actually a nonlinear, integro-differential problem. However, $\eta$ and $\chi$ become known functions of only the coordinates t, r, $\epsilon$ and $\mu$, when J and H are used as given solutions of the moment equations. This allows one to calculate a formal solution of the Boltzmann equation.

3.4.1 Model equations

Making the change of variables (cf. Yorke 1980; Mihalas & Mihalas 1984; Körner 1992; Basheck et al. 1997)

$$(r,\mu)\mapsto(s:=\mu r,p:=r\sqrt{1-\mu^2}),\qquad\mu\ge 0 ~,$$ (32)

and introducing the symmetric and antisymmetric averages of the specific intensity

  

$$j(t,s,p):= \frac{1}{2}\big({\cal I}(\mu)+{\cal I}(-\mu)\big) ~,$$ (33)

$$h(t,s,p):= \frac{1}{2}\big({\cal I}(\mu)-{\cal I}(-\mu)\big) ~,$$ (34)

Eq. (6) can be split into two equations for j and h,

 

$$\frac{1}{c}\frac{{\rm D}}{{\rm D} t}j+ \frac{\partial}{\partial s... ...artial r}~j +\mu \frac{2}{c}\frac{\partial \beta}{\partial t}~ h = s_{\rm E} ~,$$ (35)

 

$$\frac{1}{c}\frac{{\rm D}}{{\rm D}t}h+ \frac{\partial}{\partial s}... ...r}\right) h +\mu \frac{1}{c}\frac{\partial \beta}{\partial t}~ j = u_{\rm E} ~.$$ (36)

The symmetric ($s_{\rm E}$) and antisymmetric ($u_{\rm E}$) averages of the source term C are defined in analogy to Eqs. (33), (34). Following Burrows et al. (2000, with our Eq. (38) correcting a misprint in their Eq. (24)) we have made the replacements

  

$$\left(\mu-\mu^3\right)\left(\frac{\partial \beta}{\partial r}-\frac{\beta}{r}\right)~\frac{\partial j}{\partial \mu} \to \left(3\mu^2-1\right)\left(\frac{\partial \beta}{\partial r}-\frac{\beta}{r}\right)~j ~,$$ (37)

$$\left(\mu-\mu^3\right)\left(\frac{\partial \beta}{\partial r}-\frac{\beta}{r}\right)~\frac{\partial h}{\partial \mu} \to \left(4\mu^2-2\right)\left(\frac{\partial \beta}{\partial r}-\frac{\beta}{r}\right)~h ~,$$ (38)

in deriving Eq. (35) and Eq. (36), respectively, from Eq. (6). This modification can be motivated as follows: As a consequence of Eqs. (37), (38) no more partial derivatives of the intensity (and hence j and h) with respect to $\mu$ appear in the model equations (Eqs. (35), (36)) of the original Boltzmann equation (Eq. (6)). This quality allows us to construct an efficient numerical solution scheme to be described in the next section. At the same time it is ensured that the monochromatic ${\cal O}(\beta)$moment equations (Eqs. (7), (8)) can be recovered when Eq. (35) and $\mu$ times Eq. (36) are integrated over angles. Thus, Eqs. (35), (36) are consistent with the monochromatic moment equations without including the exact aberration terms (see also Burrows et al. 2000). Moreover, the effects of angular aberration on the solution of the neutrino transport are found to be small (Liebendörfer, personal communication). Therefore we do not expect significant changes, in particular because only normalized moments H/J, K/J, L/J are computed from the solution of the Boltzmann equation to be used as closure relations in the moment equations.

  
3.4.2 Finite-difference equations

For the moment we shall ignore terms which contain frequency derivatives $\partial/\partial\epsilon$ (see Eqs. (35), (36)). These Doppler terms will be included in an operator-split manner.

Equations (35), (36) are then discretized on a so-called "tangent ray grid'' (for an illustration, see Fig. 1), the geometry of which being an immediate consequence of the transformation of variables given by Eq. (32). Applying this transformation, partial derivatives of only one momentum-space coordinate s remain, whereas the second coordinate p appears only in a parametric way (cf. Baschek et al. 1997). This greatly facilitates the numerical solution of the system.

Figure 1: Eulerian radial and tangent ray grid of the computation (solid lines) and virtual Lagrangian shell at time level $t^{n}=t^{n+1}-\Delta t$ (dotted). For simplicity we have assumed a radial motion of only one shell with index i. The filled circles mark the positions where j and h are known as the solutions of the previous time step on the Eulerian grid. For $v\not \equiv 0$ these positions do not coincide with the locations $(r_{i}^{\rm old},\mu _{ik}^{n+1})$, where initial values of the radiation field are required for computing the current time step (open circles). The latter positions are subject to the requirement that the ${\rm D}/{\rm D}t$ operator has to be evaluated at fixed radial Lagrangian and angular coordinates (The radial dashed lines connect points with constant angle cosines $\mu$ of tangent rays.).

Open with DEXTER

A "tangent ray'' k is defined by its "impact parameter'' p_k=r_k at $\mu=0$. The coordinate s serves to measure the path length along the ray. On each tangent ray k, a staggered numerical mesh is introduced for the coordinate s. The zone boundaries (centers) of this mesh are given by the ray's intersections with the zone boundaries (centers) of the radial grid (cf. Fig. 1). With the "flux-like'' variable h being defined at the zone boundaries s_ik and the "density-like'' variable j being defined at the zone centers $s_{i+{\frac{1}{2}},k}$, the finite-differenced versions of Eqs. (35), (36) finally can be written down (cf. Mihalas & Klein 1982, Sect. V.2):

 

$$\frac{{j^*}^{n+1}_{i_k+{\frac{1}{2}},k}-\widetilde{j}^n_{i_k+{\fr... ...rac{1}{2}},k} +B_{i_k+{\frac{1}{2}},k}^{n+1}~h^{n+1}_{i_k+{\frac{1}{2}},k}=0 ~,$$ (39)

 

$$\frac{{h^*}^{n+1}_{i_k,k}-\widetilde{h}^n_{i_k,k}}{ct^{n+1}-ct^{n... ...}{2}},k}} +A_{i_k,k}^{n+1}~j^{n+1}_{i_k,k} +B_{i_k,k}^{n+1}~h^{n+1}_{i_k,k}=0~,$$ (40)

with the indices $k=K_0,\dots,N_r-1$ and $i_k=k,\dots,N_r-1$. If an inner boundary condition at some finite radius $r_{\rm min}>0$ is used, a number of "core rays'' $k=K_0,K_0+1,\dots,-1$ ( $-K_0\in {\rm I\!N}$), which penetrate the innermost radial shell ( $p_k \le r_{\rm min}$), can be defined in order to describe the angular distribution of the radiation at the inner boundary. For details, see e.g., Yorke (1980), Körner (1992), Dorfi (1998). The coefficients A and B are combinations of the velocity and angle-dependent terms and the source terms of Eqs. (35), (36). The definition of $\widetilde{j}$ and $\widetilde{h}$ will be given later, for the moment one may ignore the tilde symbol. We add a diffusive term to Eq. (40) analogous to the discussion in Sect. 3.3.2, replacing H r² in Eq. (27) by h.

The frequency derivatives of Eqs. (35), (36), which were ignored in the finite-difference versions, Eqs. (39), (40), are taken into account in a separate step by operator-splitting (the partially updated values for j and h were marked by asterisks in Eqs. (39), (40). The discretization of the corresponding terms is explicit in time and - provided the acceleration terms in the second line of Eqs. (35), (36) are neglected - can be performed straightforwardly using upwind differences in energy space ( $l=0,\dots,N_\epsilon-1$):

 

$$j^{n+1}_{\_,\_,l+{\frac{1}{2}}}={j^*}^{n+1}_{\_,\_,l+{\frac{1}{2}... ...a(l+1)} -\epsilon_{l} \upsilon_{\_,\_} ~{j^*}^{n+1}_{\_,\_,\iota(l)} \right) ~,$$ (41)

 

$$h^{n+1}_{\_,\_,l+{\frac{1}{2}}}={h^*}^{n+1}_{\_,\_,l+{\frac{1}{2}... ...a(l+1)} -\epsilon_{l} \upsilon_{\_,\_} ~{h^*}^{n+1}_{\_,\_,\iota(l)} \right) ~,$$ (42)

with

$$\iota(l):= \left\{ \begin{array}{ll} l-\frac{1}{2} &{\rm for... ...~, \\ [2mm] l+\frac{1}{2} &{\rm else}~, \\ \end{array}\right.$$ (43)

and

$$\upsilon_{i_k,k}:= (1-\mu_{i_k,k}^2)\left[\frac{\beta}{r}\rig... ...i_k,k}^2\left[\frac{\partial\beta}{\partial r}\right]_{i_k} ~,$$ (44)

where, for better readability, dashes replace the appropriate index pairs corresponding to radius and impact parameter, (i_k+1/2,k) in Eq. (41) and (i_k,k) in Eq. (42), respectively.

It should be possible to proceed along similar lines for including the exact aberration effects in the Boltzmann solution which, for the time being, includes aberration only in an approximate way on the basis of the replacements given by Eqs. (37), (38). In practice, however, the adopted tangent-ray geometry does not allow for a conservative discretization of the angular derivatives in a way as simple as in the case of the frequency derivatives. We plan to spend extra work on the omitted angular derivatives of the intensity in order to include them in future version of our code.

  
3.4.3 Boundary conditions

On the tangent ray grid boundary conditions must be specified for each ray k at s_k,k and at s_Nr,k. At the inner core radius ( $i=0;~ k=K_0,\dots, 0$ and $s_{-{\frac{1}{2}},k} :=s_{0,k}$) we set the boundary condition $j_{-{\frac{1}{2}},k}:={\cal I}^+_{0,k}-h_{0,k}$, with ${\cal I}^+_{0,k}:={\cal I}(t,r_{\rm min},\mu_{0,k})$. For the remaining rays ( $k=1,\dots, N_r$), symmetry and Eq. (34) imply h_k,k=0, since $\mu_{k,k}=0$.

At the outer radius ( $i=N_r;~ k=K_0,\dots,N_r$ and $s^{n+1}_{N_r+{\frac{1}{2}},k}:=s^{n+1}_{N_r,k}$) we consider Eq. (40) on the "half shell'' between $r_{N_r-{\frac{1}{2}}}$ and r_Nr (cf. Sect 3.3.3) and make the replacement $j_{N_r+{\frac{1}{2}},k}:={\cal I}^-_{N_r,k}+h_{N_r,k}$, with ${\cal I}^-_{N_r,k}:={\cal I}(t,r=r_{\rm max},-\mu_{N_r,k})$.

The physical boundary conditions are described by the functions ${\cal I}(t,r_{\rm min},\mu)$ and ${\cal I}(t,r_{\rm max},-\mu)$ with $0\le \mu \le 1$, which specify the specific intensity entering the computational volume at the inner and outer surfaces, respectively.

At the boundaries of the energy grid we use the same type of boundary conditions as described for the moment equations (cf. Sect. 3.3.3).

3.4.4 Numerical solution

By virtue of the approximations used to derive the model equations (Eqs. (35), (36)) and upon introducing the tangent ray grid, the system of Eqs. (39), (40) with suitably chosen boundary conditions can be solved independently for each impact parameter p_k, each type of neutrino (note that for simplicity we have dropped the index $\nu$ in our notation), and - because Doppler shift terms are split off - each neutrino energy bin $\epsilon_{j+{\frac{1}{2}}}$ (index also suppressed).

For the same reasons as detailed in Sect. 3.3, we have employed fully implicit ("backward Euler'') time differencing. Solving Eqs. (39), (40) therefore requires the separate solution of $N_k \times N_\epsilon \times N_\nu$ ( N_k:=N_r-K₀+1 is the number of tangent rays) pentadiagonal linear systems of dimension $\le$N_r. On vector computers, this can be done very efficiently by employing a vectorization over the index k. Once the numerical solution for j and h has been obtained, the monochromatic angular moments and thus the Eddington factors f_H=H/J, f_K=K/J and f_L=L/J can be computed using the numerical quadrature formulae

    

$$\int_0^1{\rm d}\mu~j(r_i,\mu)\approx \sum_{k=K_0}^{i} j_{ik} a_{ik} ~,$$ J(r_i) = (45)

$$\int_0^1{\rm d}\mu~\mu h(r_i,\mu)\approx \sum_{k=K_0}^{i} h_{ik} b_{ik} ~,$$ H(r_i) = (46)

$$\int_0^1{\rm d}\mu~\mu^2j(r_i,\mu)\approx \sum_{k=K_0}^{i} j_{ik} c_{ik} ~,$$ K(r_i) = (47)

$$\int_0^1{\rm d}\mu~\mu^3h(r_i,\mu)\approx \sum_{k=K_0}^{i} h_{ik} d_{ik} ~.$$ L(r_i) = (48)

Explicit expressions for the quadrature weights a_ik, b_ik and c_ik can be found in Yorke (1980), d_ik is calculated in analogy.

Unless the velocity field vanishes identically (in this case $\widetilde{j}\equiv j$ and $\widetilde{h}\equiv h$ in Eqs. (39) and (40), respectively), an additional complication arises, since the partial derivative ${\rm D}/{\rm D}t$ has to be evaluated not only at fixed Lagrangian radius (the index i in our case) but also for a fixed value of the angle cosine $\mu$. The angular grid $\{\mu_{ik}^n\}_{k=K_0\dots i}$, which is associated with the radius r_iⁿ, however, changes as r_i moves with time. As a consequence, e.g., the values of h at the old time level used in Eq. (40) are $\widetilde{h}^n_{ik}:=h(t^n,r_{i}^n,\mu_{ik}^{n+1})$ and not $h^n_{ik}:=h(t^n,r_{i}^n,\mu_{ik}^n)$, the solution of the previous time step (cf. Mihalas & Mihalas 1984, Sect. 98). At the beginning of the time step we therefore have to interpolate the radiation field for each fixed radial index i from the $\{\mu_{ik}^n\}_{k=K_0\dots i}$-grid onto the $\{\mu_{ik}^{n+1}\}_{k=K_0\dots i}$-grid. If an Eulerian grid is used for the computation, an additional interpolation in the radial direction is performed to map from the fixed r_i-grid to the coordinates given by $r_{i}^{\rm old}:=r_{i}-v_i\Delta t$ (see Fig. 1). This seemingly complicated procedure has the advantage of being computationally much more efficient than directly discretizing Eqs. (35), (36) in Eulerian coordinates. In this case the $v\partial_r$-term, which originates from making the replacement ${\rm D}_t=\partial_t+v\partial_r$, would couple grid points of different impact parameters and would therefore defy an independent treatment of the tangent rays.

3.5 Iteration procedure

A transport time step starts by adopting guess values (e.g., given by the solution of the previous time step) for J, H, ${Y_{\rm e}}$, and e. The quantities ${Y_{\rm e}}$, e and the density $\rho$ determine the thermodynamic state and thus the temperature and chemical potentials. These, together with J and H are used to evaluate the rhs of the Boltzmann equation. This allows one to compute its formal solution and the Eddington factors as detailed in Sect. 3.4. The Eddington factors are then fed into the system of moment and source-term equations, the solution of which (Sect. 3.3) yields improved values for J, H, ${Y_{\rm e}}$, and e. At this point, where required (e.g., for evaluating $\frak{C}^{(0)}$), the neutrino number density ${\cal J}$ and number flux density ${\cal H}$ are conventionally replaced by $\epsilon^{-1}J$ and $\epsilon^{-1}H$. This procedure is iterated until numerical convergence is achieved. With the Eddington factors being known from the described iteration procedure, the complete system of source-term and moment equations for neutrino energy and number is solved (once) in order to accomplish lepton number conservation. In this step ${\cal J}$ and ${\cal H}$ are treated as additional variables (cf. Sect. 3.3.5).

   
3.6 Details of coupling neutrino transport and hydrodynamics

Our experience with the operator-splitting technique has shown that considerable care is necessary in how precisely the equations of hydrodynamics are coupled with the neutrino transport and how the fractional time steps are scheduled. In the following we therefore describe the used procedures in some detail.

Since the hydrodynamics code and the transport solver in general have different requirements for numerical resolution, accuracy and stability, the discretization in both space and time of the two code components is preferably chosen to be independent of each other. In supernova simulations, for example, the time step limit (from the CFL condition) for the hydrodynamic evolution computed by a time-explicit method is typically more restrictive than in the time-implicit treatment of the transport part. Since the transport is computationally much more expensive, one wants to use a larger time step than for the hydrodynamics part. This means that a transport time step is divided into a suitable number of hydrodynamical substeps.

  
3.6.1 Schedule of updates

Starting from the hydrodynamic state $(\rho^n,e^n,{Y_{\rm e}}^n,v^n)$ at the old time level tⁿ, the time-explicit PPM algorithm first computes a solution of the hydrodynamical conservation laws without the effects of gravity and neutrinos. From the updated density, the gravitational potential can be calculated by virtue of Poisson's equation and finally the gravitational source terms are applied to the total energy equation (Eq. (3)) and the momentum equation (Eq. (2)). This completes the "hydrodynamics'' time step with size $\Delta t_{\rm Hyd}$, which is limited by the CFL-condition (see e.g., LeVeque 1992). By performing a total of $N_{\rm Hyd}^n$ substeps the hydrodynamic state is evolved from the time level tⁿ to the new time level tⁿ⁺¹ given by

$$t^{n+1}= t^n+\Delta t^n= t^n+\sum_{i=1}^{N_{\rm Hyd}^n}\Delta t^{i}_{\rm Hyd} ~,$$ (49)

where $N_{\rm Hyd}^n\in{\rm I\!N}$ is determined by the condition $\Delta t^n \le \Delta t^{\rm max}_{\rm Tr}$. Here $\Delta t^{\rm max}_{\rm Tr}$ is the maximum size for the transport time step, which is compatible with the constraints that the relative changes are at most 10% for the monochromatic neutrino energy and number density, and the relative changes of the density, temperature, internal energy density and electron fraction must not exceed values of a few percent during the time interval $\Delta t^n$. This leads to a partially updated hydrodynamic state $(\rho^{n+1},e^*,{Y_{\rm e}}^*,v^{*})$ which includes the effects due to hydrodynamic fluid motions and the acceleration by gravitational forces and neutrino pressure (see below). The corresponding state variables are then mapped onto the transport grid. During the transport time step of size $\Delta t^n$, the transport equations in combination with the evolution equations for the electron fraction and internal energy of the stellar medium in response to neutrino absorption and emission (Eqs. (28), (29)) are solved. This yields the source terms for energy and lepton number as $Q_{\rm E}=\rho^{n+1}\cdot (e^{n+1}-e^{*})/\Delta t^n$ and $Q_{\rm N}=\rho^{n+1}\cdot ({Y_{\rm e}}^{n+1}-{Y_{\rm e}}^{*})/\Delta t^n$. Since the transport grid in general does not necessarily coincide with the hydrodynamics grid, the new specific internal energy eⁿ⁺¹ and electron fraction ${Y_{\rm e}}^{n+1}$ as computed on the transport grid do not directly describe the new hydrodynamic state. Instead of interpolating the quantities eⁿ⁺¹ and ${Y_{\rm e}}^{n+1}$ back onto the hydrodynamics grid, we map the source terms by a conservative procedure and then update the electron number density and the total energy density on the hydrodynamics grid according to

  

$$[\rho \varepsilon]^{n+1}= \rho^{n+1}\varepsilon^{*}+ (Q_{\rm E}^{n+1}+v^{n+1}Q_{\rm M}^{n+1})\cdot\Delta t^n ~,$$ (50)

$$\rule{0pt}{0pt} [\rho {Y_{\rm e}}]^{n+1}= \rho^{n+1} {Y_{\rm e}}^*+ Q^{n+1}_{\rm N}\cdot\Delta t^n~.$$ (51)

Equations (50) and (51) express the effective influence of the source terms $Q_{\rm E}+vQ_{\rm M}$ and $Q_{\rm N}$ over the time of a transport step, but in fact we apply the corresponding sources (as given at the old time level) during each substep of the hydrodynamics. The quantities $\varepsilon^{*}$ and ${Y_{\rm e}}^*$ as needed for the transport part of the code, are recovered a posteriori by subtracting the accumulated effects of the neutrino sources again. The momentum equation of the stellar fluid is treated in a similar way. Accounting for the momentum transfer (acceleration) by neutrinos after each individual hydrodynamical substep (by using the old momentum source term $Q^n_{\rm M}$), however, is of crucial importance here. When a sizeable contribution of the neutrino pressure is ignored during a larger number of hydrodynamical substeps a potentially hydrostatic configuration can be severely driven out of mechanical equilibrium. Correspondingly, the fluid velocity v^* used in the transport step includes the effects due to the acceleration by neutrinos.

After the transport time step has been completed, the new neutrino stress $Q_{\rm M}^{n+1}$ is used for correcting v^* to give the new velocity vⁿ⁺¹ at time level tⁿ⁺¹:

$$[\rho v]^{n+1}=\rho^{n+1} v^{*}+ (Q_{\rm M}^{n+1}- Q_{\rm M}^{n}) \cdot \Delta t^n ~.$$ (52)

  
3.6.2 Arrangement of spatial grids

Radial discretization of the transport and hydrodynamic equations is done on independent grids. Therefore there is freedom to choose the number of zones and their coordinate values separately in both parts of the code. Only quantities which obey conservation laws have to be communicated from the hydro to the transport grid. By invoking the equation of state all other thermodynamical quantities can be derived from the density, the total energy density, momentum density and the number densities of electrons and nuclear species. In the other direction it is the neutrino source terms which have to be mapped back onto the hydro grid.

For the mapping procedure we assume the conserved quantities to be piecewise linear functions of the radius inside the grid cells, with parameters determined by the cell average values and so-called "monotonized slopes'' (for details see, e.g., Ruffert 1992 and references cited therein). In order to achieve global conservation of integral values we then simply average this function within each cell of the target grid of the mapping (cf. Dai & Woodward 1996).

Using this procedure in dynamical supernova simulations, we discovered spurious radial and temporal fluctuations of the temperature, electron fraction, entropy and related quantities inside the opaque protoneutron star unless the hydro and transport grids coincide exactly in this region. The scale of these fluctuations is sensitive to the radial cell size and the size of a time step, the relative amplitude was found to be typically of the order of a few percent, with the exact number varying between different quantities (see Rampp & Janka 2000). These fluctuations are understood from the fact that the mapping of the source terms on the one hand and the internal energy density (resp. temperature) on the other hand implies deviations from local thermodynamical equilibrium between neutrinos and stellar medium. The attempt of the neutrino transport to restore this equilibrium within a given grid cell leads to large net production or absorption rates, driving the temperature in the opposite direction and causing it to perform oscillations in time around a mean value. Despite these oscillations and fluctuations, the use of different numerical grids inside the protoneutron star did not cause noticeable problems for the accuracy of the "global'' evolution of our models, because the conservative mapping describes the exchange of energy and lepton number between neutrinos and the stellar fluid correctly on average.

Being cautious, we take advantage of the option of using different hydro and transport grids only during the early phases of the collapse and in regions of the star, where neutrinos are far from reaching equilibrium with the stellar fluid. In this case no visible fluctuations occur.

3.6.3 Conservation of total energy and lepton number

The numerical treatment of the radiation hydrodynamics problem should guarantee that the conservation laws for energy and lepton number are fulfilled.

Provided the acceleration terms ( $\propto\partial\beta/\partial t$), which are of order (v²/c²)and thus usually very small compared to the other terms, are ignored, the constancy of the total lepton number is ensured by, (a) a conservative discretization of the neutrino number equation (Eq. (30)), (b) a conservative handling of the electron number equation (Eq. (4)), and (c) the exact numerical balance of the source terms (cf. Eq. (13)) $-4\pi~ m_{\rm B}\int{\rm d}V~\int_0^\infty{\rm d}\epsilon~\frak{C}^{(0)}(\epsilon)$(defined on the transport grid) and $\int{\rm d}V~Q_{\rm N}$ (defined on the hydro grid). Point (a) requires that in Eq. (30) the flux divergence is discretized in analogy to the second line in Eq. (21) and that the $\beta~\partial{\cal J}/\partial r$ and $(2\beta/r+\partial\beta/\partial r)~{\cal J}$ terms are combined to ${\rm div} (\beta~{\cal J})$ to be discretized in analogy to the third line in Eq. (21). The energy derivative in Eq. (30) is treated in a conservative way as described in Sect. (3.3.5). Point (b) is achieved by the use of a conservative numerical integration of the electron number equation (Eq. (4)) in the spirit of the PROMETHEUS code, and requirement (c) is fulfilled by employing a conservative procedure for mapping the electron number source term from the transport grid to the hydro grid (see Sects. 3.6.1, 3.6.2). Doing so, the total lepton number remains constant in principle at the level of machine accuracy.

Different from the number transport, where the zeroth order moment equation for neutrinos by itself defines a conservation law, the derivation of a conservation law for the total energy implies a combination of the radiation energy and momentum equations. The use of a staggered radial mesh for discretizing the latter equations defies a suitable contraction of terms in analogy to the analytic case. Therefore our numerical description does not conserve neutrino energy with the same accuracy as neutrino number and the quality of total energy conservation has to be verified empirically for a given problem and numerical resolution.

For our supernova simulations, tests showed that neutrino number is conserved to an accuracy of better than 10^-11 per time step, while for neutrino energy a value below 10^-7 is achieved. With a typical number of about 50 000 transport time steps for a supernova simulation we thus find an empirical upper limit for the violation of energy conservation of $0.5\%$ of the neutrino energy. This translates to $0.05\%$ of the internal energy of the collapsed stellar core, i.e. a few times $10^{49}~{\rm erg}$ in absolute number. Errors of the same magnitude are introduced by the non-conservative treatment of the gravitational potential as a source term in the fluid-energy equation (Eq. (3)). Note that the use of different grids for the hydrodynamics and the transport does not affect the energy budget because we employ a conservative mapping of the neutrino source term between the grids (see Sects. 3.6.1, 3.6.2).

  
3.7 Approximate general relativistic treatment

We have not yet coupled our general relativistic version of the neutrino transport to a general relativistic hydrodynamics code. For the time being we work with a basically Newtonian code, which was extended to include post-Newtonian corrections of the gravitational potential. We hope that the deeper gravitational potential can account for the main effects of general relativity on stellar core collapse and the formation of neutron stars which do not approach gravitational instability to become black holes (cf. Bruenn et al. 2001). Because the general relativistic changes of the space-time metric are ignored, a consistent description of the neutrino transport requires that the fully relativistic equations are simplified such that only the effects of gravitational redshift and time dilation are retained.

3.7.1 Modified gravitational potential

By comparing the Tolman-Oppenheimer-Volkoff equation for hydrostatic equilibrium in general relativity (see, e.g., Kippenhahn & Weigert 1990, Sect. 2.6) with its Newtonian counterpart (cf. Eq. (2)) one can define a modified "gravitational potential'' which includes correction terms due to pressure and energy of the stellar medium and the neutrinos:

 

$$\Phi^{\rm GR}(r)= G \int\limits_\infty^r{\rm d}r'~ \frac{1}{{r'}^... ...\right) \frac{1}{\Gamma^2} \left(\frac{\rho_{\rm tot}c^2+p}{\rho c^2}\right) ~,$$ (53)

where $\rho_{\rm tot}c^2:=\rho (c^2 + e)$ is the total ("relativistic'') energy density and $P=4\pi/c\int_0^\infty {\rm d}\epsilon~ K$ the neutrino pressure. The calculation of the gravitational mass $m(r):=\int_0^r{\rm d}r'~4\pi {r'}^2(\rho_{\rm tot}+c^{-2}E+c^{-3}U F/\Gamma)$ takes into account contributions of neutrino energy density $E=4\pi/c\int_0^\infty {\rm d}\epsilon~ J$ and flux $F=4\pi\int_0^\infty {\rm d}\epsilon~ H$. The metric function $\Gamma$ is calculated as $\Gamma(r)=\sqrt{1+U(r)^2-2Gm(r)/rc^2}$ with the term U² accounting for the effects of fluid motion.

Equation (53) can be used in the Newtonian hydrodynamic equations (Eqs. (2), (3)) in order to approximately take into account general relativistic effects (cf. Keil 1997). The quality of this approach has to be ascertained empirically by comparison with fully general relativistic calculations. In our case such a comparison yields quite satisfactory results (see Sect. 4.3).

3.7.2 Approximate GR transport

The general relativistic moment equations describing transport of neutrino energy, momentum and neutrino number can be derived from the Lindquist-equation (cf. Eq. (5), Sect. 2.2.1). They are:

 

$$\frac{1}{c}\frac{{\rm D}}{{\rm D}t}~ J \frac{\Gamma}{R^2} \frac{\partial}{\partial R}(R^2H{\rm e}^{\Phi}) + \Gamma \partial_R {\rm e}^{\Phi}~ H$$ +

$$\frac{\partial}{\partial \epsilon}\left[ \epsilon~\Big( {\rm e}^{... ...+c^{-1} {\rm D}_t \Lambda~ K +\Gamma \partial_R {\rm e}^{\Phi}~ H \Big) \right]$$ -

$${\rm e}^{\Phi}\frac{U}{R}~(3J-K) +c^{-1}{\rm D}_t \Lambda~(J + K) ={\rm e}^{\Phi}~C^{(0)} ~,$$ + (54)

 

$$\frac{1}{c}\frac{{\rm D}}{{\rm D}t}~ H \frac{\Gamma}{R^2}\frac{\partial}{\partial R}(R^2K{\rm e}^{\Phi}) + \Gamma\partial_R {\rm e}^{\Phi}~ J +{\rm e}^{\Phi}\Gamma~\frac{K-J}{R}$$ +

$$\frac{\partial}{\partial \epsilon}\left[ \epsilon~\Big( {\rm e}^{... ...+c^{-1} {\rm D}_t \Lambda~ L +\Gamma \partial_R {\rm e}^{\Phi}~ K \Big) \right]$$ -

$$2({\rm e}^{\Phi}\frac{U}{R} +c^{-1}{\rm D}_t \Lambda)~H = {\rm e}^{\Phi}~C^{(1)} ~,$$ + (55)

for the energy transport, and

 

$$\frac{1}{c}\frac{{\rm D}}{{\rm D}t}~ {\cal J} \frac{\Gamma}{R^2} \frac{\partial}{\partial R}(R^2{\cal H}{\rm e}^{\Phi})$$ +

$$\frac{\partial}{\partial \epsilon}\left[ \epsilon~\Big( {\rm e}^{... ...D}_t \Lambda~ {\cal K} +\Gamma\partial_R {\rm e}^{\Phi}~ {\cal H} \Big) \right]$$ -

$$(2~{\rm e}^{\Phi}\frac{U}{R}+c^{-1}{\rm D}_t \Lambda)~{\cal J} ={\rm e}^{\Phi}~\frak{C}^{(0)} ~,$$ + (56)

 

$$\frac{1}{c}\frac{{\rm D}}{{\rm D}t}~ {\cal H} \frac{\Gamma}{R^2}\frac{\partial}{\partial R}(R^2{\cal K}{\rm e}^... ...rac{{\cal K}-{\cal J}}{R} +\Gamma \partial_R {\rm e}^{\Phi}~({\cal J}-{\cal K})$$ +

$$\frac{\partial}{\partial \epsilon}\left[ \epsilon~\Big( {\rm e}^{... ...}_t \Lambda~ {\cal L} +\Gamma \partial_R {\rm e}^{\Phi}~ {\cal K} \Big) \right]$$ -

$${\rm e}^{\Phi}\frac{U}{R}~({\cal H}+{\cal L}) + c^{-1} {\rm D}_t \Lambda~(2{\cal H}-{\cal L}) = {\rm e}^{\Phi}~\frak{C}^{(1)} ~,$$ + (57)

for the number transport.

In the approximate treatment we neglect the distinction between coordinate radius and proper radius ( $\partial_R \to \partial_r$, $\Gamma=1$) in Eqs. (54)-(57), and identify corresponding quantities with their Newtonian counterparts ($R\to r$, ${\rm e}^{\Phi}~U\to\beta$, $c^{-1}{\rm D}_t\Lambda \to \partial \beta/\partial r$). The same approximations are made in the "parent'' Boltzmann equation from which the moment equations of the relativistic approximation can be consistently derived. Accordingly, the approximations to Eqs. (54)-(57) contain only general relativistic redshift and time dilation effects for neutrinos. Coupling the transport with the Newtonian equations of hydrodynamics, these restrictions to a fully relativistic treatment are necessary in order to verify conservation laws for energy and lepton number of the coupled system.

Finite-difference versions of the moment equations and the corresponding parent Boltzmann equation for the approximate GR transport are obtained by applying the techniques described in Sects. 3.3 and 3.4. The lapse function is calculated by integrating the general relativistic Euler equation $\partial \ln~ {\rm e}^{\Phi}/\partial r=-(\rho_{\rm tot}c^2+p)^{-1}~(\partial p/\partial r-Q_{\rm M}/\Gamma)$ inward from the surface, where the boundary condition ${\rm e}^{\Phi}=\Gamma$ is applied (van Riper 1979).

  
3.8 Approximate multi-dimensional neutrino transport

Multi-dimensional frequency dependent neutrino transport in moving media and relativistic environments is a challenging problem for future work. Since convective phenomena were recognized to be highly important in supernovae (Herant et al. 1994; Burrows et al. 1995; Janka & Müller 1996; Keil et al. 1996 and refs. therein) one would of course like to perform simulations with Boltzmann neutrino transport also in two and three dimensions. Here we suggest an approximate approach based on a straightforward generalization of our variable Eddington factor method, which offers some advantages concerning computational efficiency. The approximation should be considered as an intermediate step between spherically symmetric and fully multi-dimensional models. Of course, the quality of the approximation for the neutrino transport which we will describe below, will finally have to be checked by a comparison with fully multi-dimensional transport calculations.

3.8.1 Basic considerations

Our approximate treatment may be a reasonably accurate and physically justifyable approach for describing neutrino transport in situations where the star does not show a global deformation (e.g., due to rotation) in layers which are opaque to neutrinos, but where inhomogeneities and anisotropies are present only on smaller scales (e.g., due to convection). Multi-dimensional hydrodynamical simulations suggest that convective processes occur in two distinct regions of the supernova core:

(a)

Convection inside the opaque protoneutron star causes deviations of the structural and thermodynamical quantities (like density or temperature) from spherical symmetry typically of the order of a few percent (Keil et al. 1996; Keil 1997). Moreover, the time scale for changes of such local fluctuations is short compared to the neutrino diffusion time scale. Hence, on the neutrino diffusion time scale no persistent local gradients in lateral and azimuthal directions are present in the dense stellar interior. This allows one to disregard the neutrino transport in these directions in a first approximation, which is at least correct in the temporal average.

(b)

Convective overturn motions occur between the neutrinospheres and the stalled hydrodynamic shock. Deviations from spherical symmetry are present there on larger angular scales (Herant et al. 1992; Herant et al. 1994; Burrows et al. 1995; Janka & Müller 1996; Müller & Janka 1997; Mezzacappa et al. 1998b; Kifonidis et al. 2000), but the neutrinos are only loosely coupled to the stellar medium in these regions. Therefore neutrino transport in the lateral directions does not play an important role and cannot have a significant influence on the local emission and absorption of neutrinos.

3.8.2 "Ray-by-ray'' transport

Under these circumstances the specific intensity ${\cal I}(t,r,\vartheta,\varphi, \epsilon,\vec{n})$can be assumed to depend mainly on r but only weakly on longitude $\varphi$ and latitude $\vartheta$ of the background medium. Hence, like in the spherically symmetric case, the dependence of the specific intensity on the direction of propagation $\vec{n}$ can be described by only one angle $\mu:=\vec{n}\cdot\vec{r}/\vert\vec{r}\vert$. The flux is thus approximated as $\vec{F}=4\pi\vec{r}/\vert\vec{r}\vert\cdot H(t,r,\vartheta,\varphi, \epsilon)$ and the scalar $P=4\pi/c\cdot K(t,r,\vartheta,\varphi, \epsilon)$ is sufficient to define the radiation stress tensor.

In the moment equations, gradients in the $\vartheta$- and the $\varphi$-direction, which describe the transport of energy and neutrino number into these lateral and azimuthal directions, are neglected, yet the parametric dependence of the (scalar) moments on the coordinates $\vartheta$ and $\varphi$ is retained and the radial transport is computed using the moment equations independently in each angular zone of the stellar model. In order to close the set of moment equations, variable Eddington factors are computed by iterating the Boltzmann equation and the corresponding moment equations on a spherically symmetric "image'' of the stellar background. The latter is defined as the angular average of physical quantities $\xi\in\{\rho,T,{Y_{\rm e}},\beta,\dots\}$ according to $\xi(t,r):=(4\pi)^{-1}\int_{-1}^{+1}{\rm d}\cos \vartheta~\int_{0}^{2\pi}{\rm d}\varphi~\xi(t,r,\vartheta,\varphi)$. Note that the variable Eddington factors are normalized moments of the neutrino phase space distribution function and thus should not show significant variation with the angular coordinates of the star. This justifies replacing them by quantities that depend only on the radial position and time.

Since for each latitude $\vartheta$ and longitude $\varphi$ the moment equations (Eqs. (7), (8)) in our approach are solved together with the evolution equations of electron fraction and internal energy due to neutrino sources (Eqs. (15), (16), local radiative equilibrium can be attained properly and conservation of energy can be fulfilled. It is not obvious to us how these fundamental requirements could be met in a more simple approximation where one uses a one-dimensional transport scheme to compute the transport on a spherically symmetric "mean star'', which is obtained at each time step by averaging the multi-dimensional hydrodynamical stellar model over angles.

3.8.3 Computational efficiency

Besides having significant advantages for easy and efficient implementation on vector and parallel computer architectures, the suggested approach also saves computer time compared to a multiple application of a one-dimensional Boltzmann solver. Let $\Delta t^{\rm CPU}_{1{\rm D}}$ be the computation time required for calculating a time step of the full one-dimensional transport problem and let $n^{\rm It}\in {\rm I\!N}$ be the number of steps that are necessary for the iteration between the moment equations and the Boltzmann equation to achieve convergence. If $n_{\vartheta,\varphi}:=n_{\vartheta}\cdot n_{\varphi}$ is the total number of angular zones that discretize the background star, the computation time for calculating one time step of the multi-dimensional problem with the described approximate method is

$$\Delta t^{\rm CPU}= \left(1+\frac{n_{\vartheta,\varphi}}{n^{\rm It}}\right)\cdot \Delta t^{\rm CPU}_{1{\rm D}} ~.$$ (58)

Note that the iteration procedure ( $n^{\rm It}> 1$) has to be performed only once, namely for calculating the Eddington factors on the angularly averaged "image'' of the background medium. With typical values of $n^{\rm It}=3\dots 10$ and large values of $n_{\vartheta,\varphi}$($\simeq$100) the computational effort for solving the multi-dimensional problem becomes about

$$\Delta t^{\rm CPU}\approx (0.1\dots 0.3)~ n_{\vartheta,\varphi}\cdot \Delta t^{\rm CPU}_{1{\rm D}} ~,$$ (59)

to be compared with $\Delta t^{\rm CPU}= n_{\vartheta,\varphi}\cdot \Delta t^{\rm CPU}_{1{\rm D}}$, which would be necessary for a straightforward multiple application of the one-dimensional method to a number of $n_{\vartheta,\varphi}$ angular zones.

  
4 Test cases

Stationary solutions of the Boltzmann equation of radiative transfer can be calculated analytically, semi-analytically or numerically to high accuracy in a few cases, where a particularly simple form of the emissivity and the opacity is assumed. These problems provide test cases for our new code for solving the equations of neutrino transport. To this end we impose suitably chosen initial and boundary conditions and solve the time-dependent neutrino transport equations numerically, until a stationary state is reached. The numerical solutions can then be compared to the reference solutions. This first class of problems, presented in Sects. 4.1 and 4.2, does neither refer to specific properties of neutrinos nor will the equations of neutrino radiation hydrodynamics, i.e. the coupling of the neutrino transport equation to the equations of hydrodynamics, be tested. Time-dependent hydrodynamical tests will be described in Sect. 4.3.

  
4.1 Radiative transfer on a static background

We first consider the simplest class of radiative transfer problems, namely those where the "background medium'' is assumed to be time-independent and static, i.e. the radial profiles of the emissivity and opacity do not change with time and the velocity and acceleration vanish everywhere and at all times. As usual, the background medium is assumed to be spherically symmetric.

Propagation of a light pulse:

Here, our method for computing the variable Eddington factor (see Sect. 3.4) by solving the Boltzmann equation with a finite difference scheme (Mihalas & Klein 1982) is compared with results obtained from a general quadrature solution along characteristics (Yorke 1980; Körner 1992; Baschek et al. 1997). For convenience, we set $c\equiv 1$ in this section.

As a test problem we consider a central point source which emits radiation with an intensity ${\cal I}_0$ into a spherical, static and homogeneous shell of matter bounded by two spheres of radius r₀ and R > r₀. The only interaction of the radiation with the medium is by absorption ( $C=-\kappa_{\rm a}~{\cal I}$). The central source is active during the time interval $-r_0 \le t < \infty$. Initially (t=0), there is no radiation inside the shell.

By symmetry, the intensity vanishes for all but the radial direction of propagation. In the absence of scattering, it is therefore sufficient to follow the propagation of the light pulse along a single radial (characteristic) ray. The dependence of the intensity on the angle can thus conveniently be suppressed in the notation, writing

$${\cal I}^+(t,s):= \delta(\mu-1)\cdot{\cal I}(t,s,\mu)~,$$

$${\cal I}^-(t,s):= \delta(\mu+1)\cdot{\cal I}(t,s,\mu) ~,$$ (60)

with s:=r-r₀ and $\delta(x)$ being the Dirac $\delta$-function. The equations to be solved are (cf. Eq. (6))

  

$$\partial_t{\cal I}^+(t,s)+\partial_s{\cal I}^+(t,s) -\kappa_{\rm a}{\cal I}^+(t,s) ~,$$ = (61)

$$\partial_t{\cal I}^-(t,s)-\partial_s{\cal I}^-(t,s) -\kappa_{\rm a}{\cal I}^-(t,s)~,$$ = (62)

subject to the boundary conditions ${\cal I}^+(t,0)={\cal I}_0$ (for $t \ge 0$) and ${\cal I}^-(t,s_{\rm max})=0$, and the initial condition ${\cal I}^{\pm}(0,s)=0$ (for $0 < s < s_{\rm max}:=R-r_0)$. The analytical (weak) solution can easily be verified to read:

  $$\begin{eqnarray*} {\cal I}^+(t,s)= \left\{ \begin{array}{ll} {\cal I}_0\cdo... ...t, \\ 0 &\quad {\rm else,} \\ \end{array} \right.\nonumber \end{eqnarray*}$$

and

$${\cal I}^-(t,s)\equiv 0.$$ (63)

With our implementation of the finite-difference method adopted from Mihalas & Klein (1982), we computed two models with a shell spanning the range $0 < s \le 200$. We use an equidistant radial grid with a resolution of $\Delta s=1$. Model "VAC'' assumes a vacuum shell ( $\kappa_{\rm a}\equiv 0$), whereas the absorbing shell of model "ABS'' is characterized by $\kappa_{\rm a}\equiv 0.01$, resulting in a total optical depth of $\tau=2$.

Figure 2 shows our results for ${\cal I}^+$ at the time t=100, together with the analytical solution (Eq. (63)). We define the position of the numerically broadened light front as the radius, where the average of the true pre and postfront value of the intensity is reached. In all cases the mean propagation speed is well reproduced, with some minor loss of accuracy for very large time steps. The shape of the light front, however, deviates from the true solution. One observes an artificial precursor together with a reduced intensity behind the front, both effects resulting in a spatial broadening of the front. A clear trend towards larger diffusive broadening with increasing time steps can be seen from Fig. 2 and Table 1. This phenomenon was also observed by Mihalas & Klein (1982), who reported very similar results for their tests.

Table 1 summarizes our results for the width of the front and compares them to simulations done with an implementation of the method of characteristics (Yorke 1980; Körner 1992). The latter method has the advantage to exactly reproduce the analytical solution without any smearing of the front, if the light front is advanced by exactly one grid cell per time step. Though the ability to reproduce the exact solution appears to be a very appealing property of the method of characteristics, the necessary condition $\Delta t=\Delta s$ can hardly ever be accomplished in realistic simulations. This is simply due to the fact that in typical astrophysical simulations the radial resolution can be varying over several orders of magnitude, whereas in most applications the global value of the time step is determined in some region of the grid. When doing radiation hydrodynamical simulations, for example, one is usually interested in resolving time scales different from those given by the speed of propagating light fronts. Accordingly, the time steps can be very different from the optimum for describing the light front propagation. Our finite difference scheme therefore seems to be preferable to the method of characteristics as far as the resolution of travelling discontinuities is concerned (cf. Table 1).

Figure 2: Outwardly directed intensity ${\cal I}^+$ (normalized to the true postfront value ${\cal I}_0$) as a function of the position s, calculated with our Boltzmann solver at time t=100 for different values of the computational time step. The bold line represents the analytical solution. Panel  a) shows Model "VAC'' with no absorption, in Panel  b) Model "ABS'' with the absorptive shell is displayed. Employing an additional diffusive term in the transport equations helps damping the spurious post-front oscillations, which are most prominent for $\Delta t/\Delta s=0.02$ (cf. the solid line compared to the dotted line in the plot inserted in Panel  a)).

Open with DEXTER

For time steps $\Delta t$ significantly smaller than the radial zone width, a case not discussed by Mihalas & Klein (1982), we observe some spurious postfront oscillations (see also Mihalas & Weaver 1982, Sect. 5) . Their maximum amplitude, however, does not grow with time, which is in agreement with a local von Neumann stability analysis for the fully implicit implementation of the equations (Mihalas & Klein 1982). If present in practical applications, the oscillations can be damped by employing an additional diffusive term in Eq. (62), which is active in nearly transparent regions ( $\tau\lesssim 1$) of the computational grid (cf. Sects. 3.3.2, 3.4.2; Blinnikov et al. 1998). Due to the correspondingly increased diffusivity of the finite-difference scheme the light front is smeared over a larger number of zones ( $\delta s=18$ for $\Delta t/\Delta s=0.02$; see also the inserted panel in Fig. 2a) as compared to the original method of Mihalas & Klein (1982). Note, however, that the representation of the front is still better than the results obtained with the method of characteristics.

   

Table 1: Width of the light front $\delta s$ (defined as the number of grid points for which ${\cal I}^+$ lies between 0.1 and 0.9 times its true postfront value) as a function of the time step $\Delta t$ (or the number of time steps # required for the light front to reach s=100). Computations with the finite difference method ("FD''; Mihalas & Klein 1982) and with the method of characteristics ("CH''; Yorke 1980) are compared. Numbers obtained with the latter method are taken directly from Yorke (1980, Tables 1, 2). They were, in addition, confirmed by an implementation of this numerical scheme by Körner (1992).

		model VAC		model ABS
$\Delta t/\Delta s$	#	$\delta s$ (FD)	$\delta s$ (CH)	$\delta s$ (FD)	$\delta s$ (CH)
0.02	5000	6	25	6	26
0.1	1000	8	-	8	-
0.5	200	18	16	16	20
1.0	100	26	1	22	1
2.0	50	37	26	29	26
10.0	10	80	-	55	-

  
Homogeneous sphere:

The so-called "homogeneous sphere'' is a test problem frequently employed for radiative transfer calculations (e.g., Bruenn 1985; Schinder & Bludman 1989; Smit et al. 1997). Physically, one can think of a static, homogeneous and isothermal sphere of radius R radiating into vacuum. Inside the sphere, the only interactions of the radiation with the background medium are isotropic absorption and thermal emission of radiation.

Despite of its simplification, the model has some important numerical and physical properties which are typically found in practical applications: Finite difference methods are challenged by the discontinuity at the surface of the sphere. Moreover, the sudden transition from radiation diffusion inside an optically thick sphere to free streaming in the ambient vacuum - a similar but less extreme situation arises in the neutrino-heating region of a core-collapse supernova - is a major source of inaccuracies for approximate radiative transfer methods like, e.g., flux-limited diffusion.

The model is defined by setting $C=\kappa_{\rm a}(b-{\cal I})$ with $b={\rm const.}$, $\kappa_{\rm a}={\rm const.}$ for $0\le r \le R$, and $\kappa_{\rm a}=b\equiv 0$ for r > R. Since for this choice of parameters, the right hand side of the Boltzmann equation does not contain terms that depend explicitly on the angular moments of the radiation field, the solution of Eq. (6), with $\beta\equiv 0$, can be obtained analytically by computing a formal solution. With the boundary conditions ${\cal I}(r=0,\mu=1) = {\cal I}(r=0,\mu=-1)$ and ${\cal I}(r=R_{\rm max},-1\le \mu \le 0) = 0$, which account for isotropy (due to spherical symmetry) of the radiation field at the center and no incoming radiation at the outer boundary, respectively, the result for the stationary state is (see e.g., Smit et al. 1997):

$${\cal I}(r,\mu)=b~(1-{\rm e}^{-\kappa_{\rm a}~s(r,\mu)}) ~.$$ (64)

Here s is defined as

  $$\begin{eqnarray*} s(r,\mu) := \left\{ \begin{array}{ll} r\mu+R g(r,\mu) & {\rm ... ...7em}x \le\mu \le +1, \\ 0 & {\rm else,} \\ \end{array}\right. \end{eqnarray*}$$

$${\rm with}\qquad g(r,\mu) :=\sqrt{1-\left(\frac{r}{R}\right)^2(1-\mu^2)}$$ (65)

$$\begin{eqnarray*}{\rm and}\qquad x :=\sqrt{1-\left(\frac{R}{r}\right)^2} ~. \end{eqnarray*}$$

The solution given by Eqs. (64), (65) depends on only three independent parameters, namely the radial position relative to the radius of the sphere, r/R, its total optical depth $\tau=\kappa_{\rm a} R$, and the equilibrium intensity b. The latter merely acts as a scale factor of the solution.

We employ a radial grid consisting of 213 radial zones to cover the range between r=0 and $r=R_{\rm max} \approx 12~~R$. Approximately 200 zones were distributed logarithmically between $r\approx0.0006~~R_{\rm max}$ and $r=R_{\rm max}$, about two thirds of which were spent to resolve the sphere. Initially we set ${\cal I}\approx 0$ everywhere and evolved the radiation field until a stationary state was reached.

Figure 3: Numerical (dotted and dash-dotted lines) and analytical (solid lines) stationary-state solutions vs. radius (normalized to the radius $R_{\rm max}$ of the outer boundary) of homogeneous sphere problems with "low'' ($\tau =4$, left column, Panels a), c)) and "high'' opacity ($\tau =26$, right column, Panels b), d)). Panels a) and b) show the zeroth moment of the intensity. In Panels c) and d) the Eddington factor K/J (dash-dotted) and the flux factor H/J (dotted) are displayed. The inserted panels show the flux factor on a logarithmic scale. Thin vertical lines mark the surface R of the homogeneous sphere.

Open with DEXTER

In Fig. 3 we display the stationary-state solutions for two models, one with $\tau=b=4$ (Figs. 3a, c) and another with $\tau=b=26$ (Figs. 3b, d), representative of spheres with low and high optical depth, respectively. In general, our numerical results agree very well with the analytical solutions. The Eddington factor K/J (as calculated from the Boltzmann equation) and the flux factor H/J (as obtained from the moment equations) both follow the exact values very closely (Figs. 3c, d), in case of the flux factor over many orders of magnitude. In both models, the transition to free streaming is reproduced very accurately (Figs. 3c, d). This region usually causes serious problems for flux-limited diffusion methods (see e.g., Bruenn 1985; Smit et al. 1997). Also the forward peaking of the radiation distribution at large distances from the surface of the sphere (cf. Eqs. (64), (65)) is described excellently by our method: at $r=R_{\rm max}$, the tangent ray grid yields approximately 150 angular grid points to resolve the radiation field from the central source, which has an opening half angle of only $\Theta\approx 5^\circ$ ( $\cos\Theta\approx 0.996$). These numbers should be compared with calculations employing the discrete angles ( $S_{\rm N}$) method: In time dependent neutrino transport calculations, typically less than 10 angular grid points can be afforded to equidistantly cover the range $0\le \cos\Theta \le 1$.

Although we have used a geometrical radial zoning for our tangent-ray grid, we do not find any systematic effects caused by a "bias'' in angular binning as described by Burrows et al. (2000). Far from a central source, they report Eddington factors and flux factors that asymptote to between 0.96 and 0.98 instead of 1.0 in this case. For example, we find values of K/J=0.996602 for the Eddington factor and H/J=0.999361 for the flux factor at $r=R_{\rm max}$ (in the model with $\tau =4$), to be compared with K/J=0.996636 and H/J=0.998626 obtained from the analytical solution (Eqs. (64), (4.1)).

In case of the mean intensity J we observe a systematic trend towards larger deviations from the true solution for spheres with larger total optical depth $\tau$. In our "low opacity'' case ($\tau =4$, see Fig. 3a), there is hardly any difference visible, whereas in the "high opacity'' case ($\tau =26$, see Fig. 3b) the numerical solution slightly overestimates the true solution in a narrow region beneath the surface of the sphere and underestimates it in the ambient vacuum region. A detailed analysis of the data yields values of $6\%$ and $10\%$ for the relative deviations from the analytical solutions in the "low opacity'' and the "high opacity'' case, respectively.

This finding is caused by the steep opacity gradient near the surface of the radiating sphere. Since we keep our radial grid to be the same in all models, the transition from high optical depth ( $\tau(r)=\kappa_{\rm a}~(R-r)\gtrsim 1$, if r<R) to transparency ($\tau(r)=0$, if r>R) occurs in a boundary layer beneath the surface of the sphere, which is less well resolved when the opacity is very large. In the "high opacity'' case, the semi-transparent layer is geometrically much thinner than in the "low opacity'' case and therefore the radial gradients of the radiation density are steeper and would require more radial zones for a proper description. In fact, in the "low opacity'' case the radial zones of our grid near the surface of the sphere are optically thin ( $\Delta\tau\approx 0.15$), whereas in the "high opacity'' case the outermost zone has already an optical depth of $\Delta\tau\approx 2.3$. A better quality of the numerical results could be achieved by increasing the spatial resolution in the vicinity of large opacity or emissivity gradients and/or by using higher-order difference schemes.

Due to the absence of scattering in the discussion of the homogeneous sphere problem, the solution of the Boltzmann equation does not require any information from the moment equations. No iteration between both parts of the code is necessary. Therefore this problem allows one to test the algorithm which solves the Boltzmann equation for the radiation intensity independently from the numerical solution of the moment equations. For the stationary state, we find that the radiation moments calculated by a numerical integration of the intensities are consistent with the moments directly obtained from the moment equations to within an accuracy of less than a percent.

  
Static scattering atmospheres:

Hummer & Rybicki (1971) published stationary-state solutions for the spherical analogue of the classical Milne problem. The model comprises a static, spherically symmetric, pure scattering atmosphere of some radius $R_{\rm max}$, with a central point source that is emitting radiation isotropically with a constant luminosity L₀. Due to the presence of scattering, the problem is of integro-differential nature and defies solution by simply computing the formal solution of the Boltzmann equation.

The opacity of the atmosphere is assumed to be solely caused by isotropic scattering with a simple power-law dependence on radius:

$$\chi(t,r)\equiv\kappa_{\rm s}(r)=r^{-n}, \quad 0 < r\le R_{\rm max}, \quad (n > 1) ~.$$ (66)

By a redefinition of the unit of length we have set a scale factor $\alpha>0$ to unity. This factor appears in the original form $\kappa_{\rm s}(r,t)=\alpha r^{-n}$ used by Hummer & Rybicki (1971, Eq. (1.1)). The emissivity vanishes within the atmosphere ( $\eta\equiv 0$) and no radiation is entering from outside ( ${\cal I}(t,R_{\rm max},\mu)=0$, for $-1\le \mu<0$).

For a number of atmospheres defined by various combinations of the parameters $n\in\{1.5,2,3\}$, $R_{\rm max}\in\{0.1,\dots,100\}$ and $L_0=(4\pi)^2$, Hummer & Rybicki (1971) computed numerically the zeroth moment J and the Eddington factor K/J of the stationary-state solution as a function of radius. They found values of better than one percent for the accuracy of their results. The stationary-state solution for the first moment H as a function of radius can easily be derived analytically from the zeroth order moment equation (Eq. (7), with $\beta\equiv 0$):

$$H(r)=\frac{L_0}{(4\pi)^2}~\frac{1}{r^2} ~\cdot$$ (67)

This means that the luminosity is conserved ( $L(r):=4\pi r^2\cdot 4\pi H(r)=L_0$). Hummer & Rybicki (1971) also gave useful asymptotic expressions applicable to regions of very low and very high optical depth, respectively.

The idealized concept of a central point source with a given luminosity L₀ is in practice modeled by imposing a suitable inner boundary condition at some finite radius $R_{\rm min}$, which bounds the atmosphere from below. Since all atmospheres with a scattering opacity according to Eq. (66) become optically thick at sufficiently small radius, it is reasonable to employ the diffusion ansatz

$${\cal I}(t,R_{\rm min},\mu)={\cal I}_0+\mu {\cal I}_1$$ (68)

for the radiation field at the inner boundary. In our program, only the quantity $h(t,R_{\rm min},\mu)=\mu {\cal I}_1$needs to be specified. Using Eq. (67) and $H=\int_0^1{\rm d}\mu~\mu h$(cf. Eq. (46)), the parameter ${\cal I}_1$ is easily verified to be ${\cal I}_1=3H(R_{\rm min})=3~L_0/(4\pi R_{\rm min})^2$.

We evolved the transport to stationary-state solutions for the two sets of parameter combinations (n=1.5, $R_{\rm min}=0.001$, $R_{\rm max}\in\{0.1,1,10,100\}$), and (n=2, $R_{\rm min}=0.01$, $R_{\rm max}\in\{0.1,1,10,20\}$). The total optical depth at $r=R_{\rm min}$ is larger than 55 for all models of the former class (n=1.5) and larger than 90 for all of the models with n=2. Hence, we verify that the inner boundary is placed at a radius which is sufficiently small to justify the use of the diffusion approximation (Eq. (68)) for the inner boundary condition. In order to test the sensitivity of the results to the corresponding angular dependence of the intensity (Eq. (68)), we have calculated a model with a different angular dependence: ${\cal I}(t,R_{\rm min},\mu)=\delta(\mu-1)\cdot{\cal I}_{\rm in}$. The parameter ${\cal I}_{\rm in}$ is chosen such that $L(R_{\rm min})=L_0$. This prescription is appropriate for an atmosphere around a central hollow sphere with radius $R_{\rm min}$, which is irradiated from below by the central point source. Figure 4d shows the Eddington factor of the stationary state for this model. The numerical solution deviates visibly from the reference solution only in the innermost radial zones.

Figure 4: Stationary state solutions for selected radiation quantities as a function of radius obtained with our radiative transfer code (solid lines). Our results are compared with the reference solutions (crosses) and asymptotic solutions (dashed lines) of Hummer & Rybicki (1971). Panel  a): mean intensity J (times r2) for the combination of model parameters $(n=2,R_{\rm min}=0.01,R_{\rm max}=0.1)$; Panel b): mean intensity J for $(n=2,R_{\rm min}=0.01, R_{\rm max}=10)$; Panel c): r2 J for models with $(n=3/2,R_{\rm min}=0.001, R_{\rm max}\in \{0.1,1,10,100\})$. The vertical arrows indicate the radial positions, where $\tau =1$ is reached in the different atmospheres; Panel d): Eddington factor K/J for the test case with a purely radial distribution of the radiation intensity at the inner boundary and the model parameters $(n=2,R_{\rm min}=0.01, R_{\rm max}=10)$. Note that the solid line connects values which are evaluated at the centers of the radial zones of our computational grid, the first one being located at $r_{{\frac{1}{2}}}=0.05$.

Open with DEXTER

In the models with n=1.5 the standard radial grid consisted of 200 logarithmically distributed zones with some additional zones giving higher resolution near the surface of the atmosphere, whereas 250 equidistant radial zones were used for the models with n=2. Additional tests with 500 radial grid points were performed for some selected models without finding significant changes. Figure 4 shows some results of our tests together with the data of Hummer & Rybicki (1971). The agreement is excellent. In all models we were able to reproduce the analytical solution for the flux (Eq. (67)) with a relative accuracy of at least 10^-4 throughout the atmosphere.

Radiative transfer in general relativity:

Figure 5: Metric functions $\Gamma$ (dashed line) and ${\rm e}^\Phi$ (solid line) of the background model (Panel  a)) and stationary-state solutions for selected radiation quantities ((Panels  b), c)) versus radius. Panel  b) shows the zeroth moment J (solid line) and first moment H (dashed line) of the intensity normalized to the non-relativistic solutions. The reference solutions of Schinder & Bludman (1989) are given by the symbols (J: triangles; H: diamonds). Error bars indicate the estimated uncertainties of reading off the reference solutions from the plots in Schinder & Bludman (1989). The Eddington factor vs. radius is shown as a solid line in Panel  c) together with the solution of Schinder & Bludman (1989, Fig. 13; diamonds).

Open with DEXTER

Schinder & Bludman (1989) considered the effects of a strong gravitational field in some of the radiative transfer problems discussed above, in particular in the cases of a homogeneous sphere and a static scattering atmosphere. They computed stationary-state solutions numerically for the general relativistic equations of radiative transfer and compared the results to the weak field limit. For a static background ( $\partial_t R=\partial_t \Lambda\equiv 0$), and by application of a change of coordinates $(R,\epsilon)\mapsto(R,\epsilon_\infty:={\rm e}^{\Phi}\epsilon)$, Schinder & Bludman (1989) simplified the general relativistic moment equations (Eqs. (54), (55)) to

  

$$\frac{1}{c}\frac{{\rm D}}{{\rm D}t} J \frac{\Gamma}{R^2}\frac{\partial}{\partial R}(R^2H{\rm e}^{\Phi}) ={\rm e}^{\Phi}~C^{(0)} ~,$$ + (69)

$$\frac{1}{c}\frac{{\rm D}}{{\rm D}t} H \frac{\Gamma}{R^2}\frac{\partial}{\partial R}(R^2K{\rm e}^{\Phi})... ...ma\frac{\partial {\rm e}^{\Phi}}{\partial R}~(J-K) = -{\rm e}^{\Phi}~C^{(1)} ~.$$ + (70)

The new coordinates have the advantage that the moment equations decouple in energy space (if energy-changing neutrino-matter interactions are ignored). Hence, all radiation quantities depend only parametrically on $\epsilon_\infty$.

Different from the moment equations, which we treat in their most general form (Eqs. (54), (55)), our implementation of the tangent ray scheme for computing the variable Eddington factors employs a number of approximations. In particular, by using straight lines for the tangent rays, general relativistic ray bending is not included. The current tests therefore serve the purpose of clarifying the influence of these approximations on the quality of the solutions of the moment equations. Schinder & Bludman (1989) found significant differences of the Eddington factors for the general relativistic and Newtonian simulations of the scattering atmosphere (while in case of the homogeneous sphere the differences were minor). The scattering atmospheres therefore seems to be an ideal test case for our purposes.

Following Schinder & Bludman (1989), we choose the parameters $R_{\rm max}=10$ and n=2. The variation of the metric functions ${\rm e}^\Phi$ and $\Gamma$ with radius is depicted in Fig. 5a. All other model parameters like boundary conditions and initial conditions are the same as in the "weak-field''-case which is defined by ${\rm e}^{\Phi}=\Gamma\equiv 1$ (see above). The stationary-state Eddington factor as computed from the Boltzmann equation in our program is displayed as a function of radius in Fig. 5c. The deviation from the relativistically correct results of Schinder & Bludman (1989) is significant. On the other hand, our result is very close to the corresponding "weak-field'' case, which is not shown in Fig. 5. This indicates that general relativistic ray bending, which is not correctly taken into account in our calculation, has a determining influence on the Eddington factor, while the contraction of rods and time dilation - both included in our calculation of the Eddington factor - are of minor importance. Indeed, the relativistically correct characteristic curves deviate considerably from straight rays, as can be seen in Schinder & Bludman (1989, Fig. 9).

The analytical stationary-state solution for the radiation flux density, which is easily verified to read $R^2{\rm e}^{\Phi}H(R,\epsilon_\infty)=A(\epsilon_\infty)$, (with $A(\epsilon_\infty)\in{\rm I\!R}$), is reproduced with an accuracy of 10^-7. This, of course, was to be expected for the atmosphere with isoenergetic scattering, since in the stationary state, the solution H(r) is solely determined by the zeroth moment equation (Eq. (69)) with no reference to the (incorrect) Eddington factor.

More remarkably, the quality of our solution for the mean intensity J(r), which is governed in the stationary state by Eq. (70) and thus is directly influenced by the Eddington factor, appears to be rather good, too. As can be seen from Fig. 5b, the slope of the general relativistic result is reproduced correctly throughout the atmosphere. The quantitative agreement is quite satisfactory as well. In judging the accuracy of our results one has to keep in mind that the model computed here is a rather extreme case in two respects: Firstly, the deviations of the metric functions ${\rm e}^\Phi$ and $\Gamma$ from unity (see Fig. 5a) are larger than those typically encountered in supernova simulations when the neutron star is not going to collapse to a black hole. Secondly, as stated in the beginning, the effect of a curved spacetime on the radiation field was observed by Schinder & Bludman (1989) to be much larger for the scattering atmosphere than in a case dominated by absorption (and emission). We expect real situations to be somewhere in between these two extremes.

  
4.2 Radiative transfer in stationary background media

Figure 6: Stationary-state angular mean J (times r2) of the intensity as measured in the comoving frame of reference. The abscissa gives the radial position. Results are shown for gray scattering atmospheres that expand according to a linear (Model "LVL''; Panel a) and the quadratic (Model "QVL''; Panel b) velocity law. The thin vertical line marks the upper boundary of the atmospheres. Different line styles of the curves correspond to different values of the parameter $\beta _{\rm max}$ which gives the maximum expansion velocity in units of c reached at the surface (r=11) of the model atmospheres. The thin solid curve corresponds to the static case. Reference solutions (symbols) were taken from Mihalas (1980, Fig. 2).

Open with DEXTER

In the following section we consider radiative transfer problems in moving, yet stationary background media. This means that in addition to time-independent radial profiles of the opacity and emissivity, a time-independent velocity field as a function of radius is prescribed. Stationary-state solutions for the radiation field are expected to exist is such cases and have been computed to high accuracy for some test problems including differentially moving atmospheres with relativistic velocities. By comparison with fully (special) relativistic calculations, we are not only able to test our implementation of the velocity dependent terms but can also judge the quality of the employed ${\cal O}(v/c)$-approximation of special relativistic effects in the equations of radiative transfer.

Upon introducing a non-vanishing velocity field a particular frame has to be specified, where physical quantities are measured. In all cases considered in this work, the latter is chosen to be the Lagrangian or comoving frame of reference. This has to be distinguished from the (numerical) coordinate grid that is used to simply label the events in spacetime. Although the transformation between different coordinate grids is trivial from an analytical point of view (e.g., the simple replacement $\partial/\partial t+v\partial/\partial r\to {\rm D}/{\rm D}t$ transforms from Eulerian to Lagrangian coordinates), the numerical treatment can be involved (cf. Sect. 3.4, where our algorithm for computing a formal solution of the radiative transfer equation in Eulerian coordinates is described). Therefore we present test results obtained with two different versions of our radiative-transfer code, one which uses Lagrangian coordinates and another one which employs an Eulerian coordinate grid.

Figure 7: Frequency integrated zeroth order angular moment of the comoving frame intensity ( $\int_0^\infty{\rm d}\xi~J(\xi)$) as a function of radius for our ${\cal O}(v/c)$ stationary-state solutions. The thin vertical line marks the outer boundary of the atmospheres. The small arrows point to the positions where an optical depth (which is calculated for $\xi <\xi _0$) of unity is reached. Reference solutions (symbols) obtained with a fully relativistic code were taken from Mihalas (1980, Figs. 3 and 5), where the read-off error is approximately given by the size of the symbols. Different maximum velocities $\beta _{\rm max}$ of the atmospheres are indicated by different line styles and symbols. In Panel  a) the results are depicted for Model (1) with the opacity step located at $\xi _0=10$ ( $\Delta =0.25$). Panel  b) shows results for Model (2) with $\xi _0=3$ ( $\Delta =0.2$).

Open with DEXTER

Differentially expanding gray scattering atmospheres:

Mihalas (1980) presents stationary-state solutions for the same type of purely scattering atmospheres considered above. In addition to the static case, he investigated differentially expanding atmospheres with relativistic velocities.

We refer to this paper and study two classes of atmospheres as test problems, which differ in the functional dependence of the expansion velocity ( $v\equiv c~\beta$) on the radius:

$$\beta(r)=\beta_{\rm max}~ \left(\frac{r-r_{\rm min}}{r_{\r... ...{\rm min}}\right)^m~, \quad {\rm with} \quad m \in \{1,2\} ~.$$ (71)

The case m=1 describes a "linear velocity law'', the case m=2 is called "quadratic velocity law''. Different models are labeled by the parameter $0 \le \beta_{\rm max} < 1$, which is the maximum expansion velocity in units of c reached at the surface of the atmosphere ( $r=r_{\max}$).

The scattering opacity of the expanding atmospheres is the same as the one used for the static cases in Sect. 4.1 except for the unit of length, which, for the ease of comparison, is adopted from Mihalas (1980):

$$\chi(t,r,\epsilon)\equiv\kappa_{\rm s}(r)=\alpha~r^{-2},\quad r_{\rm min} \le r \le r_{\rm max} ~,$$ (72)

with the parameters $r_{\rm min}=1$, $r_{\rm max}=11$ and $\alpha=10.989$. This yields the same optical depth as the atmosphere with $r_{\rm max}=1$ (and $r_{\rm min}=1/11$) in the units of Hummer & Rybicki (1971), who set $\alpha=1$ (cf. Eq. (66)). For the profile of Eq. (72), an optical depth of unity is reached at a radius of $r\approx 5.5$. Since the opacity does not depend on the frequency of the radiation, the atmospheres are referred to as "gray''. In this case, the transport solution does not depend on $\epsilon$ and the $\partial/\partial\epsilon$-terms in the comoving-frame Boltzmann equation can be dropped.

For our simulations we used a numerical grid with 200 radial points, which were initially distributed logarithmically between $r_{\rm min}$ and $r_{\rm max}$. Following Mihalas (1980) as closely as possible, we impose the boundary conditions

$${\cal I}(t,r_{\rm min},\mu) 11+3\mu,\quad 0 \le \mu \le 1,$$ =

$${\cal I}(t,R(t),\mu) 0, \qquad\quad -1 \le \mu \le 0 ,$$ = (73)

with $R(t)=r_{\rm max}+\beta_{\rm max}~ct$. The numerical values in the diffusion ansatz at the inner boundary are chosen according to the asymptotic solution of Hummer & Rybicki (1971, Eq. (2.22)). We started the simulations using the stationary-state solution of the corresponding static atmosphere (see Sect. 4.1) as an initial condition (at t=0). The results presented here were obtained by using a Lagrangian coordinate grid. Our calculations were terminated when the radiation field (as measured in the comoving frame of reference) showed no more substantial temporal variations at any given radial position. We refer to our results as "stationary-state solutions'' in this sense. We also tested the Eulerian grid version of the code with a less extensive set of models. The results were practically identical.

We compare our stationary-state solutions with results obtained by Mihalas (1980). In the latter investigation a numerical method that is accurate to all orders in v/c was employed to solve the stationary radiative transfer equation. We therefore do not only test the correct numerical implementation of the ${\cal O}(v/c)$ equations (Eqs. (6)-(8)), but we can also get a feeling for the quality of the ${\cal O}(v/c)$-approximation to the special relativistic transfer equation. All characteristic features of the solution, discussed in detail by Mihalas (1980), are reproduced correctly by our implementation of the ${\cal O}(v/c)$ equations. It is remarkable how accurately the fully relativistic solutions are reproduced (Fig. 6). The differences are $\lesssim$10%, even for $\beta_{\rm max}=0.5$.

Figure 8: Spectral distributions of the angular moment J of the comoving frame intensity for our stationary-state solutions of the ${\cal O}(v/c)$ radiative transfer equation. Spectra are displayed for a radial position corresponding to an optical depth (which is calculated for $\xi <\xi _0$) of unity and for the surface of the atmosphere. For comparison, the equilibrium intensity (Planck function) at an optical depth of unity is plotted as a thin line. Different maximum velocities ( $\beta _{\rm max}$) of the expanding atmospheres are indicated by different line styles. Panel  a) shows results for Model (1), where the opacity step (vertical arrows in the plots) is located at $\xi _0=10$ ( $\Delta =0.25$). Panel  b) shows Model (2) with $\xi _0=3$ ( $\Delta =0.2$). Our results may be compared with Mihalas (1980, Figs. 4 and 6).

Open with DEXTER

Differentially expanding, nongray, isothermal atmospheres:

Mihalas (1980) computed stationary-state solutions also for relativistically expanding nongray atmospheres, i.e., the opacity and emissivity depended explicitly on the radiation frequency. Therefore the full frequency dependence of the comoving frame Boltzmann equation had to be retained. The atmospheres were assumed to be stationary, to emit a thermal continuum of radiation, and to be isothermal, i.e., no radial or temporal variations of the source function were present. The radiation-matter interaction was solely via absorption and emission.

Introducing the dimensionless "frequency'' $\xi=\epsilon/T$ (the temperature T is measured in units of the Boltzmann constant ${k_{\rm B}~}$), the emissivity reads

$$\eta(\xi)=\kappa_{\rm a}(\xi)~b(\xi)= \kappa_{\rm a}(\xi)~\frac{\xi^3}{\exp(\xi)-1} ~\cdot$$ (74)

The absorption opacity $\kappa_{\rm a}$ was given a step-like functional dependence on frequency, characterized by the frequency $\xi_0$ where the step is located and some characteristic width $\Delta$ of the function $s(\xi):=\exp[-(\xi-\xi_0)^2/\Delta^2]$, which mediates a continuous variation of $\kappa_{\rm a}$ across the step:

$$\kappa_{\rm a}(r,\xi)= \left\{ \begin{array}{ll} \chi_1(r... ...\ \chi_1(r),& ~~{\rm for}~~ \xi > \xi_0, \end{array}\right.$$ (75)

with $\chi_1(r):=10\alpha~r^{-2}$, $\chi_2(r):= \alpha~r^{-2}$, and the value $\alpha:=10.9989$. Thus, the optical depth at a given radius is roughly 10 times smaller for "low-frequency'' radiation ( $\xi <\xi _0$) than for radiation with "high'' frequency ( $\xi > \xi_0$).

Following Mihalas (1980), we consider the two examples (1) $\xi _0=10$, $\Delta =0.25$, and (2) $\xi _0=3$, $\Delta =0.2$. In the former case the frequency $\xi_0$, where the opacity jump is located, is considerably larger than the frequency of the maximum of the Planck function ( $\xi_{\rm max}\approx 2.82$), whereas it is close to the maximum in the latter case.

The radial extent of the atmospheres was chosen to be $0 \le r \le 11$ in both cases. For all models the linear law (m=1, cf. Eq. (71)) for the expansion velocity as a function of radius was used. The results presented below were obtained using an Eulerian coordinate grid. Our code version with a Lagrangian grid yields very similar results.

We have used 200 logarithmically distributed radial zones for the radial range $1 \le r \le 11$, supplemented by a few additional zones between $0 \le r < 1$. Boundary conditions were chosen as usual: Spherical symmetry at the coordinate center requires ${\cal I}(t,r=0,\mu=1)={\cal I}(t,r=0,\mu=-1)$, and demanding that no radiation is entering through the surface of the atmosphere translates to the condition ${\cal I}(t,r=11,\mu \le 0)=0$.

The energy grid consisted of 21 bins (very similar results were obtained using 8 and 12 bins instead) of equal width covering the range $0\le \xi \le 12$ for Model (1) and $0\le \xi \le 8$ for Model (2), respectively. Details about the stationary-state solutions and their physical interpretation can be found in Mihalas (1980, Sect. 5). In our test calculations we were able to reproduce all qualitative trends discussed there. The quantitative agreement of the zeroth order angular moment J of the specific intensity is within a few percent (see Fig. 7), even for the case of $\beta_{\rm max}=0.5$. This, as previously stated, does not only demonstrate the accuracy of our numerical implementation but even more importantly also justifies the physical approximations employed in the underlying finite difference equations. Note that in applications of our code to supernova calculations, one expects velocities that are not in excess of $v\approx 0.3c$.

A rigorous comparison of the spectral distribution of the angular moment J as calculated by our radiative transfer code (cf. Fig. 8) with the results of Mihalas (1980) is difficult, since he shows only results for the static case and the rather extreme case $\beta_{\rm max}=0.8$. The latter atmosphere obviously cannot be modeled reasonably well with our ${\cal O}(\beta)$-code. Nevertheless, by inspection of Figs. 3 and 4 of Mihalas (1980) one can infer that the shapes of our spectra as displayed in Fig. 8 exhibit the correct qualitative dependence on the expansion velocity.

Frame dependence?

As an important consistency check we verified that, to order ${\cal O}(v/c)$, the numerical solution of our implementation of the comoving frame equations shows the correct transformation properties of the stress-energy tensor of radiation (see Eq. (10)) and no unphysical artifacts are caused by the choice of a moving reference frame (cf. the critique of Lowrie et al. 1999).

For this purpose results are produced for neutrino transport in a thermally and hydrodynamically frozen post-bounce model of a supernova calculation by Bruenn (1993). This model can be characterized as follows: Neutrinos are emitted from a hot and dense hydrostatic inner core with a radius of $r\approx 110$ km. This inner core is surrounded by a supersonically infalling outer core, which is effectively transparent to the radiation. Velocities in this outer atmosphere range from $v(r=140~{\rm km})\approx -0.12 c$ to $v(r=250~{\rm km})\approx -0.08 c$.

The stress-energy tensor as derived from the stationary-state solution of the comoving frame transport equation (Eq. (6)) is transformed to the Eulerian frame (i.e. the inertial frame in which the center of the star is at rest) and can then be compared with the result of an independent calculation which solves the transport equation directly in the Eulerian frame.

	Figure 9: Comparison of the stationary-state solutions of the comoving frame transfer equation (bold lines) with the solutions after transformation to the Eulerian frame (thin lines) and a reference calculation for a static stellar background ($v\equiv 0$; dashed lines). Panel a) shows the total energy density of neutrinos (times r2 c), Panel  b) displays the total energy flux density (times r2).
Open with DEXTER

Figure 9 shows the components $E=4\pi/c\int_0^\infty {\rm d}\epsilon~ J$ and $F=4\pi\int_0^\infty {\rm d}\epsilon~ H$ of the stress-energy tensor (scaled by r² c and r², respectively). Due to the large inwardly directed velocities in the outer atmosphere, the neutrinos emitted by the inner core are blueshifted for observers which are locally comoving with the fluid flow (see the bold lines in Fig. 9). When the stress-energy tensor is transformed to the Eulerian frame (Eq. (10)), this effect obviously disappears (thin lines in Fig. 9). Comparison with the stress-energy tensor obtained by solving the transport equation directly in the Eulerian frame (dashed lines in Fig. 9) reveals good agreement for both the energy density E and the energy flux density F in the rapidly moving atmosphere. This demonstrates that no unphysical frame dependence is present in the calculations. The fluctuations of the transformed energy flux density near the center are caused by the term $\beta (J + K )$ in the expression for $H^{\rm Eul}$ in Eq. (10). Because J (and K) are large compared with H in the dense interior, even small velocities lead to a significant "advective contribution'' to $H^{\rm Eul}$.

Note that an Eulerian-frame calculation in general requires complicated velocity-dependent transformations of the source terms on the rhs of the transport equation (see Mihalas & Mihalas 1984). This, however, is unnecessary in the specific situation of the discussed model: Because the region where most of the neutrinos are produced moves with low velocities, the source terms there can be evaluated in the rest frame of the fluid. In the outer atmosphere on the other hand, where the large velocities would require a careful transformation, the source terms nearly vanish. This allowed us to simply drop the velocity-dependent terms on the lhs of Eq. (6) in order to perform a calculation with results that are in agreement with the Eulerian frame solution in the low-density part of the star.

  
4.3 Core collapse and supernova simulations

   
4.3.1 Newtonian gravity

Our new neutrino hydrodynamics code has already been applied to dynamical supernova simulations (Rampp & Janka 2000). Here we report some results of a Newtonian calculation of the collapse phase of a stellar iron core with mass $M_{\rm Fe}=1.28~{M_{\odot}}$ (plus the innermost $0.1~{M_{\odot}}$ of the silicon shell). It is the core of a $15~{M_{\odot}}$ progenitor star (Model "s15s7b2'', Woosley 1999; Heger 2000). We compare the results to a calculation published by Bruenn & Mezzacappa (1997). The input physics, in particular the stellar model and the high-density equation of state, were adopted from model "B'' of Bruenn & Mezzacappa (1997), with the exception that we do not include neutrinos other than ${\nu_{\rm e}}$. This, however, does not make a difference during the collapse phase where the electron degeneracy is so high that only ${\nu_{\rm e}}$ can be produced in significant numbers.

	Figure 10: Profiles of the total lepton fraction $Y_{\rm lep}$ (dash-dotted line) and the electron neutrino fraction $Y_{\nu _{\rm e}}$ (solid line) as a function of the enclosed mass at the time when the central density reaches 1014  ${{\rm g}~{\rm cm}^{-3}}$. For comparison the results of Bruenn & Mezzacappa (1997, Fig. 5a) are plotted with diamonds.
Open with DEXTER

The hydrodynamics was solved on a grid with 400 radial zones out to 20 000 km, which were moved with the infalling matter during core collapse. For the transport we used an Eulerian grid with 213 geometrically spaced radial zones, 233 tangent rays and 21 energy bins geometrically distributed between 0 and 380 MeV the center of the first zone being at 2 MeV.

	Figure 11: Central lepton fraction $Y_{\rm lep}:={Y_{\rm e}}+Y_{\nu _{\rm e}}$ (Panel a)) and central entropy per baryon ( b)) as functions of the central density for our Newtonian core-collapse simulation of a $15~{M_{\odot}}$ star (solid lines), compared with the results (diamonds) obtained by Bruenn & Mezzacappa (1997, Figs. 4a, 8a).
Open with DEXTER

The total lepton fraction $Y_{\rm lep}:={Y_{\rm e}}+Y_{\nu _{\rm e}}$and the entropy per baryon at the center of the core as a function of the central density are displayed in Fig. 11. The agreement with results of Bruenn & Mezzacappa (1997) is excellent. Also the total lepton fraction as a function of enclosed mass matches perfectly (Fig. 10). Defining the shock formation point according to Bruenn & Mezzacappa (1997) as the mass coordinate where the entropy per baryon after core bounce first reaches a value of $s=3~{k_{\rm B}~}$ we find $M_{\rm S}\approx 0.62~M_\odot$. This is about $0.03~M_\odot$ or 5% larger than the value given by Bruenn & Mezzacappa (1997, Table V).

Judging our results, one should, of course, take into account that Bruenn & Mezzacappa (1997) employed a multi-group flux-limited diffusion (MGFLD) approximation for treating the neutrino transport, whereas our code solves the Boltzmann equation. Mezzacappa & Bruenn (1993b,c) performed a detailed comparison of core-collapse simulations with MGFLD and with their Boltzmann solver based on the $S_{\rm N}$-method. For the quantities presented here they also found good agreement of their results throughout the inner core.

In the outer regions of the collapsing core, the situation is somewhat different: Mezzacappa & Bruenn (1993b,c) demonstrated that a number of quantities reveal significant deviations between MGFLD and the Boltzmann transport ( $S_{\rm N}$-method). In particular, they obtained an electron neutrino fraction which was by up to 20% larger in the Boltzmann run. For the total lepton fraction, however, the difference was only 1%. Our results (Fig. 10) confirm the agreement of Boltzmann and MGFLD results for the latter quantity. A more detailed comparison between the variable Eddington factor method and the $S_{\rm N}$-method, however, is desirable and currently underway.

Figure 12: Electron neutrino fraction $Y_{\nu _{\rm e}}$ (bold solid line, scale given by the axis on the left side of the plots) and electron neutrino equilibrium chemical potential $\mu _{\nu _{\rm e}}^{\rm eq}$ ("Fermi surface'', dash-dotted line, scale given by the axis on the right side) as functions of enclosed mass at the time when the central density has reached 1014  ${{\rm g}~{\rm cm}^{-3}}$ (The plots show results of collapse simulations for a $18~M_\odot$ blue supergiant progenitor star; (Woosley et al. 1997). The spectral resolution of the energy grid is indicated by the spacing of the short horizontal lines (right ordinate), which correspond to the boundaries of the energy grid. Panel a) displays results obtained with "low'' spectral resolution (21 energy bins distributed geometrically between 0 and $380~{{\rm MeV}}$), Panel b) shows the run with "high'' spectral resolution (81 energy bins covering the same spectral range). For better comparison $Y_{\nu _{\rm e}}$ of the high-resolution run is drawn as a thin solid line also in Panel a).

Open with DEXTER

Spectral resolution:

A few quantities exhibit spurious oscillations as a function of the radial (mass-) coordinate and time in our core-collapse simulations. Figure 12a shows such features for the electron neutrino fraction in the mass shells $0 \le M \lesssim 0.6~{M_{\odot}}$ when the central density is $\rho_{\rm c} \simeq 10^{14}$  ${{\rm g}~{\rm cm}^{-3}}$. Similar structures are visible in the profile of the entropy per baryon and in the neutrino luminosity as a function of radius (or mass) once the central density becomes larger than $\rho_{\rm c} \simeq 10^{13}$  ${{\rm g}~{\rm cm}^{-3}}$. They can also be seen in the central entropy as a function of the central density or time (cf. Fig. 11b). Mezzacappa & Bruenn (1993c) found a similar effect in their core-collapse calculations. They identified the finite spectral resolution of the degenerate Fermi distribution of the electron neutrinos to be the origin of these oscillations (Mezzacappa & Bruenn 1993c, Sect. 3.3).

For comparison, we have computed a collapse model with a significantly improved spectral resolution: instead of 21 bins we used 81 energy zones to describe the energy spectrum between 0 and 380 MeV. A result of this simulation is displayed in Fig. 12b. The neutrino fraction (and related quantities) are now perfectly smooth functions of the mass coordinate.

A reasonable compromise between accuracy and computational load could be achieved by using approximately 20-30 (geometrically spaced) energy bins. This seems to be sufficient to adequately treat the sharp Fermi surface of the degenerate distribution of electron neutrinos that builds up during collapse. Even for less energy bins fluctuating quantities follow the correct evolutionary trends (Fig. 12a). Therefore we conclude with the remark that the presence of oscillations in some transport quantities may be considered as undesirable. The overall evolution of the collapsing core, its deleptonization or the shock formation radius, however, were found not to be sensitive to such "noise''.

Also during the post-bounce evolution (for which the transport of ${\bar\nu_{\rm e}}$ was included), we determined no significant differences between the dynamical evolution of a model with 17 energy bins and a model computed with 27 energy bins, using an otherwise identical setup: At a time of 100 ms after bounce, the structure of the models, measured by the radial positions of fluid elements with the same mass coordinate, agreed within a relative error of <1%. This is corroborated by the fact that for given snapshots from the dynamical evolution we found the neutrino source terms $Q_{\rm E}$ and also $Q_{\rm N}$ to be practically identical in both calculations.

Radial resolution and size of the numerical time step:

Varying the number of radial zones for the hydrodynamic and transport parts of the code between 100 and 200 and using a comoving hydrodynamic grid during the phase of stellar core collapse, we found no significant changes of the results. The same holds true for the post-bounce evolution where we switch to a fixed, Eulerian grid. The only exception to this insensitivity is the central density at bounce, which shows variations of up to $10\%$, depending on the size of the innermost zones of the hydrodynamic grid.

The choice of the radial grid, however, also influences the size of the numerical time step. According to our experience the latter has to be chosen very carefully. We give two examples here:

(a)

While neutrinos escape unhindered initially, the collapsing stellar core becomes opaque when the density exceeds $\simeq$10¹² g cm^-3, and "neutrino trapping'' sets in. Since neutrinos are then carried along with the stellar plasma, the subsequent, accelerating phase of the infall is adiabatic (the sum of medium and neutrino entropy per baryon is constant: $s+s_\nu={\rm const.}$) and the total lepton fraction of a fluid element is conserved. Fulfilling the corresponding conditions, ${\rm D}Y_{\rm lep}/{\rm D}t=0$ and ${\rm D}(s+s_\nu)/{\rm D}t=0$, numerically on an Eulerian grid requires a sufficiently accurate treatment of the advection of energy and leptons. We found that treating the advection terms in the neutrino transport at first order accuracy in time (setting $\zeta=1$in Eq. (25)) is not sufficient and causes deviations from $Y_{\rm lep}={\rm const.}$ and $s+s_\nu={\rm const.}$at a level of $\gtrsim$10%. This problem can be diminished when the transport time step is reduced significantly below the limits mentioned in Sect. 3.6. The time step reduction and the corresponding increase in computer time, however, can be avoided by handling the advection "second-order'' accurate in time, chosing $\zeta\approx 0.5$ in Eq. (25). As visible from Fig. 11, this yields very good results (note that the entropy of neutrinos, $s_\nu$, is small compared to the plasma entropy s, and essentially constant after neutrino trapping).

(b)

Another problem occurred during the long-time evolution of the forming neutron star after core bounce. For too large a value of the hydrodynamic time step, we discovered an artificial inflation of the mass shells that are located in the region of a steep density gradient below the neutrinosphere. This effect was found not to be caused by insufficient energy resolution of the neutrino transport nor does it directly depend on the employed number of radial zones. Varying the latter by a factor of two and testing 17 and 25 energy bins in the interval between 0 and $380~{{\rm MeV}}$, we found that the agreement of our transport results was better than $1\%$ and thus not worse than the uncertainties due to the interpolation of the stellar background. Since the post-bounce evolution is usually calculated with the same grids for the hydrodynamics and the transport, possible errors associated with the mapping of quantities between the grids can also be excluded as the source of the problem. Moreover, even if different radial grids are used as in the model published by Rampp & Janka (2000), an expansion of the outer layers of the forming neutron star is not necessarily a consequence. The quality of our calculations, however, turned out to be dependent on the size of the time step of the hydrodynamics code. Our tests showed that a reduction of the time step was necessary to prevent the artificial expansion of the quasi-hydrostatic structure. We were so far unable to exactly locate the origin of this numerical problem, but presume that it might be connected with the accuracy of the implementation of the neutrino momentum term in the hydrodynamics code. Good results, however, were obtained, when the length of the time step was guided by the typical time scale set by the neutrino momentum transfer to the stellar medium, $\Delta t_{\rm Hyd}\simeq \rho~v/Q_{\rm M}$ (cf. Eq. (2)).

4.3.2 Approximate relativistic core collapse

We also tried to assess the quality of our approximate treatment of general relativity (GR) in the equations of hydrodynamics and neutrino transport (cf. Sect. 3.7). For this purpose we performed core-collapse simulations for the $15~{M_{\odot}}$ progenitor star "s15s7b2'' and compared characteristic quantities with the information given for a MGFLD simulation of the same star in a recent paper by Bruenn et al. (2001), and for a Boltzmann simulation by Liebendörfer et al. (2001a).

For the hydrodynamics we used 400 radial zones out to 10 000 km, which were moved with the infalling matter during core collapse. The transport was solved on an Eulerian grid with 200 geometrically spaced radial zones, 220 tangent rays and 17 energy bins geometrically distributed between 0 and 380 MeV the center of the first zone being at 2 MeV.

At bounce, which occurs at $t_{\rm b}^{\rm GR}=177.7~{{\rm ms}}$ ( $t_{\rm b}^{\rm Newt}=196.1~{{\rm ms}}$) after the start of the calculation, the central density reaches a maximum of $\rho^{\rm GR}_{\rm c}=4.8\times 10^{14}$  ${{\rm g}~{\rm cm}^{-3}}$( $\rho^{\rm Newt}_{\rm c}=3.35\times 10^{14}$  ${{\rm g}~{\rm cm}^{-3}}$) in our GR model (for comparison, the corresponding values for the Newtonian runs are given in brackets). The hydrodynamic shock forms at a mass coordinate of $M_{\rm S}^{\rm GR}=0.51~{M_{\odot}}$ ( $M_{\rm S}^{\rm Newt}=0.62~{M_{\odot}}$). For the same stellar model, Bruenn et al. (2001) found a bounce density of $\rho^{\rm GR}_{\rm c}=4.23\times 10^{14}$  ${{\rm g}~{\rm cm}^{-3}}$for the general relativistic simulation ( $\rho^{\rm Newt}_{\rm c}=3.20\times 10^{14}$  ${{\rm g}~{\rm cm}^{-3}}$). The corresponding time of the bounce was $t_{\rm b}^{\rm GR}=182.9~{{\rm ms}}$ ( $t_{\rm b}^{\rm Newt}=201.3~{{\rm ms}}$). Liebendörfer et al. (2001a) located the shock formation at a mass coordinate of $M_{\rm S}^{\rm GR}=0.53~{M_{\odot}}$ ( $M_{\rm S}^{\rm Newt}=0.65~{M_{\odot}}$).

Since the collapse time is determined by the slow initial phase of the contraction, which is sensitive to the initial conditions and the details of the numerical treatment at the start of the calculation (Liebendörfer, personal communication), it may be better to use the difference between the GR and the Newtonian simulations, instead of comparing absolute times at bounce. Our value of $t_{\rm b}^{\rm Newt}-t_{\rm b}^{\rm GR}=18.4~{{\rm ms}}$ is very close to the result of Bruenn et al. (2001): $t_{\rm b}^{\rm Newt}-t_{\rm b}^{\rm GR}=21.6~{{\rm ms}}$. A part of the remaining difference may be attributed to the fact that Bruenn et al. (2001) employed a MGFLD approximation of the neutrino transport with some smaller differences also in the input physics (e.g., Bruenn et al. 2001 as well as Liebendörfer et al. 2001a did not take into account ion-ion correlations for the coherent scattering, but we did).

At the level of comparison which is possible on grounds of published numbers we conclude that the agreement between our approximate treatment of GR and the relativistic core collapse seems to be reasonably good.

  
5 Summary and conclusions

We have presented a detailed description of the numerical implementation of our new neutrino transport code and its coupling to a hydrodynamics program which allows simulations in spherical symmetry as well as in two or three dimensions. The transport code solves the energy and time dependent Boltzmann equation by a variable Eddington factor technique. To this end a model Boltzmann equation is discretized along tangent rays and integrated for determining the variable Eddington factors which provide the closure relations for the moment equations of neutrino energy and momentum. The latter yield updated values of the neutrino energy density and neutrino flux, which are fed back into the Boltzmann equation to handle the collision integral on its right hand side.

The system of Boltzmann equation and moment equations together with the operator-splitted terms of lepton number and energy exchange with the stellar background (which influence the evolution of the thermodynamical quantities and the composition of the stellar medium and thus the neutrino-matter interactions) are iterated to convergence. Lepton number conservation is ensured by solving additional moment equations for neutrino number density and number flux.

The integration is performed with implicit time stepping, which avoids rigorous limitations by the CFL condition and ensures that equilibrium between neutrinos and stellar medium can be established accurately and without oscillations despite of stiff neutrino absorption and emission terms. The coupling of energy bins due to Doppler shifts and energy changing scattering processes is implemented in a fully implicit way. When coupling the transport part to an explicit hydrodynamics program, as done in this work, computational efficiency requires to have the option of using different time step lengths for transport and hydrodynamics. In addition, we found it advantageous to have implemented the possibility of choosing different radial grids for both components of the program, and to switch between Lagrangian and Eulerian coordinates dependent on the considered physical problem (although we always measure physical quantities of the transport in the comoving frame of reference).

We have suggested an approximation for applying the code to multi-dimensional supernova simulations which can account for the fact that hydrodynamically unstable layers develop in the collapsed stellar core. The variable Eddington factor method offers the advantage that this can be done numerically rather efficiently (concerning programming effort as well as computational load). Although the neutrino moment equations for the radial transport are solved independently in all angular zones of the spherical grid of the multi-dimensional hydrodynamical simulation, the set of moment equations is closed by a variable Eddington factor which is computed only once for an angularly averaged stellar background.

This approximation saves a significant amount of computer time compared to truely multi-dimensional transport, and yields a code structure which can easily be implemented on parallel computers. The treatment is constrained to radial transport and thus neglects lateral propagation of neutrinos. Its applicability is therefore limited to physical problems where no global asphericities occur. It must also be considered just as a first, approximate step towards a really multi-dimensional transport in convective layers inside the neutrinosphere. The approach, however, should be perfectly suitable to treat anisotropies associated with convective overturn in the neutrino-heated region between the neutron star and the supernova shock. In the latter region neutrinos are only loosely coupled to the stellar medium (the optical depth is typically only 0.05-0.2). Therefore the energy and lepton number exchange between neutrinos and the stellar medium depends on the presence of anisotropies and inhomogeneities, whereas the transport of neutrinos is essentially unaffected by non-radial motions of the stellar plasma.

In preliminary multi-dimensional simulations of transport in convecting neutron stars we have convinced ourselves that this method guarantees good numerical convergence, because the variable Eddington factor as a normalized quantity depends only weakly on lateral variations in the stellar medium. The method ensures global energy conservation and enables local thermodynamical equilibrium between stellar medium and neutrinos to be established when the optical depth is sufficiently high. We therefore think that the proposed approximation is practicable for neutrino transport in multi-dimensional supernova simulations before fully multi-dimensional transport becomes feasible. The latter is certainly a major challenge for the years to come.

Our transport code exists in a ${\cal O}(v/c)$ version for nonrelativistically moving media, and in a general relativistic version. We have not yet coupled it to a hydrodynamics program in full general relativity. Instead, we performed calculations of approximate relativistic core collapse, where corrections to the Newtonian gravitational potential were included and only the redshift and time dilation effects were retained in the transport equations (this means that space-time is considered to be flat). The results for stellar core collapse compared with published fully relativistic calculations (with Boltzmann neutrino transport by Liebendörfer et al. 2001 and with multi-group flux-limited diffusion by Bruenn et al. 2001) are encouraging and suggest that this approximation works reasonably well and accounts for the most important effects of relativistic gravity as long as one does not get very close to the limit of black hole formation.

We have presented results for a variety of idealized, partly analytically solvable test calculations (in spherical symmetry) which demonstrate the numerical efficiency and accuracy of our neutrino transport code. The "neutrino radiation hydrodynamics'' was verified by core-collapse and post-bounce simulations for cases where published results were available and could reasonably well serve for a comparison (a direct and more detailed comparion with Boltzmann calculations using the $S_{\rm N}$ method by M. Liebendörfer and A. Mezzacappa is in progress).

Tests for a number of static model atmospheres showed that the treatment of the angular dependence of the neutrino transport and phase space distribution can be handled accurately by the variable Eddington factor method. This holds for moderately strong general relativistic problems, too, even if ray bending effects are neglected in the determination of the variable Eddington factor. The numerical quality of the handling of the energy dependence of the transport was also checked by considering background models with stationary fluid flow. We included a model with mildly relativistic motion (up to v/c = 0.5) and found that the code produces accurate and physically consistent results although we employ an ${\cal O}(v/c)$ approximation to the special relativistic, comoving frame radiative transfer equation. Also the omission of some velocity-dependent terms in the model Boltzmann equation for determining the variable Eddington factor turned out not to be harmful in this respect. Closing the radiation moment equations by a variable Eddington factor seems to be remarkably robust against approximations in the treatment of the Boltzmann equation.

The tests have also demonstrated that our code performs very efficiently. Only a few iterations (typically between 3 and 7) are needed for obtaining a converged solution of the system of Boltzmann equation and moment equations even in cases where scattering rates dominate neutrino absorption. The calculation of the formal solution of the Boltzmann equation, from which the variable Eddington factor is derived, turned out to consume only 10-20% of the computer time. The major computational load comes from the inversion of the set of moment equations. The latter have a reduced dimensionality relative to the Boltzmann equation, because the dependence on the angular direction of the radiation propagation was removed by carrying out the angular integration. This advantage, however, has to be payed for by the iteration procedure with the Boltzmann equation.

Using Boltzmann solvers for the neutrino transport means a new level of sophistication in hydrodynamical simulations of supernova core collapse and post-bounce evolution. It permits a technical handling of the transport which for the first time is more accurate than the standard treatment of neutrino-matter interactions. The latter includes a number of approximations and simplifications which should be replaced by a better description, for example the detailed reaction kinematics of neutrino-nucleon interactions, phase space blocking of nucleons, effects due to weak magnetism in neutrino-nucleon reactions, or nucleon correlations in the dense medium of the neutron star (Reddy et al. 1998; Reddy et al. 1999; Burrows & Sawyer 1999; Horowitz 2002).

Boltzmann calculations will not only remove imponderabilities associated with the use of approximate methods for the neutrino transport in hydrodynamical supernova models. They will also allow for much more reliable predictions of the temporal, spectral and flavour properties of the neutrino signal from supernovae and newly formed neutron stars. This is an indispensable requirement for the interpretation of a future measurement of neutrinos from a Galactic supernova.

Acknowledgements

We thank our referee, A. Mezzacappa, for insightful comments which helped us improving the original manuscript. We are grateful to K. Kifonidis, T. Plewa and E. Müller for the latest version of the PROMETHEUS code, to K. Kifonidis for contributing a matrix solver for the transport on parallel computer platforms, and to R. Buras for implementing the three-flavour version of the code and coining subroutines to calculate the neutrino pair and bremsstrahlung rates. We thank S. Bruenn, R. Fischer, W. Keil, M. Liebendörfer and G. Raffelt for helpful conversations, and M. Liebendörfer for providing us with data of his simulations. The Institute for Nuclear Theory at the University of Washington is gratefully acknowledged for its hospitality and the Department of Energy for support during a visit of the Summer Program on Neutron Stars. This work was also supported by the Sonderforschungsbereich 375 on "Astroparticle Physics'' of the Deutsche Forschungsgemeinschaft. The computations were performed on the NEC SX-5/3C of the Rechenzentrum Garching and on the CRAY T90 of the John von Neumann Institute for Computing (NIC) in Jülich.

References

 
• Arnett, W. D. 1977, ApJ, 218, 815 In the text NASA ADS  
• Baron, E., Bethe, H. A., Brown, G. E., Cooperstein, J., & Kahana, S. 1987, Phys. Rev. Lett., 59, 736 In the text NASA ADS  
• Baron, E., & Cooperstein, J. 1990, ApJ, 353, 597 In the text NASA ADS  
• Baron, E., Cooperstein, J., & Kahana, S. 1985, Phys. Rev. Lett. 55, 126 In the text NASA ADS  
• Baron, E., Myra, E. S., Cooperstein, J., & van den Horn, L. J. 1989, ApJ, 339, 978 In the text NASA ADS  
• Baschek, B., Efimov, G. V., Von Waldenfels, W., & Wehrse, R. 1997, A&A, 317, 630 In the text NASA ADS  
• Bethe, H. A., & Wilson, J. R. 1985, ApJ, 295, 14 In the text NASA ADS  
• Blinnikov, S. I., Eastman, R., Bartunov, O. S., Popolitov, V. A., & Woosley, S. E. 1998, ApJ, 496, 454 In the text NASA ADS  
• Bowers, R. L., & Wilson, J. R. 1982, ApJS, 50, 115 In the text NASA ADS  
• Bronstein, I., & Semendjajew, K. 1991, Taschenbuch der Mathematik (Thun und Frankfurt/Main: Verlag Harri Deutsch) In the text  
• Bruenn, S. W. 1985, ApJS, 58, 771 In the text NASA ADS  
• Bruenn, S. 1989a, ApJ, 340, 955 In the text NASA ADS  
• Bruenn, S. 1989b, ApJ, 341, 385 NASA ADS  
• Bruenn, S. W. 1993, in Nuclear Physics in the Universe, ed. M. Guidry & M. Strayer, First Symposium on Nuclear Physics in the Universe (Oak Ridge, Tennessee, USA: Institute of Physics Publishing, Bristol and Philadelphia), 31 In the text  
• Bruenn, S. W., De Nisco, K. R., & Mezzacappa, A. 2001, ApJ, 560, 326 In the text NASA ADS  
• Bruenn, S. W., & Dineva, T. 1996, ApJ, 458, L71 In the text NASA ADS  
• Bruenn, S. W., & Mezzacappa, A. 1997, Phys. Rev. D, 56, 7529 In the text NASA ADS  
• Bruenn, S. W., Mezzacappa, A., & Dineva, T. 1995, Phys. Rep., 256, 69 In the text NASA ADS  
• Burrows, A. 1987, ApJ, 318, L57 In the text NASA ADS  
• Burrows, A., Hayes, J., & Fryxell, B. A. 1995, ApJ, 450, 830 In the text NASA ADS  
• Burrows, A., & Lattimer, J. 1988, Phys. Rep., 163, 51 In the text  
• Burrows, A., & Sawyer, R. F. 1999, Phys. Rev. C, 59, 510 In the text NASA ADS  
• Burrows, A., Young, T., Pinto, P., Eastman, R., & Thompson, T. 2000, ApJ, 539, 865 In the text NASA ADS  
• Carlson, B. 1967, Computational Methods in Physics, 1 In the text  
• Castor, J. 1972, ApJ, 178, 779 In the text NASA ADS  
• Cernohorsky, J. 1994, ApJ, 433, 247 In the text NASA ADS  
• Colella, P., & Woodward, P. 1984, J. Comp. Phys., 54, 174 In the text  
• Colgate, S. A., & White, R. H. 1966, ApJ, 143, 626 In the text NASA ADS  
• Cooperstein, J., & Baron, E. 1992, ApJ, 398, 531 In the text NASA ADS  
• Dai, W. W., & Woodward, P. R. 1996, J. Comp. Phys., 124, 217 In the text  
• Dorfi, E. 1998, Radiation Hydrodynamics: Numerical Aspects and Applications, ed. O. Steiner, & A. Gautschy, Lecture Notes: Saas Fee Advanced Course 27 (Berlin: Springer) In the text  
• Einfeldt, B. 1988, SIAM Jour. Numer. Anal., 25, 294 In the text  
• Fryer, C. L. 1999, ApJ, 522, 413 In the text NASA ADS  
• Fryer, C. L., & Heger, A. 2000, ApJ, 541, 1033 In the text NASA ADS  
• Fryer, C. L., Heger, A., Langer, N., & Wellstein, S. 2002, ApJ, in press In the text  
• Fryxell, B., Müller, E., & Arnett, W. 1989, Hydrodynamics and Nuclear Burning, preprint MPA-449, Max Planck Institut für Astrophysik (Garching) In the text  
• Hannestad, S., & Raffelt, G. 1998, ApJ, 507, 339 In the text NASA ADS  
• Heger, A. 2000, personal communication In the text  
• Herant, M., Benz, W., & Colgate, S. 1992, ApJ, 395, 642 In the text NASA ADS  
• Herant, M., Benz, W., Hix, W. R., Fryer, C. L., & Colgate, S. A. 1994, ApJ, 435, 339 In the text NASA ADS  
• Hix, W. R., & Thielemann, F. 1999, ApJ, 511, 862 NASA ADS  
• Horowitz, C. 1997, Phys. Rev. D, 55, 4577 In the text NASA ADS  
• Horowitz, C. J. 2002, Phys. Rev. D, 65, 043001 In the text  
• Hummer, D., & Rybicki, G. 1971, MNRAS, 152, 1 In the text NASA ADS  
• Janka, H.-T. 1991, Ph.D. Thesis, Technische Universität München, preprint MPA-587 In the text  
• Janka, H.-T. 1992, A&A, 256, 452 In the text NASA ADS  
• Janka, H.-T. 1999, unpublished  
• Janka, H.-T. 2001, A&A, 368, 527 In the text NASA ADS  
• Janka, H.-T., Kifonidis, K., & Rampp, M. 2001, in Lecture Notes in Physics, vol. 578, Proc. Workshop on Physics of Neutron Star Interiors, ed. D. Blaschke, N. Glendenning, & A. Sedrakian (Springer), 333 In the text  
• Janka, H.-T., & Müller, E. 1995, ApJ, 448, L109 In the text NASA ADS  
• Janka, H.-T., & Müller, E. 1996, A&A, 306, 167 In the text NASA ADS  
• Keil, W. 1997, Ph.D. Thesis, Technische Universität München In the text  
• Keil, W., Janka, H.-T., & Müller, E. 1996, ApJ, 473, L111 In the text NASA ADS  
• Kifonidis, K. 2000, Ph.D. Thesis, Technische Universität München In the text  
• Kifonidis, K., Plewa, T., Janka, H.-T., & Müller, E. 2000, ApJ, 531, L123 In the text NASA ADS  
• Kippenhahn, R., & Weigert, A. 1990, Stellar Structure and Evolution, Astronomy and Astrophysics Library (Berlin: Springer) In the text  
• Körner, A. 1992, Diploma Thesis, Technische Universität München In the text  
• Lattimer, J., & Swesty, F. 1991, Nucl. Phys. A, 535, 331 In the text  
• LeVeque, R. J. 1992, Numerical Methods for Conservation Laws, Lectures in Mathematics (ETH Zürich: Birkhäuser) In the text  
• Liebendörfer, M., Messer, O., Mezzacappa, A., & Hix, W. 2001a, in Proc. of the 20th Texas Symp. on Relativistic Astrophysics, ed. J. Wheeler, & H. Martel (American Institute of Physics) In the text  
• Liebendörfer, M., Mezzacappa, A., Thielemann, F., et al. 2001b, Phys. Rev. D, 63, 3004 In the text  
• Lindquist, R. 1966, Annals of Physics, 37, 487 In the text  
• Lowrie, R. B., Mihalas, D., & Morel, J. E. 2001, J. Quant. Spec. Radiat. Transf., 69, 291 In the text NASA ADS  
• Lowrie, R. B., Morel, J. E., & Hittinger, J. A. 1999, ApJ, 521, 432 In the text NASA ADS  
• Mayle, R. W., Tavani, M., & Wilson, J. R. 1993, ApJ, 418, 398 In the text NASA ADS  
• Messer, O. E. B., Mezzacappa, A., Bruenn, S. W., & Guidry, M. W. 1998, ApJ, 507, 353 In the text NASA ADS  
• Mezzacappa, A., & Bruenn, S. W. 1993a, ApJ, 405, 669 In the text NASA ADS  
• Mezzacappa, A., Z Bruenn, S. W. 1993b, ApJ, 410, 740 In the text NASA ADS  
• Mezzacappa, A., Z Bruenn, S. W. 1993c, ApJ, 405, 637 In the text NASA ADS  
• Mezzacappa, A., Calder, A. C., Bruenn, S. W., et al. 1998a, ApJ, 493, 848 In the text NASA ADS  
• Mezzacappa, A., Calder, A. C., Bruenn, S. W., et al. 1998b, ApJ, 495, 911 In the text NASA ADS  
• Mezzacappa, A., Liebendörfer, M., Messer, O., et al. 2001, , 86, 1935 In the text NASA ADS  
• Mihalas, D. 1980, ApJ, 237, 574 In the text NASA ADS  
• Mihalas, D., & Klein, R. 1982, J. Comp. Phys., 46, 97 In the text  
• Mihalas, D., & Mihalas, B. 1984, Foundations of Radiation Hydrodynamics (Oxford University Press) In the text  
• Mihalas, D., & Weaver, R. 1982, J. Quant. Spec. Radiat. Transf., 28, 213 In the text NASA ADS  
• Miralles, J. A., Pons, J. A., & Urpin, V. A. 2000, ApJ, 543, 1001 In the text NASA ADS  
• Müller, E. 1991, in Supernovae, ed. J. Audouze, R. Bludman, R. Mochkovitch, & J. Zinn-Justin, Les Houches Summer School (Elsevier Science Publishers B.V.) In the text  
• Müller, E., & Janka, H.-T. 1997, A&A, 317, 140 In the text NASA ADS  
• Myra, E. S., Bludman, S. A., Hoffman, Y., et al. 1987, ApJ, 318, 744 In the text NASA ADS  
• Plewa, T., & Müller, E. 1999, A&A, 342, 179 In the text NASA ADS  
• Plewa, T., & Müller, E. 2001, Computer Physics Communications, 138, 101 In the text NASA ADS  
• Pons, J. A., Miralles, J. A., & Ibanez, J. M. A. 1998, A&AS, 129, 343 In the text NASA ADS  
• Pons, J. A., Reddy, S., Prakash, M., Lattimer, J. M., & Miralles, J. A. 1999, ApJ, 513, 780 In the text NASA ADS  
• Press, W., Teukolsky, S., Vetterling, W., & Flannery, B. 1992, Numerical Recipes in Fortran, 2nd ed. (Cambridge University Press) In the text  
• Rampp, M. 2000, Ph.D. Thesis, Technische Universität München In the text  
• Rampp, M., & Janka, H.-T. 2000, ApJ, 539, L33 In the text NASA ADS  
• Reddy, S., Prakash, M., & Lattimer, J. M. 1998, Phys. Rev. D, 58, 3009 In the text  
• Reddy, S., Prakash, M., Lattimer, J. M., & Pons, J. 1999, Phys. Rev. C, 59, 2888 In the text NASA ADS  
• Ruffert, M. 1992, A&A, 265, 82 In the text NASA ADS  
• Schinder, P. J., & Bludman, S. A. 1989, ApJ, 346, 350 In the text NASA ADS  
• Schinder, P. J., & Shapiro, S. L. 1982, ApJS, 50, 23 In the text NASA ADS  
• Shapiro, S., & Teukolsky, S. 1983, Black Holes, White Dwarfs, and Neutron Stars - the Physics of Compact Objects (New York: John Wiley and Sons)  
• Smit, J. M., & Cernohorsky, J. 1996, A&A, 311, 347 In the text NASA ADS  
• Smit, J. M., Cernohorsky, J., & Dullemond, C. P. 1997, A&A, 325, 203 In the text NASA ADS  
• Swesty, F. 1996, J. Comp. Phys., 127, 118  
• Thompson, T. A., Burrows, A., & Horvath, J. E. 2000, Phys. Rev. C, 62, 035802 In the text  
• Totani, T., Sato, K., Dalhed, H. E., & Wilson, J. R. 1998, ApJ, 496, 216 In the text NASA ADS  
• Tubbs, D., & Schramm, D. 1975, ApJ, 201, 467 In the text NASA ADS  
• van Riper, K. A. 1979, ApJ, 232, 558 In the text NASA ADS  
• Wilson, J. R. 1985, in Numerical Astrophysics, ed. J. Centrella, J. LeBlanc, & R. Bowers, Univ. Illinois, Urbana-Champaign US (Boston: Jones and Bartlett), 422 In the text  
• Wilson, J. R., & Mayle, R. 1988, Phys. Rep., 163, 63 In the text  
• Wilson, J.R., & Mayle, R. 1993, Phys. Rep., 227, 97 In the text  
• Wilson, J. R., Mayle, R., Woosley, S., & Weaver, T. 1986, in Annals of the New York Academy of Sciences, vol. 470, Twelfth Texas Symp. on Relativistic Astrophysics, ed. M. Livio & G. Shaviv, 267 In the text  
• Woosley, S. 1999, personal communication In the text  
• Woosley, S., Heger, A., Weaver, T., & Langer, N. 1997, SN 1987A - Presupernova Evolution and the Progenitor Star, preprint MPA-1024, MPI f. Astrophysik, Garching In the text  
• Woosley, S. E., Heger, A., & Weaver, T. A. 2002, Rev. Mod. Phys., in press  
• Yamada, S., Janka, H.-T., & Suzuki, H. 1999, A&A, 344, 533 In the text NASA ADS  
• Yorke, H. 1980, A&A, 86, 286 In the text NASA ADS  
• Yueh, W. R., & Buchler, J. R. 1977, ApJ, 217, 565 In the text NASA ADS

 

Online Material

  
Appendix A: Neutrino opacities

In this appendix we shall describe the neutrino-matter interactions that are included in the current version of our neutrino transport code for supernova simulations. We shall focus on aspects which are specific and important for their numerical handling and will present final rate expressions in the form used in our code and with a consistent notation.

   

Table A.1: Overview of all neutrino-matter interactions currently implemented in the code. For each process we list the sections where we summarize fundamental aspects of the calculation of the corresponding rate and give details of its numerical implementation. The references point to papers where more information can be found about the approximations employed in the rate calculations. In the first column the symbol $\nu$ represents any of the neutrinos ${\nu_{\rm e}},{\bar\nu_{\rm e}},\nu_\mu,\bar\nu_\mu,\nu_\tau,\bar\nu_\tau$, the symbols ${\rm e}^-$, ${\rm e}^+$, ${\rm n}$, ${\rm p}$ and A denote electrons, positrons, free neutrons and protons, and heavy nuclei, respectively, the symbol ${\rm N}$ means ${\rm n}$ or ${\rm p}$.

Reaction			Rate described in	Reference
$\nu {\rm e}^{\pm}$	$\rightleftharpoons$	$\nu {\rm e}^{\pm}$	Sects. A.1.2, A.2.4	Mezzacappa & Bruenn (1993b); Cernohorsky (1994)
$\nu {\rm A}$	$\rightleftharpoons$	$\nu {\rm A}$	Sects. A.1.2, A.2.3	Horowitz (1997); Bruenn & Mezzacappa (1997)
$\nu {\rm N}$	$\rightleftharpoons$	$\nu {\rm N}$	Sects. A.1.2, A.2.1	Bruenn (1985); Mezzacappa & Bruenn (1993c)
${\nu_{\rm e}}{\rm n}$	$\rightleftharpoons$	${\rm e}^- {\rm p}$	Sects. A.1.1, A.2.1	Bruenn (1985); Mezzacappa & Bruenn (1993c)
${\bar\nu_{\rm e}}{\rm p}$	$\rightleftharpoons$	${\rm e}^+ {\rm n}$	Sects. A.1.1, A.2.1	Bruenn (1985); Mezzacappa & Bruenn (1993c)
${\nu_{\rm e}}{\rm A}'$	$\rightleftharpoons$	${\rm e}^- {\rm A}$	Sects. A.1.1, A.2.2	Bruenn (1985); Mezzacappa & Bruenn (1993c)
$\nu\bar\nu$	$\rightleftharpoons$	${\rm e}^- {\rm e}^+$	Sects. A.1.3, A.2.5	Bruenn (1985); Pons et al. (1998)
$\nu\bar\nu ~{\rm N}{\rm N}$	$\rightleftharpoons$	${\rm N}{\rm N}$	Sects. A.1.3, A.2.6	Hannestad & Raffelt (1998)

A.1 Basic considerations

In the following we discuss the different neutrino-matter interaction processes which contribute to the source term ("collision integral'') on the rhs of the Boltzmann equation:

$$\frac{1}{c}\frac{\partial}{\partial t}~f + \mu\frac{\partia... ...mu^2}{r}\frac{\partial}{\partial \mu}~f = \sum_{\{I\}} B_I ~.$$ (A.1)

The sum on the rhs of this equation runs over all considered interaction processes that can change the distribution function of a particular type of neutrino. For simplicity we have written down the Boltzmann equation for static media, here. In the general case of a moving stellar fluid, velocity-dependent terms have to appear on the lhs of Eq. (A.1) when the neutrino quantities and the neutrino-matter interaction rates are evaluated in the local rest frame of the fluid. The phase-space distribution function f is related to the specific intensity ${\cal I}$ used in the main body of this paper by ${\cal I}=c/(2\pi\hbar c)^3\cdot\epsilon^3~f$. Accordingly, the source term on the rhs of the Boltzmann equation, Eq. (6), and the corresponding moment equations for neutrino number and neutrino energy or momentum (Eqs. (30), (31), (7), (8)) can be calculated from B_I as

 

$${C}(\epsilon,\mu) \frac{c}{(2\pi\hbar c)^3}\epsilon^3 \cdot \sum_{\{I\}} B_{I}(\epsilon,\mu) ~,$$ :=

$${C^{(k)}}(\epsilon) \frac{c}{(2\pi\hbar c)^3}\epsilon^3 \cdot \sum_{\{I\}} B^{(k)}_{I}(\epsilon) ~,$$ := (A.2)

$${\frak{C}^{(k)}}(\epsilon) \frac{c}{(2\pi\hbar c)^3} \epsilon^2 \cdot \sum_{\{I\}} B^{(k)}_{I}(\epsilon) ~,$$ :=

where $B^{(k)}_{I}(\epsilon):= 1/2\int_{-1}^{+1}{\rm d}\mu~\mu^k B_I(\epsilon,\mu)$. Here and in the following, the dependence of quantities on the space-time coordinates (t,r) is suppressed in the notation. Throughout Appendix A the temperature T is measured in units of energy.

  
A.1.1 Neutrino absorption and emission

The rate of change (modulo a factor of 1/c) of the neutrino distribution function due to absorption and emission processes is given by (see Bruenn 1985)

$$B_{\rm AE}(\epsilon,\mu)= j(\epsilon)[1-f(\epsilon,\mu)]-f(\epsilon,\mu)/\lambda(\epsilon) ~,$$ (A.3)

where j denotes the emissivity and $\lambda$ is the mean free path for neutrino absorption. The factor (1-f) accounts for fermion phase space blocking effects of neutrinos. Using the Kirchhoff-Planck relation ("detailed balance''), and introducing the absorption opacity corrected for stimulated absorption,

$$\kappa_{\rm a}^*:=\frac{1}{1-f^{\rm eq}}\cdot \frac{1}{\lambda}= j+\frac{1}{\lambda} ~,$$ (A.4)

Eq. (A.3) can be rewritten as

$$B_{\rm AE}(\epsilon,\mu)= \kappa_{\rm a}^*(\epsilon)[f^{\rm eq}(\epsilon)-f(\epsilon,\mu)] ~.$$ (A.5)

In Eq. (A.5) it becomes evident that the source term drives the neutrino distribution function $f(\epsilon,\mu)$towards its equilibrium value $f^{\rm eq}(\epsilon)=(1+\exp{[(\epsilon-\mu_\nu^{\rm eq})/T]})^{-1}$, where $\mu_\nu^{\rm eq}$ denotes the chemical potential for neutrinos in thermodynamic equilibrium with the stellar medium. In equilibrium, the source term vanishes in accordance with the requirement of detailed balance.

  
A.1.2 Scattering

Reduction of the collision integral:

The rate of change of the neutrino distribution function due to scattering of neutrinos off some target particles is given by the collision integral (cf. Cernohorsky 1994, Eqs. (2.1), (2.3))

 

$$B_{\rm S}(\vec{q})= \int \frac{{\rm d}^3\vec{\vec{q'}}~}{(2\pi)^3... ...},\vec{q'}) -f_\nu(1-f'_\nu)\widetilde{R}^{\rm out}(\vec{q},\vec{q'})\right] ~,$$ (A.6)

where $f'_\nu$ depends on $\vec{q'}$ and the scattering kernels $\widetilde{R}^{\rm in}$ and $\widetilde{R}^{\rm out}$ are defined as the following phase space integrals over products of the transition rate ${\cal R}$ and the Fermi distribution functions $F_{\rm T}$ of the target particles:

 

$$\widetilde{R}^{\rm in}(\vec{q},\vec{q'}) \frac{2}{c\cdot(2\pi \hbar c)^3}\cdot \int \frac{{\rm d}^3\vec{p'... ...\pi)^3} (1-F_{\rm T})F'_{\rm T}~ {\cal R}(\vec{p'},\vec{q'};\vec{p},\vec{q}) ~,$$ =

$$\widetilde{R}^{\rm out}(\vec{q},\vec{q'}) \frac{2}{c\cdot(2\pi \hbar c)^3}\cdot \int \frac{{\rm d}^3\vec{p'... ...\pi)^3} (1-F'_{\rm T})F_{\rm T}~ {\cal R}(\vec{p},\vec{q};\vec{p'},\vec{q'}) ~.$$ = (A.7)

The vectors $\vec{q}$ and $\vec{q'}$ are the momenta of the ingoing and outgoing neutrinos, respectively, and $\vec{p}$ and $\vec{p'}$ those of their scattering targets.

The scattering kernels are usually given as functions of the energies $\epsilon$ and $\epsilon'$ of the ingoing and outgoing neutrino and the cosine $\omega=\mu\mu'+\sqrt{(1-\mu^2)(1-{\mu'}^2)}\cos(\varphi-\varphi')$ of the scattering angle, where $(\mu,\varphi)$ and $(\mu',\varphi')$ are the corresponding momentum space coordinates. Expanding the scattering kernels ( $R^{\rm in/out}:=1/(2\pi c)^3\cdot\widetilde{R}^{\rm in/out}$) in a Legendre series

$$R^{\rm in/out}(\epsilon,\epsilon',\omega)= \sum_{l=0}^{\inft... ...2l+1}{2}~\phi_l^{\rm in/out}(\epsilon,\epsilon')P_l(\omega) ~,$$ (A.8)

with the Legendre coefficients

$$\phi_l^{\rm in/out}(\epsilon,\epsilon')= \int_{-1}^{+1}{\rm ... ...ega~ P_l(\omega){R}^{\rm in/out}(\epsilon,\epsilon',\omega) ~,$$ (A.9)

and applying the addition theorem for Legendre polynomials $P_l(\omega)=P_l(\mu)P_l(\mu')+2\sum_{m=1}^{l}\frac{(l-m)!}{(l+m)!}P_l^m(\mu)P_l^m(\mu')\cos[m(\varphi-\varphi')]$ (e.g., Bronstein & Semendjajew 1991, P_l^m are the associated Legendre polynomials), the integral over $\varphi'$ in the collision integral can be performed analytically (Yueh & Buchler 1977; Schinder & Shapiro 1982). For use in our Boltzmann transport code, it is necessary to truncate the Legendre expansion at some level. Note that the orthogonality relation $\int_{-1}^{+1}{\rm d}\omega~P_l(\omega)P_{l'}(\omega)=\frac{2}{2l+1}\delta_{ll'}$ (e.g., Bronstein & Semendjajew 1991) implies that any truncation of the Legendre series $R^{\rm in/out}(\epsilon,\epsilon',\omega)$ still gives the exact integral value $\int_{-1}^{+1}{\rm d}\omega~{R}^{\rm in/out}(\epsilon,\epsilon',\omega)$, independent of the level of truncation $l_{\rm max}>0$.

The collision integral and its first two angular moments finally read:

 

$$B_{\rm S}(\epsilon,\mu) 2\pi~\int_0^\infty {\rm d}\epsilon'~ {\epsilon'}^2$$ =

$$\Big\{(1-f_\nu) \cdot \sum_{l=0}^{\infty}(2l+1)P_l(\mu) \phi_l^{\rm in}(\epsilon,\epsilon')L'_l(\epsilon')$$

$$-f_\nu \cdot \sum_{l=0}^{\infty}(2l+1)P_l(\mu) \phi_l^{\rm out}(\epsilon,\epsilon')(\delta_{l0}-L'_l(\epsilon'))\Big\} ~,$$ (A.10)

 

$$B^{(0)}_{\rm S}(\epsilon) = 2\pi~\int_0^\infty {\rm d}\epsilon'~ ... ...lon) \phi_l^{\rm out}(\epsilon,\epsilon')(\delta_{l0}-L'_l(\epsilon'))\Big\} ~,$$ (A.11)

 

$$B^{(1)}_{\rm S}(\epsilon) = 2\pi~\int_0^\infty {\rm d}\epsilon'~ ... ...on)] \phi_l^{\rm out}(\epsilon,\epsilon')(\delta_{l0}-L'_l(\epsilon'))\Big\} ~,$$ (A.12)

where the "Legendre moment'' $L_l:=\frac{1}{2}\int_{-1}^{+1}{\rm d}\mu~~P_l(\mu) f(\mu)$ of order l of the distribution function has been introduced. These Legendre moments can obviously be written as linear combinations of the angular moments $M_m:=\frac{1}{2}\int_{-1}^{+1}{\rm d}\mu~\mu^m f(\mu)$, which occur in the formulation of the transport equations used in our code. Detailed balance requires

$$\phi_l^{\rm in}(\epsilon,\epsilon')= {\rm e}^{-(\epsilon-\epsilon')/T}~\phi_l^{\rm out}(\epsilon,\epsilon') ~.$$ (A.13)

Note that exploiting the in-out invariance of the transition rate (and therefore $\phi_l^{\rm out}(\epsilon,\epsilon')=\phi_l^{\rm in}(\epsilon',\epsilon)$; e.g. Cernohorsky 1994) leads to

$$\int_0^\infty{\rm d}\epsilon~\epsilon^2 B^{(0)}_{\rm S}(\eps... ...ty{\rm d}\epsilon~\epsilon^3 B^{(0)}_{\rm S}(\epsilon)\ne 0 ~,$$ (A.14)

which means that neutrino number is conserved in the scattering process, whereas there is a nonvanishing energy exchange between neutrinos and matter due to scatterings.

Isoenergetic scattering:

If for a particular scattering process the energy transfer between neutrinos and target particles can be neglected, we have $\phi_l^{\rm out}(\epsilon,\epsilon')= \phi_l^{\rm in} (\epsilon,\epsilon')=: \phi_l(\epsilon)\cdot\delta(\epsilon-\epsilon')$, and the source terms given by Eqs. (A.10), (A.11), (A.12) simplify to

   

$$B_{\rm IS}(\epsilon,\mu)= 2\pi\sum_{l=0}^{\infty}(2l+1)P_l(\mu)~ \epsilon^2\phi_l(\epsilon)L_l(\epsilon) -2\pi \epsilon^2\phi_0(\epsilon)~f_\nu ~,$$ (A.15)

$$B^{(0)}_{\rm IS}(\epsilon) = 0 ~,$$ (A.16)

$$B^{(1)}_{\rm IS}(\epsilon) = 2\pi ~\epsilon^2 L_1(\epsilon)\cdot\big(\phi_1(\epsilon)-\phi_0(\epsilon)\big)~.$$ (A.17)

  
A.1.3 Pair processes

Applying the procedures outlined in Sect. A.1.2 to the collision integral for the process of thermal emission and absorption of neutrino-antineutrino pairs by electron positron pairs and related reactions (e.g., nucleon-nucleon bremsstrahlung), the corresponding source terms read (see Bruenn 1985; Pons et al. 1998):

 

$$B_{\rm TP}(\epsilon,\mu) 2\pi\int_0^\infty {\rm d}\epsilon'~ {\epsilon'}^2$$ =

$$\Big\{ (1-f_\nu)\cdot \phi_0^{\rm p}(\epsilon,\epsilon') -\sum_{l=0}^{\infty}(2l+1)P_l(\mu) \phi_l^{\rm p}(\epsilon,\epsilon')\bar{L}_l(\epsilon')$$

$$+f_\nu \cdot\sum_{l=0}^{\infty}(2l+1)P_l(\mu) \phi_l^{\rm a}(\epsilon,\epsilon')\bar{L}_l(\epsilon')) \Big\} ~,$$ (A.18)

 

$$B^{(0)}_{\rm TP}(\epsilon) = 2\pi\int_0^\infty {\rm d}\epsilon'~ ... ...)\phi_l^{\rm a}(\epsilon,\epsilon') L_l(\epsilon)\bar{L}_l(\epsilon') \Big\} ~,$$ (A.19)

 

$$B^{(1)}_{\rm TP}(\epsilon) = 2\pi\int_0^\infty {\rm d}\epsilon'~ ... ...n') [ (l+1)L_{l+1}(\epsilon)+l L_{l-1}(\epsilon) ]\bar{L}_l(\epsilon') \Big\}~.$$ (A.20)

Here bars indicate quantities for antineutrinos. The absorption and production kernels are related by

$$\phi_l^{\rm a}(\epsilon,\epsilon')= \left(1-{\rm e}^{(\epsilon+\epsilon')/T}\right)\cdot\phi_l^{\rm p}(\epsilon,\epsilon')$$ (A.21)

in accordance with the requirement of detailed balance.

   
A.2 Numerical implementation of various interaction processes

In the following we describe the implementation of the various neutrino-matter interactions into our transport scheme. The expressions are mainly taken from the works of Tubbs & Schramm (1975), Schinder & Shapiro (1982), Bruenn (1985), and Mezzacappa & Bruenn (1993b,c).

All average values of the source terms within individual energy bins (cf. Eq. (20)) are approximated to zeroth order by assuming the integrand to be a piecewise constant function of the neutrino energy $\epsilon$, and hence:

$$B^{(k)}_{j+1/2}:= B^{(k)}(\epsilon_{j+1/2}) ~.$$ (A.22)

With the rest-mass energy of the electron, $m_{\rm e}c^2=0.511~{{\rm MeV}}$, and the characteristic cross section of weak interactions $\sigma_0:= 4(m_{\rm e}c^2 G_{\rm F})^2/(\pi(\hbar c)^4)=1.761\times 10^{-44}~{\rm cm}^2$($G_{\rm F}$ is Fermi's constant), we define:

$${\cal G}:=\frac{\sigma_0}{4~{m_{\rm e}}^2c^4} ~\cdot$$ (A.23)

  
A.2.1 Neutrino absorption, emission and scattering
by free nucleons

In the standard description of neutrino-nucleon interactions (see Tubbs & Schramm 1975; Bruenn 1985; Mezzacappa & Bruenn 1993c), many-body effects for the nucleons in the dense medium, energy transfer between leptons and nucleons as well as nucleon thermal motions are ignored. The corresponding simplifications allow for a straightforward implementation of the processes, using Eq. (A.5) for the charged-current reactions and Eqs. (A.15)-(A.17) for the neutral-current scatterings, respectively. Expressions for $\kappa_{\rm a}^*$ and the Legendre coefficients $\phi_0(\epsilon)$, $\phi_1(\epsilon)$ can be computed from the rates given by Bruenn (1985) and Mezzacappa & Bruenn (1993c).

For absorption of electron neutrinos by free neutrons ( ${\nu_{\rm e}}+ {\rm n} \rightarrow {\rm e}^- + {\rm p}$) the final result for the opacity is:

 

$$\kappa_{\rm a}^*(\epsilon)= {\cal G}\cdot(g_{\rm V}^2+3g_{\rm A}^... ...lon)} \;\eta_{\rm np} \cdot (\epsilon+Q)\sqrt{(\epsilon+Q)^2-m_{\rm e}^2c^4} ~.$$ (A.24)

The opacity of the absorption of electron antineutrinos by free protons ( ${\bar\nu_{\rm e}}+ {\rm p} \rightarrow {\rm e}^+ + {\rm n}$) is

 

$$\kappa_{\rm a}^*(\epsilon)= {\cal G}\cdot(g_{\rm V}^2+3g_{\rm A}^... ...silon)} \;\eta_{\rm pn} \cdot(\epsilon-Q)\sqrt{(\epsilon-Q)^2-m_{\rm e}^2c^4}~,$$ (A.25)

if the neutrino energy $\epsilon\ge m_{\rm e}c^2+Q$, and $\kappa_{\rm a}^*(\epsilon)=0$, else. Here, $Q:=m_{\rm n}c^2-m_{\rm p}c^2$ denotes the difference of the rest-mass energies of the neutron and the proton, and $F_{{\rm e}^\mp}$ are the Fermi distribution functions of electrons or positrons. The constants $g_{\rm V}=1$ and $g_{\rm A}=1.254$ are the nucleon form factors for the vector and axial vector currents, respectively. The quantities $\eta_{\rm np}$ and $\eta_{\rm pn}$ take into account the effects of fermion blocking of the nucleons. Making the approximations mentioned above, in particular ignoring the recoil of the nucleon, Bruenn (1985, Eq. (C14); see also Reddy et al. 1998) derived:

 

$$\eta_{\rm np} 2\int\frac{{\rm d}^3\vec{p}~}{(2\pi\hbar c)^3}F_{\rm n}(\epsilon)[1-F_{\rm p}(\epsilon)]$$ :=

$$\frac{n_{\rm p}-n_{\rm n}}{\exp{\left[\psi_{\rm p}-\psi_{\rm n}\right]}-1} ~,$$ = (A.26)

where $n_{\rm n}$, $n_{\rm p}$ are the number densities of neutrons and protons, respectively. The corresponding degeneracy parameters $\psi_{\rm n}$ and $\psi_{\rm p}$are calculated by inverting the relation (${\rm N}$ stands for ${\rm n}$ or ${\rm p}$)

$$n_{\rm N}=\frac{1}{2\pi^2}\left(\frac{2~m_{\rm N} c^2~T}{(\hbar c)^2}\right)^{3/2} F_{1/2}(\psi_{\rm N}) ~,$$ (A.27)

where F_1/2 is the Fermi integral of index 1/2 (see, e.g., Lattimer & Swesty 1991, Sect. 2.2). The quantity $\eta_{\rm pn}$ in Eq. (A.25) is defined accordingly by exchanging the subscripts n and p in Eq. (A.26). In the nondegenerate regime, where blocking is unimportant, one verifies the familiar limits $\eta_{\rm pn}=n_{\rm p}$ and $\eta_{\rm np}=n_{\rm n}$.

Scattering of neutrinos off free nucleons (${\rm n}$ or ${\rm p}$) is mediated by neutral currents only, which makes the distinction between neutrinos and antineutrinos unneccessary. Neglecting nucleon recoil and nucleon thermal motions, the isoenergetic kernel obeys ${R}_{\rm IS}(\epsilon,\omega)\propto (1+\omega)$ (e.g. Bruenn 1985; Mezzacappa & Bruenn 1993c), which implies $\phi_{l>1}(\epsilon)\equiv 0$. The non-vanishing Legendre coefficients read

  

$$\phi_0(\epsilon) \frac{{\cal G}}{8\pi}\cdot \left\{ \begin{array}{ll} \eta_{\rm nn... ..._{\rm V}-1)^2+\frac{3}{4}g_{\rm A}^2\right] &({\rm for\ p})~,\end{array}\right.$$ = (A.28)

$$\phi_1(\epsilon) \frac{{\cal G}}{24\pi}\cdot \left\{ \begin{array}{ll} \eta_{\rm n... ..._{\rm V}-1)^2-\frac{1}{4}g_{\rm A}^2\right] &({\rm for\ p})~.\end{array}\right.$$ = (A.29)

The factors

$$\eta_{\rm nn/pp}:=T\frac{\partial n_{\rm n/p}}{\partial \mu_{\rm n/p}}$$ (A.30)

approximately account for nucleon final state blocking (Bruenn 1985), and $C_{\rm V}=1/2+2\sin^2\theta_{\rm W}$ with $\sin^2\theta_{\rm W}=0.23$for the Weinberg angle. In order to accurately reproduce the limits $\eta_{\rm nn/pp}=n_{\rm n/p}$ for nondegenerate (and nonrelativistic), and $\eta_{\rm nn/pp}=n_{\rm n/p}\cdot 3T/(2 E^{\rm F}_{\rm n/p})$for degenerate (and nonrelativistic) nucleons ( $E^{\rm F}_{\rm n/p}$ is the Fermi energy) we prefer using the analytic interpolation formula proposed by Mezzacappa & Bruenn (1993c)

$$\eta_{\rm nn/pp}= n_{\rm n/p}\frac{\xi_{\rm n/p}}{\sqrt{1+\x... ...with}\quad \xi_{\rm n/p}:= \frac{3T}{2 E^{\rm F}_{\rm n/p}} ~,$$ (A.31)

instead of a direct numerical evaluation of Eq. (A.30).

  
A.2.2 Neutrino absorption and emission by heavy nuclei

As in the case of charged-current reactions with free nucleons, neutrino absorption and emission by heavy nuclei ( ${\nu_{\rm e}}+ A' \rightleftharpoons A + {\rm e}^-$) is implemented by using Eq. (A.5). For calculating the opacity of this process we adopt the description employed by Bruenn (1985) and Mezzacappa & Bruenn (1993c):

$$\kappa_{\rm a}^*(\epsilon)= {\cal G}\cdot g_{\rm A}^2~ \frac{2}{7... ...}(\epsilon)} ~ n_{A} \cdot(\epsilon+Q') \sqrt{(\epsilon+Q')^2-m_{\rm e}^2c^4}~,$$ (A.32)

if $\epsilon\ge m_{\rm e}c^2-Q'$, and $\kappa_{\rm a}^*(\epsilon)=0$ else. Here $Q':=m_{A'}c^2-m_{A}c^2\approx \mu_{\rm n}c^2-\mu_{\rm p}c^2 +\Delta$ denotes the energy difference between the excited state A'=(N+1,Z-1) and the ground state A=(N,Z). Following Bruenn (1985), we set $\Delta=3$ MeV for all nuclei. The internal structure of the nucleus with charge Z and neutron number N=A-Z is represented by the term $2/7~N_{\rm p}(Z)N_{\rm h}(N)$, with the functions $N_{\rm p}(Z)$, $N_{\rm h}(N)$ as given by Mezzacappa & Bruenn (1993c, Eqs. (31), (32)).

  
A.2.3 Coherent scattering of neutrinos off nuclei

We use reaction rates for coherent neutrino scatterings off nuclei which include corrections due to the "nuclear form factor'' (following an approximation by Mezzacappa & Bruenn 1993c; Bruenn & Mezzacappa 1997), and ion-ion correlations (as described in Horowitz 1997). For a detailed discussion of the numerical handling of both corrections, see the appendix of Bruenn & Mezzacappa (1997). The result of the latter reference can readily be used to calculate $B^{(1)}_{\rm IS,ions}$, the contribution of coherent scattering to the rhs of the first order moment equation. For the Boltzmann equation, we have to make a minor additional approximation: The correction factor $\left\langle\mathcal{S}(\epsilon)\right\rangle_{\rm ion}$ as provided by Horowitz (1997) applies for the transport opacity and thus for the combination $\phi_0-\phi_1$ (which is given by Bruenn & Mezzacappa 1997) but not necessarily for the Legendre coefficients $\phi_0$ and $\phi_1$ individually, which are needed for calculating the collision integral $B_{\rm IS,ions}$ for the Boltzmann equation. We nevertheless assume here that $\phi_{0/1}(\epsilon)\propto \left\langle{\cal S}(\epsilon)\right\rangle_{\rm ion}$. This implies

$$\phi_{0/1}(\epsilon)\approx \frac{1}{16\pi}~{\cal G}\cdot A^... ...y)\left\langle\mathcal{S}(\epsilon)\right\rangle_{\rm ion} ~,$$ (A.33)

with

$${\cal F}_0(y) \frac{2y-1+{\rm e}^{-2y}}{y^2} ~,$$ :=

$${\cal F}_1(y) \frac{2-3y+2y^2-(2+y){\rm e}^{-2y}}{y^3} ~,$$ := (A.34)

$$\frac{2}{5}\frac{(1.07A)^{2/3}\epsilon^2(10^{-13}{\rm cm})^2}{(\hbar c)^2} ~,$$ y :=

containing the effects of the nuclear form factor, and

$$\mathcal{Z}(A,Z):= \left[(C_{\rm A}-C_{\rm V})+(2-C_{\rm A}-C_{\rm V})\frac{2Z-A}{A}\right]^2 ~\cdot$$ (A.35)

$C_{\rm V}$ is defined as before (see Eqs. (A.28), (A.29)) and $C_{\rm A}=1/2$. In effect, the angular dependence of the ion-ion correlation correction is not exactly accounted for in the Boltzmann equation. The approximation, however, only influences the closure relation for the moment equations at a higher level of accuracy. It neither affects the Legendre coefficients at larger neutrino energies (where $\left\langle\mathcal{S}(\epsilon)\right\rangle_{\rm ion}$ is close to unity) nor does it destroy the consistency of the Boltzmann equation and its moment equations. Moreover, it has been shown that the overall effect of ion-ion correlations in stellar core collapse is only marginally relevant (Bruenn & Mezzacappa 1997; Rampp 2000).

  
A.2.4 Inelastic scattering of neutrinos off charged leptons

For the scattering of neutrinos off charged leptons, we approximate the dependence of the scattering kernel on the scattering angle by truncating the Legendre expansions (Eqs. (A.10)-(A.12)) at $l_{\rm max}=1$. For our purposes, this has been shown to provide a sufficiently accurate approximation (see Smit & Cernohorsky 1996; Mezzacappa & Bruenn 1993b), in particular also because the total scattering rate is correctly represented. Expressions for the Legendre coefficients are taken from the works of Yueh & Buchler (1977), Mezzacappa & Bruenn (1993b), Cernohorsky (1994), Smit & Cernohorsky (1996):

$$\phi_l^{\rm out}(\epsilon,\epsilon')= \alpha_{\rm I} A_l^{\... ...psilon') +\alpha_{\rm II} A_l^{\rm II}(\epsilon,\epsilon') ~,$$ (A.36)

with $k\in \{{\rm I},{\rm II}\}$ and

 

$$A_l^k(\epsilon,\epsilon')= \frac{{\cal G}}{\epsilon^2~\epsilon^{\... ...\rm e}(E_{\rm e}+\epsilon-\epsilon')\big)H^k_l(\epsilon,\epsilon',E_{\rm e}) ~.$$ (A.37)

The constants $\alpha_{\rm I/II}$, which depend on the neutrino species and also the type of charged lepton (i.e. electron or positron) involved in the scattering process, are combinations of the weak coupling constants (see, e.g., Cernohorsky 1994). For the $\mu$ and $\tau$ neutrinos and antineutrinos, which are treated as a single kind of neutrino in our code, we take the arithmetic mean of the corresponding coupling constants. For l=0,1 the functions $H^{\rm I/II}_l(\epsilon,\epsilon',E_{\rm e})$(in units of ${{\rm MeV}}^5$) are given by Yueh & Buchler (1977 with the corrections noted by Bruenn 1985).

Calculating the Legendre coefficients $\phi^{\rm in/out}_l(\epsilon_j,\epsilon_{j'})$ is a computationally expensive task, since for all combinations $(\epsilon_j,\epsilon_{j'})$and all radial grid points numerical integrals over the energy of the charged leptons (Eq. (A.37)) have to be carried out. It is, however, sufficient to explicitly compute $\phi^{\rm in/out}_l(\epsilon_j,\epsilon_{j'})$ for $\epsilon_j \le \epsilon_{j'}$. The missing coefficients can be obtained by exploiting a number of symmetry properties (see Cernohorsky 1994). In doing so, not only the computational work is reduced by almost a factor of two, but even more importantly, detailed balance ( $B_{\rm NES}\equiv 0$, if $f\equiv f^{\rm eq}$) can be verified for the employed approximation and discretization of the collision integral and its angular moments to within the roundoff error of the machine, provided that all integrals over neutrino energies are approximated by simple zeroth-order quadrature formulae (cf. (Eq. A.22)). Similarly, one can also verify the conservation of particle number (Eq. (A.14)) in the corresponding finite difference representation.

In our neutrino transport code the dependence of the source terms on the neutrino distribution and its angular moments is treated fully implicitly in time, i.e. the time index of the corresponding quantities in Eqs. (A.10)-(A.12) reads fⁿ⁺¹ and Lⁿ⁺¹_l. The thermodynamic quantities and the composition of the stellar matter which are needed to calculate the Legendre coefficients for inelastic scattering of neutrinos, however, are computed from $\varepsilon^{*}$ and ${Y_{\rm e}}^*$, representing the specific internal energy and electron fraction at the time level after the hydrodynamic substeps (cf. Sect. 3.6.1). This allows one to save computer time since the Legendre coefficients must be computed only once at the beginning of each transport time step. Comparing to results obtained with a completely time-implicit implementation of the source terms (i.e. using $\varepsilon^{n+1}$ and ${Y_{\rm e}}^{n+1}$ for determining the stellar conditions that enter the computation of the Legendre coefficients) we have found no differences in the solutions during the core-collapse phase, where neutrino-electron scattering plays an important role in redistributing neutrinos in energy space and thus couples the neutrino transport to the evolution of the stellar fluid.

  
A.2.5 Neutrino pair production/annihilation by thermal ${\rm e}^+~{\rm e}^-$ pairs

Following the suggestion of Pons et al. (1998) we truncate the Legendre series (A.18)-(A.20) for the pair-production kernels at $l_{\rm max}=2$. The Legendre coefficients are given by (Bruenn 1985; Pons et al. 1998)

 

$$\phi^{\rm p}_l(\epsilon_j,\epsilon_{j'})={\cal G}\cdot \frac{T^2}... ...epsilon_j,\epsilon_{j'}) +\alpha_2^2\Psi_l(\epsilon_{j'},\epsilon_j) \right) ~,$$ (A.38)

with $x:=(\epsilon+\epsilon')/T$ and ${\cal G}$ defined by Eq. (A.23). Explicit expressions for the (dimensionless) functions $\Psi_l(\epsilon_j,\epsilon_{j'})$ and details about their efficient calculation are given by Pons et al. (1998). Also the coefficients $\alpha_1,\alpha_2$ for the different neutrino flavours can be found there (Pons et al. 1998, Table 1).

  
A.2.6 Nucleon-nucleon bremsstrahlung

Legendre coefficients for the nucleon-nucleon bremsstrahlung process, assuming nonrelativistic nucleons, are calculated by using the rates given in Hannestad & Raffelt (1998). The non-vanishing coefficients read:

  

$$\phi^{\rm p}_0(\epsilon_j,\epsilon_{j'}) {\cal G}\cdot\frac{3~C_{\rm A}^2}{(2\pi)^2}~ \big({\rm e}^{-x}-1\... ...{\Gamma_\sigma}{(\epsilon+\epsilon')^2 + 0.25~(\Gamma_\sigma~g)^2}\cdot s(x) ~,$$ = (A.39)

$$\phi^{\rm p}_1(\epsilon_j,\epsilon_{j'}) -\frac{1}{9}\phi^{\rm p}_0(\epsilon_j,\epsilon_{j'})~,$$ = (A.40)

with

$$\Gamma_\sigma:=\frac{8\sqrt{2\pi}\alpha_\pi^2}{3\pi^2} \eta_... ... \quad \eta_*:=\frac{(p_{\rm F}~c)^2}{2 m_{\rm N} c^2 ~T}\cdot$$ (A.41)

In Eq. (A.40), $p_{\rm F}=\hbar\cdot(3\pi^2~n_{\rm N})^{1/3}$denotes the Fermi momentum for the nucleons, $\alpha_\pi=15$ is the pion-nucleon "fine-structure constant'' (see Hannestad & Raffelt 1998), and $m_{\rm N}$ is the nucleon mass. The analytic form of the dimensionless fit functions g and s(x), which both depend on the nucleon density $n_{\rm N}$ and temperature T of the stellar medium can be found in Hannestad & Raffelt (1998). The sum in Eq. (A.39) runs over the individual processes $\nu\bar\nu~{\rm nn}\rightleftharpoons {\rm nn}; ({\rm N}={\rm n})$, $\nu\bar\nu~{\rm pp}\rightleftharpoons {\rm pp}; ({\rm N}={\rm p})$, and $\nu\bar\nu~{\rm np}\rightleftharpoons {\rm np}; ({\rm N}={\rm np})$. The latter process is approximately taken into account by using the particle density $n_{\rm np}:=\sqrt{n_{\rm n} n_{\rm p}}$ and multiplying the corresponding contribution to the sum in Eq. (A.39) by a factor of 28/3 (cf. Thompson et al. 2000).

The production processes of neutrino-antineutrino pairs both by the annihilation of ${\rm e}^+~{\rm e}^-$ pairs and by the nucleon-nucleon bremsstrahlung are implemented into the discretized equations (Eqs. (21), (22), (28), (29)) by applying the techniques described for the inelastic scattering reactions (see Sect. A.2.4).