# A New Code for Proto-Neutron Star Evolution

###### Abstract

A new code for following the evolution and emissions of proto-neutron stars during the first minute of their lives is developed and tested. The code is one dimensional, fully implicit, and general relativistic. Multi-group, multi-flavor neutrino transport is incorporated that makes use of variable Eddington factors obtained from a formal solution of the static general relativistic Boltzmann equation with linearized scattering terms. The timescales of neutrino emission and spectral evolution obtained using the new code are broadly consistent with previous results. Unlike other recent calculations, however, the new code predicts that the neutrino-driven wind will be characterized, at least for part of its existence, by a neutron excess. This change, potentially consequential for nucleosynthesis in the wind, is due to an improved treatment of the charged-current interactions of electron flavored neutrinos and anti-neutrinos with nucleons. A comparison is also made between the results obtained using either variable Eddington factors or simple equilibrium flux-limited diffusion. The latter approximation, which has been frequently used in previous studies of proto-neutron star cooling, accurately describes the total neutrino luminosities (to within 10%) for most of the evolution, until the proto-neutron star becomes optically thin.

## 1. Introduction

A proto-neutron star (PNS) is born after the core of a massive star collapses to supra-nuclear densities, experiences core bounce due to the repulsive portion of the nuclear interaction which launches a shock wave that may eventually serve to disrupt the entire star in a supernova, and leaves behind a compact remnant. The overlying star is ejected and some portion of the mass may or may not fall back (c.f. Janka et al., 2007). In reality the mass of the PNS may increase with time due to this accretion, but a frequent assumption that is reasonable for low mass progenitors, and one adopted here, is that the PNS evolves in isolation after the shock has exited. Because of the large release of gravitational binding energy ( ergs), the PNS is initially hot and extended compared to a cold neutron star but large portions of the mass are still at supra-nuclear densities. Due to the high density and reasonably large temperature of this nuclear material, it is opaque to neutrinos of all flavors. In the outer regions of the PNS where the density is lower, the material is semi-transparent to neutrinos. This hot extended object undergoes Kelvin-Helmholtz cooling by emitting neutrinos of all flavors over a period of up to a minute (Burrows & Lattimer, 1986), at which time it transitions to a phase of optically thin neutrino cooling.

This qualitative description was confirmed when about twenty neutrinos were observed from SN 1987A (Bionta et al., 1987; Hirata et al., 1987). There were not enough events, however, to determine much detail of the cooling process (Lattimer & Yahil, 1989; Loredo & Lamb, 2002), though limits were placed on the properties of weakly interacting particles (Keil et al., 1997). If a similar core collapse supernova were to occur today, modern neutrino detectors would see thousands of events. Detailed modeling of the neutrino emission is needed if we are to learn about the central engine of core collapse supernovae from a nearby event.

There are numerous other reasons why understanding the properties of late-time supernova neutrinos is important, despite the rarity with which they are directly detected. An important one is the impact of neutrinos on supernova nucleosynthesis. Charged current neutrino interactions in the wind blown from the surface of PNSs determine the electron fraction of the ejected material and thereby constrain its nucleosynthesis. Current uncertainties in the relative energies of the electron neutrinos and anti-neutrinos are large enough to allow for both neutron-rich and proton-rich ejecta, which may be favorable for -process (Woosley et al., 1994) and -process nucleosynthesis (Fröhlich et al., 2006; Pruet et al., 2006), respectively. Recent work points to the wind ejecta being proton rich at all times (Hüdepohl et al., 2010; Fischer et al., 2010), but the reasons for this change are only beginning to be understood (Fischer et al., 2011). Additionally, the average energies of and neutrinos also significantly affect the neutrino spallation rates that determine nucleosynthetic yields of the -process (Woosley et al., 1990), which may be responsible for a number of rare isotopes.

It is also possible that current neutrino detectors with upgrades or next generation neutrino detectors will be able to observe the diffuse background of neutrinos produced by supernovae over the lifetime of the universe (Horiuchi et al., 2009). Predictions for the diffuse MeV scale neutrino background density depend significantly on the integrated spectrum of neutrinos emitted in core-collapse supernovae (Woosley et al., 1986; Ando, 2004). The integrated neutrino emission is dominated by PNS evolution, so that accurate modeling of PNSs can contribute to understanding the diffuse supernova neutrino background.

