Redshift, gravitational bending, and the evolution of angular moments

In this subsection we test free streaming in spherically symmetric geometry. This is the opposite limit to the diffusion investigated in the previous subsection. For the time being, we stay with our stationary time slice at ms after bounce, but focus on radii larger than km. Outside this radius, we switch all interactions between the radiation field and matter off. This facilitates the comparison of the neutrino density, neutrino number flux, and neutrino luminosity with the analytical behavior. Unlike in the Newtonian limit, they are subject to time lapse, gravitational frequency shift, and gravitational bending.

For a stationary neutrino flux in the free streaming limit, Eq. (24) expresses particle number conservation

$$\frac{\partial }{\partial a}\left[ 4\pi r^{2}\alpha \rho H^{N}\right] =0.$$ (115)

Figure ()

Figure: Shown are the radial variations in the locally observed neutrino number luminosity (dashed line) and the conserved number luminosity (solid line) in a stationary-state situation at ms after bounce in the M model. Eq. () is fullfilled by construction.

shows the radial dependence of the locally observed neutrino number luminosity (dashed line) and the conserved number luminosity (solid line). We find that the latter is constant to machine precision as required by Eq. (). This is the direct result of our finite differencing of the term , because we have placed the lapse function inside the space derivative. The lapse function in the bracket converts the locally advected particle number per proper time to the particle advection per global coordinate time. There is a similar equation for the conservation of the energy luminosity. Eq. (30) for the particle energy carries a second lapse function inside the spatial derivative from the gravitational frequency shift,

$$\frac{\partial }{\partial a}\left[ 4\pi r^{2}\alpha ^{2}\rho H\right] =0.$$ (116)

Figure ()

Figure: Shown are the energy luminosities (dashed line), (dash-dotted line), and (solid line) as functions of radius in a stationary-state situation at ms after bounce in the M model. Deviations from Eq. () do not exceed in the third expression for the luminosity, which is expected to be conserved. The other two expressions illustrate that the gravitational effects contribute substantially in this test.

shows the radial dependence of the locally observed luminosity (dashed line), the luminosity with the time lapse correction only, (dash-dotted line), and the conserved luminosity, , with the gravitational redshift included. We see again, that the latter is fairly constant in the interaction free region. Unlike the number luminosity, the conservation of energy luminosity is not automatically fulfilled by the finite difference representation, and small deviations can be observed. As the local mean energies of the moving particles are determined by the ratio of energy and number flux, we may conclude from the fulfillment of Eqs. () and () that, in this regime, the gravitational energy shift of the particles is accurately implemented.

In order to fully constrain the redshift, gravitational bending, and evolution of the angular moments, we need to test the radial dependence of a third quantity. Most appropriate is the particle number density. However, some special care is necessary to clearly separate the following three effects on the angular neutrino distribution: (i) gravitational bending, (ii) the changing opening angle as a function of the geometric distance from the source, (iii) the numerical bending caused by limited angular resolution in our finite difference method. We find an analytical expression for the radial dependence of the particle density by applying the operator to the interaction-free stationary state limit of the Boltzmann equation (15),

