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)).

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.

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$.

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).

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).

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$).

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).

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).

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.

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).

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).

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.