Finally, the neutrino emission from the “photosphere” of PNSs gives the initial conditions for the study of both matter-induced and neutrino-induced neutrino oscillations (Duan et al., 2006). The differences between the spectra of various neutrino flavors, especially and , can significantly affect the impact of flavor evolution in the nearly free streaming regime (Keil et al., 2003). The rate of PNS cooling also has the potential to put limits on exotic physics, such as axions (Keil et al., 1995), the presence of quark matter or a Kaon condensate in the core (Pons et al., 2001a, b), as well as possible extensions of the standard model using data already in hand from SN 1987A.

Theoretical predictions of post-bounce neutrinos have existed for more than 25 years (Burrows & Lattimer, 1986; Mayle et al., 1987; Keil & Janka, 1995; Sumiyoshi et al., 1995; Pons et al., 1999; Fischer et al., 2010; Hüdepohl et al., 2010; Roberts et al., 2012). Since the evolution of PNSs is described by the Kelvin-Helmholtz cooling of the collapsed, shock heated remnant of a core-collapse supernova, it is fundamentally a radiation hydrodynamics problem (although the regions important for neutrino emission are not very dynamic after bounce). Over time, the treatment of radiative transfer and neutrino microphysics in simulations has become increasingly sophisticated, moving from the equilibrium flux limited diffusion (EFLD) and greatly simplified neutrino physics (Burrows & Lattimer, 1986) to full solutions of the Boltzmann equation (Fischer et al., 2010) with more realistic microphysics (Hüdepohl et al., 2010).

Here, a new fully implicit code is developed for calculating the detailed evolution of PNSs in spherically symmetric general relativity within a variable Eddington factor formalism. The structure of the paper is as follows: In section 2, the equations of neutrino transport within the projected symmetric trace-free moment formalism of Thorne (1981) are described, and generic neutrino source terms for this formalism are derived. In section 3, a method for obtaining closure relations for the moment equations via a formal solution of the Boltzmann equation are described. A fully implicit numerical implementation of neutrino transport coupled to hydrodynamics/hydrostatics is described in section 4 (with code tests described in appendix A). A fiducial model of PNS cooling is detailed in section 5. These results are compared with the results of an EFLD calculation of PNS cooling in section 6.1. The implications of these new calculations of PNS cooling on the composition of the neutrino driven wind are discussed in section 6.2. In section 6.3, the properties of the integrated neutrino emission are discussed. The convention is adopted in sections 2 through 4 to avoid a plethora of factors. In section 6.2, units with are used.

## 2. The Moment Approach to General Relativistic Radiative Transfer

The equations of radiative transfer in curved space-times were first derived by Lindquist (1966), which described the evolution of the invariant distribution function along geodesics in phase space. The general form of the general relativistic Boltzmann (or Lindquist) equation in the absence of external forces is

(1) |

where is the invariant distribution function, is the neutrino four-momentum (which is constrained to be on mass shell), and are the Christoffel symbols. The collision term on the right hand side describes the destruction and production of neutrinos on a particular phase-space trajectory by capture processes, pair annihilation, scattering, and their inverses. In addition to describing the propagation of neutrinos along trajectories in physical space, this also encodes the evolution of the energy of neutrinos along geodesics of the spacetime. In three spatial dimensions, this is a seven dimensional equation that needs to be solved for each neutrino species.

A number of numerical strategies can be employed to solve the transport problem (Mihalas & Mihalas, 1984). Foremost among these are discrete ordinate methods, where the Boltzmann equation is directly discretized in momentum space as well as in physical space (e.g. Yueh & Buchler, 1977; Mezzacappa & Messer, 1999; Liebendörfer et al., 2004), and moment-based approaches, where angular integrations of the Boltzmann equation in momentum space are performed (e.g. Thorne, 1981; Burrows et al., 2000; Rampp & Janka, 2002). These two approaches give similar results in one-dimensional models, at least in the context of core-collapse supernovae (Liebendörfer et al., 2005). An additional technique that has only been employed for solving static problems in the supernova context, but is perhaps the most capable of retaining fidelity to the underlying Boltzmann equation, is Monte Carlo neutrino transport (Janka & Hillebrandt, 1989; Keil et al., 2003).