$$\frac{\partial }{\alpha \partial a}\left( 4\pi r^{2}\alpha \rho J... ... \left( \frac{1}{r}-\frac{1}{\alpha }\frac{\partial \alpha }{\partial r}\right) \times$$

$$\left( 4\int F(\mu =0)E^{2}dE-\int \left[ F-F(\mu =0)\right] E^{2}dE\frac{1-\mu ^{2}}{\mu ^{2}}d\mu \right) = 0.$$ (117)

As before, we may compare this to our finite difference representation. The application of the operator to the finite difference representation of the Boltzmann equation leads to

$$\left( 4\pi r_{i+1}^{2}\alpha _{i+1}\rho _{i+1}J^{N}_{i+1}-4\pi r_{i}^{2}\alpha _{i}\rho _{i}J_{i}^{N}\right) =$$

$$-\alpha _{i'}\Upsilon _{i+1}\sum F_{i',j,k'}da_{i'}E_{k'}^{2}dE_{k'}\frac{\zeta _{j}}{\mu _{j'}\mu _{j'-1}}\left( \mu _{j'}-\mu _{j'-1}\right) .$$ (118)

Here, is interpolated according to Eq. () and according to Eq. (). Figure ()

Figure: This figure shows radial increments of the particle number density (actually plotted is the difference in the particle density across a zone, multiplied by ). The particle number luminosity is constant in the free streaming limit and the changes in the particle density depends on the angular distribution of the radiation field. The dashed line shows the changes induced by the increasing distance to the neutrino source in Newtonian geometry. The dash-dotted line shows the contribution from gravitational bending. The sum (solid line) represents the right hand side of Eq. (). It exactly matches the left hand side (circles). This indicates a correct derivation and implementation. However, the crosses mark a more accurate solution to the right hand side of Eq. () than the finite differencing imposed by the discrete Boltzmann equation. It has been evaluated according to Eq. (). The discrepancy is moderate at flux factors and large far from the source where should converge to .

compares the left hand side of Eq. () (circles) to the right hand side (solid line), evaluated in our ms post bounce time slice. The match is to machine precision; it demonstrates that we made no errors in the derivation of Eq. () or its implementation. We split this quantity into the contribution from geometric changes according to the opening angle of the source, proportional to , (dashed line) and the gravitational bending, proportional to , (dash-dotted line). Because the latter is much smaller, we exclude the gravitational bending contribution to as a primary source of uncertainties. We can then focus on the comparison of the analytical integral in Eq. () with its finite difference representation in Eq. (). We start with a more accurate numerical evaluation of the integral in Eq. () in order to compare with Eq. () in our time slice. The application of a Gaussian quadrature in Eq. () converts the integral to the sum

$$-\alpha _{i'}\Upsilon _{i+1}\left( 4\sum F(0)_{i',k'}E_{k'}^... ...^{2}dE_{k'}\frac{1-\mu _{j'}^{2}}{\mu _{j'}^{2}}w_{j'}\right)$$ (119)

The result of Eq. () depends quite sensitively on the choice of the distribution function in the tangential direction, . We determine it here by the maximum entropy model for Maxwell-Boltzmann statistics in angle space,

$$F_{j'}=A\exp \left( B\mu _{j'}\right) ,$$ (120)

(see e.g. Smit_VanDenHorn_Bludman_00 for an overview on maximum entropy closures). The zeroth and second moments of Eq. () define the flux factor, , as a function of the parameters and ,

$$J = \frac{2A}{B}\sinh (B)$$

$$h = \coth (B)-\frac{1}{B},$$ (121)

We derive the free parameters and by the inversion of Eq. () from the zeroth and first angular moments of the numerically obtained neutrino distribution function and define . The result is shown with cross markers in Fig. (). We observe comparable values between the two finite difference representations () and () at smaller radius, where the flux factor is smaller (). Large differences are found at larger radii, where the flux factor would be expected to be close to one. I.e., the primary source of inaccuracy in the particle density stems from the fact that the sum over angles on the right hand side of the finite difference Eq. () is a poor representation of the angle integral in the analytic Eq. ().

In order to study this discrepancy in more detail, we construct a series of analytical particle distribution functions according to the maximum entropy model in Eq. () with . This corresponds to flux factors of . Then, we evaluate the integral on the right hand side of Eq. () in different ways. First, as given by the finite differencing in Eq. () itself, then, with the more natural Gauss quadrature in Eq. (), and finally by an analytic integration. The analytic integration yields

$$\frac{1}{J}\int _{-1}^{1}\frac{1}{\mu }\frac{\partial }{\partial \mu }\left[ \left( 1-\mu ^{2}\right) A\exp \left( B\mu \right) \right] d\mu =$$

