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
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 
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
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.
Copyright ESO 2002