The moment approach results in an infinite hierarchy of coupled equations which needs to be truncated at some order in practice. Generally, only the zeroth and first order moment equations are retained and a closure relation is assumed between the first two moments and the higher order moments that enter the first two moment equations. Such schemes are referred to as variable Eddington factor methods (Mihalas & Mihalas, 1984). When only the first two moments are used, the number of equations relative to discrete ordinate methods is significantly reduced, easing the computational burden (especially in an implicit scheme like the one described below). Of course, this gain in computational efficiency is useful only if reasonable closures can be obtained. The closure relations only encode information about the angular distribution of neutrinos, so that the approximations involved in solving a linearized Boltzmann equation do not severely impact the fidelity of numerical calculations to the true solution (Mihalas & Mihalas, 1984; Ensman, 1994).

Here a variable Eddington factor approach to radiative transfer is employed, with the closure relations being obtained from a formal solution of the static relativistic Boltzmann equation. This approach is similar to that used by Burrows et al. (2000) and Rampp & Janka (2002), except for being fully general relativistic, incorporating both inelastic scattering and pair production (in contrast to only the former), using energy integrated groups rather than energy “pickets”, and in the specific method of finding the closure relations. The formalism for this method is described below.

The moments of the Boltzmann equation also most naturally give the various forms of the diffusion approximation, which has been used in the majority of PNS studies (Burrows & Lattimer, 1986; Keil & Janka, 1995; Pons et al., 1999; Roberts et al., 2012) and in a significant fraction of studies of the early core-collapse and bounce phases(Bruenn, 1985; Wilson & Mayle, 1993). The formalism is connected to EFLD in appendix B.

### 2.1. General Relativistic Generalities

In spherical symmetry, it is simplest to work in a coordinate system that anticipates a Lagrangian frame for the fluid. The metric for such a space-time is given by (Misner & Sharp, 1964)

(2) |

where is the invariant interval, is the time measured at infinity, is the areal radius, is the solid angle, and and are metric potentials. Coordinate freedom can be exploited to choose this frame to be the rest-frame of the fluid, which demands (Liebendörfer et al., 2001a)

(3) |

Here, is the baryon number density and

(4) |

where and are defined below. With this choice, the orthonormal frame associated with the coordinate frame is just the rest frame of the fluid and is just the change in enclosed baryon number with the physical volume. Therefore, this formulation is working in the Lagrangian frame, as claimed.

The equations of spherically symmetric general relativistic hydrodynamics and the Einstein equation are recorded for convenience (Misner & Sharp, 1964). Most of these results are nicely presented and detailed in similar form by Liebendörfer et al. (2001a). The time evolution of the areal radius is given by,

(5) |

which defines . The evolution of is given by

(6) |

where is the viscosity and is the pressure of the fluid. This gives the equation of hydrostatic balance when the left hand side equals zero (i.e. the Tolman-Oppenheimer-Volkov equation (Oppenheimer & Volkoff, 1939)). The enclosed gravitational mass, , is defined by

(7) |

where is the internal energy per baryon, is the total neutrino energy density in the rest frame, and is the net radial energy flux from neutrinos. The constraint equation for the metric potential is

(8) |

where a small time dependent term has been neglected.

The transport equations described in the next section are formulated in a congruence corresponding to the four-velocity field of the PNS (clearly, this is not a geodesic congruence). The behavior of this congruence is best described by expanding the covariant derivative of the four-velocity as

(9) |

where is the tangent four-vector field of the congruence, is the projection tensor (which projects into the vector subspace orthogonal to ) is the acceleration, is the expansion, is the shear, and is the rotation. Using the continuity equation, the expansion of the congruence becomes

(10) |

In spherical symmetry, the acceleration four-vector is parallel to the radial orthonormal basis vector, so that only the scalar acceleration is needed

(11) |

In spherical symmetry, the shear is characterized by a single component, the scalar shear