$$-2-Bh+\frac{B^{2}}{\sinh (B)}\left( B+\frac{B^{3}}{3\cdot 3!}+\frac{B^{5}}{5\cdot 5!}+\frac{B^{7}}{7\cdot 7!}+\ldots \right) .$$ (122)

The result for different flux factors and angular resolutions is shown in Table ().

Table: This table investigates the accuracy of the integral on the right hand side of Eq. () in a parameter study based on analytic distribution functions with different flux factors. The result of the numerical integrations are given for discretizations with , , and angular bins. The uppermost block shows the evaluation according to Eq. (). It emerges from the chosen finite differencing of the term in Eq. () and converges only slowly to the exact solution. The center block shows the same evaluation without upwind differencing for the advected (i.e. with ). The convergence is much better in this case. The lowermost block shows the evaluation based on the more natural Gauss quadrature (), which converges very rapidly to the exact result.

	B	flux factor	j =6	j =12	j =48	analytical
Eq. ()	0	0	-2.00	-2.00	-2.00	-2.00
	1	0.31	-1.10	-1.23	-1.36	-1.41
	2	0.54	-0.21	-0.22	-0.28	-0.32
	5	0.80	0.28	0.47	0.68	0.77
	10	0.90	0.06	0.14	0.26	0.31
Eq. (),	0	0	-2.00	-2.00	-2.00	-2.00
but no	1	0.31	-1.48	-1.43	-1.41	-1.41
upwind	2	0.54	-0.47	-0.36	-0.32	-0.32
diff.	5	0.80	0.66	0.74	0.77	0.77
	10	0.90	0.27	0.30	0.31	0.31
Eq. ()	0	0	-2.00	-2.00	-2.00	-2.00
	1	0.31	-1.41	-1.41	-1.41	-1.41
	2	0.54	-0.32	-0.32	-0.32	-0.32
	5	0.80	0.77	0.77	0.77	0.77
	10	0.90	0.32	0.31	0.31	0.31

We find quite poor convergence for the evaluation of the integral as it emerges from the finite differencing in the Boltzmann equation. The comparison with the center block, where no upwind differencing was used for the advection terms, demonstrates that the advective terms play a major rôle in the accuracy of the representation of this integral. A forward peaked radiation field shows a steep gradient in the angular distribution. The choice of upwind differencing always underestimates the advective flux at the edge of the angular zone. This is necessary for the stability of the scheme for arbitrary large time steps, but can lead to a substantial underestimation of the integral in Eq. (). The most populated forward bin, for example, never contributes to the integral. The distribution function asymptotically approaches a highly populated forward bin while all other angular bins are emptied. The maximum achievable flux factor is smaller than one. Although higher angular resolution improves the situation, it does not eliminate this systematic effect. The choice of a higher order advection scheme in angular space might be promising. More rigorous, however, appears the implementation of an adaptive angular grid in the transparent regime Yamada_Janka_Suzuki_99. It would be designed to minimize the advective flow in angle space such that the undesired effects are eliminated at the root. This difficulty with angular advection does not occur in methods that solve a model Boltzmann equation along particle characteristics to close the system of radiation moments equations with a variable Eddington factor Burrows_et_al_00,Rampp_Janka_02. An independent and much smaller source of inaccuracy stems from the fact that the integral in Eq. () is not written as a sum over a function evaluated at Gaussian quadrature points and multiplied with Gaussian weights. Hence, it does not take profit from the fast convergence of Gaussian quadrature. This can be seen if the result is compared with the evaluation based on the Gaussian quadrature as shown in the third block in Table ().

In our supernova application, however, stationary-state convergence tests have shown that angle Gaussian quadrature produces physically reasonable results Messer_et_al_98,Messer_00. We confirm this in section in a comparison with a dynamical run discretized with angular bins. In the regions of free streaming, where numerical angular diffusion becomes apparent, the neutrino field is already decoupled from the matter and does not influence the dynamical evolution anymore.