(12) |

Additionally, such a spherically symmetric congruence possesses no rotation, so that . The quantity

(13) |

will also be required, which is related to the extrinsic curvature (Thorne, 1981). The orthonormal frame temporal and radial derivative operators are

(14) |

and

(15) |

### 2.2. Variable Eddington Factor Transport Equations

Here, the evolution equations for the neutrino number density, energy density, number flux and energy flux are derived from the zeroth and first order moments of the relativistic Boltzmann equation. The basic results are taken from the spherically symmetric version of the projected symmetric trace-free moment formalism of Thorne (1981). This formalism reduces to an expansion of the neutrino distribution function in terms of Legendre polynomials in a flat space-time.

The moments of the distribution function are defined in spherical symmetry as

(16) |

where

(17) |

are the Legendre polynomials, and is the neutrino energy in the fluids rest frame. The first two moment equations in spherical symmetry can be read off from equation 5.10 of Thorne (1981)

(18) |

and

(19) |

where are the neutrino source terms defined in section 2.3. To close this system, define the Eddington like factors

(20) |

(21) |

which both go to zero in the limit , which corresponds to the diffusion regime. Note that these differ from the standard definition of the Eddington factors (Rampp & Janka, 2002), which is due to how I have chosen to calculate the moments. For free streaming radiation in a flat background , and these equations reduce to the linear wave equation for . A method for approximating these Eddington factors is detailed in section 3.

For problems that are close to being static on the radiation timescale, it is useful to switch the independent variable , the energy in the fluid rest frame, of to the energy at infinity, (c.f. Schinder & Bludman, 1989). In the case of PNS cooling, the energy at infinity is much closer to being a constant of the motion and therefore a more natural variable. Additionally, this choice simplifies the formal solution of the Boltzmann equation. The moments of the distribution function are then

(22) |

where the energy at infinity is defined as . This means that the replacement

(23) |

needs to be made for all radial and time derivatives, resulting in

(24) |

and

To easily deal with optically thick regions where the distribution function may possess a sharp Fermi surface, energy integrated groups are used rather than discrete energy “pickets”. The group numbers, energies, number fluxes, energy fluxes, and source terms in group are defined by

(26) |

Here, is the lower energy bound of an energy group and is the upper bound. Integrating over energy at infinity within groups gives

(27) |

and similar expressions for , , and the source terms. The operators and can then be applied to the “red shifted” equations. The evolution the neutrino group number densities are described by

(28) |

The last term on the left hand side describes the red or blue shifting of neutrinos to other groups via compression and time variation of the metric potential . If the group comprises energies from zero to infinity, the red shifting terms drop out and one is left with the standard number transport equation given in Pons et al. (1999). Applying the number operator to equation 2.2 and simplifying gives

(29) |

This includes similar terms to equation 2.2, plus a term that includes the effects of compression on the total number flux. The energy group evolution equations are

Aside from the addition of a compression term and different factors of , this is identical to equation 2.2. The energy flux group evolution equations are

(31) |

The numerical implementation of the red-shifting terms is described in section 4.6.

Additionally, neutrinos have a back-reaction on the matter they are propagating through by exchanging energy, lepton number, and momentum with the background medium. Assuming that the background possesses a thermal state, the first law of thermodynamics for the medium can be combined with the sum of equations 2.2 over all groups to find an equation for the conservation of total internal energy

(32) |

where the sums are over groups and species. Obviously, the neutrino energy source terms have exactly canceled with the source terms for the medium.

In the absence of neutrinos, the electron fraction of the background medium is fixed, i.e. . When neutrinos are included, interactions of electron flavored neutrinos exchange lepton number with the background, yielding . The total lepton number of the medium is given by . Combining the evolution equation for with equations 2.2 gives the lepton number evolution equation

(33) |

This constitutes the full set of evolution equations for the state of the medium including non-thermal neutrinos of all flavors, when the Eddington factors and are specified.

### 2.3. Neutrino Source Terms

The collision term in equation 1 describes how neutrinos move from one trajectory to another via scattering and how they are created and destroyed by the underlying medium. For the PNS problem these processes include neutral current scattering off of electrons, nucleons, and nuclei (Reddy et al., 1998), neutrino pair production via nucleon-nucleon bremsstrahlung (Hannestad & Raffelt, 1998) and electron-positron annihilation (Bruenn, 1985), and charged current processes involving electron and anti-electron flavor neutrinos and neutrons and protons, respectively (Reddy et al., 1998).

The details of these microphysical processes are eschewed by assuming that the differential cross-sections for these processes are known and referring the reader to the papers cited above, as well as the review Burrows et al. (2006). The exact details of the microphysics used in the code will be reported in a future publication, although certain aspects are discussed in sections 5 and 6.2. Many of the results in this section are well known (e.g. Bruenn, 1985; Pons et al., 1999), and are included here for completeness and to make clear the details of the exact implementation within the integrated energy group formalism described above. Explicit detailed balancing (independent of the choice of underlying scattering kernels) is emphasized.

The source function for a particular moment is given by

(34) |

which includes contributions from absorption , scattering , pair-annihilation , and their inverses . The choice of metric and reference frame implies that the scattering kernels should be evaluated in the rest frame of the fluid, simplifying things compared to hybrid frame approaches (Hubeny & Burrows, 2007).

For the absorption part, using the standard detailed balance relations gives

(35) |

where and is a Fermi-Dirac distribution. The scattering contributions are

(36) |

and

(37) |

and the pair-annihilation contributions are

(38) |

and

(39) |

The outgoing cosine is given by . The functions are related to the differential cross-section by

(40) |

with no phase space blocking term for the final neutrinos in the differential cross-section. The functions obey the detailed balance relations for scattering

(41) |

and annihilation

(42) |

For use in the moment formalism, must be expanded in terms of the Legendre polynomials as

(43) |

The distribution function, , is also expanded in a similar way. In general, this results in integrals of the form

(44) |

where . For a more detailed description of such an expansion, see Mezzacappa & Bruenn (1993).

Using this expansion, the scattering contribution to the source term is given by

(45) | |||||

and the pair annihilation contribution is given by

(46) | |||||

Clearly, all of the moments are coupled to all of the other moments by the source terms in addition to the coupling present on the LHS of the moment equations. Practically, this series must be truncated at some finite order. It is standard to use only the zeroth and first moment (Burrows et al., 2006). This convention is followed, but with the caveat that this may not be a good approximation for the annihilation terms near the free streaming regime (Pons et al., 1998).

#### 2.3.1 Zeroth Order Source Function

Now the three contributions to the source functions for the number and energy group equations are considered separately. Special attention is given to assuring that the chosen forms for the source terms explicitly push the neutrinos towards equilibrium, independent of the chosen opacity functions.

The absorption part of the source function is given by

(47) |

and

(48) |

where

(49) |

and

(50) |

The average over the inverse absorption mean free path can be performed in a number of ways. For small enough energy intervals for groups, the mean free path for the central energy of the group can be taken. It can also be assumed that the energy within a group is distributed as in a blackbody, but renormalized to the total energy within the group. Then the averaged inverse mean free path is analogous to the Planck mean opacity, but using a Fermi-Dirac distribution instead of a Planck distribution. Independent of the averaging procedure chosen, this term serves to push the neutrino energy density towards equilibrium with the background medium.

To find the scattering contribution to the zero order moment equation, it is assumed that the distribution of energy within a particular energy group is proportional to the blackbody distribution. Using this ansatz in equation 45 and then averaging over group energies gives scattering term

(51) |

and

(52) |

where

The averaged scattering kernel is defined as

(53) |

where the average is taken over the energies of the incoming and outgoing groups. This has the useful property , which derives from the detailed balance criterion given above. This form of the scattering source term naturally conserves neutrino number, although it does not push the neutrino numbers toward the expected distribution for a purely scattering process. Both source terms go to zero when neutrinos are in both chemical and energy equilibrium with the medium. The energy exchange expression does not go to zero when . Although it might be naively assumed that this corresponds to elastic scattering and therefore not contribute to the evolution of the group energy, there are in fact contributions from small energy transfer scatterings to this term as well. Although these will conserve neutrino number within the group, they can result in energy exchange with the medium. This allows the formalism to somewhat naturally deal with energy transfer due to scattering off nucleons, which generally exchanges energy on a scale that is smaller than the group spacing. The weighting function for the average over the groups necessarily involves some level of approximation. The natural incorporation of equilibrium far outweighs the small error introduced due to the approximate weighting of the scattering kernel.

Using a similar procedure to the one used for the scattering source term, the energy integrated pair production/annihilation source term is given by

(54) |

where the over bar denotes the energy density and flux of a neutrinos anti-species. This term once again naturally goes to zero when thermal equilibrium is reached (i.e. when and ). Here,

(55) |

The source term can be obtained from the above equation by the replacements , , , and , while leaving the terms unchanged.

#### 2.3.2 First Order Source Function

For the first order source function terms, one does not need to be as careful about getting forms that explicitly go to zero in equilibrium as all terms end up being proportional to the first order distribution function and therefore satisfy this constraint automatically.

The absorption contribution to the first order source function is

(56) |

and

(57) |

The scattering source term in equation 2.2 is

(58) |

The source term can be obtained from the above equation by the replacements , , , and , while leaving the terms unchanged. When scattering is iso-energetic, this reduces to

(59) |

where is the iso-energetic scattering kernel.

The first order moment equation pair annihilation source term is

(60) |

and can be found the same replacement required to find from .

## 3. Formal Solution of the Boltzmann Equation

To get the factors and , an angle dependent version of the Boltzmann equation needs to be solved. First, note that the outer layers of the PNS are in tight radiative equilibrium throughout the duration of the simulation. Therefore, all time dependence can be reasonably dropped in the equation of radiative transfer if one is only interested in the ratios of various moments. Of course, such an approximation breaks down in highly dynamical situations. For such circumstances, a closure scheme like the one described in Rampp & Janka (2002) is more appropriate. This time-independent formulation of the formal solution, which makes calculation of the Eddington factors significantly easier, is similar to the approach advocated by Ensman (1994), except that it incorporates general relativistic affects, such as the bending of geodesics. Schinder & Bludman (1989) describe a similar, but time-dependent formulation.

In a spherically symmetric static spacetime, the equation of radiative transfer is (Lindquist, 1966)

(61) |

Here, the neutrino distribution function, , has been written in terms of the energy of the neutrinos at infinity, , and is the cosine of the angle of neutrino propagation relative to the radial vector.

A formal solution to equation 3 can easily be found using the method of characteristics. The characteristic equations are

(62) |

where is the physical path length. The second equality is easily integrated to find a relationship between and along a geodesic. Any geodesic can be characterized by the radius at which . First, define the quantity

(63) |

where the subscript denotes the minimum radius of propagation. This is just the impact parameter of the trajectory. Then, for a given and , the angle of propagation along a geodesic is given by

(64) |

The first equality in the characteristic equations can be integrated to find the physical path length between any two radii for a particular characteristic if and are assumed constant over this distance, giving

(65) |

where the plus sign is for and the minus sign otherwise. This form is consistent with the assumption of constant metric functions across zones (as is used in the actual code), but it can introduce difficulties when a trajectory moves from one zone to another near the radius of minimum propagation.

Clearly, equation 3 is a non-linear integro-differential equation due to the functional dependence of the scattering terms on the local distribution function. An approximate solution to the Lindquist equation is desired were the solutions along characteristics are decoupled and the formal solution can be directly integrated. The simplest approximation is to just to make the replacement in the scattering terms. This approximation will only be valid at high optical depth and is therefore suspect for use in the decoupling region. The next order approximation is to use a distribution function inferred from our knowledge of and . Assuming that the energy is distributed within a group as within a blackbody gives

(66) |

Employing the Legendre expansion of the scattering kernel, integrating over outgoing neutrino angle, and only including the elastic scattering contribution gives the scattering source and sink terms

(67) | |||||

and

(68) |

Using the last characteristic equation, the solution of the linearized Boltzmann equation is

(69) | |||||

where the optical depth is

(70) |

This has the appealing property that there is no coupling between different s and s, so the evolution of the distribution function along each path in phase space can be solved for independently.