Chimera: A Massively Parallel Code for Core-collapse Supernova Simulations

, , , , , , , , , , , and

Published 2020 April 30 © 2020. The American Astronomical Society. All rights reserved.
, , Citation Stephen W. Bruenn et al 2020 ApJS 248 11 DOI 10.3847/1538-4365/ab7aff

Download Article PDF
DownloadArticle ePub

You need an eReader or compatible software to experience the benefits of the ePub3 file format.

0067-0049/248/1/11

Abstract

We provide a detailed description of the Chimera code, a code developed to model core collapse supernovae (CCSNe) in multiple spatial dimensions. The CCSN explosion mechanism remains the subject of intense research. Progress to date demonstrates that it involves a complex interplay of neutrino production, transport, and interaction in the stellar core, three-dimensional stellar core fluid dynamics and its associated instabilities, nuclear burning, and the fundamental physics of the neutrino–stellar core weak interactions and the equations of state of all stellar core constituents—particularly, the nuclear equation of state associated with core nucleons, both free and bound in nuclei. Chimera, by incorporating detailed neutrino transport, realistic neutrino–matter interactions, three-dimensional hydrodynamics, realistic nuclear, leptonic, and photonic equations of state, and a nuclear reaction network, along with other refinements, can be used to study the role of neutrino radiation, hydrodynamic instabilities, and a variety of input physics in the explosion mechanism itself. It can also be used to compute observables such as neutrino signatures, gravitational radiation, and the products of nucleosynthesis associated with CCSNe. The code contains modules for neutrino transport, multidimensional compressible hydrodynamics, nuclear reactions, a variety of neutrino interactions, equations of state, and modules to provide data for post-processing observables such as the products of nucleosynthesis, and gravitational radiation. Chimera is an evolving code, being updated periodically with improved input physics and numerical refinements. We detail here the current version of the code, from which future improvements will stem, which can in turn be described as needed in future publications.

Export citation and abstract BibTeX RIS

1. Introduction

Modeling core collapse supernovae (CCSNe) has become an extremely demanding computational problem over the years, requiring realistic multidimensional, general relativistic, multigroup, neutrino radiation hydrodynamics; sophisticated nuclear equations of state; and extensive nuclear reaction networks to elucidate the CCSN explosion mechanism and capture some of the important observables. The sheer breadth of interdependent physics, over a vast range of density and energy scales, is one aspect of the computational challenge. This includes general relativistic gravity; matter velocities at nonnegligible fractions of the speed of light; the production, transport, and interaction of neutrinos and anti-neutrinos of all flavors across three regimes: the tight coupling of neutrinos and matter at high densities in the core, weaker coupling far from the core, and intermediate coupling in between; the evolution of the nuclear composition both in and out of nuclear statistical equilibrium (NSE) as mediated by both strong and weak nuclear interactions; and an equation of state that spans a density range that can exceed fourteen orders of magnitude. And there is a second notable aspect of the computational challenge. The CCSN explosion mechanism also appears to be marginal in the sense that rather modest changes in the numerical modeling of the neutrino transport and/or the neutrino interactions, and/or the use of Newtonian versus general relativistic gravity, for example, can change the outcome of a simulation not only quantitatively but qualitatively as well. Therefore, the physical treatment and numerical implementation of the above-described physics has to be realistic and highly accurate (in the numerical sense) if meaningful results are to be obtained. Finally, a CCSN simulation may require millions of time steps to integrate a model forward in time through a sufficiently long period in order to determine outcomes such as explosion energies and other observables. This requires that the various numerical algorithms implemented in the code be highly optimized.

Our development of the radiation hydrodynamics code Chimera to model CCSNe has drawn on previous codes that have successfully modeled one or another physical process relevant to CCSNe. The hydrodynamics module has been built on the dimensionally split, Lagrangian-plus-remap scheme with piecewise parabolic reconstruction as formulated by Colella & Woodward (1984) and implemented in VH1 as described by Hawley et al. (2012) and Blondin & Lufkin (1993), but extended to include multi-species advection, energy absorbed or released due to compositional changes, nuclei coming in or out of NSE, multidimensional gravity, momentum and energy exchange with neutrinos, and a sliding radial grid algorithm that continually adjusts the radial grid to resolve structures that arise during the course of a simulation. The neutrino transport stems from the original formulation in Bruenn (1985) but modified and improved, and with a number of neutrino source terms either refined or added to the Chimera suite of neutrino–matter interactions. The nuclear composition, when material is not under conditions appropriate for NSE, is evolved by the thermonuclear reaction network code XNet developed by Hix & Meyer (2006) and Travaglio & Hix (2013).

In this paper, we present a detailed description of our Chimera code, which has been developed with the aim of realistically modeling CCSN and the associated observables. In the sections that follow, we describe each of our algorithms in detail. Chimera is under continuous development. When it is deemed that significant improvements or refinements have been implemented and tested, the version of the code current at that time is frozen, removed from further development, and used to execute and follow new simulations in our ongoing investigations of CCSNe. These frozen code versions, along with common sets of microphysics inputs are designated as a lettered "Series" of Chimera simulations, where simulations within a "Series" share not only the code base and algorithm choices, but the default microphysics and control parameters as well.

The first test results from Chimera for low-resolution 2D (axisymmetric) simulations were reported in Bruenn et al. (2006), as was an early description of the code (Messer et al. 2008). The first attempt at production simulations in 2D and 3D, with a reasonably complete set of physics and with reasonable resolution were reported in Bruenn et al. (2009) and retroactively designated as "Series-A." Though some flaws emerged late in the Series-A simulations, the earlier portions of the 2D simulations were used for gravitational wave extraction analysis (Yakunin et al. 2010). Subsequent to the Series-A simulations, we made several improvements to and enhanced Chimera, as well as updated some input physics. The "Series-B" simulations (Bruenn et al. 2013, 2016) used the updated code to recompute the four 2D simulations of Series-A, with enhanced resolution, from which we analyzed the gravitational wave signals (Yakunin et al. 2015) and nucleosynthesis of the ejecta (J. A. Harris et al. 2018, in preparation). The methods and algorithms described in this paper largely reflect those used in the Chimera Series B. The next group of models, Series-C, focused on the 3D modeling (Lentz et al. 2015; E. J. Lentz et al. 2018, in preparation) using the same input microphysics as Series B. Direct analyses of the gravitational wave (Yakunin et al. 2017) and neutrino detector signals (O. E. B. Messer et al. 2018, in preparation) have been performed for these simulations. The primary differences between Series-B and Series-C were code consolidation and optimizations, but a few improvements are described herein, including an improved treatment of neutrino transport through the shock (Section 6.12), a more efficient interpolation of neutrino opacities (Section 8.11), and parallel IO (Section 2.3). Further series will follow, including a general "Series D" consisting of several 3D models and 2D studies, and a "Series E" focused on the nuclear equation of state and related code improvements. All of the above mentioned series have employed the so-called ray-by-ray (RbR) plus approximation to neutrino transport in which spherically symmetric multi-energy neutrino transport is performed along each radial ray with lateral advection of neutrinos with matter performed in optically thick regions. Fully multi-energy, multi-angle neutrino transport has recently been implemented in core-collapse simulations with various levels of microphysics Kuroda et al. (2016, 2018), Ott et al. (2018), Vartanyan et al. (2018), and Skinner et al. (2019) Multi-angle transport will be a future upgrade of Chimers. Descriptions of the modifications and enhancements for these, and subsequent, series will appear with the simulation results, as needed.

Section 2 gives a general overview of Chimera including the domain decomposition used in implementing parallel computing architectures, the directional splitting used, and the sequence of computational steps in a complete time cycle. Section 3 details the implementation of the nuclear reaction network and equations of state, including the technique for juxtaposing more than one equation of state in adjacent density regimes. Section 4 presents a detailed description of our hydrodynamics algorithms, with test problems in Section 5. The neutrino transport method and numerical methods are described in Section 6, with a detailed description of the neutrino transport source terms given in Section 8. Finally, Section 7 presents results of static neutrino transport tests, and Section 9 presents comparisons of 1D CCSN simulations with the Boltzmann code Agile-Boltztran (Liebendörfer et al. 2004, 2005). Additional specific test results verifying a number of our algorithms are presented at the end of the relevant sections.

2. General Overview

To drive the component pieces integrated from other codes, Chimera uses a master-data model, where data is stored in a master copy on each computing element, transposed as needed to appropriate dimensional sweeps, checked-out to constituent codes, processed, and checked back to the master copy. At the master level, general program control, data input and output, and monitoring are also managed.

2.1. Directional Splitting

For multidimensional applications, Chimera uses the method of dimensional splitting (Strang 1968). In this method, the numerical solution of the multidimensional problem proceeds by a series of one-dimensional steps or "sweeps" to build up the full solution. The updated variables after each sweep are used as initial conditions for the next sweep. The order of the sweeps is shown in Figure 1. For second-order accuracy in time, the time step, Δt, should be selected at the beginning of the sweep sequences, and for the 2D case, $\tfrac{1}{2}{\rm{\Delta }}t$ should be used in the x-y subsequence of the sweep sequence and again in the y-x subsequence. Likewise, in the 3D case, $\tfrac{1}{4}{\rm{\Delta }}t$ should be used for each of the four x-y-z transposition subsequences. In practice, the full time step computed before each subsequence is used for that subsequence. Since the time steps computed for each of the subsequences are nearly the same, approximate time-centering is maintained. Furthermore, the larger time steps permit a more refined grid to be used for the same amount of computer time.

Figure 1.

Figure 1. Panel (a): sweep sequence used in Chimera for numerically solving one-, two-, and three-dimensional problems. Panels (b) and (c): sequence of operations performed during the (b) radial sweep and (c) theta or phi sweep.

Standard image High-resolution image

2.2. Domain Decomposition

Chimera is designed to numerically evolve CCSNe on spherical polar grids and to run on multiprocessor machines. The domain decomposition in Chimera is dictated, in part, by the present implementation of neutrino transport, which is in the radial direction. The decomposition is along rays, or bundles of rays, each ray being one zone wide in two dimensions and spanning the entire set of zones from boundary to boundary along the third. Thus, an x-ray, or radial ray refers to a set of zones along a radial line from the center to the outer edge of the grid. A y-ray, or angular ray, refers to a set of zones at a given radius and azimuth which, for a 180° angular grid, spans the arc from the "north" pole to the "south" pole. Finally, a z-ray, or azimuthal ray, refers to a set of zones at a given radius and polar angle that, in the case of a 360° azimuthal grid, completely encircles the polar axis.

As an extremely simple example of the domain decomposition in Chimera, consider a 2D grid consisting of four radial rays, each radial ray consisting of 12 radial zones with one radial ray per MPI rank, four MPI ranks in all. Figure 2 shows the logical structure of this grid. The large rectangular block bordered in thick black lines encompasses the logical grid with the x (radial), y (angular), and z (azimuthal) directions shown at the upper right. Because the grid is two-dimensional, the third or z-dimension is superfluous but is shown as an elongated cell wall in the z-direction for the sake of comparison with a 3D example below. Each zone is shown in the figure as a square-sided vertical rectangular volume. The total number of radial, angular, and azimuthal zone-centers are denoted by imax, jmax, and kmax, respectively, and the total number of zone edges by ${\mathtt{imax}}+1$, ${\mathtt{jmax}}+1$, and ${\mathtt{kmax}}+1$. In this case, ${\mathtt{imax}}=12$, ${\mathtt{jmax}}=4$, and kmax = 1.

Figure 2.

Figure 2. Example Chimera domain decomposition for two-dimensional models with one x-ray per MPI rank.

Standard image High-resolution image

The computation in Chimera is directionally split, and the (x-, y-, z-) sweeps—i.e., the (radial-, angular-, azimuthal-) sweeps—refer to the direction of computation. During the x-sweep, a set number, or bundle, of radial rays (in this example just one) is assigned to each MPI rank. In Figure 2, a radial ray is shown in green for a particular, but otherwise arbitrary, MPI rank. The dimensions ij_ray_dim and ik_ray_dim denote the dimensions of a bundle of x-rays in the y- and z-directions, respectively, so that ij_ray_dim × ik_ray_dim is the number of x-rays in the bundle. In this example, where the bundle of rays per MPI rank consists of just one ray, both ij_ray_dim and ik_ray_dim are unity. The local indices ij_ray and ik_ray locate a particular ray within a bundle relative to the upper left corner of the bundle, so that ij_ray ∈ [1, ij_ray_dim] and ik_ray ∈ [1, ik_ray_dim]. In this example, both of these indices are unity as there is only one ray per bundle. During the x-sweep, the radial hydrodynamics, radial RbR transport, and nuclear reactions are evolved along with global gravity solves.

Following the x-sweep, a transpose to the y- or z-oriented rays is performed. In this 2D example, the transpose is just to the y-oriented rays. Because there are 12 radial zones and only four angular zones in this example, each MPI rank now consists of a bundle of three y-rays, so that the total number of MPI ranks remain the same. This is delineated in Figure 2 by the rays enclosed by red for a particular, but arbitrary, bundle. Now, (ji_ray, jk_ray) takes the place of (ij_ray, ik_ray) and locates a particular angular ray in the xz plane of the bundle. The widths in the xz plane of each bundle are given by j_ray_dim and ik_ray_dim, which in this example are equal to 3 and 1, respectively. j_ray_dim and k_ray_dim are the number of radial zones on each MPI rank after transposing to the y-oriented rays and z-oriented rays, respectively.

Figure 3 illustrates a slight variation of the previous domain decomposition example using the same logical grid. The total number of radial, angular, and azimuthal zones are the same as before, but in this example, bundles of two x-rays are associated with each MPI rank, two MPI ranks in all. In Figure 3, a particular, but arbitrary, bundle of x-rays is shown bounded by thick green lines. In this case, the yz dimensions of the bundle are given by ij_ray_dim = 2 and ik_ray_dim = 1, respectively, and the particular, but arbitrary, radial ray designated by the narrow-lined green X at the front and top has the local indices in the bundle (ij_ray = 2, jk_ray = 1). Following the x-sweep, a transpose to the y-oriented rays is made, and because there are only two MPI ranks in this case, the number of y-rays associated with each MPI rank is six, as there are 12 y-rays in all. This is delineated in Figure 3 by the rays enclosed by the thick red lines for a particular, but arbitrary, y-bundle, The xz dimensions of the y-bundles are j_ray_dim = 6 and ik_ray_dim = 1, respectively, and the particular, but arbitrary, y-ray designated by the narrow-lined red X at the front and top has the local indices in the bundle (ji_ray = 5, jk_ray = 1).

Figure 3.

Figure 3. Example Chimera domain decomposition for two-dimensional models with two x-rays per MPI rank. Some lines have been slightly offset to render them visible.

Standard image High-resolution image

Figure 4 illustrates a general domain decomposition example using the same logical grid for the x- and y-rays, but adding a z-ray consisting of 6 zones. Now there are ${\mathtt{imax}}\times {\mathtt{jmax}}\times {\mathtt{kmax}}=12\times 4\times 6=288$ zones in all, and ${\mathtt{jmax}}\times {\mathtt{kmax}}=4\times 6=24$ x-rays. Let us suppose we wish to compute with 4 MPI ranks. We must then use 4 bundles of x-rays, each bundle assigned to one MPI rank and consisting of 6 x-rays. One such particular, but arbitrary, bundle is shown in Figure 4, outlined in the thick green lines. The y-z dimensions of the bundle are ij_ray_dim = 2 and ik_ray_dim = 3, respectively. A particular, but arbitrary, x-ray is shown in this bundle by a thin-lined green X on its front face, having the local indices (ij_ray = 1, ik_ray = 2); another x-ray located in an adjacent bundle is delineated by another thin-lined green X on its front face having the local indices (ij_ray = 2, ik_ray = 3). On executing the x-sweep, each MPI rank performs the computation required to complete the individual x-sweep for each x-ray in its bundle.

Figure 4.

Figure 4. Example Chimera domain decomposition for a three-dimensional model. Some of the lines have been slightly offset to render them visible.

Standard image High-resolution image

Following the x-sweep, a transpose to the y- and z-oriented rays is performed, either in that order, or reversed order, as shown above in Figure 1. Consider the transpose to the y-oriented rays. There are ${\mathtt{imax}}\times {\mathtt{kmax}}=12\times 6=72$ y-rays. With four MPI ranks, we assign to each MPI rank four bundles of y-rays with 18 y-rays each, having the xz dimensions j_ray_dim = 6 and ik_ray_dim = 3, respectively. A particular, but arbitrary, bundle of y-rays is shown bounded by thick red lines in Figure 4, and a particular, but arbitrary, y-ray in that bundle is singled out by a thin-lined red X, having the local indices (ji_ray = 4, jk_ray = 3).

Following or preceding the y-sweep, a transpose to the z-oriented rays is performed. There are ${\mathtt{imax}}\times {\mathtt{jmax}}\,=12\times 4=48$ z-rays. With four MPI ranks, we assign to each MPI rank four bundles of z-rays with 12 z-rays each, having the xy dimensions k_ray_dim = 4 and ij_ray_dim = 2, respectively. A particular, but arbitrary, bundle of z-rays is shown bounded by thick blue lines in Figure 4, and a particular, but arbitrary, z-ray in that bundle is singled out by a thin-lined blue X, having the local indices (${\mathtt{ki}}\_{\mathtt{ray}}=3,{\mathtt{kj}}\_{\mathtt{ray}}=2$).

In general, for a 2D or 3D grid with a domain decomposition consisting of x-rays, y-rays, and z-rays, to ensure load balancing, the number number of rays in x-bundles, y-bundles, and z-bundles must satisfy

Equation (1)

where NMPI is the number of MPI ranks.

2.3. IO Subsystem

The input and output (IO) of data by simulation codes presents a significant challenge, and good IO is required to achieve acceptable computational performance. Chimera has several components to its IO subsystem. A typical full explosion model requires ${ \mathcal O }(1000)$ hours to complete. This requires writing checkpoint files for restarting the simulation and data files for later analysis, plotting, and visualization. We have implemented two schemes for general restart and analysis IO: (1) the original, serial method used for much of our early 2D work and (2) a parallel method required for effective computation in 3D.

2.3.1. Serial IO Method

The original serial IO scheme had separate components for restart and analysis output. The restart IO consisted of one file written per MPI rank, containing only the necessary data to restart the simulation and to maintain tracking of conservation, including the positions of the Lagrangian tracer particles associated with the processors domain (see Section 2.4). These sets of files were written on typically 100-cycle intervals to alternating files in a pair. At predetermined points in the computation, "permanent" restart files were saved, to be used to provide the initial data if a simulation had to be "rewound" to fix a problem.

These restart files were supplemented with plot files written on fixed time intervals, typically every 0.2 ms after core bounce. Each file consisted of single variable (e.g., entropy, radial velocity, ${\bar{\nu }}_{e}$ luminosity, etc.) data for the entire grid and was assembled by gathering the data to the root processor before writing it out. For 3D runs with more than ∼10,000 radial rays, the memory available on the root node was typically insufficient. These plot files did not contain all information needed for restart, but they did contain derived quantities (like luminosity) useful for visualization and analysis.

The above binary files were also supplemented with in situ analysis output of global properties (shock radius information, explosion energies, radial traces along fixed angles, etc.) written to plain text files. To generate the "trace files" for each Lagrangian tracer, the thermodynamic, abundance, and neutrino quantities were interpolated to the particle position and recorded in a binary tracer file.

2.3.2. Parallel IO Method

To achieve scaling to larger process counts in 3D and improved file performance, we implemented parallel IO with the HDF5 library.9 The HDF5 library permits complex file structures, metadata, and file portability. Initially, the HDF5 implementation was a replacement for the restart IO, without the alternating file scheme, with the analysis and plotting data added from the binary one-variable files. To improve the time resolution of the gravitational wave analysis from the Δt = 0.2 ms resolution used in Series-A (Yakunin et al. 2010) and Series-B (Yakunin et al. 2015) analyses, we added a finer-resolution sampling of the quantities needed to compute the matter contribution to the gravitational wave signal—i.e., density and velocity—to the HDF5 "Restart" files.

For Series-D, we are moving toward a fully HDF5-based system for large IO. We have implemented fixed-time-interval "Frame" HDF5 files, without the extra data for gravitational wave analysis, to replace the single-variable equivalents. The density and velocity information at finer intervals is retained in "GW" HDF5 at intervals matching those of the "Restart" files from which the "GW" data have now been removed. This separation of the data renders the "Restart" HDF5 files unnecessary after serving their primary role as checkpoint files for restarting a simulation. (The "Frames" files contain the exact same data as the "Restart" files and can be used to restart a simulation, as well, if needed.) The combination of "Frames"+"GW" HDF5 files retained for analysis is smaller than the set of "Restart" HDF5 files retained for the C-series.

2.4. Tracer Particles

In addition to the IO of grid-based data, Chimera also outputs data for passive Lagrangian tracer particles. These are used to record the thermodynamic and neutrino exposure histories of individual mass elements that are then used for post-processing analysis. Following each directional sweep of the hydrodynamics, the position of a tracer particle in that direction at time tn and position $({r}^{n},{\theta }^{n},{\phi }^{n})$ is advanced to ${t}^{n+1}$ according to the simple Euler method:

Equation (2)

Equation (3)

Equation (4)

assuming constant velocity $({u}_{r}^{n},{u}_{\theta }^{n},{u}_{\phi }^{n})$ through the time interval ${\rm{\Delta }}{t}^{n}={t}^{n+1}-{t}^{n}$. For the B-series models, which were all 2D and included 4000–8000 particles, individual output files were maintained for each tracer. These contained the physical quantities of interest, linearly interpolated in radius to the tracer particle positions from the zone-center (cell-averaged) values of the computational grid; the lone exception being the interpolation of differential neutrino number fluxes, which are defined at radial zone edges. Output for individual particles was performed whenever the physical quantities of interest changed by 10% (see Harris et al. 2017, for more details).

For the C-series and later models, which included 3D models with as much as 400,000 tracers, this approach proved impractical; thus, Chimera simply records the tracer positions in "Frames" and "Restart" files, with interpolation of quantities of interest relegated to a post-processing step. Nucleosynthesis tests by Harris et al. (2017) show that output intervals of Δt ∼ 1 ms are sufficient in CCSN to resolve the features in the thermodynamic profiles needed for accurate post-processing of the Lagrangian tracers. Thus, the computational time expended on the temporally detailed individual particle traces written out in the B-series models can now be better spent elsewhere.

2.5. Scalability

To test the scalability of Chimera, we computed 3D models with all of the standard physics of the B- and C-series, started from the same 15 ${M}_{\odot }$ progenitor used in the B15-WH07 run (Bruenn et al. 2013, 2016), with 512 radial zones. Models were run from the onset of collapse for a couple of hours for nθ of 32, 64, 128, and 256, with nϕ = 2nθ, on the Cray XK7 ("Titan") at OLCF, utilizing "maximal decomposition" with one radial ray per process (MPI rank). The results plotted in Figure 5 use the version of the code used for the C15-3D model (Lentz et al. 2015) and show slightly more than doubled wall time for the largest model relative to the smallest, but the typical size of our actual production 3D runs are closer to the 32768-ray (128 × 256) model, which takes about 25% longer than the smallest model to complete 100 steps.

Figure 5.

Figure 5. Chimera weak scaling for 100 cycles of full calculation on the Cray XK7 (OLCF Titan) for 32 × 64 (2048) rays to 256 × 512 (131072) rays, with one ray per core

Standard image High-resolution image

2.6. Chimera Series

The lettered series of Chimera simulations are not in the strict sense of software engineering and development practices "versions" of the code, but rather represent a combination of default physics, parameter choices, and a common code base used in those simulations. As it has occurred, all of the series discussed in this paper have utilized different code bases, with progressive improvement with each subsequent series, but in the future, we may use the series designation to separate groups of simulations with different defaults run with the same code. The series label combined with a basic description of the progenitor and any deviations from the series defaults are used to construct unique simulation designations to identify a model across publications.

Initially, there were separate 2D and 3D branches of the code and both were used for the simulations retroactively designated "Series A" (Bruenn et al. 2009) with the 2D branch using only the binary IO for restart and plotting and the 3D branch using HDF5 for restarts. The Series-B runs included many improvements to the code and simulation design, reflecting the acquired experience of Series A. While the 2D-only Series B code was running, the 3D code was optimized and became the only branch of Chimera for all subsequent runs of any dimension and included changes to the code structure needed for the Yin-Yang grid. The Yin-Yang mesh was not yet satisfactory when Series C was started and the major differences between that and Series D are the finalized Yin-Yang grid and an improved radial-mesh motion routine. Major improvements were made to the handling of the NSE state and the boundary with the non-NSE material handled by the network before Series-B, with geometric upgrades before Series D permitting free islands of NSE or non-NSE material to embed, as necessary, in the other. Some operational choices in the dual EoS scheme were made for Series B, while Series E added pre-tabulated equations of state (R. E. Landfield et al. 2020, in preparation). The temporal sub-cycling of the non-radial sweeps deep in the core generated entropy and was replaced for Series B by a frozen core, wherein the non-radial hydrodynamics was skipped, which was in turn replaced by a spherical averaging scheme late in Series B and used in all later series.

In the list that follows, we describe major differences between the code revisions used in the series indicated. For descriptions of the utilized physics and default choices for resolution, the associated simulation papers should be referenced.

  • 1.  
    Series A
    • (a)  
      EoS:LS180, ρ > $1.7\times {10}^{8}\,{{\rm{g\; cm}}}^{-3}$
    • (b)  
      Hydrodynamics: subcycle lateral sweep
    • (c)  
      NSE: spherical boundary
    • (d)  
      IO: binary restart and plot files
  • 2.  
    Series B
    • (a)  
      EoS: LS220, ρ > $1\times {10}^{11}\,{{\rm{g\; cm}}}^{-3}$
    • (b)  
      Hydrodynamics: frozen core, carbuncle update
    • (c)  
      NSE: 17-species (α, n, p, ${}^{56}\mathrm{Fe}$), boundary a function of latitude
  • 3.  
    Series C
    • (a)  
      Opacities: shared points (Section 8.11.5)
    • (b)  
      Hydrodynamics: averaged core
    • (c)  
      Transport: shock fix (Section 6.12), Minerbo closure replaces geometric closure
    • (d)  
      IO: merged restart and plot quantities into single HDF5 file
    • (e)  
      Mesh: Yin-Yang 3D grid-related changes

3. Equation of State

In order for the hydrodynamics, nuclear transmutations, and neutrino transport to be tied closely to the thermodynamics, the equation of state (EoS) must be invoked several times each cycle (See Figure 1). Furthermore, the EoS must provide not only the quantities needed for the hydrodynamics—e.g., the pressure, internal energy, and entropy as a function of density, temperature, and electron fraction—but the element composition and chemical potentials, as well, as these are needed for the computation of the opacities which, in turn, are needed in the neutrino transport.

Furthermore, the derivatives of a number of these thermodynamic quantities are needed to compute, by a Newton–Raphson iteration, updates in one or more independent variables given updates in other independent and dependent variables. For example, we need the derivative of the internal energy with respect to temperature to update the temperature given updates in the internal energy, density, and electron fraction. And we need the derivatives of the internal energy with respect to the temperature and electron fraction (as well as derivatives of the neutrino opacities with respect to these variables) to build the Jacobian for the neutrino transport solve.

3.1. General EoS Methods

Two general cases must be distinguished when considering the thermodynamic state of the fluid, NSE, and non-NSE. In NSE, the thermodynamic state of the fluid in a given zone is specified by the values of its density, ρ, temperature, T, and electron fraction, ${Y}_{{\rm{e}}}$. To accommodate the demand for frequent EoS interrogations and the need for derivatives of some of the dependent thermodynamic quantities, Chimera constructs a thermodynamic grid in ($\mathrm{log}\rho ,\mathrm{log}T,{Y}_{{\rm{e}}}$)-space defined by a user specified number, ${\mathtt{dgrid}}({\mathtt{m}})$, of evenly spaced points of $\mathrm{log}\rho $ per decade change in ρ, a user specified number, ${\mathtt{tgrid}}({\mathtt{m}})$, of evenly spaced points of $\mathrm{log}T$ per decade change in T, and a user specified number, ${\mathtt{ygrid}}({\mathtt{m}})$, of points in ${Y}_{{\rm{e}}}$ over the range [0, 1]. The index ${\mathtt{m}}=1,2,3$ allows the user to select three different thermodynamic grid resolutions for the three density ranges $\rho \lt {\rho }_{\mathrm{es}}(1)$, ${\rho }_{\mathrm{es}}(1)\lt \rho \lt {\rho }_{\mathrm{es}}(2)$, and $\rho \gt {\rho }_{\mathrm{es}}(2)$, where ${\rho }_{\mathrm{es}}(1)$ and ρes(2) are user selected densities. A particular dependent thermodynamic function, corresponding to the thermodynamic state ($\mathrm{log}\rho ,\mathrm{log}T,{Y}_{{\rm{e}}}$), is computed by linear interpolation from its values at the eight surrounding grid points, which satisfy

Equation (5)

where the ρi, Tj, and ${Y}_{{\rm{e}},k}$ are the values of ρ, T, and ${Y}_{{\rm{e}}}$ at the grid points. Figure 6 shows a cell of the thermodynamic grid within which the thermodynamic state, ($\mathrm{log}\rho ,\mathrm{log}T,{Y}_{{\rm{e}}}$), is located. We will refer to the eight grid points surrounding a given mass zone as the "surrounding grid points," and the cell itself simply as the "EoS cell." We emphasize that not all grid points of the thermodynamic grid have thermodynamic quantities evaluated and stored there, but only those grid points surrounding the thermodynamic states of mass zones, with a cell of grid points tied to each zone. Initially, the thermodynamic state of each mass zone gets a suite of thermodynamic quantities computed and stored at the eight surrounding grid points. During a simulation, when the changing thermodynamic state of a mass zone causes the state to enter a different EoS cell, the needed thermodynamic quantities are in turn computed and stored on the grid points of that cell. Thus, thermodynamic quantities (and neutrino opacities) are computed on an "as needed" basis, keeping the thermodynamic state of each mass zone surrounded by the needed thermodynamic quantities on the eight nearest ($\mathrm{log}\rho ,\mathrm{log}T,{Y}_{{\rm{e}}}$)-grid points. Lastly, to avoid involving an excessive number of quantities in internode communication when transposing from one set of rays (radial, angular, or azimuthal) to another, the EoS grids along these rays are maintained independently.

Figure 6.

Figure 6. The thermodynamic state, $(\mathrm{log}\rho ,\mathrm{log}T,{Y}_{{\rm{e}}})$, of a zone inside a grid element in $\mathrm{log}\rho ,\mathrm{log}T$, and ${Y}_{{\rm{e}}}$ space.

Standard image High-resolution image

A total of 14 dependent EoS variables comprise the thermodynamic vector that is computed and stored at each of the grid points of a cell surrounding the thermodynamic state of a mass zone. These quantities are the pressure, p; specific internal energy, e; specific entropy, s; neutron chemical potential, μn; proton chemical potential, μp; electron chemical potential, μe; neutron mass fraction, Xn; proton mass fraction, Xp; representative heavy nucleus mass fraction XH (nuclei with mass numbers greater than helium), along with the mass number, A, charge number Z, and mean binding energy per particle bA of the representative heavy nucleus; the adiabatic exponent ${{\rm{\Gamma }}}_{s}={\left(\partial p/\partial \rho \right)}_{s,{Y}_{{\rm{e}}}}$, and the specific internal energy, eint, with the particle rest masses and arbitrary constants subtracted out. The helium mass fraction, Xα, is not stored but is computed as ${X}_{\alpha }=1-{X}_{{\rm{n}}}-{X}_{{\rm{p}}}-{X}_{{\rm{a}}}$. The specific internal energy, eint, is used to compute the quantity ${{\rm{\Gamma }}}_{e}=p/{e}_{\mathrm{int}}\rho +1$, utilized by the Riemann solver to prevent unphysical ${{\rm{\Gamma }}}_{e}\lt 4/3$, which can cause post-shock oscillations in the some of the thermodynamic variables (Buras et al. 2006b). The formula used to interpolate a thermodynamic quantity from its values at the eight surrounding grid points is differentiated (exactly) to obtain partial derivatives of that quantity with respect to either ρ, T, or ${Y}_{{\rm{e}}}$, as needed.

To provide a sense of the accuracy of our EoS interpolation scheme, Figures 7 and 8 show the relative deviation of the pressure, specific internal energy, and the neutron and proton chemical potentials obtained by direct output from our stellar EoS (described below) versus interpolation in the EoS grid, for the grid resolution listed in the figures. The ρ, s, and ${Y}_{{\rm{e}}}$ profiles used for generating Figures 7 and 8 are representative of models near bounce but before the formation of the shock, and many tens of ms after bounce during shock stagnation. We will refer to these profiles as "Near Bounce" (Figure 7) and "Shock Stagnation" (Figure 8). It is clear from the figures that increasing the density resolution from 10 to 20 grid points per decade in density does not decrease the relative deviation of these thermodynamic quantities except for densities above nuclear saturation. (The black and red curves lie on top of each other below the saturation density.) Increasing the temperature resolution from 20 to 40 points per decade in temperature (orange curve) reduces the relative deviation for all of the graphed quantities throughout most of the density range displayed. Increasing the electron fraction resolution from 50 to 100 over the range [0, 1] in ${Y}_{{\rm{e}}}$ decreases the relative deviation for quantities (namely, abundances and chemical potentials) in regimes (low entropy) where there is a strong dependence on ${Y}_{{\rm{e}}}$—e.g., where there is partial dissociation. The B-series models were performed with the EoS grid resolution of (d-, t-, ygrid) = 20, 20, 50, which, except for a few exceptions described below, typically obtains values of interpolated thermodynamic quantities within a percent or so of the values obtained directly from the EoS (in most cases less than a percent). The D-series models are being performed with the higher grid resolution of (d-, t-, ygrid) = 20, 40, 100, which typically gives a relative deviation about five times smaller.

Figure 7.

Figure 7. Relative deviation of representative thermodynamic quantities obtained by direct output from the EoS vs. by interpolation, for the listed EoS grid resolutions and for a "Near Bounce" profile typical of a stellar core just prior to the formation of the bounce shock. Panels (a), (b), (c), and (d) show the relative deviation for the pressure, internal energy, neutron chemical potential, and proton chemical potential, respectively.

Standard image High-resolution image
Figure 8.

Figure 8. As in Figure 7 but for a "Shock Stagnation" profile representative of a stellar core during the epoch of shock stagnation.

Standard image High-resolution image

A few features of the graphs deserve comment. One feature is the slight kink in the temperature at ρ = 1011 ${{\rm{g\; cm}}}^{-3}$ in the "Near Bounce" profile (Figure 7(a)), and is is due to the LS EoS—C EoS (see Section 3.2) transition at that density. Another is the substantial relative error in the neutron chemical potential at a density of 2 × 1014 ${{\rm{g\; cm}}}^{-3}$ for the "Near Bounce" profile, and at a density of 2.5 × 1014 ${{\rm{g\; cm}}}^{-3}$ for the "Shock Stagnation" profile. These are the densities for the respective profiles at which the neutron chemical potentials pass through zero, so any slight deviation in the interpolated versus the directly obtained values for these quantities will be amplified by their small absolute values when computing their relative deviations. The region where the neutron chemical potentials change sign is shown in Figure 9 for the "Near Bounce" profile but is representative of the "Shock Stagnation" profile as well. The spikes in the relative deviation of all quantities at 1.3 × 1014 ${{\rm{g\; cm}}}^{-3}$ for the "Near Bounce" profile deserve mention, as well. These are caused by the nuclei–nuclear matter phase transition at that density, and the abrupt change in composition there. Table interpolation smooths this transition across the width of the density grid, while direct calls to the EoS see this transition as a discontinuity. These spikes do not appear in the "Shock Stagnation" profile because the higher entropy results in matter being completely dissociated at the above density, leading to a smoother transition across the nuclei–nuclear matter phase transition.

Figure 9.

Figure 9. Neutron chemical potential for the "Near Bounce" profile, obtained directly from the EoS (black line) vs. by interpolation for (d-, t-, ygrid) = 20, 40, 100 (red line). Slight differences in the chemical potential while it passes through zero at densities 2–$3\times {10}^{14}\,{{\rm{g\; cm}}}^{-3}$ are the cause of the large relative errors there. The "Shock Stagnation" profile yields similar results and discrepancies.

Standard image High-resolution image

3.2. NSE EoSs

For CCSN simulations, the pressure, specific internal energy, and specific entropy are taken as the sum of contributions from different species, namely,

Equation (6)

where the subscripts "ion," "e + e+," and "rad" denote contributions from nuclei, electrons and positrons, and photons, respectively. For the B-, C-, and D-series runs, Chimera employs, for densities above 1011 ${{\rm{g\; cm}}}^{-3}$, the K = 220 MeV incompressibility version of the Lattimer & Swesty (1991) (LS) EoS for the ion and photon components. (The retroactively named A-series used the K = 180 MeV version of LS EoS.) The LS EoS utilizes a compressible liquid drop model for nuclei modeled after the Lamb et al. (1978) (LLPR) formalism, and considers an ion composition of free neutrons and protons, helium, and a representative heavy nucleus. For matter in NSE at densities below 1011 ${{\rm{g\; cm}}}^{-3}$ with ${Y}_{{\rm{e}}}\lt {}^{26}{/}_{56}$, the Cooperstein (1985) (C) EoS is used. The C EoS does not treat the high-density parameters of the liquid drop model (nuclear incompressibility modulus, surface energy, symmetry energy) as consistently as the LS EoS, but it computes the mass fraction of helium more accurately than the LS EoS in the regime below 1011 ${{\rm{g\; cm}}}^{-3}$ where it is being employed. Advection of material across this EoS boundary requires consistent tracking of the specific internal energy. (See Section 4.4.3 for details.) To improve the fidelity of the composition of matter that may eventually become part of the ejecta, in regions of NSE where ${Y}_{{\rm{e}}}\geqslant {}^{26}{/}_{56}$ (the value of Z/A for ${}^{56}\mathrm{Fe}$) the NSE calculation in C EoS has been upgraded to a 17-species representation of the composition, including free neutrons, free protons, the 14 even-Z and even-A nuclei between ${}^{4}\mathrm{He}$ and ${}^{60}\mathrm{Zn}$ plus ${}^{56}\mathrm{Fe}$. This NSE calculation functions similarly to others in the literature (e.g., Clifford & Tayler 1965; Hartmann et al. 1985) but is limited by the small isotope set.

3.3. Nuclear Network and Non-NSE Region

In zones where the timescale required to reach NSE is larger than other physical timescales (e.g., those associated with the evolution of stellar core fluid), the nuclear composition is evolved using the XNet thermonuclear reaction network code. In these regions, the thermodynamic state depends on the isotopic composition, as well as ρ and T, and the electron fraction (${Y}_{{\rm{e}}}$) is calculated from ${Y}_{e}={\sum }_{i}{Z}_{i}{Y}_{i}$, where Zi is the proton number of an isotope and Yi is the molar abundance of that isotope.

The initial value problem presented by a nuclear reaction network for an isolated region (individual zone) can, in principle, be solved by a wide range of methods discussed in the literature. However the physical nature of the problem, reflected in the wide range of reaction timescales, renders these numerical systems stiff. The challenges of solving such stiff astrophysical systems are detailed in a number of review articles on the subject (see, e.g., Hix & Meyer 2006; Travaglio & Hix 2013). In the A-, B- and C-series of Chimera models, XNet utilizes the fully implicit Backward-Euler method, introduced to nuclear astrophysics by Arnett & Truran (1969). Data for these reactions is drawn from the REACLIB compilation (Rauscher & Thielemann 2000). Unfortunately, the A- and B-series Chimera models neglect screening of nuclear reactions. The nuclear state is updated for each non-NSE zone in each time step. As needed in each zone, the network is sub-cycled until the hydrodynamic time step is reached.

Several pre-built networks available in Chimera are shown in Table 1. The A- and B-series models all utilize the simple 14-species $\alpha $-network alpha. The active nuclear material evolved in XNet excludes free protons, free neutrons, and an auxiliary heavy nucleus that are advected with the nuclear composition. In the C-series models, we switched to the alpnp network that adds protons and neutrons to the network, though the free nucleons are effectively inert, as their are no reactions included that connect them to other species. The properties (mass and charge number, binding energy, and mass fraction) of the auxiliary heavy nucleus are taken initially from the part of the composition of the progenitor that cannot be mapped onto the network species. For material that has come out of NSE, the properties are taken from the sum of nuclei not included in the network composition vector. For networks alpha and alpnp, this consists of ${}^{56}\mathrm{Fe}$ for ${Y}_{{\rm{e}}}\geqslant {}^{26}{/}_{56}$, or the representative heavy nucleus from the nuclear EoS for ${Y}_{{\rm{e}}}\leqslant {}^{26}{/}_{56}$. For the various D-series models underway, the base network has been updated to anp56, which adds ${}^{56}\mathrm{Fe}$ as an additional unconnected, inert species and permits the network to map directly to the 17-species NSE used by the extended C-EoS and also reduces the mass fraction traced by the auxiliary heavy nucleus. These modifications to the network infrastructure primarily serve the development of even larger networks (sn150 and sn160) for CCSN simulations. The D-series includes 2D and 3D simulations utilizing the sn160 network.

Table 1.  Available Nuclear Networks

Network Species
alpha ${}^{4}\mathrm{He}$, ${}^{12}{\rm{C}}$, ${}^{16}{\rm{O}}$, ${}^{20}\mathrm{Ne}$, ${}^{24}\mathrm{Mg}$, ${}^{28}\mathrm{Si}$, ${}^{32}{\rm{S}}$, ${}^{36}\mathrm{Ar}$, ${}^{40}\mathrm{Ca}$, ${}^{44}\mathrm{Ti}$, ${}^{48}\mathrm{Cr}$, ${}^{52}\mathrm{Fe}$, ${}^{56}\mathrm{Ni}$, ${}^{60}\mathrm{Zn}$
alpnp ${}^{}{\rm{n}}$ a, ${}^{1}{\rm{H}}$ a, ${}^{4}\mathrm{He}$, ${}^{12}{\rm{C}}$, ${}^{16}{\rm{O}}$, ${}^{20}\mathrm{Ne}$, ${}^{24}\mathrm{Mg}$, ${}^{28}\mathrm{Si}$, ${}^{32}{\rm{S}}$, ${}^{36}\mathrm{Ar}$, ${}^{40}\mathrm{Ca}$, ${}^{44}\mathrm{Ti}$, ${}^{48}\mathrm{Cr}$, ${}^{52}\mathrm{Fe}$, ${}^{56}\mathrm{Ni}$, ${}^{60}\mathrm{Zn}$
anp56 ${}^{}{\rm{n}}$ a, ${}^{1}{\rm{H}}$ a, ${}^{4}\mathrm{He}$, ${}^{12}{\rm{C}}$, ${}^{16}{\rm{O}}$, ${}^{20}\mathrm{Ne}$, ${}^{24}\mathrm{Mg}$, ${}^{28}\mathrm{Si}$, ${}^{32}{\rm{S}}$, ${}^{36}\mathrm{Ar}$, ${}^{40}\mathrm{Ca}$, ${}^{44}\mathrm{Ti}$, ${}^{48}\mathrm{Cr}$, ${}^{52}\mathrm{Fe}$, ${}^{56}\mathrm{Fe}$ a, ${}^{56}\mathrm{Ni}$, ${}^{60}\mathrm{Zn}$
 
  ${}^{}{\rm{n}}$, ${}^{1-2}{\rm{H}}$, ${}^{3-4}\mathrm{He}$, ${}^{6-7}\mathrm{Li}$, ${}^{\mathrm{7,9}}\mathrm{Be}$, ${}^{\mathrm{8,10,11}}{\rm{B}}$, ${}^{12-14}{\rm{C}}$, ${}^{13-15}{\rm{N}}$, ${}^{14-18}{\rm{O}}$, ${}^{17-19}{\rm{F}}$, ${}^{18-22}\mathrm{Ne}$, ${}^{21-23}\mathrm{Na}$,
sn150 ${}^{23-26}\mathrm{Mg}$, ${}^{24-27}\mathrm{Al}$, ${}^{28-32}\mathrm{Si}$, ${}^{29-33}{\rm{P}}$, ${}^{32-36}{\rm{S}}$, ${}^{33-37}\mathrm{Cl}$, ${}^{36-40}\mathrm{Ar}$, ${}^{37-41}{\rm{K}}$, ${}^{40-48}\mathrm{Ca}$, ${}^{43-49}\mathrm{Sc}$,
  ${}^{44-50}\mathrm{Ti}$, ${}^{46-51}{\rm{V}}$, ${}^{48-54}\mathrm{Cr}$, ${}^{50-55}\mathrm{Mn}$, ${}^{52-58}\mathrm{Fe}$, ${}^{53-59}\mathrm{Co}$, ${}^{56-62}\mathrm{Ni}$, ${}^{57-63}\mathrm{Cu}$, ${}^{59-66}\mathrm{Zn}$
 
  ${}^{}{\rm{n}}$, ${}^{1-2}{\rm{H}}$, ${}^{3-4}\mathrm{He}$, ${}^{6-7}\mathrm{Li}$, ${}^{\mathrm{7,9}}\mathrm{Be}$, ${}^{\mathrm{8,10,11}}{\rm{B}}$, ${}^{12-14}{\rm{C}}$, ${}^{13-15}{\rm{N}}$, ${}^{14-18}{\rm{O}}$, ${}^{17-19}{\rm{F}}$, ${}^{18-22}\mathrm{Ne}$, ${}^{21-23}\mathrm{Na}$,
sn160 ${}^{23-26}\mathrm{Mg}$, ${}^{25-27}\mathrm{Al}$, ${}^{28-32}\mathrm{Si}$, ${}^{29-33}{\rm{P}}$, ${}^{32-36}{\rm{S}}$, ${}^{33-37}\mathrm{Cl}$, ${}^{36-40}\mathrm{Ar}$, ${}^{37-41}{\rm{K}}$, ${}^{40-48}\mathrm{Ca}$, ${}^{43-49}\mathrm{Sc}$, ${}^{44-51}\mathrm{Ti}$,
  ${}^{46-52}{\rm{V}}$, ${}^{48-54}\mathrm{Cr}$, ${}^{50-55}\mathrm{Mn}$, ${}^{52-58}\mathrm{Fe}$, ${}^{53-59}\mathrm{Co}$, ${}^{56-64}\mathrm{Ni}$, ${}^{57-65}\mathrm{Cu}$, ${}^{59-66}\mathrm{Zn}$, ${}^{62-64}\mathrm{Ga}$, ${}^{63-64}\mathrm{Ge}$

Note.

aInert species, which are advected but not reactive.

Download table as:  ASCIITypeset image

When not in NSE, the thermodynamic quantities for the EoS cell are computed assuming the same composition of nuclei on all eight vertices. This gives rise to an apparent inconsistency, however, as the electron fraction, ${Y}_{{\rm{e}}}$, at some or all of the EoS grid points will not correspond to the ${Y}_{{\rm{e}}}$ of the composition of nuclei. Furthermore, the electron–positron contribution to the thermodynamic vector at each EoS grid point is computed from the values of ρ, T, and ${Y}_{{\rm{e}}}$ at that grid point, the result being that the ${Y}_{{\rm{e}}}$ of the electron–positron gas is not consistent with the ${Y}_{{\rm{e}}}$ of the composition of nuclei at a grid point. This procedure, however, allows us to take finite derivatives with respect to the electron fraction, albeit with some approximation, but accurate enough to stabilize the neutrino transport in the non-NSE regions. At the same time, when the thermodynamic state of a mass zone is interpolated from the EoS cell, the contribution of the ions will be based on a ${Y}_{{\rm{e}}}$ common to all EoS cell points, while the contributions of the electron–positron gas will be interpolated to the same ${Y}_{{\rm{e}}}$, and the two contributions will reflect the same value of ${Y}_{{\rm{e}}}$. The ions used to compute the thermodynamic properties are those in the network, the heavy nucleus advected with the active composition, and, in the case of the original alpha network used in Series A and B, the free nucleon mass fractions whose values are also stored with the thermodynamic vector.

3.4. NSE Transition

Chimera's treatment of the transition of matter into NSE is comparable to that used in other CCSN codes of similar capability (Buras et al. 2006a; Müller et al. 2012; Nakamura et al. 2014; Skinner et al. 2019). The transition condition is motivated by the temperatures and densities at which complete silicon burning would occur within the current global time step. For temperatures above this threshold, the use of the nuclear network is superfluous, as the network will achieve NSE every time step. For transitioning into NSE, there will be a slight change in the nuclear binding energy, due to the change of the representative nucleus, the ensemble of nuclei, etc., and the extent of this change is one metric for gauging the accuracy of our assumption of NSE. In order to maintain hydrodynamic stability across this transition, we adjust temperature to maintain constant pressure for a given density and electron fraction. This may result in a small change in the specific internal energy (including binding energy), but the dynamical impact that results from the transition between inconsistent thermodynamic states is minimized. The advection of material across an NSE/non-NSE interface in either direction, as well as the transition into NSE and freeze-out from NSE of entire zones, includes the appropriate gain or loss of nuclear binding energy (see Section 4.4.4 for details). In the A-series simulations, the NSE interface was a sphere of fixed radius, while in the B-series simulations, the NSE boundary was independent for each radial ray. The C- and later series allow multiple NSE/non-NSE interfaces along any coordinate direction.

To determine whether a zone is in NSE and may, therefore, be omitted from nuclear burning, Chimera applies an empirically determined linear relationship between the NSE transition temperature, ${T}_{\mathrm{NSE}}$, and the density:

Equation (7)

where C1 ≡ 5.333 ${\rm{K}}\,{{\rm{g}}}^{-1}\,{\mathrm{cm}}^{3}$ and ${C}_{2}\equiv 5.433\times {10}^{9}\,{\rm{K}}$. At the beginning of a global timestep, any non-NSE zone for which $T\geqslant {T}_{\mathrm{NSE}}$ is transitioned to NSE. A zone which is in NSE at the beginning of a timestep will be transitioned out of NSE if $T\lt ({T}_{\mathrm{NSE}}-2\times {10}^{8}\,{\rm{K}})$ and if the representative heavy nucleus, split into ${}^{56}\mathrm{Ni}$ and ${}^{56}\mathrm{Fe}$, will result in less than half of the mass fraction being ${}^{56}\mathrm{Fe}$. If this latter condition is not met, NSE is maintained as long as is practical, until the condition is met or until $T\lt 4.9\times {10}^{9}\,{\rm{K}}$, as the NSE representation of neutron-rich composition is better than that allowed by this limited network.

For simplicity, the transition out of NSE occurs when the temperature drops below this NSE condition (Equation (7) for Chimera). However, for the rapidly changing conditions in expanding CCSN matter, the assumption of NSE has been shown to break down when the temperature falls below 6 GK (Meyer et al. 1998). For transitioning out of NSE using the alpha network (but before evolving the network), the nuclear binding energy does not change. The NSE composition is evaluated using some analytically calculated state, but the temperature is then adjusted so that, for a given density and electron fraction, we will get the same specific internal energy (including binding energy) using a local EoS cube interpolation.

3.5. Electron–Positron EoS

The computation of the electron–positron component of the EoS is divided into five regimes (Figure 10) as described below. The first major division is based on whether the electrons are relativistic or nonrelativistic. The electron–positron gas is regarded as nonrelativistic if

Equation (8)

where the nonrelativistic electron Fermi energy, ${e}_{\mathrm{nr}\mathrm{Fermi}}$, is given by

Equation (9)

where ${n}_{{\rm{e}}}={n}_{{{\rm{e}}}^{-}}-{n}_{{{\rm{e}}}^{+}}$ is the net electron number density, ${n}_{{{\rm{e}}}^{-}}$ is the electron density, and ${n}_{{{\rm{e}}}^{+}}$ is the positron density. Otherwise, the electron–positron gas is treated as relativistic.

Figure 10.

Figure 10. Regimes in the Tρ plane, showing the different schemes for computing the electron–positron EoS quantities at ${Y}_{{\rm{e}}}=0.5$.

Standard image High-resolution image

If the electron–positron gas is relativistic, there are three regimes where different approximations are used: high temperature, very degenerate, and intermediate. In the intermediate regime, the Fermi integrals for the electron–positron thermodynamic functions are integrated numerically. Since this is the most computationally intensive procedure, it is first ascertained whether the thermodynamic state can be considered to be in the high-temperature or the very degenerate regime. If the thermodynamic state is found to be in neither of those regimes, the calculation of the thermodynamic functions begins with net electron density, given by

Equation (10)

where in the last expression $\beta ={m}_{{\rm{e}}}{c}^{2}/{kT}$, ${\eta }_{{\rm{e}}}={\mu }_{{\rm{e}}}/{kT}$, and the following substitutions have been made: $y=p/{m}_{{\rm{e}}}c$, ${z}^{2}={y}^{2}+1$, and $x=\beta (z-1)$. To obtain ${\eta }_{{\rm{e}}}$, the right-hand side of Equation (10) is integrated numerically by means of a 48-point Gauss–Laguerre scheme using a guess for ηe, and then iterated for the ηe until the right-hand side equals ne to within one part in 106. Once ηe is obtained, ${p}_{{\rm{e}}+{\rm{p}}}$ and ${e}_{{\rm{e}}+{\rm{p}}}$, the latter of which includes the electron–positron rest mass energy, are obtained by 48-point Gauss–Laguerre numerical integration of the following Fermi integrals

Equation (11)

Equation (12)

The entropy is obtained from the pressure, specific internal energy, and chemical potential from the thermodynamic relation

Equation (13)

To determine whether the high-temperature approximation may be applied, Equation (10) is rewritten with the substitution z = x + β to get

Equation (14)

The high-temperature approximation consists of setting β = 0 in Equation (14) and performing the integrations analytically

Equation (15)

with the solution for ηe

Equation (16)

where $\chi ={\pi }^{2}{\left({\hslash }c/{kT}\right)}^{3}{n}_{{\rm{e}}}$. The high-temperature approximation is then applied if β < 2/3, or if β < 2 and ηe given by Equation (16) satisfies ηe > 2β. The high-temperature approximation is applied if values of ρ and T lie above and to the right of the solid red line in Figure 10. The solid red horizontal line at the left is given by the first condition above, the rest of the red line is given by the second condition. The solid green line terminating in the curved red segment is given by the condition ηe > 2β. Once ηe is obtained, analytic expressions for ${p}_{{\rm{e}}+{\rm{p}}}$ and ${e}_{{\rm{e}}+{\rm{p}}}$ in the high-temperature approximation are given by

Equation (17)

Equation (18)

so that

Equation (19)

where the F's are the usual Fermi integrals. Substituting Equation (19) into Equation (13) for the entropy gives

Equation (20)

After computing the pressure from Equation (17), the entropy is computed from Equation (20), and then the specific internal energy is computed from the thermodynamic relation

Equation (21)

The thermodynamic state is considered to be in the very degenerate regime if, as computed by Equation (16) or determined by iterating Equation (10), ηe > 35. In this case, the relativistic Sommerfeld approximation is used for the thermodynamic functions, starting with

Equation (22)

which is iterated on ${x}_{\eta }={pc}/{m}_{{\rm{e}}}{c}^{2}$, where p is the electron momentum, until the right-hand side is equal to ${n}_{{{\rm{e}}}^{-}}$ to one part in 106. Once xη is determined, the other thermodynamic functions are computed as follows:

Equation (23)

Equation (24)

Equation (25)

and

Equation (26)

If the nonrelativistic criteria in Equations (8) are satisfied, then there are two regimes in which different approximations are applied. If ${\epsilon }_{{\rm{F}}}/{kT}\gt 35$, the nonrelativistic Sommerfeld approximation is used, otherwise the nonrelativistic Fermi integrals for the electrons are integrated numerically. In the latter case, the Fermi integral for the electron number density

Equation (27)

is integrated by a 48-point Gauss–Laguerre scheme and iterated on ηe until the right and left sides of Equation (27) are equal to one part in 106. Once ηe has been obtained, μe is computed by ${\mu }_{{\rm{e}}}={kT}{\eta }_{{\rm{e}}}+{m}_{{\rm{e}}}{c}^{2}$, and the electron pressure, specific internal energy, and entropy are obtained by

Equation (28)

and

Equation (29)

In the case that ${\epsilon }_{{\rm{F}}}/{kT}\gt 35$, the Sommerfeld approximations for the nonrelativistic Fermi expressions are used. Thus, the electron density is given by

Equation (30)

Equation (30) is iterated on ηe as described above, μe is then computed also as described above, and the electron pressure, specific internal energy, and entropy are obtained by

Equation (31)

Equation (32)

The accuracies of the various approximations are indicated by the relative deviations of the electron–positron pressure and specific internal energy as computed by pairs of these approximations at the boundaries separating their respective regimes of applicability. This is shown in Figure 11, which indicates that relative deviations are at most a few tenths of a percent. The discontinuity between approximation regimes is smoothed by the finite resolution of the EoS table.

Figure 11.

Figure 11. Panel (a): for the T profile of the Gauss–Laguerre—Sommerfeld boundary (blue line; right axis), the normalized deviation ($\parallel {X}_{{\rm{A}}}-{X}_{{\rm{B}}}\parallel /{X}_{{\rm{B}}};$ left axis) for electron pressure (${p}_{{\rm{e}}};$ black solid line) and electron specific kinetic energy (${e}_{{\rm{e}},\mathrm{kin}};$ red line) between the Gauss–Laguerre integration and Sommerfeld approximation. The short segments on the right side show the relative deviation of the electron–positron pressure (${p}_{{{\rm{e}}}^{-}+{{\rm{e}}}^{+}};$ orange) and the electron–positron specific kinetic energy (${e}_{{{\rm{e}}}^{-}+{{\rm{e}}}^{+},\mathrm{kin}};$ green) as computed with the relativistic and nonrelativistic formalisms at the relativistic–nonrelativistic transition density. Panel (b): for the T profile (blue) at the boundaries of (i) the high-temperature approximation and Gauss–Laguerre integration regions (solid lines) and (ii) the Relativistic Sommerfeld approximation and Gauss–Laguerre integration regions (dashed lines), we plot the normalized deviation of electron–positron pressure (${p}_{{{\rm{e}}}^{-}+{{\rm{e}}}^{+}};$ black) and specific internal energy (${e}_{\mathrm{int},{{\rm{e}}}^{-}+{{\rm{e}}}^{+}};$ red).

Standard image High-resolution image

3.6. Double-γ EoS

A "2γ," thermodynamically consistent EoS has been developed, simple enough to be completely analytic, yet rich enough to be used for testing the hydrodynamics and transport modules of Chimera. A system of completely degenerate and relativistic free electrons, partially degenerate interacting neutrons, and nondegenerate free protons, is modeled by the free energy

Equation (33)

where mB is the baryon mass, gn = gp = 2 are the statistical weights of the neutrons and protons, k and h have their usual meanings, and γ1, γ2, F1, F2, Ecoef, and ${X}_{{\rm{p}}}\equiv {N}_{{\rm{p}}}$/$(\rho /{m}_{{\rm{B}}})\,=$ ${Y}_{{\rm{e}}}\equiv ({N}_{{{\rm{e}}}^{-}}-{N}_{{{\rm{e}}}^{+}})/(\rho /{m}_{{\rm{B}}})=1-{X}_{{\rm{n}}}\,=$ ${N}_{{\rm{n}}}/(\rho /{m}_{{\rm{B}}})$ are free parameters. The mass fraction of free protons, Xp, is typically taken to be 0.5, and the parameters γ1 and F1 are typically chosen to be

Equation (34)

Equation (35)

so that the first term on the right-had side of Equation (33) represents completely degenerate and extremely relativistic free electrons. Given the above expression for the free energy, the pressure is

Equation (36)

We choose F2 so that the contributions of the first two terms for p in Equation (36) become equal at ρnuc, which, for the case in which Yp is 0.5, requires that

Equation (37)

or

Equation (38)

The entropy of the system is given by

Equation (39)

where ${S}_{{\rm{i}}}$, ${\rm{i}}={\rm{n}},{\rm{p}}$ is given by

Equation (40)

The internal energy of the system is given by

Equation (41)

Lastly, the chemical potentials are given by

Equation (42)

Equation (43)

and

Equation (44)

That all of the thermodynamic functions are derived from an analytic thermodynamic potential, namely, the free energy, insures that the EoS is thermodynamically consistent. Furthermore, being completely analytic and simple, the interpolation scheme described above is not necessary when using this EoS. Derivatives of thermodynamic functions can be obtained analytically, and expressions for the independent variables, ρ, T, and ${Y}_{{\rm{e}}}$ can be solved for directly without having to invert a complicated expression by resorting to an iteration scheme.

4. Hydrodynamics

Chimera's hydrodynamics are evolved using a dimensionally split, Lagrangian-plus-remap version of the Piecewise Parabolic Method (PPMLR; Colella & Woodward 1984) as implemented in VH1 (Hawley et al. 2012) but extended to include multi-species advection, multidimensional gravity, neutrino coupling (energy and momentum) to the hydrodynamics, and radial grid movement. PPM is a high-order Godunov-type scheme (Godunov 1959; van Leer 1979) in which fluxes at zone interfaces are calculated from solutions of the Riemann problem and, therefore, are well suited for capturing shocks and contact discontinuities within one or two zones. The PPMLR evolves the zone averages of the density, ρ, fluid velocity, ${\boldsymbol{u}}$, specific internal energy, ${e}_{\mathrm{int}}$, electron fraction, ${Y}_{{\rm{e}}}$, mass fractions, Xn, of nuclear ion species in regions that are not in NSE (non-NSE regions), and zero-angular moments of the neutrino distribution functions ${\psi }_{\nu }^{(0)}$.

There are several advantages to a Lagrangian-plus-remap for CCSN simulations. For explicit differencing, the time step constraints are less severe, as they are only applied during the Lagrangian step and depend only on the Lagrangian wave speeds rather than the sum of the Lagrangian wave plus advection speeds. Furthermore, in, for example, a quasi-Eulerian approach, the remap does not have to map the grid back to its placement prior to the Lagrangian step but can allow the grid to evolve to accommodate itself to changing physical situations, such as moving with the fluid during core infall, thereby keeping the fluid well resolved or ensuring adequate zoning in the vicinity of the neutrinosphere during the formation of the density cliff or tracking the interface between two compositionally different fluids.

While the hydrodynamics can be performed in either Cartesian, cylindrical or spherical coordinates, the neutrino radiation transport is performed in the "RbR" approximation (Section 6) and requires a spherical coordinate system. We will therefore limit our discussion to spherical coordinate systems, as these are the coordinate systems that have been used for all Chimera CCSN simulations.

4.1. General Overview

The method of solution of the hydrodynamics equations is a finite-volume method, wherein conserved quantities are represented as integrated over computational volumes, or zones, and changes to these variables occur by means of sources or sinks or by means of fluxes through zone boundaries due to the relative velocity between the fluid and these boundaries. Chimera uses a Lagrangian-plus-remap version of the PPM method, which can be described by considering the zone-integrated conserved quantity, $q({\boldsymbol{x}},t)$, integrated over a computational element, ΔV(t), whose boundaries may be time-dependent

Equation (45)

The time rate of change of Q is

Equation (46)

where the second line results from the use of the Reynolds transport theorem, with ${{\boldsymbol{u}}}_{g}$, which we will refer to as the grid velocity, being the velocity of the bounding surface ${\rm{\Delta }}S(t)$ of the zone volume ${\rm{\Delta }}V(t)$. The first term in the third line is the flux of $q({\boldsymbol{x}},t)$ across the surface ΔS(t) due to the fluid velocity ${\boldsymbol{u}}$, assuming linear flux and the surface ΔS(t) to be fixed, and the second term represents a possible volume source ${{ \mathcal S }}_{q}({\boldsymbol{x}},t)$ of $q({\boldsymbol{x}},t)$. To calculate ${dQ}/{dt}$, Chimera splits the calculation into two steps,

Equation (47)

In "step 1," the grid velocity is set equal to the fluid velocity and the resulting equation

Equation (48)

is solved. This is just the Lagrangian form of the equation for $\partial Q({\boldsymbol{x}},t)/\partial t$. In "step 2," the fluid velocity is set to zero and the initial configuration of ${{ \mathcal S }}_{q}({\boldsymbol{x}},t)$ is that of its final configuration after "step 1." The grid is then given a prescribed velocity so that

Equation (49)

If the grid velocity, ${{\boldsymbol{u}}}_{g}$, is chosen to be the negative of the fluid velocity, ${\boldsymbol{u}}$, then "step 1" plus "step 2" would be equivalent to

Equation (50)

which is just the Eulerian form of the equation for $\partial Q({\boldsymbol{x}},t)/\partial t$. For the θ- and ϕ-grid, we choose to set ${{\boldsymbol{u}}}_{g}=-{\boldsymbol{u}}$ to keep these grids stationary. For the radial grid, however, we use the freedom of choice for ${{\boldsymbol{u}}}_{g}$ to make the grid dynamically adaptive, allowing it to move in such a way as to maintain good resolution during such epochs as core collapse or the formation of a density cliff in the vicinity of the neutrinospheres.

4.2. PPM Interpolation Scheme

The solution of the hydrodynamics equations proceeds in a dimensionally split manner. We will describe the solution as it proceeds in a particular, but otherwise arbitrary, dimension and will refer to a specific dimension (e.g., r, θ, or ϕ) only when expressions specific to this dimension arise. In order to construct the finite difference representations of the underlying partial differential equations that Chimera solves, a discrete grid is set up dividing the interval ξmin to ξmax into a total of I zones, where ξ is the parameter (r, θ, or ϕ) specifying the coordinate distance along a ray in a given dimension. Figure 12 illustrates our indexing convection. At each end of the grid, six ghost zones are appended to hold boundary values of the quantities stored in the real zones 1... I. In a given dimension, both the Lagrangian and the remap steps in PPMLR hydrodynamics begin by constructing zone interface values of primitive quantities, such as ρ, p, and the components of ${\boldsymbol{u}}$, from zone-average values of these quantities. Lufkin & Hawley (1993) have shown that a differencing scheme that uses zone averages to represent fluid variables will not converge to the continuity equation (essential for conserving quantities during advection) unless the interpolation scheme, a(ξ), for the zone-averaged quantities, ai, satisfies

Equation (51)

where Vi is the volume of zone i, and ${dV}/d\xi $ is the one-dimensional Jacobian determinant for V.

The interpolation scheme used here is the same as that described in Colella & Woodward (1984), as modified by Blondin & Lufkin (1993) to accommodate curvilinear coordinates while at the same time representing a linear velocity field accurately near coordinate boundaries—e.g., r = 0 in the radial dimension. The procedure for determining the zone interface value at i − $\tfrac{1}{2}$ of the zone-averaged variables, ai, is to construct a quartic polynomial, $A(\xi )$, for each zone interface i − $\tfrac{1}{2}$ such that it takes on the respective values, ${A}_{i-\tfrac{5}{2}}\cdots {A}_{i+\tfrac{3}{2}}$ at the five points ${\xi }_{i-\tfrac{5}{2}}\cdots {\xi }_{i+\tfrac{3}{2}}$, where

Equation (52)

The desired interpolating polynomial, a(ξ), is given by the integrand of the indefinite integral

Equation (53)

that is,

Equation (54)

By construction, each cubic polynomial a(ξ) so obtained has the desired property

Equation (55)

The interface value, ${a}_{i-\tfrac{1}{2}}$, is obtained from Equation (54) by evaluating it at $\xi ={\xi }_{i-\tfrac{1}{2}}$:

Equation (56)

where the explicit expression for A(ξ) is given by Equations (12) and (13) in Blondin & Lufkin (1993).

The interpolation for ${a}_{i-\tfrac{1}{2}}$ given by Blondin & Lufkin (1993) differs from Equations (1.6) and (1.7) of Colella & Woodward (1984). In the former publication, the zone-averaged quantities ai are multiplied by the geometry-dependent correction factors ${\rm{\Delta }}{V}_{i}/{\rm{\Delta }}{\xi }_{i}$. Denoting these geometry corrected quantities by ${a}_{i}^{* }$, we have

Equation (57)

where the geometry-dependent correction factor ${\rm{\Delta }}{V}_{i}/{\rm{\Delta }}{\xi }_{i}$ is given by

Equation (58)

and where

Equation (59)

The average slope, $\delta {a}_{i}^{* }$, of the parabolas are then computed from Equation (1.7) of Colella & Woodward (1984) but using the quantities ${a}_{i}^{* };$ e.g.,

Equation (60)

where

Equation (61)

The $\delta {a}_{i}^{* }$ are then modified as follows (see; Colella & Woodward 1984, Equation (1.8)):

Equation (62)

The modifications for both cases implement monotonicity constraints, ensuring that no new maxima or minima appear (i.e., that ${a}_{i-\tfrac{1}{2}}$ lies in the range of ${a}_{i}^{* }$ and ${a}_{i-1}^{* }$), and in the case of ${\rm{\Delta }}{a}_{i+\tfrac{1}{2}}^{* }{\rm{\Delta }}{a}_{i-\tfrac{1}{2}}^{* }\gt 0$ lead to somewhat steeper representations of discontinuities.

Given ${a}_{i}^{* }$ and ${\delta }_{m}{a}_{i}^{* }$, the interface value ${a}_{i+\tfrac{1}{2}}$ is now obtained from the cubic interpolating polynomial, Equation (54), which is Equation (1.6) of Colella & Woodward (1984) with the geometry-dependent corrections applied as above in Equation (57) and below in Equation (64):

Equation (63)

Equation (64)

where

Equation (65)

where the values at ${\xi }_{\tfrac{1}{2}}=0$ are to avoid singularities in Equation (64). Finally, the range of ${a}_{R,i}$ is limited to be within the range of ai and ${a}_{i+1}$:

Equation (66)

Equation (67)

and

Equation (68)

where ${a}_{L,i+1}$ is the value of a(ξ) at the left interface of zone i + 1, and ${a}_{R,i}$ is its value at the right interface of zone i. At zone boundaries, ${a}_{L,i}$ and ${a}_{R,i}$ must be modified in certain cases as follows. In the r-direction, if ${\xi }_{\tfrac{1}{2}}=0$, then ${a}_{L,i}$ must be modified in a manner depending on whether the variable is odd (e.g., velocity) or even (e.g., density, specific energy) at the origin:

Equation (69)

In the θ-direction, if ${\xi }_{\tfrac{1}{2}}=0(\theta =0)$ and reflecting boundary conditions are imposed,

Equation (70)

and if ${\xi }_{I+\tfrac{1}{2}}=0(\theta =\pi )$,

Equation (71)

In the presence of shocks, post-shock oscillations sometimes occur in some of the fluid variables, e.g., entropy (Colella & Woodward 1984, Section 4). One method of suppressing these oscillations is to introduce some additional dissipation in the vicinity of a shock. The method used here is to lower the order of the interpolation (i.e., flatten the interpolation profile) in the vicinity of a shock. Thus, in the vicinity of a shock, aL,I and aR,I are modified as follows:

Equation (72)

where the "flattening parameter," 0 ≤ f ≤ 1, is zero away from shocks and approaches a user preset value fset ≤ 1 in the limit of a strong shock with a steep profile.

The flattening parameter fi is computed similarly to that described in Fryxell et al. (2000) and begins with the computation of the pressure profile ${p}_{\mathrm{profile}i}$ given by

Equation (73)

The presence of a shock is indicated if this quantity approaches 1, i.e., when ${p}_{i+2}-{p}_{i-2}\approx {p}_{i+1}-{p}_{i-1}$. For smooth flows, ${p}_{i+2}-{p}_{i-2}\approx 2\times ({p}_{i+1}-{p}_{i-1})$ and ${p}_{\mathrm{profile}i}\approx 1/2$. A pre-flattening parameter ${f}_{\mathrm{pre}i}$ is then computed from

Equation (74)

where ω1 and ω1 are user selected parameters. We have found that 0.75 and 10 work well for ω1 and ω2 during the early stages of a simulation when the shock is nearly spherically symmetric, completely eliminating post-shock oscillations, and 0.6 and 10 work well thereafter better capturing oblique shocks. The pre-flattening parameter is set to zero if

Equation (75)

which indicates an insufficient pressure jump, and if

Equation (76)

which indicates a velocity jump in the wrong direction. Finally, the flattening parameter is given by

Equation (77)

This limits the value of fi to $0\leqslant {f}_{\mathrm{set}}$ and sets the value fi based on the maximum value of ${f}_{\mathrm{pre}}$ in zone i and the neighboring zones. Experimentation has indicated that fset = 1 works well when the shock is approximately spherically symmetric, and fset = 0.5 works well thereafter.

With the values of ${a}_{L,i}$ and ${a}_{R,i}$ for each zone i determined, a piecewise parabolic interpolation function, a(ξ), is constructed with a(ξ) given by a parabolic profile in each zone:

Equation (78)

where

Equation (79)

and where ${a}_{6,i}$. The parameter ${a}_{6,i}$ must now be determined so that Equation (51) is satisfied. The expressions for ${a}_{6,i}$ are given by

Equation (80)

where

Equation (81)

Equation (82)

Equation (83)

and

Equation (84)

Equation (85)

Equation (86)

In addition to the constraints imposed by Equations (62) and (66)–(68), the following monotonicity constraints are imposed to avoid the possibility of the interpolating function taking on values not between ${a}_{L,i}$ and ${a}_{R,i}$, which could otherwise lead to its developing spurious oscillations:

Equation (87)

Equation (88)

where

Equation (89)

The need to use the monotonized expression for the left and right states arises if the parabola exceeds either of these states. With Δai given by Equation (79), we thus have

Equation (90)

Equation (91)

With these now monotonized values of ${a}_{L,i}$ and ${a}_{R,i}$, the quantities ${\rm{\Delta }}{a}_{i}$, and ${a}_{6,i}$ are recomputed from Equations (79) and (80). This completes the piecewise parabolic interpolation for the profile, a(ξ), of zone-averaged variables ai.

4.3. Lagrangian Step

The equations describing the change in the hydrodynamic variables during the Lagrangian step are:

Equation (92)

Equation (93)

Equation (94)

Equation (95)

Equation (96)

where V(t) and S(t) are the volume and surface of a grid element or mass zone as it comoves with the fluid, $\hat{{\boldsymbol{n}}}$ is a unit vector normal to dS and pointing out of the mass zone, ${e}_{\mathrm{grav}}$ is the gravitational potential, ${{\boldsymbol{f}}}_{\nu }$ is the specific neutrino stress, and ${{\boldsymbol{f}}}_{\mathrm{ccor}}$ denotes the centrifugal and Coriolis forces. The time derivatives are taken at constant mass (i.e., are Lagrangian time derivatives). The expressions in the brackets on the right-hand sides of Equations (95) and (96) denote the changes in the electron fraction, ${Y}_{{\rm{e}}}$, due to neutrino transport and in the composition mass fractions, Xn, due to both neutrino transport and nuclear reactions, and are calculated elsewhere in the computational sweep, as described in Sections 6.8 and 3.3. During the Lagrangian hydrodynamics step, these expressions are set to zero.

Equation (94) is a common formulation of energy conservation on which difference schemes are subsequently constructed (e.g., Stone & Norman 1992; Bryan et al. 1995; Fryxell et al. 2000; Sutherland 2010). Chimera uses two alternative formulations of energy conservation, depending on the circumstances. Using Equations (92) and (93) in Equation (94), the following expression for ${e}_{\mathrm{int}}$ can be derived:

Equation (97)

which is just the first law in the absence of changes in entropy and electron fraction. This equation is applicable for updating the specific internal energy during the Lagrangian step in regions well away from shocks. It is more accurate numerically, in some cases, than updating the specific total energy as it does not ultimately involve the subtraction of potentially large values of the specific kinetic energy and the gravitational potential (if the latter is also included as a component in the specific total energy). In the vicinity of shocks, however, the specific total energy must be updated from the solution of the Riemann shock tube problem for the pressure and velocity at the zone interfaces. To do this, the term involving $\rho {\boldsymbol{u}}\cdot {\rm{\nabla }}{e}_{\mathrm{grav}}$ on the right-hand side of the energy Equation (94) is transformed as follows:

Equation (98)

Substituting Equation (98) in Equation (94) results in the following equation for the total energy (including gravitational energy)

Equation (99)

Equation (99) is used in the radial-sweep to update the specific total energy in the vicinity of a shock. For the θ- and ϕ-sweeps, changes in the gravitational potential are very small, and Equation (94) is used with the gradient of the gravitational potential on the right-hand side treated as a force.

To perform the Lagrangian update, the first step is to compute the displacement of the zone interfaces, after which Equations (92), (93), (95), (96), and (97) or (99) are used to update ρ, the components of ${\boldsymbol{u}}$, ${Y}_{{\rm{e}}}$, and ${e}_{\mathrm{int}}$. The displacement of each zone interface during the Lagrangian step is determined by solving a Riemann problem for the velocity of the contact discontinuity at the zone interface. This requires averages of the needed quantities over the domains of dependences of the left and right states. Rather than solving the exact Riemann problem, which is time consuming, as it is complicated and involves multiple calls to the EoS, Chimera uses the approximate but very accurate method developed by Colella & Glaz (1985). This method parameterizes the EoS by the slowly varying quantity γ, given by

Equation (100)

and the adiabatic exponent, Γ, defined by

Equation (101)

A solution for the approximate Riemann problem requires the values of the quantities ρ, u (the component of velocity in the direction of the directional splitting), p, γ, and Γ to the left and right of each zone interface. To maintain high-order accuracy, the values of each of these quantities are averaged over their domain of dependence of the zone interface as determined by the time step. This is accomplished for each zone interface by tracing the two characteristics from the interface at time t + Δt to the ξ axis at time t. Having speeds of ±cs, where cs is the local sound speed, the two characteristics intersect the ξ axis on either side of the interface at the points ${\xi }_{i+\tfrac{1}{2}}+{c}_{{\rm{s}}}{\rm{\Delta }}t$ and ${\xi }_{i+\tfrac{1}{2}}-{c}_{{\rm{s}}}{\rm{\Delta }}t$. Letting ai represent any one of the above quantities, the average, $\langle a{\rangle }_{L,i+\tfrac{1}{2}}^{n+\tfrac{1}{2}}$, of ai over the domain of dependence to the left of the zone interface, ${\xi }_{i+\tfrac{1}{2}}$, is obtained by integrating the parabolic profile $a(\xi )$ of ai over the interval ${\xi }_{i+\tfrac{1}{2}}-{c}_{{\rm{s}},i}{\rm{\Delta }}t$ to ${\xi }_{i+\tfrac{1}{2}}$ and averaging, and is given by

Equation (102)

where Δai is given by Equation (79), Δξi is given by Equation (59), and ${c}_{{\rm{s}},i}^{2}={{\rm{\Gamma }}}_{i}{p}_{i}/{\rho }_{i}$. Likewise, the average, $\langle a{\rangle }_{R,i+\tfrac{1}{2}}^{n+\tfrac{1}{2}}$, of a over the domain of dependence to the right of the zone interface, ${\xi }_{i+\tfrac{1}{2}}$, is obtained by integrating the parabolic profile a(ξ) of ai over the interval ${\xi }_{i+\tfrac{1}{2}}$ to ${\xi }_{i+\tfrac{1}{2}}+{c}_{{\rm{s}},i+1}{\rm{\Delta }}t$ and averaging, and is given by

Equation (103)

The time-averaged left and right states, $\langle {p}^{+}{\rangle }_{L,i+\tfrac{1}{2}}$, $\langle {p}^{+}{\rangle }_{R,i+\tfrac{1}{2}}$, $\langle u{\rangle }_{L,i+\tfrac{1}{2}}$, $\langle u{\rangle }_{R,i+\tfrac{1}{2}}$ of p and u are obtained from Equations (102) and (103) above, and the time-averaged pressure is then corrected for the presence of gravitational, neutrino, centrifugal, and Coriolis forces by

Equation (104)

Equation (105)

Given the time-averaged states of p and u to the left and right of each zone interface, ${\xi }_{i+\tfrac{1}{2}}$, the time-averaged values, ${p}_{i+\tfrac{1}{2}}^{n+\tfrac{1}{2}}$ and ${u}_{i+\tfrac{1}{2}}^{n+\tfrac{1}{2}}$, of the pressure and velocity of the zone interface itself are computed by connecting ${p}_{i+\tfrac{1}{2}}^{n+\tfrac{1}{2}}$ and ${u}_{i+\tfrac{1}{2}}^{n+\tfrac{1}{2}}$ to the time-averaged left and right states by the Rankine–Hugoniot relations; that is:

Equation (106)

Equation (107)

Equation (108)

Equation (109)

Equations (106)–(109) are iterated for ${p}_{i+\tfrac{1}{2}}^{n+\tfrac{1}{2}}$ and ${u}_{i+\tfrac{1}{2}}^{n+\tfrac{1}{2}}$ by the secant method.

Having determined ${p}_{i+\tfrac{1}{2}}^{n+\tfrac{1}{2}}$ and ${u}_{i+\tfrac{1}{2}}^{n+\tfrac{1}{2}}$ for each of the zone interfaces, the Lagrangian update proceeds as follows. With the values of ${u}_{i+\tfrac{1}{2}}^{n+\tfrac{1}{2}}$ determined, the zone interfaces are considered impenetrable and their positions are updated by

Equation (110)

where the superscripts n and $n+1^{\prime} $ denote the value of a variable at time t and at the end of the Lagrangian step at time $t+{\rm{\Delta }}t$, respectively. A superscript n + $\tfrac{1}{2}$ denotes a time-centered value of a variable. We reserve the superscript $n+1$ for the value of a variable after both the Lagrangian and the remap step have been completed. From Equation (92), the density is then updated by

Equation (111)

Because of the conservation of mass in each zone during the Lagrangian step, as expressed by Equation (92), and in differenced form by Equation (111), Equation (95) for the change in ${Y}_{{\rm{e}}}$ and Equation (96) for the change in the composition mass fractions, Xi, with their right-hand sides set to zero, state the obvious: ${Y}_{{\rm{e}}}$ and Xi are unchanged during the Lagrangian step.

Equation (93) with the help of Equation (92) can be rewritten in differential form as

Equation (112)

where, again, the time derivative is a Lagrangian time derivative. In component form, Equation (112) becomes

Equation (113)

For the radial sweep, the neutrino stress term, ${f}_{\nu ,r}$, is computed as described by the right-most term on the right-hand side of Equation (250). Because of the RbR approximation adopted by Chimera for neutrino transport, the values for ${f}_{\nu ,\theta }$ and ${f}_{\nu ,\phi }$ appearing in the θ- and ϕ-sweeps cannot be obtained directly as an outcome of the transport as can ${f}_{\nu ,r}$ for the radial sweep. To include ${f}_{\nu ,\theta }$ and ${f}_{\nu ,\phi }$ in an approximate way, we regard the matter as being completely neutrino-opaque at densities above 1012 ${{\rm{g\; cm}}}^{-3}$ and completely transparent at lower densities. The neutrino distribution in each zone is thus assumed to behave like an isotropic, completely relativistic gas at densities above 1012 ${{\rm{g\; cm}}}^{-3}$, whose effect on the hydrodynamics is computed by means of their corresponding pressure, ${p}_{\nu }$, and specific energy, ${e}_{\nu }$. For densities below 1012 ${{\rm{g\; cm}}}^{-3}$, ${p}_{\nu }$ and ${e}_{\nu }$ are set to zero. The neutrino stress for the θ- and ϕ-sweeps is either that of an isotropic gas entrained with the matter (above 1012 ${{\rm{g\; cm}}}^{-3}$) or zero (below 1012 ${{\rm{g\; cm}}}^{-3}$). As a result, we set ${f}_{\nu ,\theta }$ and ${f}_{\nu ,\phi }$ to zero and include ${p}_{\nu }$ (nonzero above 1012 ${{\rm{g\; cm}}}^{-3}$) in the material pressure, and the sum is PPM interpolated and incorporated into ${p}_{i+\tfrac{1}{2}}^{n+\tfrac{1}{2}}$ and is then used to update the specific internal energy, with ${e}_{\nu ,i}$ incorporated into and later extracted from ${e}_{\mathrm{int},i}$.

The finite difference approximations to Equations (113) are

Equation (114)

Equation (115)

where ${A}_{i+\tfrac{1}{2}}^{n+\tfrac{1}{2}}$ is defined by $\delta {V}_{i+\tfrac{1}{2}}={A}_{i+\tfrac{1}{2}}^{n+\tfrac{1}{2}}{u}_{i+\tfrac{1}{2}}^{n+\tfrac{1}{2}}{\rm{\Delta }}t$. $\delta {V}_{i+\tfrac{1}{2}}$ is the volume swept out by the change in the position of the ith + 1 interface in the time interval ${\rm{\Delta }}t$ and is given by

Equation (116)

${f}_{\mathrm{cenfugal},i}^{n+\tfrac{1}{2}}$ includes the last term in Equations (113) for the r-direction sweep and the second-to-last term in Equations (113) for the θ-direction sweep. The last term in Equation (113) for the θ-direction sweep is an expression for the change in ${u}_{\theta }$ due to a change in the radial position by virtue of angular momentum conservation. Rather than including this term in Equation (113), its contribution to the change in uθ is updated in the r-sweep by the equivalent expression

Equation (117)

Similarly, the last two terms in Equation (113) for the ϕ-direction sweeps are expressions for angular momentum conservation about the z-axis due to a change in the θ and radial positions, respectively. Rather than including these terms in Equation (113), their contributions to the change in uϕ are updated in the θ- and r-direction sweeps by the equivalent expressions

Equation (118)

The time-centered value of the gravitational potential, ${\phi }_{{\rm{g}}}$, is obtained by extrapolation, as described in Section 4.6. The two centrifugal force terms are differenced as follows:

Equation (119)

Equation (120)

where the time-centering of ri in the r-direction and ${\theta }_{i}$ in the θ-direction is performed explicitly by averaging their values at time n and $n+1$, and the time-centering of the other variables is approximately accomplished by the symmetric way in which the directional splitting is performed, as described in Section 2.1. Finally, the neutrino stress term ${f}_{\nu ,i}^{n}$ is not time centered. A second execution of the neutrino transport would be required to center it. This term is small and slowly varying in time, so we include it to the first-order only and avoid the significant additional cost that would be paid were we to include it with second-order accuracy.

The specific internal energy is updated differently depending on whether the zone is in the vicinity of a shock or away from shocks. In the vicinity of a shock, the results of solving the Rankine–Hugoniot equations must be used, as the flow there is non-isentropic. In this case, the specific total energy is updated as given, in general form, by Equation (99) for the radial sweep, and Equation (94) for the θ- and ϕ-sweeps. Specifically, for the radial-sweep, the specific total energy, ${e}_{\mathrm{tot}}={e}_{\mathrm{int}}+{e}_{\mathrm{kin}}+{e}_{\mathrm{grav}}$, is updated by

Equation (121)

where ${u}_{i}^{n+1^{\prime} }$ is given by Equation (114). The $\partial {e}_{\mathrm{grav}}/\partial t$ term in Equation (99) has been omitted at this stage but is included later in the radial sweep. For the θ- and ϕ-sweeps, the specific energy, ${e}_{\mathrm{tot}^{\prime} }={e}_{\mathrm{int}}+{e}_{\mathrm{kin}}$, is updated by

Equation (122)

followed by

Equation (123)

Away from shocks, the specific energy could still be updated as above, but errors might then arise during the remap step when subtracting the specific kinetic energy from the specific total energy. The problem arises from the use of $\langle {\boldsymbol{u}}{\rangle }^{2}$ in the expression for the specific kinetic energy rather than $\langle {{\boldsymbol{u}}}^{2}\rangle $. The two expressions can differ importantly in supersonic flow and near reflection boundaries where the gradient of ${\boldsymbol{u}}$ can be large (see Blondin & Lufkin 1993, for a discussion of this point). Computing $\langle {{\boldsymbol{u}}}^{2}\rangle $ would be costly in a multidimensional simulation, as it would involve a multidimensional integration over the components of ${\boldsymbol{u}}$ and may not be well defined. Instead, Chimera updates the specific internal energy using the first law of thermodynamics, assuming isentropic and constant-composition flow, as non-isentropic changes due to nuclear and neutrino sources are computed elsewhere by operator splitting.

The update of the specific internal energy thus takes the form

Equation (124)

where the time-centering of the pressure p is accomplished by a predictor-corrector loop, completing the Lagrangian step.

4.4. Remap Step

Following the Lagrangian step, in the case of the θ- and ϕ-sweeps, the grid is remapped back to the configuration that prevailed before the Lagrangian step, thus, making the combination of a Lagrangian step and a remap step effectively an Eulerian step. In the case of the radial sweep, the grid is remapped back to a configuration specified by the regridder, which will be described in Section 4.5. Via the regridder options, the user can specify that the grid following the Lagrangian step be left as is, remapped back to the configuration that prevailed before the Lagrangian step, or, by invoking one of the regridder algorithms, remapped to a configuration that tends to optimize the resolution of structures that develop during a simulation.

4.4.1. Remapping Mass, Momenta, and Angular Momenta

For quantities like the mass, specific momenta (momenta per gram), and specific angular momenta, the remapping procedure is straightforward. Denoting, as before, the values of the grid and other variables after the Lagrangian step by the superscript $n+1^{\prime} $, and after the remap step by the superscript $n+1$, the difference $\delta \xi $, between the values of the grid variables at $n+1^{\prime} $ and $n+1$, given by

Equation (125)

is computed, as is the volume $\delta {V}_{i+\tfrac{1}{2}}$ contained within $\delta {\xi }_{i+\tfrac{1}{2}}$. The latter is given by

Equation (126)

If $\delta {\xi }_{i+\tfrac{1}{2}}\gt 0$, the grid interface ${\xi }_{i+\tfrac{1}{2}}^{n+1}$ is placed by the remap step between ${\xi }_{i-\tfrac{1}{2}}^{n+1^{\prime} }$ and ${\xi }_{i+\tfrac{1}{2}}^{n+1^{\prime} }$. The mass advected across the zone interface ${\xi }_{i+\tfrac{1}{2}}$ is next computed by first constructing a PPM profile, $\rho (\xi )$, of the density ${\rho }_{i}^{n+1^{\prime} }$. The value, $\langle \rho {\rangle }_{L,i+\tfrac{1}{2}}^{n+1^{\prime} }$, of the interpolated density $\rho (\xi )$ averaged over the interval $\delta {\xi }_{i+\tfrac{1}{2}}$ to the left of interface ${\xi }_{i+\tfrac{1}{2}}^{n+1^{\prime} }$ is then computed by Equation (102), with $\delta {\xi }_{i+\tfrac{1}{2}}$ replacing ${c}_{{\rm{s}},i}{\rm{\Delta }}t$. The mass advected is then given by

Equation (127)

The zone mass after the remap step, ${{ \mathcal M }}_{i}^{n+1}$, is then computed from the the zone mass before the remap step, ${{ \mathcal M }}_{i}^{n+1^{\prime} }$, and the advected mass $\delta {{ \mathcal M }}_{i+\tfrac{1}{2}}^{n+1^{\prime} }$ by

Equation (128)

and the density ${\rho }_{i}^{n+1}$ is then computed by

Equation (129)

Having determined the advected mass $\delta {{ \mathcal M }}_{i+\tfrac{1}{2}}^{n+1^{\prime} }$, the remapping of any specific (i.e., per gram) quantity ${a}_{i}^{n+1^{\prime} }$ proceeds by determining the amount, $\delta {{ \mathcal A }}_{i+\tfrac{1}{2}}^{n+1^{\prime} }$, of that quantity advected. As in the case of the density, a piecewise parabolic interpolation profile, $a(\xi )$, of ${a}_{i}^{n+1^{\prime} }$ is constructed and the average $\langle a{\rangle }_{L,i+\tfrac{1}{2}}^{n+1^{\prime} }$ of that quantity over the interval $\delta {\xi }_{i+\tfrac{1}{2}}$ is computed using Equation (102). The mass of $a(\xi )$ contained within $\delta {\xi }_{i+\tfrac{1}{2}}$ is finally computed as the penultimate step by

Equation (130)

If $\delta {\xi }_{i+\tfrac{1}{2}}\lt 0$, the grid interface ${\xi }_{i+\tfrac{1}{2}}^{n+1}$ is placed by the remap step between ${\xi }_{i+\tfrac{1}{2}}^{n+1^{\prime} }$ and ${\xi }_{i+\tfrac{3}{2}}^{n+1^{\prime} }$, and the quantity of $a(\xi )$ within $\delta {\xi }_{i+\tfrac{1}{2}}$ is advected from the right to the left of interface ${\xi }_{i+\tfrac{1}{2}}^{n+1}$. The procedure is the same as described above for the $\delta {\xi }_{i+\tfrac{1}{2}}\gt 0$ case, except that $\langle \rho {\rangle }_{R,i+\tfrac{1}{2}}^{n+1^{\prime} }$ and $\langle a{\rangle }_{R,i+\tfrac{1}{2}}^{n+1^{\prime} }$ are computed from ${\rho }_{i+1}^{n+1^{\prime} }$ and ${a}_{i+1}^{n+1^{\prime} }$ by Equation (103), giving $\delta {{ \mathcal M }}_{i+\tfrac{1}{2}}^{n+1^{\prime} }$ and $\delta {{ \mathcal A }}_{i+\tfrac{1}{2}}^{n+1^{\prime} }$ to the right of interface ${\xi }_{i+\tfrac{1}{2}}^{n+1^{\prime} }$, with $\delta {\xi }_{i+\tfrac{1}{2}}$ replacing ${c}_{{\rm{s}},i}{\rm{\Delta }}t$, as before, but now with a negative value of $\delta {V}_{i+\tfrac{1}{2}}$. Therefore, the quantity $\delta {{ \mathcal M }}_{i+\tfrac{1}{2}}^{n+1^{\prime} }$ computed from $\langle \rho {\rangle }_{R,i+\tfrac{1}{2}}^{n+1^{\prime} }$ by an equation analogous to Equation (127), and the quantity $\delta {{ \mathcal A }}_{i+\tfrac{1}{2}}^{n+1^{\prime} }$ computed from $\langle a{\rangle }_{R,i+\tfrac{1}{2}}^{n+1^{\prime} }$ by an equations analogous to Equation (130), are both negative.

The remap step is finally completed by performing the advection:

Equation (131)

Negative values of $\delta {{ \mathcal A }}_{i+\tfrac{1}{2}}^{n+1^{\prime} }$ or $\delta {{ \mathcal A }}_{i-\tfrac{1}{2}}^{n+1^{\prime} }$ mean simply that the advection proceeds from right to left rather than in the other direction. The advection step is conservative since the same amount of a given quantity enters a zone as leaves the adjacent zone. Again, the sign of $\delta {{ \mathcal A }}_{i+\tfrac{1}{2}}^{n+1^{\prime} }$ determines whether the advection is from zone i to zone $i+1$, or vice versa.

4.4.2. Remapping Composition and Electron Fraction

The algorithm for the advection of the composition mass fractions, Xn, and net electron fraction, ${Y}_{{\rm{e}}}$, depends on whether the matter on either side of the zone interface is in NSE or not (non-NSE). If the the advection is between zones in NSE, the composition is given by the EoS as a function of the values of ρ, ${e}_{\mathrm{int}}$, and ${Y}_{{\rm{e}}}$ of the material being advected. In this case, no explicit advection of composition mass fractions needs to be performed. The advection of ${Y}_{{\rm{e}}}$ in this case proceeds as described in Section 4.4.1 above. If the advection is between zones in non-NSE, then along with ρ and ${e}_{\mathrm{int}}$, the advection of the composition mass fractions, Xn, is carried out explicitly as described in Section 4.4.1 and in accordance with the consistent multi-fluid advection method of Plewa & Müller (1999). That is, the average value of each composition mass fraction Xn to be advected is computed by Equations (102) or (103) depending on whether $\delta \xi $ is positive or negative. The resulting composition in the mass to be advected is then normalized to unity before performing the advection. In order for the net ${Y}_{{\rm{e}}}$ advected in this case to be consistent with the net proton fraction of the advected composition, the ${Y}_{{\rm{e}}}$ advected is taken to be that of the net ${Y}_{{\rm{e}}}$ of the advected composition.

If the advection is from a zone in non-NSE to one in NSE, the composition of the material to be advected is computed as described above for the non-NSE to non-NSE case. The material advected is then assumed to become part of the NSE material in the acceptor zone, and only the independent thermodynamic variables ρ, ${e}_{\mathrm{int}}$, and ${Y}_{{\rm{e}}}$ of the material to be advected are advected. Finally, if material is advected from a zone in NSE to one in non-NSE, the material advected is first "deflashed," that is, its NSE composition mass fractions are extracted from the EoS, stored in a temporary mass fraction array, and advected as described above for the case of two adjacent non-NSE zones.

4.4.3. Multiple EoSs and the Energy Remap

Chimera is designed to accommodate multiple EoSs that are applied in contiguous density ranges. As explained above in Section 3.2, for example, Chimera has used the LS EoS at densities above 1011 ${{\rm{g\; cm}}}^{-3}$, and a different EoS below that density. Since Chimera updates the specific energy directly rather than the temperature and uses the updated specific energy and other needed thermodynamic variables to update the temperature, slight differences in the energy zeros at the boundary between two EoS's could result in unphysical temperature updates. Specific energy differences between two EoSs could also arise from peculiarities or approximations peculiar to each EoS. In either case, we will refer to the potential difference in specific energies given by two EoS's for the same thermodynamic state as a zero-energy offset. Unphysical temperature updates could happen, for example, during the remap step if matter from a zone linked to one EoS is advected into an adjacent zone linked to a different EoS. The quantity of energy advected would contain the difference in the energy zeros as well as the physically relevant energy. This problem could affect the radial sweep but not the θ- or ϕ-sweeps as the same EoS is always used along a θ- or ϕ-directed ray.

Figure 12.

Figure 12. Schematic representation of the grid used to construct the finite-difference equations for the hydrodynamics in any given dimension. Integer zone indices represent zones, and half-integer zone indices represent zone interfaces.

Standard image High-resolution image

To avoid this problem, Chimera overlaps by four radial zones in either direction the specific internal energy at the boundary between EoSs as shown in Figure 13, where radial zones i − 1 and lower are tied to a particular EoS A, while zones i and higher are tied to a different EoS B. The two EoSs have different zero-energy offsets, as indicated by their vertical displacement, and an overlap of four zones in either direction. This overlap enables PPM profiles of the specific internal energy for zones i − 1 and below to be constructed using EoS A up to the interface i − $\tfrac{1}{2}$. Likewise, PPM profiles of the specific internal energy for zones i and above can be constructed using EoS B. The example in Figure 13 is one in which the zone interface ${\xi }_{i-\tfrac{1}{2}}$ is remapped a distance $\delta {\xi }_{i-\tfrac{1}{2}}$ to the left of its original position, its new position being indicated by the vertical red dashed line. This entails that a quantity of specific internal energy contained within $\delta {\xi }_{i-\tfrac{1}{2}}$ be advected from the left to the right of zone interface ${\xi }_{i-\tfrac{1}{2}}$. Chimera performs this advection by advecting the energy within $\delta {\xi }_{i-\tfrac{1}{2}}$ as given by EoS A out of zone i − 1, and advecting the energy within the same $\delta {\xi }_{i-\tfrac{1}{2}}$ but as given by EoS B into zone i. Since the specific internal energy advected out of a zone and into the adjacent zone is consistent with the different EoSs tied to each of the two zones, the unphysical temperature jump that would occur if the zero-energy offset had not been accounted for in the advection is avoided.

Figure 13.

Figure 13. Schematic representation of the advection of energy from the left to the right of zone interface $i-\tfrac{1}{2}$. Zones $i-1$ and below are tied to EoS A while zones i and above are tied to EoS B. EoS A and EoS B have a zero energy offset indicated by their vertical displacement from each other. In remapping zone interface i − $\tfrac{1}{2}$ a distance $\delta {\xi }_{i-\tfrac{1}{2}}$ from its initial location to the location indicated by the red vertical dashed line, the energy within $\delta {\xi }_{i-\tfrac{1}{2}}$ given by EoS A must be advected from zone i − 1 to zone i. This is accomplished by advecting the energy within $\delta {\xi }_{i-\tfrac{1}{2}}$ given by EoS A out of zone i − 1, and advecting the energy within the same interval $\delta {\xi }_{i-\tfrac{1}{2}}$ but given by EoS B into zone i. The four-zone overlap on each side exists so that a PPM profile of the energy in zone i − 1 can be constructed from both EoSs.

Standard image High-resolution image

4.4.4. Nuclear Binding Energy

In advecting the specific energy between two adjacent zones during a remap, the specific internal energy is split into a nuclear binding energy component, ${e}_{\mathrm{bind}}$, and the rest of the energy, and the two components are advected separately. For example, the specific internal energy, ${e}_{\mathrm{int}}$, is split into ${e}_{\mathrm{bind}}$ and ${e}_{\mathrm{th}}={e}_{\mathrm{int}}-{e}_{\mathrm{bind}}$, and remapped as follows:

Equation (132)

where ${e}_{\mathrm{th},i}^{n+1^{\prime} }$ and ${e}_{\mathrm{bind},i}^{n+1^{\prime} }$ are the specific internal energy minus the specific binding energy, and the specific binding energy, respectively, of zone i, and $\delta {E}_{\mathrm{th},i\pm \tfrac{1}{2}}^{n+1^{\prime} }$ and $\delta {E}_{\mathrm{bind},i\pm \tfrac{1}{2}}^{n+1^{\prime} }$ are the internal energy minus the binding energy and the binding energy, respectively, transferred through outer zone edge $\left(i+\tfrac{1}{2}\right)$ and inner zone edge $\left(i-\tfrac{1}{2}\right)$ of zone i. The masses $\delta {{ \mathcal M }}_{i}^{n+1^{\prime} }$ and $\delta {{ \mathcal M }}_{i}^{n+1}$ are the masses of zone i before and after the remap, respectively, as given by Equation (127). This mode of energy advection is appropriate for advecting non-NSE material during a remap, as the non-NSE material being advected does not necessarily have the same composition and, therefore, binding energy, as the original material in the zone from which it is being advected. The composition of the material being advected is obtained by integrating over the PPM profile constructed for each ionic mass fraction, and then normalizing the sum to unity. Since the PPM profiles of different ionic mass fractions may be differently shaped, the composition that results after integrating over the portion of the profiles being advected and normalizing can result in a composition and binding energy different from the original. In the energy advection procedure described above, the binding energy of the advecting material is computed once its composition is ascertained, and the rest of the advected energy, $\delta {e}_{\mathrm{th}}$, is also computed by integrating the PPM profile of ${e}_{\mathrm{th}}={e}_{\mathrm{int}}-{e}_{\mathrm{bind}}$ over the advecting mass. The net energy advected should thus reflect both the correct nuclear binding energy of the advecting material, as well as its thermal component.

4.4.5. Energy Remap for the θ- and ϕ-sweeps and the Preliminary Remap for the Radial Sweep

The energy remap described above is the ultimate step in the θ- and ϕ-sweep hydrodynamics, and the penultimate step in the radial sweep hydrodynamics. Away from shocks, the specific internal energy is remapped as given by Equation (132), as opposed to remapping the sum of the specific internal plus specific kinetic energy, ${e}_{\mathrm{tot}^{\prime} ,i}={e}_{\mathrm{int}}+{e}_{\mathrm{kin}}$. This is permissible as the flow is isentropic apart from the contributions of nuclear and neutrino source terms, which have been included elsewhere in the radial sweep (see the text above Equation (124) for the motivation for this approach). In the vicinity of shocks, the specific energy, ${e}_{\mathrm{tot}^{\prime} }$, which has been evolved during the Lagrangian step using the Rankine–Hugoniot equations, must be remapped, and the specific internal energy is extracted afterwards. We define ${e}_{\mathrm{tot}^{\prime\prime} }={e}_{\mathrm{tot}^{\prime} }-{e}_{\mathrm{bind}}={e}_{\mathrm{th}}+{e}_{\mathrm{kin}}$ and remap ${e}_{\mathrm{tot}^{\prime\prime} }$ and ${e}_{\mathrm{bind}}$ separately, as described above. To use consistent values for the left and right states for calculating a PPM profile of ${e}_{\mathrm{tot}^{\prime\prime} }$ (i.e., consistent with the velocity remap), we define

Equation (133)

Having calculated the PPM profile of ${e}_{\mathrm{tot}^{\prime\prime} }$, the quantity of ${e}_{\mathrm{tot}^{\prime\prime} }$ advected across a given zone interface is computed by Equations (127) or (130), and the remapping of ${e}_{\mathrm{tot}^{\prime\prime} }$ is performed by an equation analogous to Equation (132). After the remapping of ${e}_{\mathrm{tot}^{\prime\prime} }$ and ${e}_{\mathrm{bind}}$, the specific internal energy is extracted from their sum, ${e}_{\mathrm{tot}^{\prime} }$.

4.4.6. Recomputation of the Gravitational Potential and the Computation of $\partial {e}_{\mathrm{grav}}/\partial t$

Following the remap in the radial-sweep of the mass, momenta, angular momenta, the independent thermodynamic variables, and the preliminary remap of the specific energy, the specific gravitational potential, ${e}_{\mathrm{grav}}^{n+1}$, is computed. In the case in which the spherically symmetric component of the gravitational potential is computed by means of a general relativistic approximation, the remapped pressure, specific energy, and neutrino contributions are used as sources of gravity as well as the density, which necessitated the preliminary remap of the specific energy. It is at this stage of the radial-sweep that the quantity $\partial {e}_{\mathrm{grav}}/\partial t$ is computed and added, in accordance with Equation (99), to the specific total energy, ${e}_{\mathrm{tot}}^{n+1^{\prime} }$, that was computed by Equation (121) during the Lagrangian step. The most direct procedure for calculating $\partial {e}_{\mathrm{grav}}/\partial t$ would be to calculate the gravitational potential, ${e}_{\mathrm{grav}}^{n+1^{\prime} }$, after the Lagrangian step, interpolate the initial gravitational potential, ${e}_{\mathrm{grav}}^{n}$, to the Lagrangian grid to get the quantity ${e}_{\mathrm{grav},I-L}^{n}$, and approximate $\partial {e}_{\mathrm{grav}}/\partial t$ by $\left({e}_{\mathrm{grav}}^{n+1^{\prime} }-{e}_{\mathrm{grav},I-L}^{n}\right)/{\rm{\Delta }}t$. This would work well for one-dimensional simulations, but for multidimensional simulations, a given radial grid edge, ${\xi }_{i+\tfrac{1}{2}}$, after the Lagrangian step is a function of θ and/or ϕ, making the gravitational potential difficult to compute at this point. Instead, the gravitational potential is computed after the remap of the radial grid and interpolated as a function of θ and/or ϕ to the Lagrangian grid, thereby obtaining ${e}_{\mathrm{grav},F-L}^{n+1}$. The time derivative of the gravitational potential added to ${e}_{\mathrm{tot}}^{n+1^{\prime} }$ is thus given, as a function of θ and/or ϕ, by

Equation (134)

4.4.7. Final Radial-sweep Remap of the Total Energy

The final remapping of the specific total energy $\left({e}_{\mathrm{kin}}+{e}_{\mathrm{int}}+{e}_{\mathrm{grav}}\right)$ in the radial sweep begins with the specific energy, ${e}_{\mathrm{tot}^{\prime\prime} ,i}^{n+1^{\prime} }$, given by

Equation (135)

where ${e}_{\mathrm{tot},i}$ has been updated during the Lagrangian step, as given by Equation (121). Consistent left and right states of ${e}_{\mathrm{tot}^{\prime\prime} ,i}^{n+1^{\prime} }$ are determined as specified by Equations (133), PPM profiles of ${e}_{\mathrm{tot}^{\prime} ,i}^{n+1^{\prime} }={e}_{\mathrm{tot}^{\prime\prime} ,i}^{n+1^{\prime} }+{e}_{\mathrm{kin},i}^{n+1^{\prime} }$ are then obtained, and the amount of ${e}_{\mathrm{tot}^{\prime} \,i}^{n+1^{\prime} }$ to be advected across the zone interface are given by Equation (130) and the procedure outlined in the discussion below it. Remapping then proceeds in accordance with an equation analogous to Equation (132), and the final specific internal energy is extracted from ${e}_{\mathrm{tot},i}^{n+1}={e}_{\mathrm{tot}^{\prime} ,i}^{n+1}+{e}_{\mathrm{bind},i}^{n+1}$ by

Equation (136)

This completes the remap step of the radial-sweep hydrodynamics.

4.4.8. Suppression of Carbuncles

When shocks are aligned with one of the coordinate directions in multidimensional simulations, they are susceptible to an "odd–even decoupling" or "carbuncle" instability (Quirk 1994; Liou 2000; Sutherland et al. 2003). This can lead to a strong rippling of the shock front, which could, in turn, excite hydrodynamic instabilities in the post-shock region. Our scheme for suppressing this instability is the use of a "local oscillation filter" similar in philosophy to that described by Sutherland et al. (2003). This approach is local and does not affect the well-resolved features of the flow elsewhere. To suppress the carbuncle instability in the radial direction, which is where this instability typically arises in a supernova simulation, the angular (θ) and azimuthal (ϕ) remap steps are each followed by an examination of the radial velocities along the angular and azimuthal rays to search for radial velocity extrema. If there are at least three radial velocity extrema in any group of five adjacent zones, and if a shock is present, these zones are marked for "smoothing." For the particular case in which zones m and $m+1$ are marked for smoothing, the flux $\delta {{ \mathcal A }}_{m+\tfrac{1}{2}}$ of the quantity a, defined by

Equation (137)

is computed, where ${{ \mathcal M }}_{m}$ is the mass of zone m, and ${c}_{\mathrm{smooth}}$ is an empirical parameter. Experimentation has shown that a value of 0.075 for ${c}_{\mathrm{smooth}}$ works well. The final step in the procedure is to sweep across the angular and azimuthal rays and exchange the flux $\delta {{ \mathcal A }}_{m+\tfrac{1}{2}}$ between the zones marked for smoothing in a step analogous to that described by Equation (131). This reduces the difference in the values of the quantity a between adjacent zones, thereby inhibiting the growth of this difference. Applying this procedure to the quantities ur and ${u}_{\phi }$, and to ${L}_{\phi }$ (to smooth ${u}_{\theta }$), proved sufficiently robust to suppress the carbuncle instability.

4.5. Radial Regridder

The PPM Lagrangian-Remap format permits the grid after the Lagrangian step to be remapped to a grid other than the initial grid from which the Lagrangian step originated. While the θ- and ϕ-grids are remapped back to their initial grids after the Lagrangian step, making them effectively Eulerian, Chimera uses the remapping freedom to provide the user with a number of remapping options for the radial grid to ensure that the grid continues to resolve important structures that arise during the course of a simulation. One option is for the radial grid to be remapped back to the initial grid after the Lagrangian step, making the grid effectively Eulerian as it is for the θ and ϕ grids. Another option is for the remapped grid to follow the mean motion of the fluid, referred to here and below as pseudo-Lagrangian, making the grid purely Lagrangian in the case of spherically symmetric fluid flow.

Currently, a number of more sophisticated options are available specific to the pre-bounce or post-bounce phase of a CCSN simulation. For both the pre-bounce phase and the post-bounce phase, an inner-outer boundary dividing the radial grid into an inner and an outer section is determined based on a number of user selected criteria. These criteria can differ between the pre-bounce and the post-bounce phases and can differ at user-selected time intervals during the post-bounce phase. During both the pre-bounce and the post-bounce phases, the outer grid can be selected to be pseudo-Lagrangian, which is useful if there are sharp chemical discontinuities in non-NSE material that need to be preserved, or Eulerian if advection through the outer boundary of a prescribed distribution of material is important. During the pre-bounce phase, the inner grid starts out as pseudo-Lagrangian, but blends into another grid between two user-selected densities. This second grid is constructed so that adjacent zones satisfy ${\rm{\Delta }}{r}_{i+1}={\rm{constant}}\times \ {\rm{\Delta }}{r}_{i}$, referred to here as a "zoomed grid," with the properties that the width of the outer zone of this zoomed grid is equal to the zone width of the first zone of the outer grid, and the inner zone tends to a user selected width when that zone reaches $3\times {10}^{14}$ ${{\rm{g\; cm}}}^{-3}$. The result is a smooth and smoothly evolving grid that can be tuned to provide the desired grid resolution at the proto-neutron star surface when it forms. During the post-bounce phase, the inner grid remains a zoomed grid from the core center to a density of 1014 ${{\rm{g\; cm}}}^{-3}$, with the central zone width such that it would attain a user-selected zone width at a density of $3\times {10}^{14}$ ${{\rm{g\; cm}}}^{-3}$. A second zoomed grid covers the density range from 1014 to 1012 ${{\rm{g\; cm}}}^{-3}$, and a third covers the density range from 1012 to 1010 ${{\rm{g\; cm}}}^{-3}$. Both of these latter two grids have the same number of zones, which are equal to a user-selected value. This ensures there are a sufficient number of zones to resolve the neutrinosphere as the proto-neutron star shrinks and as the density cliff forms near its edge. Finally, a fourth zoomed grid covers the region from the density of 1010 ${{\rm{g\; cm}}}^{-3}$ to the outer edge of the inner grid. The result of this regridding is again a smooth and smoothly evolving user-controlled grid designed to resolve the critical features that arise during the course of a CCSN simulation.

4.6. Gravity Solver

Self gravity can be chosen to be either one- or multi-dimensional, with a further choice of either a Newtonian or an approximate general relativistic monopole component. The approximate general relativistic gravitational potential used for the monopole component in the latter case is a modified Tolman–Oppenheimer–Volkoff (TOV) potential suggested by Marek et al. (2006, Case A) and described briefly below in Section 4.6.1. Multidimensional gravity is obtained by expanding the Newtonian gravitational potential in a multipole expansion as described by Müller & Steinmetz (1995) and below in Section 4.6.2. Approximate general relativistic multidimensional gravity is obtained by replacing the Newtonian monopole in the multipole expansion by the approximate general relativistic monopole.

4.6.1. One-dimensional Gravitational Potential

Newtonian monopole gravity is trivial. The radial zone-edged and zone-centered gravitational accelerations, ${g}_{i+\tfrac{1}{2}}$ and gi, respectively, are given by

Equation (138)

where ${M}_{i+\tfrac{1}{2}}$ is the rest mass enclosed in a volume of radius ${R}_{i+\tfrac{1}{2}}$, Ri the mass-averaged zone-centered radius, and G the gravitational constant. The radial zone-edged and zone-centered gravitational potentials, ${e}_{\mathrm{grav},i+\tfrac{1}{2}}$ and ${e}_{\mathrm{grav},i}$, respectively, are given by

Equation (139)

where at the outer edge of the radial grid

Equation (140)

Approximate GR monopole gravity is computed by first iterating the following two equations (Marek et al. 2006, Case A) for ${M}_{\mathrm{TOV}}$:

Equation (141)

Equation (142)

where ${e}_{\mathrm{int},i}$ and ${e}_{\nu ,i}$ are the specific energy densities of matter and neutrinos, respectively, and ${u}_{r,i}$ is the radial velocity. With ${M}_{\mathrm{TOV},I}$ computed as $\tfrac{1}{2}\left({M}_{\mathrm{TOV},I-\tfrac{1}{2}}+{M}_{\mathrm{TOV},I+\tfrac{1}{2}}\right)$, the zone-centered gravitational force is computed by

Equation (143)

where ${p}_{\mathrm{gas},i}$ and ${p}_{\nu ,i}$ are the matter pressure and spherically averaged neutrino pressure, respectively. Once gi is computed, ${g}_{i+\tfrac{1}{2}}$ is computed as $\left({g}_{i}+{g}_{i+1}\right)/2$, with ${g}_{0+\tfrac{1}{2}}=0$ and ${g}_{I+\tfrac{1}{2}}$ extrapolated from gI and ${g}_{I-\tfrac{1}{2}}$. The zone-edged gravitational potential, ${e}_{\mathrm{grav},I+\tfrac{1}{2}}$, is then computed by Equation (139), and the zone-centered gravitational potential is given by ${g}_{i}\,=\left({g}_{i-\tfrac{1}{2}}+{g}_{i+\tfrac{1}{2}}\right)/2$.

4.6.2. Multipole Expansion of the Gravitational Potential—Axisymmetry

To incorporate nonspherical gravity, Chimera uses a scheme based on the method described by Müller & Steinmetz (1995) of expanding the integral Newtonian Poisson equation in a multipole expansion. When implementing approximate general relativistic gravity, the Newtonian monopole is replaced with the approximate general relativistic monopole (Marek et al. 2006, Case A) described above. Multipole gravity is implemented in both axisymmetric and three-dimensional simulations.

For axisymmetric simulations, this scheme utilizes the identity

Equation (144)

where

Equation (145)

${\rm{\Theta }}(x)$ is the Heaviside function, and P is the Legendre polynomial of order ${\ell }$, to expand the Poisson integral in a Legendre series:

Equation (146)

Here, the ${P}_{{\ell }}^{m}$ are the associated Legendre functions and the integral over $\phi ^{\prime} $ has been performed utilizing the assumption of axisymmetry. In differencing Equation (146), the spherical coordinate system utilized by Chimera enables the radius $r^{\prime} $ to be integrated over each spherical shell and the potential to be computed at zone interfaces. Given the singular nature of the Poisson equation, this avoids the problem of the gravitational self-interaction, which can lead to a nonconvergence of the multipole expansion, as pointed out by Couch et al. (2013). Equation (146) in differenced form then becomes

Equation (147)

Equation (148)

where

Equation (149)

The Legendre polynomials up to the specified order and the angular integrations of these polynomials over the angular zone widths are generated as an initialization step. The Legendre polynomials are first computed at the angular zone edges using

Equation (150)

and the recurrence relation for  > 1:

Equation (151)

The Legendre polynomials are then integrated over the zone widths, using

Equation (152)

and then

Equation (153)

To assess the accuracy of the gravitational potential expansion as a function of the number of multipoles used, and for code verification, the gravitational potentials computed by Equations (147) and (148) for a Maclaurin spheroid are compared to those given by an exact analytic solution in Section 5.7.

4.6.3. Multipole Expansion of the Gravitational Potential—Non-axisymmetry

For non-axisymmetric simulations, this scheme utilizes the identity

Equation (154)

to expand the gravitational potential in spherical harmonics, where ${Y}_{{\ell }}^{m* }$ denotes the complex conjugate of the spherical harmonic ${Y}_{{\ell }}^{m}$, and ${r}_{\lt }^{{\ell }}/{r}_{\gt }^{{\ell }+1}$ is defined by Equation (145). The expansion is given by

Equation (155)

where

Equation (156)

Equation (157)

and where

Equation (158)

and

Equation (159)

The occasional use of $i=\sqrt{-1}$ here should not cause confusion with the radial index i. The spherical harmonics ${Y}_{{\ell }}^{m}(\theta ,\phi )$ have been written as functions of the normalized associated Legendre functions ${\tilde{P}}_{{\ell }}^{m}(\cos \theta )$:

Equation (160)

and the latter are given in terms of the unnormalized associated Legendre polynomials ${P}_{{\ell }}^{m}(\cos \theta )$ by

Equation (161)

where ${\tilde{P}}_{{\ell }}^{m}(\cos \theta )$ satisfies the relation

Equation (162)

The normalized associated Legendre functions, ${\tilde{P}}_{{\ell }}^{m}(\cos \theta )$, and their integrals, ${P}_{\mathrm{int}}({\ell },m,j)$, defined below by Equation (172) are calculated during an initial setup step by using the subroutines developed by NGA10 based on the algorithms of Paul (1978) and Gerstl (1980).

To transform $A(r^{\prime} ,{\ell },m)$ into a form suitable for calculation, we first change variables from $\theta ^{\prime} $ to $y^{\prime} \equiv \cos \theta ^{\prime} $, to obtain

Equation (163)

which, separating real and imaginary parts, can be written as

Equation (164)

where

Equation (165)

Equation (166)

Next, we define the variables

Equation (167)

and substitute them into Equations (165) and (166) to obtain

Equation (168)

Equation (169)

Note the symmetry conditions

Equation (170)

To discretize the ${A}_{r}(r^{\prime} ,{\ell },m)$, we introduce a j index, $j=1,\cdots ,J$, and a k index, k = 1, ⋯, K, to denote the zone centers of the angular variable, θ, and the azimuthal variable, ϕ, respectively, analogous to the index i = 1, ⋯, I, which is illustrated in Figure 12 and used in this section for the radial index. Half-integer values of i, j, and k refer to zone edges. Primed and unprimed indices will refer to source and field quantities, respectively. We have

Equation (171)

where

Equation (172)

and where

Equation (173)

Similarly,

Equation (174)

where

Equation (175)

To construct the gravitational potential, we have, from Equations (155)–(157) and (160),

Equation (176)

Now, the expression ${\tilde{P}}_{{\ell }}^{m}(\cos \theta ){e}^{{im}\phi }A(r^{\prime} ,{\ell },m)$ in Equation (176) can be written

Equation (177)

where the real and imaginary parts of $B(r^{\prime} ,{\ell },m)$ are given by

Equation (178)

Equation (179)

Using the symmetry conditions expressed by Equations (162) and (170), we have

Equation (180)

and

Equation (181)

which gives

Equation (182)

and

Equation (183)

Thus, the imaginary part of Equation (177) cancels out when summed over negative and positive m, and Equation (182) for the real part of $B(r^{\prime} ,{\ell },-m,\phi )$ can be used to replace the summation of m from $-{\ell }$ to ${\ell }$ to a summation from zero to .

The construction of the gravitational potential in Chimera begins with the computation of $A(r^{\prime} ,{\ell },m)$, defined in Equation (158) and discretized as described by Equations (171)–(175). With ${P}_{\mathrm{int}}$, Sint, and Cint given by Equations (172), (173), and (175), respectively, and computed in the initial setup, $A(r^{\prime} ,{\ell },m)$ is computed from

Equation (184)

Using Equations (182) and (183), the quantities ϕin(r) and ϕout(r), given, respectively, by ${e}^{{im}\phi }{C}_{(r)}^{{\ell }m}$ and ${e}^{{im}\phi }{D}_{(r)}^{{\ell }m}$ (Equations (156) and (157)), are then computed recursively by

Equation (185)

Equation (186)

if ${\ell }\ne 2$, and

Equation (187)

if ${\ell }=2$, and

Equation (188)

if ${\ell }=0$. The gravitational potential is finally calculated as

Equation (189)

To assess the accuracy of the gravitational potential expansion as a function of multipole number and to verify the code, the gravitational potentials computed by Equation (189) for a Maclaurin spheroid are compared to those given by an exact analytic solution in Section 5.7.

5. Hydrodynamics Test Problems

We have subjected the Chimera hydrodynamics code modules to a number of test problems. The key problems are described here.

5.1. Point-blast Explosion

The ability of a supernova code to simulate a spherical outgoing shock is an important first test. An analytic solution for a spherical outgoing shock is available for the "point-blast explosion" problem, which consists of the instantaneous deposition of an amount of energy, E0, at a point in a zero-gravity, stationary, uniform medium of constant density, ${\rho }_{0}$. The energy E0 is required to be very large in comparison with the initial energy of the medium. The analytic solution for this problem was found by Taylor (1950) and Sedov (1959), and various parts of the solution can be found in a number of text books, such as Mihalas & Mihalas (1984) or Zel'dovich & Raizer (1967), with the most complete account being given in Landau & Lifshitz (1959). Because of typos in the publications cited above, we have rederived the solution and have found that Equation (99.8) of Landau & Lifshitz (1959) should be multiplied by the square of the normalized radius (see below) and that the expression for ${\nu }_{5}$ in Equation (99.10) of Landau & Lifshitz (1959) should be replaced by

Equation (190)

To set up this problem, a uniform, $\gamma =5/3$ gas maintained in hydrostatic equilibrium within a spherical volume by an external pressure boundary condition and having a constant density ${\rho }_{0}=0.1$ ${{\rm{g\; cm}}}^{-3}$ was divided into 200 zones of equal, 1 cm width. The gas was given a constant ambient temperature of 10−8 MeV. A point explosion at the center of the spherical mass was simulated by instantaneously depositing an energy, ${E}_{0}\simeq 6.06\times {10}^{17}$ erg, by increasing the temperature of the first zone to 1 MeV. Rather than using a simple gamma-law EoS, we used the Chimera non-NSE EoS, consisting of half neutrons and half protons. The electron and photon contributions were set to zero. In this way, the EoS used was equivalent to a gamma-law EoS with $\gamma =5/3$, but the Chimera EoS machinery (e.g., composition remapping, EoS interpolation) described in Section 3 was tested as well.

Denoting the distance of the shock from the origin by the term rs, the time-dependence of rs is given by

Equation (191)

where ${\xi }_{s}$ is determined by

Equation (192)

where u, ρ, and p are the fluid velocity, density, and pressure, respectively, between the origin and the shock. The subscript "s" denotes their immediate post-shock values. Numerically integrating Equation (192), using Equations (99.10) of Landau & Lifshitz (1959), gives ${\xi }_{s}=1.17$.

A comparison of rs as a function of time, computed by Chimera versus the analytic result given by Equation (191) is shown in Figure 14(a) and demonstrates agreement to within a few percent. Figure 14(b) shows that the velocity, density, and pressure of the fluid behind the shock normalized by their immediate post-shock values also agree well with the analytic solutions.

Figure 14.

Figure 14. Panel (a): shock radii vs. time for the point-blast explosion problem for Chimera (black plus signs) and the analytic expression (red line) given by Equations (191) and (192). Panel (b): velocity (red), density (green), and pressure (blue) behind the shock as a function of radius, normalized by their immediate post-shock values for Chimera (pluses) relative to the analytic solution given by Equations (99.8) and (99.10) of Landau & Lifshitz (1959) with the corrections noted in the text.

Standard image High-resolution image

5.2. Sod Shock Tube Problem

A standard hydrodynamics problem admitting an analytic solution is the Sod shock tube problem (Sod 1978). As Sod's formulation was done in the context of a plane-parallel geometry, we approximate such a geometry with Chimera's spherical grid by working with a variable $x=R-r$, where $R={10}^{5}$ cm and −2 ≤ r ≤ 2 cm. We cover the $-2\leqslant r\leqslant 2$ range of r with a grid consisting of 200 zones, initially equally spaced. Following Sod's original formulation (Sod 1978), we set up a Riemann problem with pressure, p = 1 erg cm−3, density, $\rho =1$ ${{\rm{g\; cm}}}^{-3}$, velocity, u = 0 ${{\rm{cm\; s}}}^{-1}$, for $x\lt 0$, and p = 0.1 erg cm−3, $\rho =0.125$ ${{\rm{g\; cm}}}^{-3}$, and u = 0 ${{\rm{cm\; s}}}^{-1}$, for $x\gt 0$. The same EoS was used as described above for the point-blast problem, resulting in a constant γ of 5/3. The results of the test at time $t=0.5\,{\rm{s}}$ are shown in Figure 15.

Figure 15.

Figure 15. Comparison of the analytic solution for $\gamma =5/3$ (solid lines) and Chimera (crosses) solutions to the Sod shock tube problem for velocity (red), density (green), and pressure (blue), as a function of distance from the initial discontinuity at 0.5 s. Chimera simulations were performed (a) in Lagrangian mode with the values of the flatten parameters suggested by Colella & Woodward (1984), (b) in Lagrangian mode with a different set of flatten parameters that suppress the low-amplitude post-shock oscillations, and (c) in Eulerian mode with the values of the flatten parameters suggested by Colella & Woodward (1984).

Standard image High-resolution image

Numerical hydrodynamics schemes employing Riemann solvers, such as the scheme employed in Chimera, can introduce low-amplitude, post-shock oscillations in flows involving shocks unless extra dissipation is added. To suppress these, Chimera introduces a small amount of dissipation by reducing locally the order of the interpolation scheme in the neighborhood of sufficiently strong shocks (see Colella & Woodward 1984, Section 4 and Appendix). In particular, the left- and right-hand states are modified by Equations (4.1) of Colella & Woodward (1984), namely

Equation (193)

where the "flattening parameter" $(0\leqslant {f}_{i}\leqslant 1)$ determines the mixture of first-order and higher-order PPM interpolations in constructing the left and right states. For the test results shown in Figure 15(a), Chimera was run in Lagrangian mode with the Colella & Woodward (1984) suggested value of 0.75 for ω(1) in the equation following their Equation (A.2), and a maximum value of the flattening parameter, fi = 0.5. As can be seen, the agreement between the analytical and the numerical results is very good, but very-low-amplitude, post-shock oscillations are evident, particularly in the velocity. In Figure 15(b), we show the results of the same test except that the parameter ω(1) was set to 0.6, and the maximum value of the flattening parameter was set to 1. This introduces more dissipation, and there is now no evidence of any post-shock oscillations in the numerical solutions. Finally, Figure 15(c) shows the results of the test with Chimera now run in Eulerian mode with the Colella & Woodward (1984) suggested value of ω(1) and the maximum value of fi set to 0.5, as in the test shown in Figure 15(a). The agreement between the numerical and analytical results is again very good, and there is no sign of post-shock oscillations.

5.3. Shu–Osher Shock Tube Problem

A test suggested by Shu & Osher (1989) involves structure, testing the resolution of the numerical hydrodynamics scheme: A moving shock interacts with sine waves in density. Initially, ρ = 3.85713 ${{\rm{g\; cm}}}^{-3}$, v = 2.639369 ${{\rm{cm\; s}}}^{-1}$, and P = 10.33333 erg cm−3, for x < −0.8, and ρ = 1 + 0.2sin 5x ${{\rm{g\; cm}}}^{-3}$, v = 0 ${{\rm{cm\; s}}}^{-1}$, and P = 1 erg cm−3, for x > −0.8. The results are plotted in Figure 16. The black line shows the solution obtained with a grid of 3200 evenly spaced zones, which is taken as the reference solution. The dashed green line and the red "X"s show the solution obtained with a grid of 800 and 200 zones, respectively. Clearly, the solution has converged with a grid of 800 zones. With 200 zones, the solution still shows the detailed structure, albeit with somewhat reduced amplitude.

Figure 16.

Figure 16. Evolved density for the shock tube test suggested by Shu & Osher (1989), for a grid of 3200 (black solid), 800 (dashed green), and 200 (red with X) zones.

Standard image High-resolution image

5.4. Radial Advection Test

A test of the radial advection and, in particular, the geometry corrections for a spherical geometry, given the choice of such a geometry in Chimera, consists of solving the Euler equations with an adiabatic Γ-law EoS for an initial self-similar radial outflow problem (Mignone 2014) described by

Equation (194)

where ρ(ξ,0) is an arbitrary function and α0 is a constant. An exact analytic solution of this problem is given, for spherical symmetry, by

Equation (195)

A test that focuses on the remapping procedure is that of a self-similar outflow at constant density with the imposition of constant pressure, as described by Blondin & Lundqvist (1993) and Mignone (2014). Following their example, we set α0 = 100 and utilize a grid of 100 evenly spaced zones. The solution is characterized by constant velocity/radius and constant density. The solution of this problem for the first 10 zones, where the geometry-dependent corrections in the remap are most important, is shown in Figure 17. The numerical solution exhibits a constant velocity/radius that matches the analytic solution.

Figure 17.

Figure 17. Numerical solution for velocity/radius for the constant density radial outflow problem. The solid line is the analytic solution with velocity/radius = 99.24, and the numerical solution is shown by the square symbols for the first 10 zones, where the geometry-dependent interpolations are most important.

Standard image High-resolution image

5.5. Angular Advection Test

As a test of the angular advection algorithms and the logic involved in deflashing a zone (transitioning it from NSE to non-NSE) and flashing a zone (transitioning it from non-NSE to NSE), an electron fraction ${Y}_{{\rm{e}}}$ pattern was advected in angle (θ) across most of 256 angular zones, beginning at zone 18 (θ = 0.22) and continuing across the grid. The ${Y}_{{\rm{e}}}$ pattern consisted of a linear rise over 10 zones, from a value of 0.3472 to a value of 0.4960. A density of $2.47\times {10}^{7}\,{{\rm{g\; cm}}}^{-3}$ and a temperature of 0.45 MeV were chosen, which are representative of conditions where material flows in or out of NSE. Between zones 80 and 81 (θ = 0.98), a transition from NSE to non-NSE was imposed, and between zones 168 and 169 (θ = 2.1), a transition from non-NSE to NSE was imposed. As can be seen in Figure 18(a), the pattern in ${Y}_{{\rm{e}}}$ is nicely preserved as it is advected across the angular grid. A similar test with similar results performed for the azimuthal advection algorithms is shown in Figure 18(b).

Figure 18.

Figure 18. Advection of an electron fraction pattern across an NSE to non-NSE transition (solid vertical line) followed by a non-NSE to NSE transition (dashed vertical line) at several times (colored curves) in the (a) θ-direction and (b) ϕ-direction.

Standard image High-resolution image

5.6. Energy Conservation Test

When run in normal mode, the Chimera hydrodynamics does not impose total energy conservation, where the total energy is defined by the integral on the left-hand side of Equation (99); therefore, its ability to conserve total energy is a rigorous test of the hydrodynamics algorithms. As a test of Chimera's ability to conserve energy with the realistic EoS described in Section 3.2 and its numerical implementation, we performed two Newtonian hydrodynamics simulations initiated from a 15 ${M}_{\odot }$ Woosley & Heger (2007) progenitor and carried out for 2 s post-bounce, at which point the bounce shock had traversed 38,000 of our 43,000 km extent of the radial grid.

In the first simulation, referred to as NHpar, Chimera was run in its normal mode. Away from shocks, the specific internal energy was updated during the Lagrangian step by the first law, Equation (124), and remapped as described by Equation (132). In the vicinity of a shock, the specific total energy, defined immediately above Equation (121), was evolved during the Lagrangian step by Equation (121), Equation (134), and Equation (135), and remapped by an equation analogous to Equation (132), with the specific internal energy then extracted by Equation (136). In the second simulation, referred to as NHtot, the specific total energy was updated during the Lagrangian step and remapped for all zones, whether they were in the vicinity of a shock or not. Run in this mode, Chimera automatically conserves the total energy apart from the source term given by Equation (134), which is added to the total energy after the Lagrangian step. Total energy conservation during the Lagrangian step can be ascertained by multiplying Equation (121) by the mass ${V}_{i}^{n+\tfrac{1}{2}}{\rho }_{i}^{n+\tfrac{1}{2}}$ and noting that the terms involving the pressure cancel in pairs when summed over the zones, leaving only the surface terms, and noting also that the term ${f}_{\nu i}^{n}$ is zero for these hydrodynamics-only runs. During the remap step, the terms on the right-hand side of Equation (132) cancel in pairs after multiplying by the zone mass $\delta {{ \mathcal M }}_{i}^{n+1}$. Careful bookkeeping kept track of all nonphysical changes in the specific energy having no dynamical effect, such as occurs most importantly when material is advected between adjacent spatial zones characterized by different EoSs with slightly different energy zeros, as described in Section 4.4.3.

The results of this test are shown in Figure 19. The Lagrangian trajectories, at 0.025 ${M}_{\odot }$ intervals, for the NHtot simulation shown in Figure 19(a) were almost identical to those of the NHpar simulation and, thus, are not shown. The shock trajectories (red and dashed green lines) for the two simulations, shown in both panels of Figure 19 are essentially on top of each other. The total energy (blue and violet lines) is shown in Figure 19(b). The same arbitrary constant has been added to the energy of each simulation to bring the energies within the range of the energy ordinate. They are gratifyingly flat, with both showing a small blip at bounce, and the NHpar simulation (blue line) shows a small rise of about 1049 erg from 700 ms to the end of the simulation. In that simulation, there is also a small decline in total energy from bounce to about 700 ms, which we traced to the remapping of the velocity, thus conserving momentum by design, rather than remapping the square of the velocity, which would conserve kinetic energy by design. This slight decline in kinetic energy does not appear in the evolved total energy of the NHtot simulation (violet line). In this case, the specific total energy is the remapped quantity and therefore automatically conserved, modulo the small source term on the right-hand side of the specific total energy equation, Equation (134), as noted above. Figures 20(a) and 20(b) plot the evolution of the gravitational, internal, and kinetic energy components of the total energy for the simulations NHpar (dashed lines) and NHtot (solid lines). Figure 20(a) plots the evolution of the energy components of the two simulations but at a scale for which differences are not readily discernible. Figure 20(b) plots the last 200 ms of the two simulations at a much finer scale. The same arbitrary constant has been added to each pair of energy components of the two simulations so they fit within the energy range of the ordinate. It is seen that the energy differences are small, with most of the energy differences arising from the ∼3 × 1049 erg difference in the internal energy and the ∼1 × 1049 erg difference in the kinetic energy.

Figure 19.

Figure 19. Panel (a): Newtonian hydrodynamics core collapse simulations, NH${}_{\mathrm{par}}$ and NH${}_{\mathrm{tot}}$, as described in the text. The black lines show the Lagrangian trajectories of zones enclosing increments of 0.025 ${M}_{\odot }$ for the NH${}_{\mathrm{par}}$ simulation. Shock trajectories for simulations NH${}_{\mathrm{par}}$ and NH${}_{\mathrm{tot}}$ are shown by the solid red and dashed green lines, respectively. Panel (b): shock trajectories for the two simulations are plotted as in panel (a), and the total energy with arbitrary offset (but the same offset for the two simulations; right scale) for simulations NH${}_{\mathrm{par}}$ (blue) and NH${}_{\mathrm{tot}}$ (violet).

Standard image High-resolution image
Figure 20.

Figure 20. Panel (a): gravitational, internal, and kinetic energy components of the total energy for NH${}_{\mathrm{tot}}$ (solid lines) and NH${}_{\mathrm{par}}$ (dashed lines). The difference between the two simulations is too small to discern at the scale of a Bethe. Panel (b): gravitational, internal, and kinetic energy components of the final 200 ms for NH${}_{\mathrm{tot}}$ (solid lines) and NH${}_{\mathrm{par}}$ (dashed lines), with an arbitrary offset (but the same for corresponding pairs of energy components), shown at a scale for which differences are apparent.

Standard image High-resolution image

Finally, an important quantity relating to the explosion energy obtained in CCSNe simulations and frequently used in comparing simulation outcomes between groups is the "diagnostic energy" Ediag, which is the sum of the gravitational, thermal, and kinetic energy in each zone in turn summed over all zones for which the sum in that zone is positive. At 2500 ms from the initiation of core collapse, the diagnostic energy for both of the above models had essentially become constant and was 0.477 B for the NHtot simulation and 0.475 B for the NHpar simulation.

Both simulations, NHpar and NHtot, were performed on an adaptive radial grid of 720 zones, as described in Section 4.5. To ascertain whether this grid resolution results in a radially converged numerical solution, we have performed the NHpar simulation with 240 and 480 zones. The results are shown in Figure 21. It is clear that the numerical solutions have essentially converged at a radial grid resolution of 480 zones but not with a radial grid resolution of 240 zones.

Figure 21.

Figure 21. Panel (a): shock trajectory as a function of post-bounce time for the simulation of model NH${}_{\mathrm{par}}$ with 240, 480, and 720 radial zones. Panel (b): diagnostic energy (sum of the total kinetic, gravitational, and internal minus rest-mass energies) as a function of post-bounce time.

Standard image High-resolution image

5.7. Gravitational Potential Tests

To verify the accuracy of the gravitational potential expansions given in Sections 4.6.2 and 4.6.3, and to ascertain an appropriate maximum number of multipoles to use in simulations, we have compared with an analytic solution the results of the axisymmetric, Equations (147) and (148), and non-axisymmetric, Equations (187)–(189), expansions for the gravitational potential of a Maclaurin spheroid. We consider the spheroid described by Couch et al. (2013), viz., a spheroid of uniform density of ρ = 1 ${{\rm{g\; cm}}}^{-3}$ embedded in a background of vanishing density. The spheroid has a semimajor axis of 1 m, an eccentricity of 0.9, and is located in a spherical volume of radius 2 m. So that Equations (147), (148), and (187)–(189) can be used to compute the gravitational potential, we set up a spherical computational grid with the spheroid located at the center. A 720 × 240 grid is used for the axisymmetric multipole expansion, reflecting the grid resolution in some of our recent 2D simulations, and a grid of 540 × 180 × 180 is used for the non-axisymmetric multipole expansion, reflecting the angular resolution used in some of our 3D simulations. Where the boundary of the spheroid cuts through a given zone, the density of that zone is adjusted in proportion to the percentage of the zone interior to the boundary. The analytic solution that we use for a point interior to the spheroid is given by Equation (21) of Couch et al. (2013), and for a point exterior to the spheroid, we use the solution given by Equation (1) of Hofmeister et al. (2018).

The results of the comparisons are shown in Figure 22. The zone-weighted mean of the deviation from the analytic solutions of our multipole expansions decreases with max, the maximum multipole used, and is below 0.02 percent for max ≥ 10. The maximum deviation also decreases with max and is below 1% for max ≥ 10. The choice of max in a simulation is obviously a compromise between accuracy and computational time. Chimera simulations typically use a value of max = 10.

Figure 22.

Figure 22. Panel (a): maximum percentage deviation of the gravitational potentials computed by the axisymmetric multipole expansion (red line) and the non-axisymmetric multipole expansion (blue line) from the analytic solution, as a function of ${{\ell }}_{\max }$, the maximum multipole used. Panel (b): similar to Panel (a) but for the zone-weighted mean deviation of the multipole gravitational potentials from the analytic solution, as a function of ${{\ell }}_{\max }$.

Standard image High-resolution image

6. Neutrino Transport

Neutrino transport is a key process that must be modeled accurately in the simulation of CCSNe. The neutrino-driven explosion mechanism depends sensitively on the coupling to matter of a small fraction of the enormous neutrino luminosity that ensues upon the collapse of the stellar core. Additionally, accurate neutrino transport modeling is important for computing the neutrino emission expected from a given nearby CCSNe to best enable us to work backward from the sequence of neutrino detections accompanying such an event to establishing details of the explosion mechanism in the deep interior of the stellar core. In this section, we describe our algorithms for modeling neutrino transport. Sections 6.16.8 provide (i) the derivation of the neutrino Boltzmann equation, which establishes the metric and independent variables we use and forms the basis of our transport scheme, (ii) the angular moment equations obtained from the Boltzmann equation, (iii) our method of flux limiting, (iv) of operating splitting the resultant transport equations into a transport piece and an energy advection piece, and (v) the derivation of the terms required for coupling neutrinos to the matter hydrodynamics. Finally, in Section 6.9 by means of Equations (273), (274), (278) followed by Equations (277), (291), and (292), we present the full differencing scheme used to advance the neutrino transport through a Lagrangian step. A number of velocity-dependent terms are dropped in closing the angular moment equations by replacing the first angular moment equation by a diffusion-like equation with a flux-limiter. This resulted in transport solutions near shocks that were less than satisfactory due to the neutrino distribution moments being defined relative to the fluid frame and, therefore, subject to the effects of the large velocity discontinuities there. Section 6.12 describes a modification of our transport scheme that more accurately models these discontinuities in the neutrino distribution moments in the presence of shocks. This modification replaces Equations (273) and (274) for the Lagrangian step by Equations (306) and (300), respectively.

Neutrino energy advection during transport is operator split from spatial advection and is described in Section 6.13. Our scheme for updating the neutrino zeroth angular moment during a Lagrangian step is detailed by Equations (322)–(325). Because our energy grid is tied to a lapse function (Equation (216)), Section 6.14 describes the use of the neutrino energy advection machinery developed in Section 6.13 to update the neutrino distribution due to changes in the lapse function resulting from changes in the configuration of the core during a hydrodynamics step. As our neutrino transport scheme is based on a Lagrangian-Remap formulation of numerical hydrodynamics, a transport remap step must follow the Lagrangian update. Section 6.15 describes our scheme for spatially remapping the neutrino distributions that is associated with the remapping of the grid following a Lagrangian hydrodynamics step. The scheme is summarized by Equations (338) and (339). Section 6.16 specifies the scalar Eddington factors used to represent higher angular moments of the neutrino distributions in terms of lower moments, the former of which appear in the equations for the lowest moments.

The derivation of the neutrino transport equations and the energy–momentum transfer between neutrinos and the background matter have been carried out, with an eye toward future developments of Chimera, in full spherically symmetric general relativity. Several approximations have been made, however, in the current versions of Chimera. Some velocity terms have been dropped in deriving the diffusion-like equation and the flux-limiter relating the first moment of the neutrino distribution function to the zeroth moment (Equations (224) and (225)). Additionally, the proper length parameter "Γ," Equation (249), has been set equal to unity. It was found that retaining Γ as computed by Equation (249) does not significantly affect the neutrino transport during a CCSN simulation (Bruenn et al. 2001). What we do retain is the redshift parameter "a" as given by Equations (262) and (263), as this parameter does significantly affect the luminosity and mean energy of neutrinos emerging from the neutrinosphere (Bruenn et al. 2001).

Chimera employs multi-neutrino energy transport, and to delineate the neutrino energy grid structure, we use indices $k+\tfrac{1}{2},\ k=0,1,2,\cdots ,{N}_{\epsilon }$—i.e., half-integer values to denote energy-grid zone edges. Energy-grid zone centers are denoted by indices $k,\ k=1,2,\cdots ,{N}_{\epsilon }$—i.e., with integer values. With ${\epsilon }_{k+\tfrac{1}{2}}$ denoting the value of the neutrino energy at zone edges, the value of the neutrino energy at zone centers is defined by

Equation (196)

so that

Equation (197)

where ${\rm{\Delta }}{\epsilon }_{k}={\epsilon }_{k+\tfrac{1}{2}}-{\epsilon }_{k-\tfrac{1}{2}}$. With this definition of epsilonk, $4\pi {\epsilon }_{k+\tfrac{1}{2}}^{2}{\rm{\Delta }}{\epsilon }_{k}$ is the energy-space volume between zone edges ${\epsilon }_{k-\tfrac{1}{2}}$ and ${\epsilon }_{k+\tfrac{1}{2}}$.

6.1. Boltzmann Equation

The time evolution of the invariant occupation function $f=f(x,p)$ [number of neutrinos per state at the phase point $(x,p)]$ is given by the coordinate invariant Boltzmann equation (e.g., Lindquist 1966)

Equation (198)

where is the affine path-length, which we choose to be defined such that

Equation (199)

and the right-hand side of Equation (198) denotes the change in f due to sources—i.e., emission, absorption, and scattering. In the spirit of the RbR approximation in which transport along each radial ray is assumed to be spherically symmetric, we will consider the spherically symmetric form of the transport equation. In addition, we will express the Boltzmann equation in the comoving or fluid frame with the intention of executing our neutrino transport along with the Lagrangian radial hydrodynamics step, to be followed by a remap of both the matter and the neutrino quantities to the original grid or to the grid displaced according to a regridding algorithm. We will hereafter denote all quantities defined with respect to the fluid frame by a subscript "0."

The derivation, leading to Equation (215), of the Boltzmann transport equation for spherically symmetric spacetimes, beginning with Equation (198), can be found with varying detail in a number of references (e.g., Lindquist 1966; Castor 1972; Mihalas & Mihalas 1984; Baron et al. 1989; Mezzacappa & Matzner 1989), and we include enough detail so that important quantities used subsequent to Equation (215) are clearly defined. We assume that neutrinos follow null geodesics between localized interactions, so that their paths between interactions are given by

Equation (200)

where $\left\{\begin{array}{c}\alpha \\ \beta \ \ \gamma \end{array}\right\}$ is a Christoffel symbol of the second kind (connection coefficients in the coordinate basis), usually denoted by ${{\rm{\Gamma }}}_{\beta \gamma }^{\alpha }$, but we reserve the latter symbol to denote the Ricci rotation coefficients (connection coefficients in the orthonormal basis).

To evaluate the source functions, we choose a local, comoving, orthonormal set of basis vectors (${{\boldsymbol{e}}}_{t},{{\boldsymbol{e}}}_{m},{{\boldsymbol{e}}}_{\theta },{{\boldsymbol{e}}}_{\phi }$) parallel to a spherical polar coordinate basis to resolve the components of the neutrino four-momentum. In terms of this orthonormal set of basis vectors, the components of four-vectors, such as the four-momentum, ${p}_{0}^{\hat{{\rm{a}}}}$, are denoted by characters with hatted Latin indices. In terms of the original components, ${p}_{0}^{\alpha }$, the components with respect to the orthonormal basis are given by

Equation (201)

Substituting the second part of Equation (201) into Equation (200) and rearranging the indices, ${{dp}}^{\hat{{\rm{a}}}}/d{\ell }$ is given by

Equation (202)

where the Ricci rotation coefficients, ${{\rm{\Gamma }}}_{\hat{{\rm{a}}}\hat{{\rm{c}}}}^{\hat{{\rm{b}}}}$, are defined by Equation (202). Using Equations (199), (201), and (202), Equation (198), written in terms of the ${p}^{\hat{{\rm{a}}}}$, becomes

Equation (203)

The components of ${p}_{0}^{\hat{{\rm{a}}}}$ with respect to the above orthonormal basis are given by

Equation (204)

where ${\epsilon }_{0}$ and ${\mu }_{0}$ are the neutrino energy and the direction cosine of the neutrino three-momentum with respect to the radial direction, respectively, both measured in the comoving frame. Because of the mass shell condition

Equation (205)

only three of the four ${p}_{0}^{\hat{a}}$'s are independent. We choose these independent components to be ${p}_{0}^{\hat{{\imath }}},\hat{{\imath }}=1,2,3$. In terms of these independent four-momentum components, the derivatives, $\partial /\partial {p}_{0}^{\hat{{\imath }}}$ are given by

Equation (206)

We choose a synchronous gauge with a spherically symmetric, orthogonal metric given by

Equation (207)

where ${x}^{0}={ct}$, ${x}^{1}=m$, ${x}^{2}=\theta $, and ${x}^{3}=\phi $. With this metric, the transformation functions relating coordinate and orthonormal bases are given by

Equation (208)

We choose this form of the metric so that the relativistic equations will closely parallel the Newtonian fluid equations, and various levels of Newtonian approximations can be easily made. The metric function $a=a(t,m)$ is the lapse function and relates the interval of proper time of an observer attached to the motion of a fluid element to an interval of coordinate time, and defined so that coordinate and proper time are equal at infinity. The metric function $b=b(t,m)$ will be chosen so that the coordinate m can be identified with the enclosed rest mass. The metric function R is the areal radius (i.e., the two-sphere area $=4\pi {R}^{2}$). The condition that ties the coordinate system to the comoving frame is that the four-velocity of the fluid, ${u}^{\nu }$, be given by ${u}^{\nu }\equiv {{dx}}^{\nu }/d\tau \equiv c\,{{dx}}^{\nu }/{ds}=[{u}^{0},0,0,0]$, which with the metric given by Equation (207) requires that ${u}^{0}=c/a$. This definition of ${u}^{\nu }$ implies that

Equation (209)

where u is the first component of the four-velocity as observed from a frame of constant areal radius R (May & White 1967). To specify b so that m can be identified with the enclosed rest mass, we note that the rest mass density ρ satisfies the local conservation condition

Equation (210)

where the semicolon denotes covariant differentiation and where we have used the expression for ${u}^{\nu }$ given immediately above Equation (209).

With the metric (Equation (207)), the proper volume, dV, is given by ${dV}=4\pi {R}^{2}b\,{dm}$, or, in terms of the rest mass dm contained in dV, by ${dV}={dm}/\rho $. It follows from these two expressions for dV that the requisite choice of b is

Equation (211)

With the above choices for a spherical, comoving coordinate system and a comoving, orthonormal four-vector basis, the coordinate invariant volume elements become

Equation (212)

Equation (213)

where g is the determinant of the metric, and ${\epsilon }_{\delta \alpha \beta \gamma }$ and ${\epsilon }_{{ijk}}$ are the Levi-Civita alternating symbols. The invariant distribution function, ${f}_{0}$, introduced at the beginning of this section, is defined so that the number of world lines crossing the volume element, dV, with four-momenta in the range dP about ${{\boldsymbol{p}}}_{0}$ is given by (Lindquist 1966)

Equation (214)

where the right-hand sides of Equations (212)–(214) are the expressions for dV, dP, and dN in our choice of coordinate system and four-vector basis. Finally, evaluating the transformation coefficients from the metric (Equation (207)) and using them for the Ricci coefficients in Equation (203) and using Equations (204) and (206), the Boltzmann equation becomes

Equation (215)

Note that with the substitutions $a={e}^{{\rm{\Phi }}}$, $b={e}^{{\rm{\Lambda }}}$, and ${\rm{\Gamma }}={(1/b)(\partial R/\partial m)}_{t}$, Equation (215) reduces to Equation (3.7) of Lindquist (1966). For economy of notation and for clarity of presentation, we have suppressed the dependency of ${f}_{0}$ on t, m, ${\epsilon }_{0}$, and ${\mu }_{0}$, and the dependencies of the metric variables a, b, Γ, and R on t and m, and will do so with new dependent variables as they are introduced as long as this brings no ambiguity.

We now modify Equation (215) by transforming from independent variables $(t,m,{\epsilon }_{0},{\mu }_{0})$ to $(t,m,{E}_{0},{\mu }_{0})$ where

Equation (216)

Therefore, with this variable transformation

Equation (217)

Though mathematically imprecise, we will use the same symbol to describe the neutrino distribution functions, originally functions of $(t,m,{\epsilon }_{0},{\mu }_{0})$, as functions of $(t,m,{E}_{0},{\mu }_{0})$. We also define Γ by ${\rm{\Gamma }}={(1/b)(\partial R/\partial m)}_{t}$ so that derivatives with respect to m at constant time can be replaced by derivatives with respect to R, which is directly related to our grid, by the identity

Equation (218)

Applying the above definition of Γ and variable transformation to Equation (215) gives

Equation (219)

The transformation given by Equation (216) affords several advantages. First, with constant values of ${E}_{0k+\tfrac{1}{2}}$ anchoring the energy grid, the energy grid ${\epsilon }_{0k+\tfrac{1}{2}}$ will be scaled to higher values as the lapse dips below unity at high densities—i.e., for constant ${E}_{0k+\tfrac{1}{2}}$, the neutrino energy will scale as ${\epsilon }_{0k+\tfrac{1}{2}}\propto 1/a$. This will permit a smaller upper bound to the neutrino energy grid farther out from the center where $a\simeq 1$, resulting in the grid energies being more closely spaced there for a given number of grid points, while still permitting sufficient energy headroom at high densities to accommodate the high-energy neutrinos that are produced at the high-matter densities prevailing near the core center. Second, the radial derivative of a in the factor multiplying the energy derivative of ${f}_{0}$ in Equation (215) has been replaced by a time derivative in Equation (219). This factor now contains terms involving only time derivatives and, therefore, vanishes for a static spacetime. Thus, apart from energy-changing interactions, for a static spacetime, there will be no flow of neutrinos through the neutrino energy grid as they propagate outward. The gravitational redshifting will consequently be accomplished automatically. In non-static spacetimes, the advection of neutrinos through the energy grid can be performed algebraically (Sections 6.13 and 6.14), another advantage of this scheme. Using this choice of energy gridding, we will assume that the spacetime is constant over a time step and make a small correction to the neutrino distribution at the end of the time step to correct for the change in a during a time step and other processes that shift the neutrinos in energy, such as their advection across spatial zones with differing lapses.

6.2. Moment Equations

Let the nth angular moment of ${f}_{0}$ be denoted by ${\psi }_{0}^{(n)}$, that is

Equation (220)

and let

Equation (221)

Then the first two angular moments of Equation (219) are given by

Equation (222)

and

Equation (223)

6.3. Flux Limiting

To close this system of moment equations, we employ flux limiting to derive a relation between ${\psi }_{0}^{(1)}$ and ${\psi }_{0}^{(0)}$. In analogy with the procedure described in Levermore & Pomraning (1981), let the scalar Eddington factors ${\eta }^{(n)}$ be defined as the ratio of ${\psi }_{0}^{(n)}$ to ${\psi }_{0}^{(0)}$ so that ${\psi }_{0}^{(n)}={\eta }^{(n)}{\psi }_{0}^{(0)}$. Substituting ${\psi }_{0}^{(n)}={\eta }^{(n)}{\psi }_{0}^{(0)}$ in Equations (222) and (223), solving Equation (222) for ${(\partial {\psi }_{0}^{(0)}/\partial t)}_{m,{E}_{0}}$ and substituting the result into Equation (223) gives

Equation (224)

The factors multiplying ${\psi }_{0}^{(0)}$ and ${(\partial {\psi }_{0}^{(0)}/\partial {E}_{0})}_{t,m}$ in Equation (224) in both the diffusion limit (${\eta }^{(1)}\to 0;$ ${\eta }^{(2)}\to {}^{1}{/}_{3};$ ${\eta }^{(3)}\to 0$) and the free streaming limit (${\eta }^{(n)}\to 1$) are zero. We therefore make the approximation that these two factors are zero everywhere. Equation (224) then becomes

Equation (225)

Multiplying the right-hand side of Equation (225) by ${\psi }_{0}^{(1)}$/${\psi }_{0}^{(1)}$, i.e., unity, and solving for ${\psi }_{0}^{(1)}$, we get

Equation (226)

In the diffusion limit, ${\eta }^{(1)}\to 0$ and ${\eta }^{(2)}\to {}^{1}{/}_{3}$, and Equation (226) reduces to the standard diffusion equation if we neglect the second term in the numerator and identify the transport mean free path, ${\lambda }^{(t)}$, as

Equation (227)

In regimes other than the diffusion regime, we regard $\left[{\eta }^{(2)}-{\left({\eta }^{(1)}\right)}^{2}\right]$ as a free parameter, which we write as

Equation (228)

Using Equations (227) and (228) in Equation (226), we get a diffusion-like equation for ${\psi }_{0}^{(1)}$:

Equation (229)

Equations (222) and (229) for each energy zone and for each neutrino species (${\nu }_{e}$, ${\bar{\nu }}_{e}$, ${\nu }_{\mu \tau }$, ${\bar{\nu }}_{\mu \tau }$) with a prescription for ${ \mathcal F }$ are the MGFLD equations that are solved in Chimera for the B-series simulations. The modification of this scheme in the vicinity of shocks, used for our C-series and later simulations, is described in Section 6.12. The parameter ${ \mathcal F }$ is referred to as the "flux-limiter," and should be unity in the diffusion limit and tend to zero in such a way that ${\psi }_{0}^{(1)}$  = ${\psi }_{0}^{(0)}$ in the limit of free streaming.

The derivation of Equations (229) and (227) differ from that of the corresponding equations (Equations (A25) and (A26)) in Bruenn (1985). To derive a diffusion-like equation, there the derivative $\partial {\psi }_{0}^{(1)}/\partial t$ in the nonrelativistic version of Equation (223) and all velocity dependent terms were set to zero. However, the results are similar. Apart from the relativistic time dilation factor a, the expressions for $1/{\lambda }_{i}^{(t)}$, given here by Equation (227) and given in Bruenn (1985) by Equation (A26), are the same. Equation (226) for ${\psi }_{0}^{(1)}$ here is similar to Equation (A25) ${\psi }_{0}^{(1)}$ in Bruenn (1985), with the exception here of the relativistic proper distance correction Γ and the factor of $\left[{\eta }^{(2)}-{\left({\eta }^{(1)}\right)}^{2}\right]$, which we take as our flux limiter. The flux limiter in Bruenn (1985) is introduced as a modification of ${\lambda }_{i}^{(t)}$. Finally, the terms in the numerator of Equations (226) here and in Equation (A25) of Bruenn (1985) are neglected, other than the term involving $\partial {\psi }_{0}^{(0)}/\partial R$. They involve the redistribution of neutrinos in angle due to anisotropies in the source terms but are small given that they depend on the product of ${\psi }_{0}^{(0)}$ and ${\psi }_{0}^{(1)}$. Equation (229) is given in differenced form by Equation (274).

6.4. Flux Limiter

The flux limiter, ${ \mathcal F }$, constructed to satisfy the diffusion and free-streaming limits, consists of two parts. The first part is a specific implementation of the usual scheme for interpolating between these two limits, namely

Equation (230)

From Equations (229) and (230), we see that in the optically thick, diffusion regime, for which ${\lambda }^{(t)}({\epsilon }_{0})\to 0$, ${{ \mathcal F }}_{\mathrm{intrp}}\to 1$ and Equation (229) limits to the diffusion equation

Equation (231)

while in the optically thin, free-streaming regime, where ${\lambda }^{(t)}({\epsilon }_{0})\to \infty $, Equation (229) limits to the free-streaming condition

Equation (232)

As it stands, this scheme suffers from the generic problem of an overly rapid transition to the free-streaming limit (i.e., the angular distribution becomes too forwardly peaked) when matter goes from optically thick to optically thin abruptly, such as when the "density cliff" forms in the post-bounce core of a supernova progenitor. To avoid this problem, a second piece of the flux limiter is constructed. It prevents the neutrino angular distribution from becoming more forwardly peaked than the geometrical limit. This geometrical limit can be expressed for $R\gt {R}_{\nu }$, where ${R}_{\nu }\equiv {R}_{\nu }({\epsilon }_{0})$ is the radius of the neutrinosphere for a particular flavor and energy, by ${\psi }_{0}^{(1)}=\tfrac{1}{2}[1+{\mu }_{0\nu }({\epsilon }_{0})]{\psi }_{0}^{(0)}$. This expression assumes that the neutrino distribution function ${f}_{0}$ is constant for rays satisfying ${\mu }_{0}\lt {\mu }_{0\nu }$, and zero otherwise, where ${\mu }_{0\nu }$ is the cosine of the angle subtended between a line extending from a point at R to the neutrinosphere limb and a line extending from the point at R to the core center, as observed by a comoving observer. This geometrical piece of the flux-limiter is then given by

Equation (233)

and ${{ \mathcal F }}_{\mathrm{geom}}(R,t,{\epsilon }_{0})=1$ interior to the neutrinosphere, $R\leqslant {R}_{\nu }({\epsilon }_{0})$. The comoving angle ${\mu }_{0\nu }$ is given in terms of the fixed frame angle ${\mu }_{\nu }$ by

Equation (234)

where $\beta =v/c$, and the fixed-frame angle ${\mu }_{\nu }$ is given by

Equation (235)

where ${ \mathcal G }$ is the correction for gravitational ray bending given by

Equation (236)

where Mg is the gravitational mass. The geometrical flux limiter, ${{ \mathcal F }}_{\mathrm{geom}}$, by itself would relate ${\psi }_{0}^{(1)}$ to ${\psi }_{0}^{(0)}$ outside the neutrinosphere by

Equation (237)

The final flux limiter is given by

Equation (238)

6.5. Operator Splitting

To solve Equation (222) along with Equation (229), we operator split Equation (222) into a transport equation and an energy advection equation

Equation (239)

where $(\partial {\psi }_{0}^{(0)}/\partial t{)}_{m,{E}_{0}}^{T}$, the rate of change of ${\psi }_{0}^{(0)}$ due to sources and transport, is given by

Equation (240)

and $(\partial {\psi }_{0}^{(0)}/\partial t{)}_{m,{E}_{0}}^{E}$, the rate of change of ${\psi }_{0}^{(0)}$ due to energy advection, is given such that

Equation (241)

In advancing ${\psi }_{0}^{(0)}$ and ${\psi }_{0}^{(1)}$ over a time step, the pair of Equations (240) and (229) are solved in one step along with the associated energy and lepton conservation equations introduced below, then Equation (241) is solved in a second step. We will refer to these two separate steps as a transport step and an energy advection step, respectively, and describe in more detail below how each is performed.

6.6. Einstein's Equations

In order to derive the Einstein equations, which are needed to obtain expressions for the gravitational mass, Mg, and the metric parameters Γ and a, we need the stress-energy tensor, ${ \mathcal T }\,=$ ${}^{(m)}{ \mathcal T }\,+$ ${}^{(\nu )}{ \mathcal T }$, where ${}^{(m)}{ \mathcal T }$ and ${}^{(\nu )}{ \mathcal T }$ are the matter and neutrino contributions, respectively. Following Mihalas & Mihalas (1984), we begin with the definition of the radiation (i.e., neutrino) stress-energy tensor

Equation (242)

where the sum q is over all neutrino species. In the local comoving orthonormal frame, the components of the stress-energy tensor are

Equation (243)

where, as before, we use Latin indices with hats here to distinguish components with respect to the local comoving orthonormal frame from those with respect to the coordinate basis (distinguished using Greek letters).

Using Equations (204) for ${p}^{\hat{{\rm{a}}}}$ and ${p}^{\hat{{\rm{b}}}}$ in Equation (243) for ${}^{(\nu )}{{ \mathcal T }}^{\hat{{\rm{a}}}\hat{{\rm{b}}}}$, the radiation stress-energy tensor ${}^{(\nu )}{{ \mathcal T }}^{\hat{{\rm{a}}}\hat{{\rm{b}}}}$ can be written in terms of the local neutrino energy density, ${E}_{\nu }$, flux ${F}_{\nu }$, and pressure, ${P}_{\nu }$, as

Equation (244)

Transforming ${}^{(\nu )}{{ \mathcal T }}^{\hat{{\rm{a}}}\hat{{\rm{b}}}}$ written in terms of ${E}_{\nu },{F}_{\nu }/c,{P}_{\nu }$ back to the coordinate basis using Equations (208), and adding the stress-energy components of a perfect fluid, given by

Equation (245)

with ${u}^{\alpha }$ defined immediately above Equation (209) and ${g}_{\alpha \beta }$ given by Equation (207), the combined matter–neutrino stress energy tensor is given by

Equation (246)

where ρ is the proper rest mass density, and Em and Pm are the matter internal energy density and pressure, respectively.

The Einstein field equations are given by

Equation (247)

where ${{ \mathcal R }}_{\nu }^{\ \ \mu }$ is the Ricci tensor. After some algebra (e.g., May & White 1966, 1967), we find from the Einstein equations that the gravitational mass is given by

Equation (248)

and the metric parameter Γ is given by

Equation (249)

where u is the velocity, $u={a}^{-1}{({dR}/{dt})}_{m}$. Einstein equations can also be used to derive the radial equation of motion, which is given by

Equation (250)

where w, the specific enthalpy, is given by

Equation (251)

and where the sum over ν is a sum over all neutrino and antineutrino species. Equation (250), with $w={\rm{\Gamma }}=a=1$ (Newtonian approximation), and with the centrifugal term added, is the radial Equation (113) and, in differenced form, Equation (114). Here, we consider the neutrino component of these equations; i.e., the last term of Equation (250), which corresponds to the term ${f}_{\nu r}$ in the radial Equation (113). The neutrino contributions to the θ- and ϕ-components of the velocity and energy hydrodynamics equations are described in the text following Equations (113) and in Equations (115), (122), and (123).

Using Equation (199) and the first of Equation (208), the last term of Equation (250) in the Newtonian approximation can be written as

Equation (252)

The differenced form of the expression given by the last term in Equation (252) for ${f}_{\nu ,r}$ is given in Section 6.9.

6.7. Matter–Neutrino Energy–Momentum Exchange

To determine the energy–momentum exchange between the matter and neutrinos, we begin with the hydrodynamics equation

Equation (253)

where, as before, the semicolon denotes covariant differentiation. The left-hand side is the divergence of the stress-energy tensor of the matter, and ${G}^{\alpha }$ is the four-force density—i.e., the negative of the matter to neutrino energy–momentum transfer rate per unit volume. To determine the latter, we operate on Equation (215), or Equation (219), by the negative of the four-momentum density operator (the integral of the product of the invariant momentum volume element and the neutrino four-momentum), which, in the orthonormal basis, is given by

Equation (254)

and use Equations (244) and (212) to get

Equation (255)

Transforming Equation (255) to the coordinate basis, we have for ${G}^{\alpha }$,

Equation (256)

The energy equation is obtained by projecting Equation (253) along the fluid four-velocity:

Equation (257)

where we have used Equation (256) for ${G}^{\alpha }$. Using Equation (245) for ${}^{(m)}{T}^{\alpha \beta }$ on the left-hand side of Equation (253), invoking conservation of rest mass $[{u}_{;\beta }^{\beta }=-{u}^{\beta }{\rho }_{,\beta }/\rho $], and denoting by ${e}_{{\rm{m}}}$ the internal energy per unit rest mass (${e}_{{\rm{m}}}={E}_{{\rm{m}}}/\rho $), ${u}^{\alpha (m)}{T}_{\alpha ;\beta }^{\beta }$ is given by

Equation (258)

and Equation (257) becomes

Equation (259)

which states that changes in the matter internal energy arise from work by local compression or expansion and from energy exchange with neutrinos.

The other nontrivial equation in (253) is the radial equation

Equation (260)

With ${}^{(m)}{T}_{\ \ \ ;\beta }^{1\beta }$ given by

Equation (261)

Equation (260) gives

Equation (262)

where w is given by Equation (251).

Equation (262) relates the spatial gradient of the matter pressure, ${P}_{{\rm{m}}}$, and the rate of neutrino–matter momentum exchange to the spatial gradient of the lapse function. Upon integrating a from the surface using Equation (262), we start with the surface (subscript s) boundary condition, which we take to be the Schwarzschild solution on the exterior, and neglect the small deviations from this induced by the radiation; namely,

Equation (263)

so that the coordinate time is that of a distant clock. In practice, our models are extended enough that as is extremely close to unity. The differencing of the lapse is given below by Equations (294) and (295).

6.8. Matter–Neutrino Lepton Exchange

The local lepton number density of the matter, ${}^{(m)}{n}_{{\ell }}$, is given from charge conservation by

Equation (264)

where ${n}_{{\rm{p}}}$, ${n}_{{\rm{B}}}$, ${n}_{{{\rm{e}}}^{-}}$, and ${n}_{{{\rm{e}}}^{+}}$ are the proper proton, baryon, electron, and positron number densities, respectively, ${Y}_{{\rm{e}}}$ is the proton fraction, and ${m}_{{\rm{B}}}$ is the mean baryon mass. Referring to Equation (214), the local lepton number density in neutrinos, ${}^{(\nu )}{n}_{{\ell }}$, is given by

Equation (265)

Transport produces global changes in ${}^{(\nu )}{n}_{{\ell }}$, which, by itself, does not change ${}^{(m)}{n}_{{\ell }}$ (or ${Y}_{{\rm{e}}}$). Ignoring the effects of transport, the total lepton number, ${}^{\mathrm{total}}{n}_{{\ell }}{=}^{(m)}{n}_{{\ell }}{+}^{(\nu )}{n}_{{\ell }}$, is a locally conserved quantity, and this allows us to relate the change in ${Y}_{{\rm{e}}}$ to the change in ${}^{(\nu )}{n}_{{\ell }}$. In particular,

Equation (266)

Expanding the covariant derivative of ${Y}_{{\rm{e}}}$ with the use of Equation (210), which expresses rest mass conservation, Equation (266) can be written

Equation (267)

Now expanding the covariant derivative of ${}^{(\nu )}{n}_{{\ell }}$, Equation (267) becomes

Equation (268)

With the use of Equations (215) (without the transport terms), (199), (211), and (214), the term on the right-hand side of Equation (268) becomes

Equation (269)

where the last expression in Equation (269) arises from the fact that the only terms that do not vanish upon integrating over solid angle and neutrino number are emission and absorption, and the latter are isotropic. Here, ${j}_{{\nu }_{{\rm{e}}}}({\epsilon }_{0})$ and $1/{\lambda }_{{\nu }_{e}}^{a}({\epsilon }_{0})$ are the electron-neutrino emission and absorption per state per unit length, and ${j}_{{\bar{\nu }}_{e}}({\epsilon }_{0})$ and $1/{\lambda }_{{\bar{\nu }}_{e}}^{a}({\epsilon }_{0})$ are the electron–antineutrino counterparts. Equations (266)–(269) give the following equation for the evolution of the electron fraction due to sources:

Equation (270)

6.9. Lagrangian Transport Step

In the Lagrangian transport step, the pair of transport Equations (240) and (229) (or, rather, Equation (272) in place of Equation (240), where Equation (272) is Equation (240) written in conservative form) are solved along with the neutrino-specific term of the energy Equation (259), namely,

Equation (271)

and Equation (270) for the change in the matter proton number. The pair of Equations (240) and (229) for each neutrino species and Equation (270) are solved together fully implicitly, while Equation (271) is solved immediately afterwards, as described in more detail below. It was found that executing Equation (271) subsequent to the others leads to a substantially better conditioned set of equations and no loss of stability and that the internal energy update could then be easily omitted for those radial zones for which the internal energy update is derived from a total energy update. The role of neutrinos in the latter case will be described later. In either case, the updated temperature is obtained from the updated internal energy and electron fraction by numerically inverting the equation of state.

Equation (240), written in conservative form, is given by

Equation (272)

where we have omitted the superscript T. Consider a radial mesh where, as before, the zone edges are labeled by $i+{}^{1}{/}_{2}$ and the zone centers by i. Consider also an energy mesh labeled similarly by k. Let the subscript q denote, as before, the neutrino species, and let the superscript n denote the nth time slice. Then, for each neutrino species, Equation (272) is differenced conservatively as follows:

Equation (273)

with ${\psi }_{0\,i+\displaystyle \frac{1}{2},k,q}^{(1)\,n+1}$ given by

Equation (274)

(or by its modification given by Equation (300) below), where ${\mathrm{Area}}_{i+\tfrac{1}{2}}=4\pi {R}_{i+\tfrac{1}{2}}^{2}$ is the proper area of the outer boundary of zone i, ${R}_{i}=({R}_{i-\tfrac{1}{2}}+{R}_{i+\tfrac{1}{2}})/2$, ${\mathrm{Vol}}_{i}=4\pi \left({R}_{i+\tfrac{1}{2}}^{3}-{R}_{i-\tfrac{1}{2}}^{3}\right)/3{\rm{\Gamma }}$ is the proper volume enclosed between zones edges $i-{}^{1}{/}_{2}$ and $i+{}^{1}{/}_{2}$, ${\rm{\Delta }}{t}^{n+\tfrac{1}{2}}={t}^{n+1}-{t}^{n}$, and ${\lambda }_{i+\tfrac{1}{2},k,m}^{(t)\,n+1}=\sqrt{{\lambda }_{i,k,m}^{(t)\,n+1}{\lambda }_{i+1,k,m}^{(t)\,n+1}}$. To demonstrate that Equation (273) is conservative in neutrino number, first multiply it by $c\,{\mathrm{Vol}}_{i}\,{\rm{\Delta }}{t}^{n+\tfrac{1}{2}}/{a}_{i}^{2}$ to rewrite it as

Equation (275)

According to the discussion following Equation (219) and our choice of the lapse function limit, $a(R,t)\to 1$ as $R\to \infty $, a neutrino propagating along a geodesic passing through ${R}_{i-\tfrac{1}{2}}$, Ri, and ${R}_{i+\tfrac{1}{2}}$ in a static spacetime will have locally measured energies at these points related by

Equation (276)

Thus, multiplying Equation (275) by ${\left({hc}\right)}^{-3}{\left({\epsilon }_{0k\infty }\right)}^{2}{\rm{\Delta }}{\epsilon }_{0k\infty }$ and using Equation (276), we get

Equation (277)

The left-hand side of Equation (277) is the change in the number of neutrinos in ${\mathrm{Vol}}_{i}$ between energies ${\epsilon }_{0i,k}$ and ${\epsilon }_{0i,k}+{\rm{\Delta }}{\epsilon }_{0i,k}$ during the proper time interval ${a}_{i}{\rm{\Delta }}{t}^{n+\tfrac{1}{2}}$. The first term on the right-hand side subtracts the number of neutrinos that leave through the outer surface of ${\mathrm{Vol}}_{i}$ between energies ${\epsilon }_{0i+\tfrac{1}{2},k}$ and ${\epsilon }_{0i+\tfrac{1}{2},k}+{\rm{\Delta }}{\epsilon }_{0i+\tfrac{1}{2},k}$ in proper time ${a}_{i+\tfrac{1}{2}}{\rm{\Delta }}{t}^{n+\tfrac{1}{2}}$. The second term on the right-hand side adds the number of neutrinos that enter through the inner surface of ${\mathrm{Vol}}_{i}$ between energies ${\epsilon }_{0i-\tfrac{1}{2},k}$ and ${\epsilon }_{0i-\tfrac{1}{2},k}+{\rm{\Delta }}{\epsilon }_{0i-\tfrac{1}{2},k}$ in proper time ${a}_{i-\tfrac{1}{2}}{\rm{\Delta }}{t}^{n+\tfrac{1}{2}}$. The last term is the net number of neutrinos emitted, absorbed, or scattered in/out of ${\mathrm{Vol}}_{i}$ between energies ${\epsilon }_{0i,k}$ and ${\epsilon }_{0i,k}+{\rm{\Delta }}{\epsilon }_{0i,k}$ in proper time ${a}_{i}{\rm{\Delta }}{t}^{n+\tfrac{1}{2}}$.

Equation (270) is differenced straightforwardly as

Equation (278)

During the Lagrangian transport step, the set of Equations (273) and (274), modified in the presence of shocks as given by Equations (300), (306), and (278) are solved implicitly for ${Y}_{{\rm{e}},i}^{n+1}$ and ${\psi }_{0\,i,k,q}^{(0)\,n+1}$ for each neutrino species q.

The temperature change accompanying the Lagrangian transport step is computed following the solutions of Equations (273), (274), and (278), from the change in the internal energy as given by the following differenced version of Equation (271):

Equation (279)

and the change in ${Y}_{{\rm{e}}}$ given by Equation (278), where ${N}_{\nu }$ is the number of neutrino species (typically four), and ${{ \mathcal N }}_{\mathrm{stwt},q}$ is the statistical weight of each neutrino species, q, (typically 1 for ${\nu }_{e}$ and ${\bar{\nu }}_{e}$, and 2 for ${\nu }_{\mu \tau }$ and ${\bar{\nu }}_{\mu \tau }$). Specifically,

Equation (280)

It was found that the temperature change during a time step always had very little effect on the terms in the transport equation, and therefore on the stability of the difference scheme. On the other hand, including the temperature change as part of the implicit solution of Equations (273), (274), and (278) caused the system of equations to be rather ill-conditioned. It was therefore deemed numerically expedient to solve for the temperature change after, rather than simultaneously with, the solution of the transport equations. The solution of Equations (273), (274), and (278) followed by Equation (279) and (280) completes the Lagrangian transport step.

6.10. Solution of the Transport Equations

The numerical solution of the set of transport equations from time step n to time step $n+1$ proceeds by an outer iteration for the corrections $\delta {\psi }_{0\,i,k,q}^{(0)\,n,i+1}$ and $\delta {Y}_{{\rm{e}},i}^{n,i+1}$ to the quantities ${\psi }_{0\,i,k,q}^{(0)\,n,i}$ and $\delta {Y}_{{\rm{e}},i}^{n,i};$ i.e.,

Equation (281)

where the superscript "i" denotes the quantities after the ith iteration. When the corrections have become sufficiently small, the quantities ${\psi }_{0\,i,k,q}^{(0)\,n,i+1}$ and ${Y}_{{\rm{e}},i}^{n,i+1}$ are regarded as having converged to ${\psi }_{0\,i,k,q}^{(0)\,n+1}$ and ${Y}_{{\rm{e}},i}^{n+1}$, respectively. The criterion for convergence is that

Equation (282)

for all ${\psi }_{0\,i,k,q}^{(0)\,n,i+1}$ and ${Y}_{{\rm{e}},i}^{n,i+1}$, where typically ${\mathrm{tol}}_{\psi }={10}^{-6}$, ${\mathrm{tol}}_{{Y}_{{\rm{e}}}}={10}^{-6}$, ${\mathrm{tolmin}}_{\psi }=1$, and ${\mathrm{tolmin}}_{{Y}_{{\rm{e}}}}=0.1$.

The transport equations, Equations (273) and (274), when linearized comprise a set ${N}_{\nu }\times {N}_{k}$ coupled linear equations of the form:

Equation (283)

where the coefficients "LHS0" and "RHS0" are independent of the corrections, and where ${n}_{q}=({N}_{\nu }-1)\times q$, and the notation "LHS0" and "RHS0" denote the left-hand side (transport) and right-hand side (sources) of the transport equations, respectively. Equation (278) provides an additional equation, which is coupled to the preceding set of equations, and when linearized is of the form

Equation (284)

where q = 1 and q = 2 refer to ${\nu }_{e}$ and ${\bar{\nu }}_{e}$, respectively. Equations (283) and (284) are combined to form a matrix equation for the corrections, of the form

Equation (285)

where the corrections $\delta {\psi }_{0\,i,k,q}^{(0)\,n,i+1}$ and $\delta {Y}_{{\rm{e}},i}^{n,\ i+1}$ are incorporated in the array $\delta {u}_{i,\alpha },\alpha =1,\cdots ,\,{N}_{\nu }\times {N}_{k}+1$.

The boundary condition at the center is simply that the neutrino flux is zero; i.e., ${\psi }_{0\,i=1,k,q}^{(1)\,}=0$. Constant luminosity is assumed at the surface, which requires that ${\psi }_{0\,I+1,k,q}^{(1)\,}{R}_{I+1}^{2}={\psi }_{0\,I,k,q}^{(1)\,}{R}_{I}^{2}$. Following the discussion immediately above Equation (233), a geometric relation between ${\psi }_{0\,I,k,q}^{(1)\,}$ and ${\psi }_{0\,I,k,q}^{(0)\,}$ of the form ${\psi }_{0\,I,k,q}^{(1)\,}=(1/2)(1+{\mu }_{0I,k,q}){\psi }_{0\,I,k,q}^{(0)\,}$ is assumed. Then, without velocity or GR corrections to ${\mu }_{0I,k,q}$, the surface boundary condition for ${\psi }_{0\,I,k,q}^{(0)\,}$ becomes

Equation (286)

where ${R}_{\nu k,q}$ is the neutrinosphere radius for neutrinos of energy k and type q.

6.11. Neutrino Stress and the Lapse

Once ${\psi }_{0\,i,k}^{(0)\,n+1}$ has been obtained, the neutrino–matter stress ${f}_{\nu ,r}$ given by Equation (252) and used in Equations (114) and (121) for the radial hydrodynamics solve is computed as follows. Using Equation (227) and then (229) in Equation (252), we have

Equation (287)

where in the last equation, we have set Γ and w to their nonrelativistic values of unity. Equation (287) is differenced as

Equation (288)

The replacement

Equation (289)

has been made so the final expression for the neutrino stress will have the correct limiting behavior in the short neutrino mean-free-path limit (${{ \mathcal F }}_{i+\tfrac{1}{2},k,q}^{n+1}\to 1$)

Equation (290)

which is the force per unit mass exerted by an isotropic relativistic gas, and in the long neutrino mean-free-path limit which, with the help of Equation (229), is given by

Equation (291)

where we assume that at large R, where the neutrino mean free paths are large, ${\epsilon }_{0i+1,k}\simeq {\epsilon }_{0i,k}$. As before, ${\rm{\Delta }}{M}_{i}$ denotes the zone mass. Finally, the zone-centered neutrino stress is computed as the average of the edge values, viz.,

Equation (292)

To update the lapse given by Equation (262), note first that a comparison of Equation (262) with Equation (252) shows that

Equation (293)

Setting ${\rm{\Gamma }}=1$, ${M}_{{\rm{g}}}$ [given by Equation (248)] = $M\equiv {M}_{\mathrm{rest}\mathrm{mass}}$, and w = 1, the update of the lapse begins at the outer edge with

Equation (294)

then inward to the center with

Equation (295)

where the zone-centered quantities defined at zone edges are arithmetic averages of their zone-centered values on either side of the zone edge, and ${\rm{\Delta }}{M}_{i+\tfrac{1}{2}}/4\pi {R}_{i+\tfrac{1}{2}}^{2}{\rho }_{i+\tfrac{1}{2}}$ has been used for ${R}_{i+1}-{R}_{i}$. Zone-edged values of a are defined by ${a}_{i+\tfrac{1}{2}}=({a}_{i}+{a}_{1+1})/2$.

6.12. Transport Through a Shock

Comparisons of transport tests with a Boltzmann solver revealed a shortcoming of our flux-limited diffusion scheme when encountering a discontinuity in the fluid velocity—namely, a shock. The problem arises from our neglect of the velocity-dependent terms in going from Equation (224) to Equation (225), which should be important for computing the transport at the shock. The reason was traced ostensibly to the fact that both ${\psi }^{(0)}$ and ${\psi }^{(1)}$, being defined with respect to the fluid frame, are therefore physically both discontinuous at a discontinuity in the fluid velocity. However, as a result of the fact ${\psi }^{(1)}$ is given by the gradient of ${\psi }^{(0)}$ (e.g., Equation (229) and its differenced version, Equation (274)) ${\psi }^{(0)}$ is constrained to be continuous across the shock. In the case, in which the fluid velocity ahead of the shock has a large negative radial velocity that discontinuously transitions to a much smaller negative post-shock radial velocity, which characterizes the change in the fluid velocity across the bounce shock as it stagnates and then revives, both ${\psi }^{(0)}$ and ${\psi }^{(1)}$ should exhibit, outwardly, a positive discontinuous change through the shock. However, ${\psi }^{(0)}$, rather than exhibiting this discontinuous behavior, ramps up over many zones to its final pre-shock value ahead of the shock. This behavior of ${\psi }^{(0)}$ is required in order for the required flux through the shock to be computed via Equation (274). This rise of ${\psi }^{(0)}$ through the post-shock region, occurring in the neutrino heating region in the case of the bounce shock, has the unfortunate effect of causing a rise in the mean neutrino energy

Equation (296)

The result is an unphysical additional neutrino heating in the region behind the bounce shock, potentially causing the shock to revive too soon.

We illustrate this problem by comparing 1D simulations performed using Chimera and the relativistic Boltzmann code Agile-Boltztran, which is described in Section 9. The progenitor used is the 15 ${M}_{\odot }$ model evolved to the point of core collapse by Woosley & Heger (2007), and the neutrino physics is that described in Bruenn (1985). The rather straightforward collection of neutrino physics used in this comparison is chosen so that its implementation in both codes can be easily made identical. The unphysical rise in the ${\nu }_{e}$-rms and ${\bar{\nu }}_{e}$-rms energies as computed by Chimera is illustrated by the green lines in Figure 23(a). Compared with the rms energies computed by Agile-Boltztran (plotted by the red lines), the rms energies computed by Chimera are 1–2 MeV higher by the time the neutrinos reach the shock. A consequence of this neutrino rms energy rise is illustrated in Figure 23(b) by the increase in the shock radius due to increased neutrino heating in the Chimera simulation (shown in green) as compared with that given by the Agile-Boltztran simulation (shown in red). The shock position in the Chimera simulations has been pushed about 15 km farther out in radius compared with its position as given by Agile-Boltztran.

Figure 23.

Figure 23. Panel (a): ${\nu }_{e}$-rms and ${\bar{\nu }}_{e}$-rms energies at ${t}_{\mathrm{bounce}}=100$ ms; Panel (b): shock radius as a function of post-bounce time using Chimera with (blue) and without (green) the shock transport algorithm, compared to a calculation using Agile-Boltztran (red); Panel (c): total energy as a function of time for the two Chimera simulations redone using Newtonian gravity; and Panel (d): ${E}_{\mathrm{grav}}$, total gravitational potential energy, ${E}_{\mathrm{internal}}$, total internal energy, ${E}_{\mathrm{ke}}$, total kinetic energy, ${E}_{\nu }$, total on-grid neutrino energy, and ${E}_{\mathrm{rad}}$, neutrino energy lost by neutrinos exiting the grid for the Chimera runs plotted in Panel (c). The sum of these energies is the total energy plotted in Panel (c).

Standard image High-resolution image

To describe the causes of this problem in more detail and our shock transport algorithm designed to eliminate this problem, consider Equation (274) at a velocity discontinuity situated at zone interface $i+{}^{1}{/}_{2}$. In this circumstance, ${\psi }_{0\,i+1,k}^{(0)\,n+1}$ and ${\psi }_{0\,i,k}^{(0)\,n+1}$ are defined in two different reference frames. In the case of the outwardly directed bounce shock, the former would be defined in the pre-shock frame, the latter in the post-shock frame. To avoid the unphysical ramping up of ${\psi }^{(0)}$ toward the shock in the post-shock frame as we approach the shock from below, ${\psi }_{0\,i,k,q}^{(0)\,n+1}$ is replaced by ${\psi }_{0\,i,k,q}^{(0)\,n+1,\uparrow }$ in Equation (274), where ${\psi }_{0\,i,k,q}^{(0)\,n+1,\uparrow }$ is ${\psi }_{0\,i,k,q}^{(0)\,n+1}$ transformed to the same reference frame as ${\psi }_{0\,i+1,k,q}^{(0)\,n+1}\ ({\epsilon }_{0i});$ i.e., the reference frame of radial zone $i+1$. Note that the energy argument ${\epsilon }_{0i}$, referenced by the subscript k, is the same in the transformed and untransformed function ${\psi }^{(0)}$, reflecting the velocity independence of the Chimera energy grid.

To obtain an expression for ${\psi }_{0\,i,k,q}^{(0)\,n+1,\uparrow }$, let the quantities with subscripts i and i + 1 denote their values defined with respect to the frame corresponding to radial coordinate i and i + 1, respectively. The transformed neutrino functions will be denoted by ${\psi }_{0\,i,k,q}^{(0)\,n+1,\uparrow }$, as above, to emphasize that they are transformed from frame i to frame i + 1 rather than originally residing in radial zone i + 1. Suppressing the indices k, n, and q, we use the invariance of ${f}_{0}({\epsilon }_{0},{\mu }_{0})$—e.g., ${f}_{0i+1}({\epsilon }_{0i+1},{\mu }_{0i+1})={f}_{0i}({\epsilon }_{0i},{\mu }_{0i})$—and the transformation properties of the independent variables ${\epsilon }_{0}$ and ${\mu }_{0}$,

Equation (297)

to derive Equation (298) below, where we have kept terms to order ${ \mathcal O }(v/c)$. The use of ${\psi }_{0\,i,k,q}^{(0)\,n+1,\uparrow }$ in place of ${\psi }_{0\,i,k,q}^{(0)\,n+1}$ in Equation (274) enables the correct flux to be computed through zone interface $i+{}^{1}{/}_{2}$, while the difference between ${\psi }_{0\,i,k,q}^{(0)\,n+1}$ and ${\psi }_{0\,i,k,q}^{(0)\,n+1,\uparrow }$ constitutes the discontinuity in ${\psi }_{0}^{(0)}$ at the shock. It must be noted that the replacement of ${\psi }_{0\,i,k,q}^{(0)\,n+1}$ in Equation (274) by ${\psi }_{0\,i,k,q}^{(0)\,n+1,\uparrow }$ is only performed when the flux is being computed at interface $i+{}^{1}{/}_{2}$. When computed at interface $i-{}^{1}{/}_{2}$, ${\psi }_{0\,i,k,q}^{(0)\,n+1}$ is not modified.

To implement this part of the shock transport algorithm, Chimera computes $d{\psi }_{0\,i,k,q}^{(0)\,\mathrm{shock}}$, the difference between ${\psi }_{0\,i,k,q}^{(0)\,n+1,\uparrow }$ and ${\psi }_{0\,i,k,q}^{(0)\,n+1}$,

Equation (298)

in the presence of a shock, and $d{\psi }_{0\,i,k,q}^{(0)\,\mathrm{shock}}=0$ otherwise, where

Equation (299)

and the zone-centered value ${\psi }_{0\,i,k,q}^{(1)}={\eta }_{i,k,q}^{(1)}{\psi }_{0\,i,k,q}^{(0)}$, where η is defined in Section 6.16. With $d{\psi }_{0\,i,k,q}^{(0)\,\mathrm{shock}}$ computed, Equation (274) is modified to read

Equation (300)

A modification is needed in Equation (273) to complete this scheme. In Equation (273), the change in ${\psi }_{0\,i,k,q}^{(0)\,n+1}$ due to transport is determined by the fluxes through the outer zone edge, $i+{}^{1}{/}_{2}$, proportional to ${\psi }_{0\,i+\displaystyle \frac{1}{2},k,q}^{(1)\,n+1}$, and by the fluxes through the inner zone edge, $i-{}^{1}{/}_{2}$, proportional to ${\psi }_{0\,i-\displaystyle \frac{1}{2},k,q}^{(1)\,n+1}$. However, ${\psi }_{0\,i+\displaystyle \frac{1}{2},k,q}^{(1)\,n+1}$ at the outer zone edge is computed in the $i+1$ frame by Equation (300) with ${\psi }_{0\,i,k,q}^{(0)\,n+1\,\uparrow }={\psi }_{0\,i,k,q}^{(0)\,n+1}+d{\psi }_{0\,i,k,q}^{(0)\,n+1,\mathrm{shock}}$ replacing ${\psi }_{0\,i,k,q}^{(0)\,n+1}$, while ${\psi }_{0\,i-\displaystyle \frac{1}{2},k,q}^{(1)\,n+1}$ at the inner zone edge is computed in the ith frame. Consistency is restored by replacing ${\psi }_{0\,i+\displaystyle \frac{1}{2},k,q}^{(1)\,n+1}$ in Equation (273) by ${\psi }_{0\,i+\displaystyle \frac{1}{2},k,q}^{(1)\,n+1,\downarrow }$, where the latter is the former transformed from the frame at $i+1$ to the frame at i. This will ensure that the fluxes computed through both the inner and outer edges of zone i are with respect to the same frame.

To derive an expression for ${\psi }_{0\,i+\displaystyle \frac{1}{2},k,q}^{(1)\,n+1,\downarrow }$, we will be transforming quantities from the frame at radial zone $i+1$ to the frame at radial zone i. To keep the notation manageable, we will again suppress the indices k, n, and q, and denote quantities defined at frame i and $i+1$, as before, by subscript i and $i+1$, respectively, and neutrino functions transformed from frame $i+1$ to frame i by a down-arrow to distinguish them from neutrino functions actually residing at radial zone i. In particular, the quantity ${\psi }_{0\,i+\displaystyle \frac{1}{2},k,q}^{(1)\,n+1}$ is computed in the $i+1$ frame by Equation (300) above, and we will continue to denote its value when transformed to the i frame by ${\psi }_{0\,i+\displaystyle \frac{1}{2},k,q}^{(1)\,n+1,\downarrow }$. The transformation is obtained by equating the number of neutrinos passing through an area dS perpendicular to the radial direction in a time dt as seen from the two frames (Mihalas & Mihalas 1984, p.413); namely,

Equation (301)

Using ${\epsilon }_{0i+1}\,{\mu }_{0i+1}\,d{\mu }_{0i+1}={\epsilon }_{0i}\,{\mu }_{0i}\,d{\mu }_{0i}$ in the expression for ${{dN}}_{0i+1}$, and transforming the remaining variables to the i-frame, we get

Equation (302)

Dividing by $c\,{\epsilon }_{0\,i}^{2}\,d{\epsilon }_{0i}\,{dS}\,{{dt}}_{0i}$ and considering terms to order ${ \mathcal O }(v/c)$,

Equation (303)

To implement this part of the shock transport algorithm, Chimera computes $d{\psi }_{0\,i+\tfrac{1}{2},q}^{(1)\,\mathrm{shock}}({\epsilon }_{0i})$, the difference between ${\psi }_{0\,i+1,q}^{(1)\,\downarrow }\,({\epsilon }_{0i})$ and ${\psi }_{0\,i+1,q}^{(1)}({\epsilon }_{0i})$, both defined at zone edge $i+{}^{1}{/}_{2}$, by

Equation (304)

in the presence of a shock, and by $d{\psi }_{0\,i+\tfrac{1}{2},k,q}^{(1)\,\mathrm{shock}}=0$, otherwise, where

Equation (305)

and where ${\psi }_{0\,i,k,q}^{(2)}={\eta }_{i,k,q}^{(2)}{\psi }_{0\,i,k,q}^{(0)}$. With $d{\psi }_{0\,i+\tfrac{1}{2},k,q}^{(1)\,\mathrm{shock}}$ computed, Equation (273) is modified to read

Equation (306)

Numerical experiments have shown that rolling in this algorithm from $| {v}_{i+1}-{v}_{i}| $ = 0.01–0.02 c with a corresponding rollout of the energy advection algorithm—i.e., at the shock, this algorithm replaces the energy-advection algorithm—gave excellent results. In the above comparisons of the neutrino rms energies and shock trajectories computed with the Chimera and Agile-Boltztrancodes, the Chimera results with the use of the above shock transport algorithm, plotted in Figure 23 by the lines in blue, show much better agreement with the results from Agile-Boltztran. The discontinuities in ${\psi }_{0}^{(0)}$ are clearly evident, as is the absence of the unphysical rise in ${\psi }_{0}^{(0)}$ as the shock is approached from below. Further examples of Chimera test results without the use of this scheme (Series B) and with the use of this scheme (Series C) are given in Sections 7 and 9 and demonstrate the scheme's ability to give accurate transport solutions across a shock.

As a test of the ability of Chimera to conserve total energy, the simulations described above were redone with the Newtonian gravitational potential substituted for the general relativistic monopole potential. Substitution of the general relativistic monopole in place of the Newtonian monopole in the multipole expansion of the self-gravitational potential is not conservative and, therefore, not suitable for a test of total energy conservation. We also set the lapse function to unity. The results are shown in Figure 23, Panels (c) and (d). In both simulations, energy is conserved with high accuracy up to bounce, where a discontinuous change in energy of about −0.35 B occurs in both simulations. Thereafter, the total energy of the simulation with the shock transport algorithm turned off is essentially flat for the next ∼140 ms, after which it decreases by about −0.1 B during the subsequent ∼100 ms. The total energy of the simulation with the shock transport algorithm turned on increases by about 0.1 B for the ∼140 ms following bounce and then remains essentially flat for the next ∼100 ms. The components of the total energy of the two Chimera simulations are plotted in Panel (d). There are small differences discernible in the components of the total energy in the two simulations, particularly between 300 and 400 ms.

6.13. Neutrino Energy Advection Step

Following the Lagrangian transport step, Chimera executes the Lagrangian energy advection step, which consists of solving Equation (241), which has been operator split along with the "transport" Equation (240) from Equation (222). Omitting the superscript E from ${\left(\partial {\psi }_{0}^{(0)}/\partial t\right)}_{m,{E}_{0}}$ with the understanding that this time derivative will hereafter refer to the change in ${\psi }_{0}^{(0)}$ due to the energy advection step, and considering for the moment only the terms in Equation (241) involving ${\psi }_{0}^{(0)}$ (terms involving ${\psi }_{0}^{(2)}$ will be considered shortly), we rearrange these terms as follows:

Equation (307)

or, multiplying by ${E}_{0}^{2}({{bR}}^{2}/{a}^{3})$,

Equation (308)

Note that we could have multiplied Equation (307) by any power of ${E}_{0}$ in deriving the final form of Equation (308), but the second term would always have been left with one more power of ${E}_{0}$ than the first term, and this feature will be important below. The motivation for our choice of the quantity ${\epsilon }_{0}^{2}{\psi }_{0}^{(0)}$ is two fold: (1) this choice enables us to absorb the inhomogeneous terms involving ${\psi }_{0}^{(0)}$ into the energy derivative operator, and (2) it facilitates a physical interpretation, described below, which provides a useful guide for differencing, and is based on the fact that ${E}_{0}^{2}{\psi }_{0}^{(0)}$ is proportional to the neutrino number density per unit energy.

Consider now the terms in Equation (241) involving ${\psi }_{0}^{(2)}$. To derive an expression which, after multiplying by ${E}_{0}^{2}$, is similar to the second term in Equation (308), we rearrange the terms involving ${\psi }_{0}^{(2)}$ as

Equation (309)

Substituting Equation (308) and (309) (the latter multiplied by ${E}_{0}^{2}$) in Equation (241), the energy advection equation becomes

Equation (310)

where we have again introduced the scalar Eddington factor, ${\eta }^{(2)}={\eta }^{(2)}(R,t,{\epsilon }_{0})$, defined by

Equation (311)

as originally introduced at the beginning of Section 6.3 and defined mathematically in Section 6.16. Finally, defining the quantities ${{\rm{\Psi }}}_{0}^{(0)}$ and ${ \mathcal V }$ as

Equation (312)

and

Equation (313)

Equation (310) can be written in the compact form:

Equation (314)

If we take into account the factor of ${E}_{0}$ in ${ \mathcal V }$ and neglect the possible dependence of ${\eta }^{(2)}$ on ${E}_{0}$ over a time step, then writing ${ \mathcal V }={E}_{0}{{ \mathcal V }}_{0}$, we use the fact that ${{ \mathcal V }}_{0}$ is now independent of ${E}_{0}$ to write Equation (314) as

Equation (315)

which has the solution

Equation (316)

where the general solution has been particularized by requiring that ${{\rm{\Psi }}}_{0}^{(0)}(m,{t}^{n+1},{E}_{0}^{n+1})\to {{\rm{\Psi }}}_{0}^{(0)}(m,{t}_{i}^{n},{E}_{0}^{n})$ as ${t}^{n+1}\to {t}^{n}$, where ${E}_{0}^{n+1}$ and ${E}_{0}^{n}$ are related by

Equation (317)

and where Equation (315) has been integrated over the time step ${t}^{n}\to {t}^{n+1}$.

A convenient expression from which the finite differencing of Equation (316) can be developed, which additionally lends itself to a physical interpretation, is obtained by (i) using Equation (312) to transform ${{\rm{\Psi }}}_{0}^{(0)}$ back to ${\psi }_{0}^{(0)}$ in Equation (316), (ii) multiplying the result by ${\left(4\pi \right)}^{2}{\left({hc}\right)}^{-3}m{\rm{\Delta }}{E}_{0}^{n+1}$, where m is a given comoving mass, and (iii) recalling that $b=1/(4\pi {R}^{2}\rho )$. The result is

Equation (318)

where $V=m/\rho $, and Equation (317) (first factor on the right-hand side of Equation (316)) has been used to replace ${\rm{\Delta }}{E}_{0}^{n+1}$ by ${\rm{\Delta }}{E}_{0}^{n}$ on the right-hand side of Equation (318). The left-hand side of Equation (318) is ${\rm{\Delta }}{N}_{\nu }(m,{t}^{n+1},{E}_{0}^{n+1})$, the number of neutrinos of energy ${\epsilon }_{0}^{n+1}={E}_{0}^{n+1}/{a}^{n+1}$ within the energy width ${\rm{\Delta }}{\epsilon }_{0}^{n+1}={\rm{\Delta }}{E}_{0}^{n+1}/{a}^{n+1}$ in a comoving spatial volume V of mass m after the energy advection step, and the right-hand side is ${\rm{\Delta }}{N}_{\nu }(m,{t}^{n},{E}_{0}^{n})$, the number of the same set of neutrinos before the energy advection step. Equation (318) therefore states that the neutrinos, numbering ${\rm{\Delta }}{N}_{\nu }(m,{t}^{n},{E}_{0}^{n})$, in a comoving volume Vn with energies between ${\epsilon }_{0}^{n}={E}_{0}^{n}/{a}^{n}$ and ${\epsilon }_{0}^{n}+{\rm{\Delta }}{\epsilon }_{0}^{n}=({E}_{0}^{n}+{\rm{\Delta }}{E}_{0}^{n})/{a}^{n}$ are shifted in energy, while conserving number, to energies between ${\epsilon }_{0}^{n+1}={E}_{0}^{n+1}/{a}^{n+1}$ and ${\epsilon }_{0}^{n+1}+{\rm{\Delta }}{\epsilon }_{0}^{n+1}=({E}_{0}^{n+1}+{\rm{\Delta }}{E}_{0}^{n+1})/{a}^{n+1}$, where ${E}_{0}^{n}$ and ${E}_{0}^{n+1}$ are related by Equation (317).

As a particular example of this scheme, consider the neutrino diffusion limit in which ${\eta }^{(2)}={}^{1}{/}_{3}$. Recalling that $b=1/(4\pi {R}^{2}\rho )$, Equation (313) with ${E}_{0}$ factored out becomes

Equation (319)

It follows that

Equation (320)

and using Equation (320) in Equation (317), we finally get

Equation (321)

which states that in the diffusion limit, the neutrino energy scales as ρ1/3 under compression or expansion, which is the expected property of a relativistic gas.

Given the solutions, Equations (317) and (318), of the neutrino energy advection Equation (241), the Chimera numerical scheme for updating the neutrinos in energy proceeds in three steps and is implemented during the radial sweep (Figure 1) after the Lagrangian hydrodynamics step in which the density ρ and the areal radius R are updated. The steps are as follows:

  • 1.  
    We begin with an energy Lagrangian step arising from the changes in ρ, R, or a. This shifts the neutrino energies to their new values and modifies ${\psi }_{0}^{(0)}$ through Equations (317) and (318) in a way that conserves the total number of neutrinos in a given comoving fluid volume. Labeling the values of quantities after this Lagrangian step with the superscript $n^{\prime} +1$, Equations (317) and (318), in differenced form, give for this step
    Equation (322)
    where the energy grid is displaced by
    Equation (323)
    and
    Equation (324)
    where da, dR, and db are the changes in a, R, and b resulting from the Lagrangian hydrodynamics step or the update of the gravitational potential.
  • 2.  
    Next, we perform a remap of the energy grid, which has been displaced during step 1 from ${E}_{0\,k}^{n}$ to ${E}_{0\,k}^{n^{\prime} +1}$, back to the initial set of values ${E}_{0\,k}^{n}$. That is, letting quantities with the superscript $n+1$ denote their final values at the completion of the remap, this step entails
    Equation (325)
    This step uses the PPM advection technology and is described in more detail in Section 6.14. Briefly, for a given radial zone i, the quantities ${\epsilon }_{0\,i,k}^{2}{\psi }_{0\,i,k,q}^{(0)}$, which are proportional to the neutrino number density per unit energy, are given piecewise parabolic profiles as a function of energy, which are then averaged over the displacement of the energy grid occasioned by the remap. These averages provide a high-order representation of the neutrino fluxes at the energy grid edges, which are used to remap the energy grid to its final values. While the neutrino number densities are automatically conserved in this step, the fluxes are scaled by a constant overall factor of the order unity to ensure that the total neutrino energy in a given spatial volume before and after this remap is conserved as well. (This remap step just shifts the energy grid and should therefore conserve both neutrino number and energy.) For more details of this latter procedure, refer to Equations (335)–(337) and the accompanying discussion.
  • 3.  
    In the third step, the neutrino occupation probabilities ${\psi }_{0}^{(0)}$ are checked to ensure that none exceed unity. If ${\psi }_{0\,i,k,q}^{(0)}\gt 1$ for one or more values of k for a given i, each such ${\psi }_{0\,i,k,q}^{(0)}$ is then subject to the following variant of the algorithm described in Bruenn (1985, Equation (B8)), which conserves number and energy but limits ${\psi }_{0}^{(0)}\leqslant 1$. To describe the algorithm used here, let ${{\rm{\Delta }}}^{(-)}{\psi }_{0\,i,k,q}^{(0)}$ and ${{\rm{\Delta }}}^{(-)}{\psi }_{0\,i,k+2,q}^{(0)}$ correspond to the number of neutrinos removed from zone (i, k) and zone $(i,k+2)$, respectively, and let ${{\rm{\Delta }}}^{(+)}{\psi }_{0\,i,k+1,q}^{(0)}$ correspond to the number of neutrinos added to zone $(i,k+1)$. Then, specifying ${{\rm{\Delta }}}^{(-)}{\psi }_{0\,i,k,q}^{(0)}$ and using the conservation of number and energy are sufficient to fix the values of ${{\rm{\Delta }}}^{(+)}{\psi }_{0\,i,k+1,q}^{(0)}$ and ${{\rm{\Delta }}}^{(-)}{\psi }_{0\,i,k+2,q}^{(0)}$ so that number and energy is conserved. In particular, neutrino number conservation requires that
    Equation (326)
    while energy conservation requires that
    Equation (327)
    For a given ${{\rm{\Delta }}}^{(-)}{\psi }_{0\,i,k,q}^{(0)}$, Equations (326) and (327) give
    Equation (328)
    and
    Equation (329)
    Proceeding from k = 0 to ${N}_{k}-2$, ${\psi }_{0\,i,k,q}^{(0)}$ is checked, and if found to exceed unity, then ${{\rm{\Delta }}}^{(-)}{\psi }_{0\,i,k,q}^{(0)}={\psi }_{0\,i,k,q}^{(0)}-1$ is applied to Equations (328) and (329) and the computed values of ${{\rm{\Delta }}}^{(+)}{\psi }_{0\,i,k+1,q}^{(0)}$ and ${{\rm{\Delta }}}^{(-)}{\psi }_{0\,i,k+2,q}^{(0)}$ are added and removed from zones $(i,k+1)$ and $(i,k+2)$, respectively. The neutrino energies ${\epsilon }_{0i,{N}_{k}-1}$ and ${\epsilon }_{0i,{N}_{k}}$ are always chosen to be large enough so that they are never overfilled. Finally, if ${{\rm{\Delta }}}^{(-)}{\psi }_{0\,i,k+2,q}^{(0)}$ is found to exceed ${\psi }_{0\,i,k+2,q}^{(0)}$, ${{\rm{\Delta }}}^{(-)}{\psi }_{0\,i,k,q}^{(0)}$ is adjusted to prevent ${\psi }_{0\,i,k+2,q}^{(0)}$ from becoming negative. In this case, the algorithm only works approximately. In practice, conservation of leptons is always found to be conserved to machine accuracy, so if the above algorithm is ever employed, it never encounters a value of ${{\rm{\Delta }}}^{(-)}{\psi }_{0\,i,k,q}^{(0)}$ that would cause ${\psi }_{0\,i,k+2,q}^{(0)}$ to become negative.

6.14. Neutrino Energy Advection Due to Changes in the Lapse

The lapse function, a, is given by Equations (262) and (263) and in differenced form by Equations (294) and (295). During the evolution of the core, the lapse function ai of mass shell i may change as a result of a hydrodynamics step due to changes in the configuration of the core. New values of the lapse are computed in Chimera after the gravitational potential update, which follows the Lagrangian hydrodynamics step and remap, and after the gravitational potential update, which immediately precedes the transport solve. A change in lapse also affects neutrinos radially advected between adjacent zones, as described in Section 6.15, as these zones may have different lapses. Because our energy grid is tied to the lapse, these changes in the lapse change the energy grid, viz., ${\epsilon }_{0i,k+\tfrac{1}{2}}={E}_{0k+\tfrac{1}{2}}/{a}_{i}$. As before, the ${E}_{0k+\tfrac{1}{2}}={\epsilon }_{0k+\tfrac{1}{2},\infty }$ are the energy grid edges at radius infinity. Let ${a}_{i}^{n}\to {a}_{i}^{n+1}$ be the change in the lapse as a result of the evolution preceding a gravitational potential update. Let the superscript n and $n+1$ now denote the values of quantities given after the prior energy advection update and after the current energy advection update, respectively. Then the change in the energy grid from n to $n+1$ is given by

Equation (330)

and ${\psi }_{0\,i,k}^{(0)}$ must change in such a way that the total number and energy of neutrinos in each radial zone i remain the same. Ignoring factors common to both sides, we must have

Equation (331)

which is just Equation (318) with ${V}^{n+1}={V}^{n}$ and ${E}_{0}^{n+1}={E}_{0}^{n}$. As in the energy remap step following the energy advection step described above, we also must have

Equation (332)

where

Equation (333)

The quantity ${\rm{\Delta }}{\left({\epsilon }_{0i,k}^{n+1}\right)}^{3}$ is defined similarly, and the equivalence sign is a consequence of our definition of ${\epsilon }_{0i,k}$, Equation (196).

To compute ${\psi }_{0\,i,k,q}^{(0)\,n+1}$ from ${\psi }_{0\,i,k,q}^{(0)\,n}$ due to a change in the lapse, the quantities ${\left({\epsilon }_{0i,k}^{n}\right)}^{2}{\psi }_{0\,i,k,q}^{(0)\,n}$, related to the neutrino number density per unit energy, are PPM interpolated to the grid edges ${\epsilon }_{0\,i,k+\tfrac{1}{2}}^{n}$. From the interpolated ${\left({\epsilon }_{0i,k}^{n}\right)}^{2}{\psi }_{0\,i,k+\tfrac{1}{2},q}^{(0)\,n}$, the flux ${\mathrm{Flx}}_{i+\displaystyle \frac{1}{2},k+\displaystyle \frac{1}{2},q}$ across the grid boundary—i.e., the number of neutrinos in the overlapping shells between the old and new grid—is computed from the grid displacement given by Equation (330). With ${\mathrm{Flx}}_{i,k+\displaystyle \frac{1}{2},q}$ in hand, ${\psi }_{0\,i,k,q}^{(0)\,n+1}$ is then computed from ${\psi }_{0\,i,k,q}^{(0)\,n}$ by performing the advection step

Equation (334)

Equation (334), when summed from 1 to ${N}_{\epsilon }$, automatically conserves neutrino number if the fluxes at the boundaries are zero; i.e., ${\mathrm{Flx}}_{i,1-\tfrac{1}{2},q}={\mathrm{Flx}}_{i,{N}_{\epsilon }+\tfrac{1}{2},q}=0$, as the two terms on the right-hand side of this equation will then sum to zero, leaving us with Equation (331). However, neutrino energy is not necessarily conserved. To ensure energy conservation, we modify Equation (334) by multiplying the two flux terms by a constant factor ${\xi }_{i,q};$ i.e.,

Equation (335)

Equation (335) with any constant value for the parameter ${\xi }_{i,q}$ automatically conserves total neutrino number, as does Equation (334), but the parameter is now adjusted so that Equation (335) conserves energy as well. To determine ${\xi }_{i,q}$, we equate the initial neutrino energy, ${E}_{\nu i,q}$, to the final energy, ${E}_{\nu ,f,q}$, the latter being given by the right-hand side of Equation (335) summed over k. We get, substituting Equation (335) into Equation (332),

Equation (336)

Solving Equation (335) for ${\xi }_{i,q}$ gives

Equation (337)

Equation (335), with ${\xi }_{i,q}$ given by Equation (337), conserves total neutrino number and energy and is the method by which Chimera advances ${\psi }_{0\,i,k,q}^{(0)\,n}$ due to a change in the lapse.

6.15. Neutrino Spatial Advection Step

As described in Section 4.5, following the Lagrangian hydrodynamics step in the case of the θ- and ϕ-sweeps, the grid is remapped back to the configuration that prevailed before the Lagrangian step. In the case of the radial sweep, the grid is remapped back to a configuration specified by the regridder. In either case, the zero moments of the neutrino distribution functions must be advected through the grid in the remapping step. The procedure is similar to that described in Section 4.4.1 for the mass-specific quantities like the mass-specific momenta, and mass-specific angular momenta but simpler, since the quantities to be remapped are volume specific rather than mass specific, thus omitting the necessity of first computing the advected mass.

If $\delta {\xi }_{i+\tfrac{1}{2}}\lt 0$, where $\delta {\xi }_{i+\tfrac{1}{2}}$ is the displacement of the grid element ${\xi }_{i+\tfrac{1}{2}}$ during the remap step, the advection proceeds from the left of ${\xi }_{i+\tfrac{1}{2}}$ to the right; i.e., from zone i to $i+1$. If ${a}_{i,k,q}^{n+1^{\prime} }$ represents ${\psi }_{0\,i,k,q}^{(0)\,n+1^{\prime} }$, the zeroth angular moment of neutrinos of species q, radial/angular/azimuthal index i, and energy zone k, a PPM interpolation profile, ${a}_{i,k,q}(\xi )$, of ${a}_{i,k,q}^{n+1^{\prime} }$ is constructed, as described in Section 4.2. The average, $\langle a{\rangle }_{L,i+\tfrac{1}{2},k,q}^{n+1^{\prime} }$, of ${a}_{i,k,q}^{n+1^{\prime} }$ ("L" denoting "left") over the zone interface displacement $\delta {\xi }_{i+\tfrac{1}{2}}$ is then computed, and the volume displacement $\delta {V}_{i+\tfrac{1}{2}}^{n+1^{\prime} }$ is computed from Equation (126). The quantity of ${a}_{i,k,q}^{n+1^{\prime} }$ advected is then finally computed by

Equation (338)

If $\delta {\xi }_{i+\tfrac{1}{2}}\gt 0$, during the remap, the advection proceeds from the right to the left of interface ${\xi }_{i+\tfrac{1}{2}}^{n+1^{\prime} };$ i.e., from $i+1$ to i. In this case, the quantity advected, $\delta {{ \mathcal A }}_{i+\tfrac{1}{2},k,q}^{n+1^{\prime} }$, is computed from an equation similar to Equation (338) but using $\langle a{\rangle }_{R,i+\tfrac{1}{2},k,q}^{n+1^{\prime} }$ ("R" denoting "right") computed from ${a}_{i+1,k,q}^{n+1^{\prime} }$.

In either case, the advection then proceeds as

Equation (339)

where again we denote the variables after the Lagrangian step but before the remap step by the superscript $n+1^{\prime} $, and after the remap step by the superscript $n+1$.

When Chimera is run in relativistic mode, the lapse function becomes a function of the location in the core, and the neutrino energy zones also become a function of location through the lapse, as given by Equation (216). In this case, Chimera adds a couple of steps to the procedure described above for remapping ${\psi }_{0}^{(0)}$. If $\delta {\xi }_{i+\tfrac{1}{2}}\lt 0$, neutrinos are advected from zone i to zone $i+1$, and the different lapse functions in these two zones must be accounted for. This is accomplished by appropriately advecting in energy $\langle a{\rangle }_{L,i+\tfrac{1}{2},k,q}^{n+1^{\prime} }$ before computing $\delta {{ \mathcal A }}_{i+\tfrac{1}{2},k,q}^{n+1^{\prime} }$ by Equation (338). This energy advection is performed by the algorithm described by Equations (322)–(325), with ai replacing an, ${a}_{i+1}$ replacing ${a}^{n+1^{\prime} }$, and with ${\rho }^{n}={\rho }^{n+1^{\prime} }$, dR = 0, and db = 0. If $\delta {\xi }_{i+\tfrac{1}{2}}\gt 0$, the neutrinos are advected from zone $i+1$ to i, and a similar modification of $\langle a{\rangle }_{R,i+\tfrac{1}{2},k,q}^{n+1^{\prime} }$ is performed by the neutrino energy advection algorithms just described but with ${a}_{i+1}$ replacing an, and ai replacing ${a}^{n+1^{\prime} }$.

6.16. Scalar Eddington Factors

The scalar Eddington factor, ${\eta }^{(2)}={\psi }_{0}^{(2)}/{\psi }_{0}^{(0)}$, appears in the neutrino energy advection Equations (310), (313), and (323). Consider first the case for $R\gt {R}_{\nu }$, where ${R}_{\nu }$ is the radius of the neutrinosphere and is a function of both neutrino energy and flavor. Analogous to the derivation of the geometric piece of the flux limiter, Equation (233), we use Equation (220) for ${\psi }_{0}^{(2)}$ to write

Equation (340)

where this assumes that the neutrino distribution function ${f}_{0}$ is constant for rays satisfying ${\mu }_{0}\lt {\mu }_{0\nu }$, and zero otherwise, where ${\mu }_{0\nu }$ is defined in the text above Equation (233). The angle ${\mu }_{0}$ is given by Equations (234)–(236).

Where $R\lt {R}_{\nu }$, the diffusion limit applies, ${\psi }_{0}^{(1)}$ is well defined, and we solve the equation ${\psi }_{0}^{(1)}=\tfrac{1}{2}(1+{\mu }_{0\nu }){\psi }_{0}^{(0)}$ for ${\mu }_{0}$ to get

Equation (341)

where the flux factor ${\eta }^{(1)}$ is the ratio of ${\psi }_{0}^{(1)}$ to ${\psi }_{0}^{(0)}$. Using Equation (341) in Equation (340), we get the correct behavior in the strong diffusion limit, where ${\psi }_{0}^{(1)}\to 0$, as this implies that ${\eta }^{(1)}\to 0$ and ${\eta }^{(2)}\to {}^{1}{/}_{3}$.

Versions of Chimera subsequent to version C) have used the Minerbo closure (Minerbo 1978), constructed on the maximum entropy principle. This closure is the classical limit of the of closure constructed by Cernohorsky & Bludman (1994) on the same principle but using the Fermi–Dirac distribution, and it is much easier to implement and gives very similar results (Murchikova et al. 2017). For this closure, the quantity p is defined by

Equation (342)

and ${\eta }^{(2)}$ is given by

Equation (343)

where ${\eta }_{\mathrm{thick}}^{(2)}=1/3$, ${\eta }_{\mathrm{thin}}^{(2)}=1$, and as before, ${\eta }^{(1)}={\psi }^{(1)}/{\psi }^{(0)}$. The closure given by Equations (342) and (343) gives results slightly closer to those given by a Boltzmann solver than those given by Equations (340) and (341).

7. Stationary State Transport Tests

In the spirit of the tests suggested by Müller et al. (2010), involving the stationary solution of the transport equations, we set up several analogous test problems. As in Müller et al. (2010), we consider a central source consisting of a homogeneous isothermal sphere of radius R of energy-independent absorption opacity and zero scattering opacity radiating into a medium of negligible absorption and scattering opacities.

7.1. Gravitational Redshift

In this test, we assume that the medium external to the central source is static and impose a lapse profile as a function of R given by

Equation (344)

where R0 is chosen so that $a(4{\rm{km}})=0.70$. Equation (344) for a(R) is suggested by the behavior of a(R) outside a spherically symmetric source of gravity in the post-Newtonian approximation and numerically approximates the behavior of a(R) in realistic core collapse models ∼50–100 ms after bounce. An analytic expression for ${\psi }_{0}^{(1)}(R,{E}_{0})$ can be derived for this problem, from which follows an analytic expression for the neutrino luminosity L(R). With the assumption of a static medium, the derivatives with respect to t all vanish, and Equation (222) becomes

Equation (345)

which gives

Equation (346)

where Rp is some fiducial radius and ${a}_{p}=a({R}_{p})$. Using Equation (346), the neutrino luminosity, L(R), is then given by

Equation (347)

The luminosity from the central source thus goes as ${a}^{-2}(R)$ in the surrounding static medium for this problem. Figure 24(a) compares the analytic expression (Equation (347)) with the solutions for the ${\nu }_{e}$ and ${\bar{\nu }}_{e}$ luminosities given by Chimera. The numerical solutions agree extremely well with the analytic solution, which is not surprising, given that the variable transformation (Equation (216)) and the subsequent differencing of the transport equation essentially guarantee this agreement.

Figure 24.

Figure 24. Panel (a): comparison of ${\nu }_{e}$ (red) and ${\bar{\nu }}_{e}$ (blue) luminosities with the analytic solution of Equation (347). Panel (b): same as (a), but for the mean energies, compared to the analytic solution given by Equation (350). The solutions have been normalized to unity at R = 1000 km.

Standard image High-resolution image

An analytic expression for the mean neutrino energy can also be derived. With the assumption of a static medium, Equation (223) becomes

Equation (348)

where the scalar Eddington factor ${\eta }^{2}$ is defined below Equation (223). In the limit of a small angular diameter source, ${\eta }^{(2)}\to 1$ and Equation (348) reduces to an equation for ${\psi }_{0}^{(0)}$ similar to Equation (345) for ${\psi }_{0}^{(1)}$, with a solution similar to Equation (346). More generally, we note that the condition of the energy independence of the absorption and emission opacities, in addition to the static condition, implies that ${\eta }^{(2)}(R,{E}_{0})={\eta }^{(2)}(R)$. Therefore, Equation (348) for ${\psi }_{0}^{(0)}$ depends only on R, and the solution can be written as

Equation (349)

where again ${\psi }_{0}^{(0)}({R}_{p},{E}_{0})$ is the value of ${\psi }_{0}^{(0)}(R,{E}_{0})$ at some fiducial radius Rp, and $\chi (R)$ is the solution of Equation (348) for $\chi (R)$ with ${\psi }_{0}^{(0)}({R}_{p},{E}_{0})$ substituted for ${\psi }_{0}^{(0)}(R,{E}_{0})$, under the condition that $\chi ({R}_{p})=1$. The mean energy is therefore given by

Equation (350)

For this problem, the mean energy from the central source thus goes as ${a}^{-1}(R)$ in the surrounding static medium. Figure 24(b) shows that the analytic expression (Equation (350)) agrees extremely well with the solutions for the ${\nu }_{e}$ and ${\bar{\nu }}_{e}$ mean energies given by Chimera.

7.2. Imposed Shock Velocity Profile

In this test, a radial velocity field is imposed that mimics the velocity profile encountered during the accretion phase of a CCSN. The velocity profile suggested by Müller et al. (2010) is

Equation (351)

Inside the homogeneous, isotropic central spherical source of 4 km radius, we again turn off all scattering opacities and, rather than employ a frequency-independent absorption opacity as in the preceding test, employ instead the Bruenn (1985) free nucleon absorption opacity corresponding to the state ($\rho ,T,{Y}_{{\rm{e}}}$) = (${10}^{11}\,{{\rm{g\; cm}}}^{-3},4\,{\rm{MeV}},0.5$). Outside the central source, the absorption and scattering opacities vanish, so the luminosity L and the mean neutrino energy $\langle E\rangle $ are constant with R in the lab frame. The nonzero velocity regime is at a large enough radius, compared with the source radius, that the neutrino flow can be well approximated as radially free streaming. In this case, both $f(\epsilon ,\mu )$ and ${f}_{0}({\epsilon }_{0},{\mu }_{0})$, the invariant neutrino occupation probabilities, vanish except at $\mu =\pi $, ${\mu }_{0}=\pi $, respectively. The comoving neutrino energy density ${E}_{0\nu }$ is then related to the lab frame energy density ${E}_{\nu }$ by

Equation (352)

Likewise, the comoving neutrino number density, ${N}_{0\nu }$, is related to the lab frame number density, ${N}_{\nu }$, by

Equation (353)

The comoving neutrino luminosity is related to the (constant) lab frame luminosity by

Equation (354)

and the comoving mean neutrino energy, $\langle {E}_{0\nu }\rangle $, is related to the (constant) lab frame mean neutrino energy, $\langle {E}_{\nu }\rangle $, by

Equation (355)

Figure 25 shows the results of the B-series and the C-series Chimera transport, versus the analytical expressions given by Equations (354) and (355) for the neutrino luminosity and mean energies, respectively. Measured from zero, the luminosity and the mean energy as given by the B-series transport overshoot the correct luminosity and mean energy at the shock by ∼7% and ∼13%, respectively, and these solutions remain above the analytic solution at large R by ∼6%. Moreover, the B-series solution for the comoving mean energy fails to resolve the steep rise of the mean energy at the shock, as given by the analytic solution. The C-series solutions for both the luminosity and the mean energy, however, agree quite well with the analytic solutions, only deviating below the analytic solution at large R by ∼1.5%, demonstrating the impact of the added special treatment of transport across the shock, discussed in Section 6.12.

Figure 25.

Figure 25. Comoving ${\nu }_{e}$ (a) luminosities and (b) mean energies computed with B-series Chimera (green) and C-series Chimera (red; including the shock transport treatment in Section 6.12) compared with the analytic solution for luminosity (Equation (354), Panel (a)) and mean energy (Equation (355), Panel (b)) for the velocity profile specified in Equation (351). Solutions have been normalized to unity at R = 10 km. Normalized results for ${\bar{\nu }}_{e}$ are similar, but the actual ${\bar{\nu }}_{e}$ luminosities are considerably reduced in magnitude, and the actual mean energies are slightly reduced in magnitude.

Standard image High-resolution image

8. Transport Sources

Solving the neutrino radiation hydrodynamics also requires detailed neutrino–matter interactions to couple the radiation field to the fluid, which drives heating, cooling, and changes in ${Y}_{{\rm{e}}}$. Table 2 lists the scattering, absorption-emission, and pair-production opacities currently incorporated in Chimera. For each zone, the logarithms of the opacities are stored at each of the eight corners of a cell in ρT${Y}_{{\rm{e}}}$ lattice space that surrounds the $(\rho ,T,{Y}_{{\rm{e}}})$ value of each zone. The values at $(\rho ,T,{Y}_{{\rm{e}}})$ and the ${Y}_{{\rm{e}}}$-derivatives needed by the Jacobian for the solution of the transport equation are then obtained by a three-dimensional interpolation from the eight corner values in the same manner as for the thermodynamic variables of the EoS. The spacing of the ρT${Y}_{{\rm{e}}}$ grid matches the user-selected spacing of the EoS grid but is typically 20 points per decade in $\mathrm{log}\rho $ and $\mathrm{log}T$, with 0.01 intervals in ${Y}_{{\rm{e}}}$. This resolution ensures that interpolated opacities match those that are computed directly, to within 1%.

Table 2.  Summary of Neutrino Opacities

Process Description References
$\nu +{e}^{\pm }\leftrightharpoons \nu +{e}^{\pm }$ Sections 8.2, 8.3 Bruenn (1985), Mezzacappa & Bruenn (1993b)
$\nu +A\leftrightharpoons \nu +A$ Section 8.4 Bruenn & Mezzacappa (1997), Horowitz (1997)
$\nu +n,p\leftrightharpoons \nu +n,p$ Section 8.5 Reddy et al. (1998), Buras et al. (2006b)
${\nu }_{e}+n\leftrightharpoons {e}^{-}+p$ Section 8.6 Bruenn (1985), Mezzacappa & Bruenn (1993b)
${\bar{\nu }}_{e}+p\leftrightharpoons {e}^{+}+n$ Section 8.6 Bruenn (1985), Mezzacappa & Bruenn (1993b)
${\nu }_{e}+A^{\prime} \leftrightharpoons {e}^{-}+A$ Section 8.7 Langanke et al. (2003), Hix et al. (2003)
${\nu }_{e}+{\bar{\nu }}_{e}\leftrightharpoons {e}^{-}+{e}^{+}$ Sections 8.8, 8.9 Bruenn (1985), Mezzacappa & Bruenn (1993b)
$n,p+n,p\leftrightharpoons n,p+n,p+{\nu }_{e}+{\bar{\nu }}_{e}$ Sections 8.8, 8.10 Hannestad & Raffelt (1998)

Download table as:  ASCIITypeset image

8.1. Scattering: General

In spherical symmetry, the rate of change of the neutrino occupation probability, $f({\mu }_{0},{\epsilon }_{0})$, (we still suppress the dependence of f on t and R), due to scattering process "XX" is given by

Equation (356)

where ${R}_{\mathrm{XX}}^{\mathrm{out}}({\epsilon }_{0},{\epsilon }_{0}^{{\prime} },\cos \theta )$ is the "out-scattering kernel," (i.e., the neutrino-unblocked rate per final and initial neutrino state for scattering from energy ${\epsilon }_{0}$ to energy ${\epsilon }_{0}^{{\prime} }$ through angle θ), ${R}_{\mathrm{XX}}^{\mathrm{in}}({\epsilon }_{0},{\epsilon }_{0}^{{\prime} },\cos \theta )$ is the "in-scattering kernel," (i.e., the neutrino-unblocked rate per final and initial neutrino state for scattering from energy ${\epsilon }_{0}^{{\prime} }$ to energy ${\epsilon }_{0}$ through angle θ), and θ is given in terms of the individual neutrino propagation directions, $({\mu }_{0},\varphi )$ and $({\mu }_{0}^{{\prime} },\varphi ^{\prime} )$, by

Equation (357)

where $\delta \varphi =\varphi ^{\prime} -\varphi $ is the azimuthal direction between the incident and scattered neutrino. The Fermi blocking of neutrino states is incorporated explicitly in Equation (356) by the $(1-f)$ factors rather than in the scattering kernels. The scattering kernels ${R}_{\mathrm{XX}}^{\mathrm{out}/\mathrm{in}}$ have the symmetry

Equation (358)

which follows simply from the fact that an in-scattering from ${\epsilon }_{0}^{{\prime} }$ to ${\epsilon }_{0}$ is the same as an out-scattering from ${\epsilon }_{0}^{{\prime} }$ to ${\epsilon }_{0}$. Additionally, the kernels are related to each other by detailed balance at $\beta =1/({kT})$:

Equation (359)

which follows by substituting equilibrium Fermi–Dirac distributions for f in Equation (356) and setting the left-hand side of that equation to zero. Both Equations (358) and (359) should be respected in any approximation scheme.

For neutrino scattering, the complication that the occupation functions f and $f^{\prime} $ are expressed in terms of ${\mu }_{0}$ and ${\mu }_{0}^{{\prime} }$, respectively, while the scattering kernels are expressed in terms of θ, where ${\mu }_{0}$, ${\mu }_{0}^{{\prime} }$, and θ are related by Equation (357), is overcome by Legendre expanding the scattering kernels and keeping only the terms to the first order; i.e.,

Equation (360)

where the Legendre coefficients are given by

Equation (361)

Note that Equations (358) and (359) imply

Equation (362)

Equation (363)

Applying the moment operators ${\left(4\pi \right)}^{-1}\int d{\rm{\Omega }}$ and ${\left(4\pi \right)}^{-1}\int {\mu }_{0}d{\rm{\Omega }}$ to Equation (356) and using the definitions in Equations (357), (360), and (220), we get the moments of the scattering terms of the collision integral:

Equation (364)

8.2. Neutrino–Electron Scattering

For neutrino–electron scattering (NES), the scattering functions ${{\rm{\Phi }}}_{0,\mathrm{NES}}^{\mathrm{in}}({\epsilon }_{0},{\epsilon }_{0}^{{\prime} })$ and ${{\rm{\Phi }}}_{1,\mathrm{XX}}^{\mathrm{out}}({\epsilon }_{0},{\epsilon }_{0}^{{\prime} })$ are taken from Equation (C50) of Bruenn (1985). The scattering functions for down-scattering and isoenergetic scattering are computed directly, using Equation (362) where possible, and up-scattering is computed by use of the Equation (363). Figure 26 shows some representative inverse mean free paths and fractional energy transfers computed with thermodynamic states chosen to be the same as those of Buras et al. (2006b), for comparison, assuming empty neutrino final states for both a typical Chimera geometrically spaced energy grid of 20 zone-centered energies from 4–250 MeV and a much finer evenly spaced grid from 1–300 MeV.

Figure 26.

Figure 26. Inverse mean free path (left scale, black) and fractional energy transfer, $({\epsilon }_{0\mathrm{final}}-{\epsilon }_{0\mathrm{initial}})/{\epsilon }_{0\mathrm{initial}}$ (right scale, red), as a function of the incident neutrino energy for ${\nu }_{e}$–electron scattering (solid) and ${\bar{\nu }}_{e}$–electron scattering (dashed). The symbols (${\nu }_{e}$, circles; ${\bar{\nu }}_{e}$, triangles) show the results for the typical Chimera grid of 20 energy zones geometrically spaced from 4–250 MeV. The solid lines are the results of a reference calculation using 600 zones evenly spaced from 1–300 MeV. The thermodynamic conditions are the same as those in Buras et al. (2006b), for comparison.

Standard image High-resolution image

8.3. Neutrino–Positron Scattering

Neutrino–positron scattering is computed using the same scattering functions as NES but with the coefficient function of CV and CA interchanged, ${({C}_{V}+{C}_{A})}^{2}\leftrightharpoons {({C}_{V}-{C}_{A})}^{2}$, and the chemical potentials in the Fermi functions replaced by their negatives, ${\mu }_{{e}^{+}}=-{\mu }_{{e}^{-}}$. Figure 27 shows the inverse mean free paths and fractional neutrino energy transfers due to neutrino–positron scattering using the same reference states as for NES (Figure 26).

Figure 27.

Figure 27. Same as Figure 26 but for ${\nu }_{e}$–positron scattering and ${\bar{\nu }}_{e}$–positron scattering.

Standard image High-resolution image

8.4. Neutrino–Nucleus Scattering

Neutrino–nucleus ($\nu {\rm{A}}$) scattering is treated as isoenergetic such that

Equation (365)

with the moments

Equation (366)

and therefore,

Equation (367)

The calculations of ${{\rm{\Phi }}}_{0,\nu {\rm{A}}}^{\mathrm{out}}({\epsilon }_{0})$ and ${{\rm{\Phi }}}_{1,\nu {\rm{A}}}^{\mathrm{out}}({\epsilon }_{0})$ are performed using the formalism of Bruenn & Mezzacappa (1997), and results for the inverse mean free path are plotted in Figure 28 for the indicated thermodynamic conditions. The symbols show the results using the typical Chimera neutrino energy grid of 20 energy zones geometrically spaced from 4–250 MeV, the solid lines show the results obtained using a 600-zone energy grid evenly spaced from 1 to 300 MeV. These thermodynamic conditions were chosen as representative of the infall epoch when material entropy is low and nuclei dominate the composition. As expected for an isoenergetic scattering process, the fineness of the neutrino energy grid does not affect the inverse mean free paths, as they depend only on the incident neutrino energy. The nuclear form factor correction of the scattering rates that reduces the inverse mean free paths at high energies, where the incident neutrino wavelengths become comparable to or smaller than the inter-nucleon distances in the nuclei, is included. The black lines and symbols show the inverse mean free paths uncorrected for the liquid structure function (ion–ion correlations). The red lines and symbols show inverse mean free paths with the liquid structure function included, computed as described by Horowitz (1997). The liquid structure function has the effect of substantially reducing the inverse mean free paths at low energies, as can be seen by comparing the black and red lines for ${\epsilon }_{0\nu }\leqslant 10\,\mathrm{MeV}$.

Figure 28.

Figure 28. Neutrino–nucleus neutral-current scattering inverse mean free paths for neutrinos (solid lines) and anti-neutrinos (dashed lines), with the nuclear form factors included for the indicated thermodynamic conditions. Black lines and symbols show neutrino–nucleus scattering without ion–ion correlation corrections. Red lines and symbols show the scattering with the ion–ion correlations applied. Filled circles (neutrinos) and triangles (anti-neutrinos) show the results for the typical Chimera neutrino energy grid of 20 geometrically spaced energy zones. The solid lines show the results using a 600-zone energy grid evenly spaced from 1–300 MeV. The ordinate is log scaled to show the effect of the ion–ion correlation corrections, which for these conditions, mainly affects low-energy neutrinos, which have small inverse mean free paths.

Standard image High-resolution image

8.5. Neutrino–Nucleon Scattering

Differential rates for neutrino elastic scattering on free nucleons are calculated using the formalism of Reddy et al. (1998). These rates include nucleon recoil and degeneracy effects, as well as special relativity, and the implementation of these rates in Chimera is similar to that of NES but differs in several important respects. First, analytic expressions for the Legendre moments of these rates are not available, necessitating their calculation by numerical integration. Second, energy transfer between a neutrino and a nucleon is small but is an important component of the energy transfer between neutrinos and matter, as the scattering rates are relatively large. The smallness of the energy transfer is illustrated in Figure 29, which shows the maximum relative energy transfer of a neutrino scattering from a stationary nucleon. In the important energy range 5–30 MeV, typical of the energies of neutrinos emerging from the neutrinosphere, the maximum relative change in neutrino energy is only 0.01 to 0.05. The actual change in energy depends on the angle between the scattered and incident neutrino and tends to zero as this angle tends to zero. As the relative width of the neutrino energy zoning typically used in Chimera (viz. 20 energy zones geometrically spaced from 4 to 250 MeV) is ∼0.26, and therefore considerably larger than the neutrino–nucleon energy exchange in a scattering except for neutrino energies larger than ∼200 MeV, a much finer energy grid is needed to adequately resolve the energetics of this process.

Figure 29.

Figure 29. Maximum relative change in the energy of a neutrino, scattering on a stationary nucleon, as a function of the incident neutrino energy. The insert shows details for incident energies typical of neutrinos emerging from the neutrinosphere.

Standard image High-resolution image

In the case of neutrino–electron or neutrino–positron scattering, for example, where the energy transfer is not small in comparison with the widths of our energy grid, Equation (356) is differenced using zone-centered values of energy in both the neutrino distribution function and the scattering kernels; i.e.,

Equation (368)

However, for small energy transfers compared with the energy grid width, the scattering kernel ${R}_{\mathrm{XX}}^{{\rm{a}}}({\epsilon }_{0k},{\epsilon }_{0k^{\prime} },\cos \theta )$ will be effectively zero if ${\epsilon }_{0k}\ne {\epsilon }_{0k^{\prime} }$, and the scattering will become essentially isoenergetic, with negligible energy transfer. To develop a better approximation for the energy transfer via smaller energy transfer scatterings, we continue to evaluate the neutrino distribution function at the energy zone centers, but evaluate the scattering kernels on a refined energy grid. To accomplish this, we regard the scattering kernels as functions of ${\epsilon }_{0}$ and ${\epsilon }_{0}^{{\prime\prime} }$ and operate on Equation (368) by the unity operator

Equation (369)

The left-hand side, being a constant operand, is unchanged by this operation. However, by integrating over ${\epsilon }_{0}^{3}$, the right-hand side of Equation (368) becomes proportional to ${dE}({\mu }_{0},{\epsilon }_{0k})/{dt}$, which we desire to compute accurately. Operating on the right-hand side of Equation (368) by Equation (369) gives

Equation (370)

The integrals in Equation (370) over ${\epsilon }_{0}$ and ${\epsilon }_{0}^{{\prime\prime} }$ can now be replaced by summations over a refined energy grid.

The Legendre coefficients ${{\rm{\Phi }}}_{{\ell },\mathrm{XX}}^{\mathrm{in}/\mathrm{out}}({\epsilon }_{0},{\epsilon }_{0}^{{\prime} })$ of the scattering kernels are given by Equation (361) and used in the collision integrals of the transport equations as given by Equation (364). In the computation of the Legendre coefficients for neutrino–nucleon scattering, the summations of the scattering kernels over refined energy grids in Equation (370) are directly incorporated by defining

Equation (371)

where ${{\rm{\Phi }}}_{{\ell },{\rm{n}},{\rm{p}}}^{\mathrm{in}}({\epsilon }_{0k},{\epsilon }_{0k^{\prime} })$ is obtained from ${{\rm{\Phi }}}_{{\ell },{\rm{n}},{\rm{p}}}^{\mathrm{out}}({\epsilon }_{0k},{\epsilon }_{0k^{\prime} })$ by Equation (359), and where, as in Equation (370), the integrals in Equation (371) can be replaced by summations over a refined energy grid.

Technically, the lower and upper limits of the energy integrals in Equation (371) should be the zone-edged energies ${\epsilon }_{0k-\tfrac{1}{2}}$ and ${\epsilon }_{0k+\tfrac{1}{2}}$ in the initial energy integral, and ${\epsilon }_{0k^{\prime} -\tfrac{1}{2}}$ and ${\epsilon }_{0k^{\prime} +\tfrac{1}{2}}$ in the final. However, the scattering kernels can become so narrow in energy that restricting the integration range to a narrower energy width does not change the result and allows a given energy grid to be spread more finely over the smaller energy integration width. A useful integration range for energy as given by the widths of the scattering kernels is angle dependent and within the interval ${\epsilon }_{0{\rm{L}}}$ to ${\epsilon }_{0{\rm{U}}}$ where

Equation (372)

where

Equation (373)

and where α is a free parameter. Experiments found that the choice of ${\alpha }^{2}=5$ gives good results in the sense that the scattering kernel is essentially negligible outside these limits. The derivation of Equations (372) and (373) is given by Equations (378)–(382) below. Whether the zone-edged energies of the Chimera energy grid or the energies given by Equation (372) are used as integration limits is case-dependent and described below.

Three cases are considered depending on the relative values of ${\epsilon }_{0}^{\mathrm{in}}={\epsilon }_{0k}$ and ${\epsilon }_{0}^{\mathrm{out}}={\epsilon }_{0k^{\prime} }$, where ${\epsilon }_{0k}$ and ${\epsilon }_{0k^{\prime} }$ are the zone-centered, incoming and outgoing neutrino energy group.

  • 1.  
    If $k=k^{\prime} $, the initial energy integration over ${\epsilon }_{0}$ as given by the right-hand side of Equation (371) is omitted, the zone centered energy ${\epsilon }_{0k}$ is used as input, and the final energy integration is performed with a 32-point Gauss–Legendre integration from ${\epsilon }_{0{\rm{L}},{\rm{f}}}$ to ${\epsilon }_{0{\rm{U}},{\rm{f}}}$, where
    Equation (374)
    The first two Legendre moments ${{\rm{\Phi }}}_{{\ell },{\rm{n}},{\rm{p}}}^{\mathrm{out}}({\epsilon }_{0k},{\epsilon }_{0k})$ are then obtained by a 32-point Gauss–Legendre angular quadrature appropriately weighted by the Legendre polynomials ${P}_{{\ell }}(\cos \theta )$.
  • 2.  
    If $k^{\prime} =k-1$, the zero and first Legendre moments are then obtained from Equation (371) by first executing eight-point Gauss–Legendre energy quadratures for each of the energy integrations. The integration limits for the initial energy integration are from ${\epsilon }_{0{\rm{L}},{\rm{i}}}$ to ${\epsilon }_{0{\rm{U}},{\rm{i}}}$, given by
    Equation (375)
    where ${\epsilon }_{0{\rm{U}}}$ is given by Equation (372) evaluated at the boundary value ${\epsilon }_{0k+\displaystyle \frac{1}{2}}$. The limits of integration for the final energy are from ${\epsilon }_{0{\rm{L}},{\rm{f}}}$ to ${\epsilon }_{0{\rm{U}},{\rm{f}}}$, where
    Equation (376)
    where ${\epsilon }_{0{\rm{L}}}$ is given by Equation (372) evaluated at ${\epsilon }_{0k+\displaystyle \frac{1}{2}}$. A 32-point Gauss–Legendre angular quadrature appropriately weighted by the Legendre polynomials ${P}_{{\ell }}(\cos \theta )$ is then executed to obtain the Legendre coefficients.
  • 3.  
    If $k^{\prime} \lt k-1$, the differences between initial and final energy falls outside the small energy scattering range, so the necessity of refining the energy grids is less critical. The zero and first Legendre moments are obtained from Equation (371) by first integrating the scattering functions over the initial energy from ${\epsilon }_{0k-\tfrac{1}{2}}$ to ${\epsilon }_{0k+\tfrac{1}{2}}$ and the final energy grid from ${\epsilon }_{0k^{\prime} -\tfrac{1}{2}}$ to ${\epsilon }_{0k^{\prime} +\tfrac{1}{2}}$ using four-point Gauss–Legendre energy quadratures. The Legendre moments are then obtained by a four-point Gauss–Legendre angular quadrature.

We have also included corrections for weak magnetism to the scattering rates by the procedure outlined by Buras et al. (2006b), which disentangles the weak magnetism correction from the corrections given by Horowitz (2002) for both weak magnetism and recoil, ${\chi }_{\mathrm{WM},\mathrm{Rec}}^{\mathrm{nc},{\rm{n}},{\rm{p}}}$, and the corrections given for recoil only, ${\chi }_{\mathrm{Rec}}^{\mathrm{nc},{\rm{n}},{\rm{p}}}$. The resulting weak magnetism correction factor, ${\xi }^{\mathrm{nc}}({\epsilon }_{0})$, is given by the ratio

Equation (377)

where ${n}_{{\rm{n}}}$ and ${n}_{{\rm{p}}}$ are the number densities of free neutrons and protons, respectively. The Reddy rates, which do not include weak magnetism, are corrected for weak magnetism by multiplying the scattering Legendre moments by ${\xi }^{\mathrm{nc}}({\epsilon }_{0})$.

The inverse mean free paths and mean relative energy transfers are plotted for ${\nu }_{e}$ in Figure 30 and for ${\bar{\nu }}_{e}$ in Figure 31. The figures compare the scattering opacities computed with the procedure described above (solid lines) with those computed without including weak magnetism (dotted lines) and with those computed without including weak magnetism and with the isoenergetic approximation described by Bruenn (1985) (dashed lines). The effect of including recoil, nucleon final-state blocking, and special relativity (Reddy et al. 1998) is to reduce the rates compared with the isoenergetic approximation. This is mainly because the isoenergetic approximation is equivalent to assuming that the nucleon has an infinite mass. Including recoil takes into account the finite nucleon mass and, therefore, reduces the center of mass energy of the colliding system, thus reducing the opacity for a given incident neutrino energy. Final-state blocking and special relativity work in the same direction but are less important for the conditions considered here. The effect of weak magnetism is to increase the opacities for neutrinos and decrease them for antineutrinos. The net effect of including recoil, blocking, and relativistic corrections and weak magnetism corrections on the magnitude of the opacities is, thus, considerably more pronounced for antineutrinos than for neutrinos. The effect of including the above corrections is essential for including the neutrino–nucleon energy transfer. The mean relative energy transfers are shown in Figures 30 and 31. They are obviously zero when computed using the isoenergetic approximation, but they are nonzero and can play a significant role when recoil is taken into account. The effect of weak magnetism is to modify the magnitude of the cross sections but has little if any effect on the relative energy transfer.

Figure 30.

Figure 30. Inverse mean free paths for scattering on free nucleons by ${\nu }_{e}$-neutrinos (black, scale on left) and the relative energy transfer to the neutrino (red, scale on right), for the thermodynamic conditions listed on the upper left of the plots. The dashed lines give the inverse mean free paths using the isoenergetic approximation of Bruenn (1985) uncorrected for weak magnetism, the dotted lines give the inverse mean free paths as given by the formalism of Reddy et al. (1998), which include recoil, nucleon final-state blocking, and special relativity but not weak magnetism corrections, and the solid lines are the latter inverse mean free paths corrected for weak magnetism. The inverse mean free paths plotted by the solid, dashed, and dotted lines were computed using a 1500-zone neutrino energy grid evenly spaced between 1 and 300 MeV. The symbols show the inverse mean free paths computed with recoil, etc. and weak magnetism corrections, using the typical Chimera energy grid of 20 zones geometrically spaced between 4 and 250 MeV.

Standard image High-resolution image
Figure 31.

Figure 31. Same as Figure 30 but for ${\bar{\nu }}_{e}$–nucleon scattering.

Standard image High-resolution image

As a reference, the dotted, dashed, and solid lines showing the inverse mean free paths and relative energy transfers were all obtained using an energy grid of 1500 zones, evenly spaced from 1 to 300 MeV. An energy grid of the same number of zones but geometrically spaced in the same energy range, rather than evenly spaced, gave very similar results except for very low incident neutrino energies ($\leqslant 2$ MeV). The filled circles show the inverse mean free paths and relative energy transfers computed with the above corrections using the typical Chimera energy grid of 20 energy zones (with sub-grids as described above) increasing geometrically between 4 and 250 MeV. The inverse mean free paths computed with the Chimera energy grid reproduce nicely those computed with the much more refined grid. The relative energy transfers, however, are somewhat overestimated at low energies. The reason, most likely, is that Chimera stores the lowest two moments of the scattering functions, which are computed using the sub-grid and from which the inverse mean free paths are directly related, while the relative energy transfers are not stored but computed from the relatively coarse 20-energy-zone Chimera grid.

To derive Equations (372) and (373), we begin with the dynamic structure function $S({q}_{0},\bar{q})$ for neutrino–nucleon scattering by setting $\hat{\mu }={\mu }_{2}-{\mu }_{4}=0$ and $z={q}_{0}/T$ in the expression for neutrino–nucleon absorption (Equations (21) and (22) of Reddy et al. 1998) to obtain

Equation (378)

where

Equation (379)

and where ${q}_{0}={\epsilon }_{0}^{\mathrm{in}}-{\epsilon }_{0}^{\mathrm{out}}$, $\bar{q}={\bar{p}}^{\mathrm{in}}-{\bar{p}}^{\mathrm{out}}$ are, respectively, the energy and momentum transferred by the neutrino, the subscripts "2" and "4" refer to the incident and final nucleon, respectively, μ is the nucleon chemical potential, M is the nucleon mass, and θ is the angle between the incident and scattered neutrino directions. In the nondegenerate limit $(\mu \ll 0)$

Equation (380)

where we have set ${\mu }_{2}={\mu }_{4}=\mu $. Using Equation (380) in Equation (378), we get

Equation (381)

We wish to find upper and lower values, ${\epsilon }_{0{\rm{U}}}$ and ${\epsilon }_{0{\rm{L}}}$, such that

Equation (382)

Neglecting the third term in the numerator of the exponential, Equation (382) is a quadratic in ${\epsilon }_{\mathrm{out}}$ with the two solutions given by Equation (372) above.

8.6. Neutrino Absorption and Emission on Free Nucleons

Emission and absorption of neutrinos is an important process, for which we plot in Figure 32, for ${\nu }_{e}$ the inverse mean free paths (black lines), given by

Equation (383)

where ${{\chi }_{{\rm{a}}}}^{(0)}$ and ${{\chi }_{{\rm{e}}}}^{(0)}$ are the absorption and emission inverse mean free paths, respectively, and ${\epsilon }_{0{\rm{e}}}-{\epsilon }_{0\nu }$ (red lines) is the difference between the energies of the emitted electron and the absorbed neutrino. Figure 33 shows the same for ${\bar{\nu }}_{e}$. The dashed lines give the inverse mean free paths for the isoenergetic approximation of Bruenn (1985), uncorrected for weak magnetism. The dotted lines give the inverse mean free paths as given by Reddy et al. (1998), which include recoil, nucleon final-state blocking, and special relativity but not weak magnetism corrections. The solid lines give the inverse mean free paths for the latter but corrected for weak magnetism. Like neutrino scattering on free nucleons, the inverse mean free path for absorption and emission on free nucleons is reduced when recoil is taken into account by the reduced center of mass energy of the collision. At very low incident energies the emitted electron energy, ${\epsilon }_{0{\rm{e}}}$, tends to be greater than the incident ${\nu }_{e}$ energy, ${\epsilon }_{0\nu }$, because of the neutron–proton mass difference and the thermal motions of the nucleons. At high energies, ${\epsilon }_{0{\rm{e}}}$ decreases below ${\epsilon }_{0\nu }$ due to part of the incident collision energy being taken up by the final nucleon. The same is true for ${\bar{\nu }}_{e}$ except that the neutron–proton mass difference is negative in this case. In the isoenergetic approximation, ${\epsilon }_{0{\rm{e}}}-{\epsilon }_{0\nu }$ is just the neutron–proton mass difference, positive for ${\nu }_{e}$, and negative for ${\bar{\nu }}_{e}$.

Figure 32.

Figure 32. Inverse mean free paths for the absorption/emission on free nucleons of ${\nu }_{e}$-neutrinos (black, scale on left) and the difference in energy between the absorbed neutrino and emitted electron (red, scale on right), for the thermodynamic conditions listed on the upper left of the plots. The dashed lines give the inverse mean free paths for the isoenergetic approximation of Bruenn (1985), uncorrected for weak magnetism, the dotted lines give the inverse mean free paths as given by (Reddy et al. 1998), which include recoil, nucleon final-state blocking, and special relativity but not weak magnetism corrections, and the solid lines are the latter inverse mean free paths corrected for weak magnetism. The inverse mean free paths plotted by the solid, dashed, and dotted lines were computed using a neutrino energy grid of 1500 zones spaced between 1 and 300 MeV. The symbols (circles for ${\nu }_{e}$, triangle for ${\bar{\nu }}_{e}$) show the inverse mean free paths computed with the typical Chimera energy grid of 20 zones geometrically spaced between 4 and 250 MeV.

Standard image High-resolution image
Figure 33.

Figure 33. Same as Figure 32 but for ${\bar{\nu }}_{e}$ absorption/emission on free nucleons.

Standard image High-resolution image

As a reference, the $1/\lambda $ and ${\epsilon }_{0{\rm{e}}}-{\epsilon }_{0\nu }$ shown by the dashed, dotted, and solid lines were computed for a given angle θ between the incident neutrino and emitted electron, by integrating the final electron energy using an energy grid of 1500 zones of equal width spaced between 1 and 300 MeV. The first two Legendre moments of the absorption and emission kernels were then computed with a 64-point Gauss–Legendre angular quadrature. The filled circles show $1/\lambda $ and ${\epsilon }_{0{\rm{e}}}-{\epsilon }_{0\nu }$ as computed by Chimera with the typical energy grid of 20 zones geometrically spaced between 4 and 250 MeV. To compute $1/\lambda $, Chimera uses a 64-point Gauss–Legendre quadrature to integrate over the final electron energy, with limits given by Equation (372) above with the term ${dm}\,{c}^{2}=({M}_{1}-{M}_{2}){c}^{2}$ added, where M1 and M2 are the initial and final nucleon masses, respectively. A 64-point Gauss–Legendre angular quadrature is then used to compute the first two Legendre moments of the absorption and emission inverse mean free paths.

8.7. Neutrino Absorption and Emission on Nuclei

Calculation of the rate of electron capture on heavy nuclei and the resulting neutrino emission in the collapsing core requires three components: the appropriate electron capture reaction rates, the spectra of emitted neutrinos, and knowledge of the nuclear composition. In simulations of the collapsing stellar iron core, the composition is calculated by the equation of state assuming NSE, instead of being tracked in detail via a reaction network. As discussed in Section 3.2, the information on the nuclear composition typically provided by the equation of state is limited to the mass fractions of free neutrons and protons, α-particles, and the sum of all heavy nuclei, as well as the identity of an average heavy nucleus. In Chimera, we use a prescription for nuclear electron capture first utilized by Langanke et al. (2003) and Hix et al. (2003). This treatment is based on shell model electron capture rates from Langanke & Martínez-Pinedo (2000, LMP) for $45\lt A\leqslant 65$ and 80 reaction rates from a hybrid Shell-Model-Monte-Carlo (SMMC)—Random-Phase-Approximation (RPA) calculation (Langanke et al. 2001a, 2003, LMS) for a sample of nuclei with $66\leqslant A\leqslant 112$. The approximation of Langanke et al. (2001b) is used for the distribution of emitted neutrinos. To calculate the needed abundances of the heavy nuclei, a Saha-like NSE is assumed, including Coulomb corrections to the nuclear binding energy (Hix & Thielemann 1996; Bravo & García-Senz 1999), but neglecting the effects of degenerate nucleons (El Eid & Hillebrandt 1980). The combined set of LMP and LMS model rates is used to calculate an average neutrino emissivity per heavy nucleus. The full neutrino emissivity is then the product of this average and the number density of heavy nuclei calculated by the equation of state. With the limited coverage of rates for $A\gt 65$, this approach provides a reasonable estimate of what the total electron capture rate would be if rates for all nuclei were available. This averaging approach also makes the rate of electron capture consistent with the composition returned by the equation of state, while minimizing the impact of the limitations of our NSE treatment. A public version of an updated version of this rate tabulation is available from Juodagalvis et al. (2010). This tabulation, which is planned for inclusion in a future version of Chimera, includes more extensive coverage of heavier, more neutron-rich nuclei via a Fermi–Dirac (FD) parameterization of level occupation in place of the more costly SMMC approach.

In Figure 34, we plot the ${\nu }_{e}$ inverse mean free paths for absorption and emission on nuclei for both the IPM as formulated by Bruenn (1985) and the more sophisticated LMP–LMS formulations. Both the solid and dashed lines are the results of computations using a linear energy grid of 200 zones aligned with the energy grid of the electron capture table on which the LMP–LMS inverse mean free paths were tabulated. The filled circles are the results obtained using the typical Chimera energy grid. At densities below a few times 1010 ${{\rm{g\; cm}}}^{-3}$, the inverse mean free paths given by the IPM dominate (top frame of Figure 34), but at higher densities, electron capture reduces ${Y}_{{\rm{e}}}$ and drives the nuclear abundances toward neutron richness, including nuclei with neutron numbers $N\gt 40$. The IPM results in the vanishing of the electron capture inverse mean free paths for $N\geqslant 40$ nuclei due to the filling of the neutron ${}^{1}{f}_{5/3}$ orbital. In the LMS treatment of electron captures on nuclei, it was shown that Pauli blocking due to the filling of the neutron ${}^{1}{f}_{5/3}$ orbital is overcome by correlations and temperature effects. Consequently, at densities above a few times 1010 ${{\rm{g\; cm}}}^{-3}$, the LMS–LMP inverse mean free paths remain finite, whereas those given by the IPM vanish (see Figure 34(b) and (c)).

Figure 34.

Figure 34. Inverse mean free paths for the ${\nu }_{e}$ emission (black lines and circles) and ${\nu }_{e}$ absorption (green lines and circles) on nuclei, for the thermodynamic conditions listed on the upper left of the plots. The dashed lines give the inverse mean free paths calculated from the independent particle model (IPM) as formulated by Bruenn (1985). The solid lines give the inverse mean free paths computed from tables based on the LMS-LMP formulation. Data plotted by both the solid and dashed lines were computed using a linear energy grid of 200 zones from 0 to 100 MeV, exactly matching the energy grid of the LMS-LMP electron capture table. Filled circles show the inverse mean free paths computed with the typical Chimera energy grid of 20 zones geometrically spaced between 4 and 250 MeV.

Standard image High-resolution image

8.8. Pair Production: General

In spherical symmetry, the rate of change of the neutrino (antineutrino) occupation probability, $f({\mu }_{0},{\epsilon }_{0})$, due to pair production process "XX" is given by

Equation (384)

where $\bar{f}({\mu }_{0}^{{\prime} },{\epsilon }_{0}^{{\prime} })$ is the antineutrino (neutrino) occupation probability, ${R}_{\mathrm{XX}}^{{\rm{a}}}({\epsilon }_{0},{\epsilon }_{0}^{{\prime} },\cos \theta )$ and ${R}_{\mathrm{XX}}^{{\rm{p}}}({\epsilon }_{0},{\epsilon }_{0}^{{\prime} },\cos \theta )$ are, respectively, the neutrino–antineutrino annihilation and neutrino–antineutrino unblocked-creation rates per neutrino–antineutrino states for process XX, ${\epsilon }_{0}$ is the neutrino energy, ${\mu }_{0}$ is the cosine of the neutrino propagation direction with respect to the radial direction, ϕ is the azimuthal propagation directions, and $\cos \theta $ is defined in Equation (357). Unprimed quantities refer to neutrinos (antineutrinos), and primed quantities refer to antineutrinos (neutrinos). ${R}_{\mathrm{XX}}^{{\rm{a}}}$ and ${R}_{\mathrm{XX}}^{{\rm{p}}}$ are related by

Equation (385)

where ${E}_{{{\rm{e}}}^{-}}$ and ${E}_{{{\rm{e}}}^{+}}$ are the associated electron and positron energies, respectively. Expanding the annihilation and creation kernels in a Legendre expansion and keeping the first two terms, as was done in Equation (360) for the scattering kernels, gives

Equation (386)

Applying the moment operators ${\left(4\pi \right)}^{-1}\int d{\rm{\Omega }}$ and ${\left(4\pi \right)}^{-1}\int {\mu }_{0}d{\rm{\Omega }}$ to Equation (384) and using the definitions in Equations (357), (386), (384), and (220), gives the moments of the pair interaction terms:

Equation (387)

The pair annihilation kernels in Equation (387) have been multiplied by angular cutoff correction factors, $\chi ({\epsilon }_{0},{\epsilon }_{0}^{{\prime} })$, to account for the non-isotropic nature of the neutrino flow (e.g., Goodman et al. 1987; Cooperstein et al. 1987). For Chimera, these angular correction factors are derived as follows. The pair neutrino–antineutrino annihilation rate is proportional to the square of the center-of-mass energy, or to ${(1-{{\boldsymbol{n}}}_{\nu }\cdot {{\boldsymbol{n}}}_{\bar{\nu }})}^{2}$, where ${{\boldsymbol{n}}}_{\nu }$ and ${{\boldsymbol{n}}}_{\bar{\nu }}$ are unit vectors in the propagation direction of ν and $\bar{\nu }$, respectively. Inside the neutrinosphere, the anisotropy of the neutrino distributions are small, and the occupation distributions $f$ and $\bar{f}$ can be approximated by

Equation (388)

which follows from the definitions in Equation (220). Using Equation (357) for ${{\boldsymbol{n}}}_{\nu }\cdot {{\boldsymbol{n}}}_{\bar{\nu }}$, we find that the ratio, ${\chi }_{1}$, of ${(1-{{\boldsymbol{n}}}_{\nu }\cdot {{\boldsymbol{n}}}_{\bar{\nu }})}^{2}$ evaluated with an anisotropic neutrino distribution in the numerator and an isotropic neutrino distribution in the denominator, having the same number and energy spectrum, is given by

Equation (389)

Outside the neutrinospheres, the neutrino occupation distributions become increasingly anisotropic and the neutrino–antineutrino center of mass collision energy becomes increasingly smaller as the beaming becomes radial and more nearly collinear. Assuming that the neutrinosphere emits isotropically and that f is constant along rays leading back to the neutrinosphere, outside the neutrinosphere, f is given by

Equation (390)

Equation (391)

where ${R}_{\nu }({\epsilon }_{0})$ and ${R}_{\bar{\nu }}({\epsilon }_{0}^{{\prime} })$ are the radii of the neutrinospheres of the neutrinos and antineutrinos, respectively, and R is the radial coordinate at which the correction is being applied. Using Equations  (390) and (391) for the anisotropic neutrino radiation, we find that the ratio, ${\chi }_{2}$, of ${(1-{{\boldsymbol{n}}}_{\nu }\cdot {{\boldsymbol{n}}}_{\bar{\nu }})}^{2}$ evaluated, as before, with an anisotropic neutrino distribution in the numerator and an isotropic neutrino distribution in the denominator, having the same number and energy spectrum, is given by

Equation (392)

We note that in the case ${\mu }_{0\,i}^{{\prime} }({\epsilon }_{0})={\mu }_{0i}({\epsilon }_{0})$, ${\chi }_{2}({\epsilon }_{0},{\epsilon }_{0})$ reduces to

Equation (393)

To smoothly transition between ${\chi }_{1}$ and ${\chi }_{2}$, we define the neutrino–antineutrino pair annihilation angular corrections, χ, to be given by

Equation (394)

and extend the definitions of ${\mu }_{0i}$ (and ${\mu }_{0\,i}^{{\prime} }$) in Equation (390) (and (391)) to

Equation (395)

(similarly for ${\mu }_{0\,i}^{{\prime} })$. Equation (395) will cause ${\mu }_{0i}$ to rapidly approach −1 and ${\chi }_{2}$ to approach 1 as R decreases below ${R}_{\nu }$, and Equation (394) will then cause χ to be governed by ${\chi }_{1}$.

Figure 35.

Figure 35. Pair annihilation rates computed by Chimera without (green) and with (red) angular cutoff correction factors, compared with the rates computed by Boltztran with eight-point (black), 12-point (orange), and 16-point (blue) angular quadrature.

Standard image High-resolution image

To test the above angular cutoff correction factor, we compare in Figure 35 the pair annihilation rates computed by Chimera without and with the angular cutoff correction factor to the rates computed by the Boltzman solver Boltztran with three different beam resolutions. Because the angular dependence of the transport has not been integrated over in Boltztran, the angular resolution is limited only by the number of beams employed, and the Boltztran rates can therefore be regarded as a reference modulo the number of beams employed. The core configuration used in this comparison was obtained by evolving the progenitor s15s7b2 (Woosley & Weaver 1995) in 1D to 100 ms post-bounce, at which point the core configuration was frozen and the neutrino transport was continued until a steady state was achieved. Without the angular cutoff correction factor, the pair annihilation rates computed by Chimera become significantly larger than the Boltztran rates beyond 70 km, which is 10 km beyond the mean radius of the neutrinosphere. They become several orders of magnitude larger beyond 170 km. With the angular cutoff correction factor included, the behavior of the pair annihilation rates computed by Chimera more closely approximates those of Boltztran, though differences of the order of five are seen beyond 120 km.

Figure 36.

Figure 36. Inverse mean free path (black) for electron–positron annihilation and pair production for ${\nu }_{e}$ (solid) and ${\bar{\nu }}_{e}$ (dashed), computed from Equation (396) on a 600-point evenly spaced energy grid from 0–300 MeV, for listed thermodynamic conditions. Inverse mean free paths were computed with the ${\nu }_{e}$ (${\bar{\nu }}_{e}$) occupation distribution set to zero and the ${\bar{\nu }}_{e}$(${\nu }_{e}$) distribution equilibrated with the matter and energy production rates per baryon for ${\nu }_{e}$ (solid) and ${\bar{\nu }}_{e}$ (dashed) computed from Equation (397) with both the ${\nu }_{e}$ and ${\bar{\nu }}_{e}$ occupation distributions set to zero. Filled circles (${\nu }_{e}$) and triangles (${\bar{\nu }}_{e}$) are for values computed with a typical Chimera energy grid of 20 logarithmic energy zones.

Standard image High-resolution image

8.9. Neutrino–Antineutrino Pair Annihilation and Production from Electron–Positron Pairs

The zero and first moments of the kernels for neutrino–antineutrino pair annihilation into electron–positron pairs and the inverse process are taken from the analytic expressions of Bruenn (1985, Equations (C62)–(C74), with a typo corrected by removing the term $-{a}_{0}$ from the bracket in Equation (C68)). The integration over the electron energy, ${E}_{{\rm{e}}}$, was performed for the case that ${\epsilon }_{0}\gt {\epsilon }_{0}^{{\prime} }$ by a 24-point Gauss–Legendre energy quadrature for each of the intervals from zero to ${\epsilon }_{0}^{{\prime} }$, from ${\epsilon }_{0}^{{\prime} }$ to ${\epsilon }_{0}$, and from ${\epsilon }_{0}$ to ${\epsilon }_{0}+{\epsilon }_{0}^{{\prime} }$, with a similar set of integrations for the case that ${\epsilon }_{0}^{{\prime} }\gt {\epsilon }_{0}$.

In Figure 35, we plot the inverse mean free path,

Equation (396)

where for neutrinos (antineutrinos) we have taken ${\psi }_{0\,\nu }^{(0)}({\epsilon }_{0})$ for the neutrinos (antineutrinos) to be zero and ${\bar{\psi }}_{0}^{(0)}({\epsilon }_{0}^{{\prime} })$ for the antineutrinos (neutrinos) to be in thermal equilibrium with the matter at the stated thermodynamic conditions, and plotted is the neutrino (antineutrino) spectral energy production rate per baryon,

Equation (397)

where we have set the neutrino–antineutrino blocking factors to zero.

8.10. Neutrino–Antineutrino Pair Annihilation and Production from Nucleon–Nucleon Bremsstrahlung

Neutrino–antineutrino pair annihilation and production from nucleon–nucleon bremsstrahlung has been computed at various levels of approximation (Friman & Maxwell 1979; Raffelt & Seckel 1995; Hannestad & Raffelt 1998). The last paper examined the process and provided an interpolation scheme that interpolates the nucleon–nucleon bremsstrahlung kernel between the nucleon nondegenerate limit treated by Raffelt & Seckel (1995) and the degenerate limit treated (for the case of axion emission) by Ishizuka & Yoshimura (1990). We use Equation (35) in Equation (23) of Hannestad & Raffelt (1998), which naturally breaks up into ${{\rm{\Phi }}}_{0,\mathrm{brem}}^{{\rm{a}}/{\rm{p}}}({\epsilon }_{0},{\epsilon }_{0}^{{\prime} })$ and ${{\rm{\Phi }}}_{1,\mathrm{brem}}^{{\rm{a}}/{\rm{p}}}({\epsilon }_{0},{\epsilon }_{0}^{{\prime} })$ terms. We also use Equation (394) for the neutrino–antineutrino pair annihilation correction. In Figure 37, we plot the inverse mean free paths, given by Equation (396), and the spectral energy production rates per baryon, given by Equation (397), with the nucleon–nucleon bremsstrahlung kernels substituted for the electron–positron pair annihilation kernels for select thermodynamic states.

Figure 37.

Figure 37. Inverse mean free paths for ${\nu }_{e}$ (black solid lines) and ${\bar{\nu }}_{e}$ (black dashed lines) computed from Equation (396) for nucleon–nucleon-bremsstrahlung kernels on an energy grid of 600 evenly spaced zones from 0 to 300 MeV, for the thermodynamic conditions listed. The ${\nu }_{e}$ (${\bar{\nu }}_{e}$) inverse mean free paths were computed with the ${\nu }_{e}$ (${\bar{\nu }}_{e}$) occupation distribution set to zero and the ${\bar{\nu }}_{e}$(${\nu }_{e}$) distribution equilibrated with the matter. Inverse mean free paths were computed with the Chimera energy grid of 20 energy zones for ${\nu }_{e}$ (black circles) and ${\bar{\nu }}_{e}$ (black triangles). Energy production rates per baryon from ${\nu }_{e}$ (solid red lines) and ${\bar{\nu }}_{e}$ (dashed red lines) were computed from Equation (397) for the nucleon–nucleon-bremsstrahlung kernels on the above-described 600 zone energy grid, with both the ${\nu }_{e}$ and ${\bar{\nu }}_{e}$ occupation distributions set to zero. The spectral energy production from ${\nu }_{e}$ (red circles) and ${\bar{\nu }}_{e}$ (red triangles) was computed with the 20 energy zone Chimera energy grid.

Standard image High-resolution image

Mathematically, inelastic neutrino scattering on nucleons is the bremsstrahlung process with a final-state neutrino crossed into the initial state and could be included with the same kernel. Moreover, this process should reduce to elastic neutrino–nucleon scattering when the nucleon–nucleon interaction tends to zero. However, the formulation of Hannestad & Raffelt (1998) neglects nucleon recoil, and we have used instead of the Hannestad & Raffelt (1998) formalism for inelastic neutrino–nucleon scattering the formalism of Reddy et al. (1998), described earlier, for elastic neutrino scattering on nucleons, which includes recoil effects as well as nucleon degeneracy and relativistic effects. At very high densities, where inelastic neutrino scattering on nucleons becomes important, we have therefore neglected this mode of energy exchange between the neutrinos and the matter, with the hope of recovering a significant part of it through the inclusion of nucleon recoil.

8.11. Opacity Interpolation

As noted earlier in this section, the Chimera method for building and interpolating neutrino interactions in $(\mathrm{log}\rho ,\mathrm{log}T,{Y}_{{\rm{e}}})$ uses a local cube of points in $(\rho ,T,{Y}_{e})$ that bracket the $(\rho ,T,{Y}_{e})$-value of each spatial cell and that lie on a regular mesh in $(\mathrm{log}\rho ,\mathrm{log}T,{Y}_{{\rm{e}}})$-space (cf. Section 3.1; Mezzacappa & Bruenn 1993a). The derivatives computed within the local $(\rho ,T,{Y}_{e})$-cube are fully consistent with the interpolated value. This method has a disadvantage with the on-the-fly re-computation of opacities to fill the local $(\rho ,T,{Y}_{e})$-cubes when the values no longer bracket the cell values. Computation of some of the opacities is expensive, notably the opacities with the sub-grid integration described in Section 8.5, and could consume up to ∼25%–40% of the simulation time. This cost was compounded by computational load imbalance from the irregular number of local $(\rho ,T,{Y}_{e})$-cubes updates required for each transport time step, which typically numbered between zero and a few new $(\rho ,T,{Y}_{e})$-cubes.

Starting with the C-series models, we implemented the scheme described below, which greatly reduced the on-the-fly computations and is better suited to future use of pre-computed tables. Specifically, we chose to implement a sparse "local pool" of $(\rho ,T,{Y}_{e})$-tuples. Using such a "local pool," rather than retaining the existing individual $(\rho ,T,{Y}_{e})$-cubes for each cell with the addition of a "reuse" algorithm, has lower memory requirements, as many cells use the same logical cube of $(\rho ,T,{Y}_{e})$-points. The memory savings is potentially larger when computing transport on multiple adjacent "rays" that are similar in $(\rho ,T,{Y}_{e})$-space. The reduced memory usage was particularly helpful when using machines with smaller memory footprints (∼1–2 GB of memory per core).

8.11.1. Energy Interpolation

The energy grid used in the Chimera neutrino transport solver is not the comoving observer's energy, ${\epsilon }_{0}$, but the energy of a comoving observer outside the gravitational well of the supernova, $a{\epsilon }_{0}\equiv {E}_{0}$ (Equation (216)), which depends on radius and, thus, complicates the sharing of opacity $(\rho ,T,{Y}_{e})$-points. To account for this difference, we must interpolate the energy grid of the opacity table into the specific points needed for the $\alpha \epsilon $ grid in each group, from a reference grid, without changing the total interaction rate.

The energy grid, both ${\epsilon }_{0k}$ and ${E}_{0k}$, are logarithmically spaced at the group centers

Equation (398)

(and group edges, with k replaced by $k+1/2$ above). Likewise, the difference between the two energy grids for a non-unity lapse can be written as a logarithmic shift

Equation (399)

for any value of the energy.

8.11.2. A Simple Integral Scheme

Informed by the requirements of the Chimera energy grid, we can construct a more general interpolation scheme for any function, f, tabulated on a logarithmically spaced grid, ${\varepsilon }_{\bar{{\imath }}}$, to another logarithmically spaced grid, ${\varepsilon }_{i}$, that differs by a multiplicative constant that changes with time, ${\varepsilon }_{\bar{{\imath }}}={\varepsilon }_{i}/\alpha $, where both grids have the same logarithmic spacing

Equation (400)

Here and in the following, we will omit the subscript naught and all subscripts will refer to the energy grid. The integral scheme comes from the need to evaluate integrals in the form $\int f({\varepsilon }_{})d{\varepsilon }_{}$ numerically by summation of the terms ${\rm{\Sigma }}{f}_{i}\delta {\varepsilon }_{i}$, where $\delta {\varepsilon }_{i}={\varepsilon }_{i+1}-{\varepsilon }_{i}$ and $\delta {\varepsilon }_{i}{f}_{i}$ is the integration of f over that same interval, $\delta {\varepsilon }_{i}$. When the shift, $\mathrm{log}(1/\alpha )$, is smaller than the grid spacing, $\mathrm{log}\zeta $, we can define the overlap between the $\delta {\varepsilon }_{i}$ and $\delta {\varepsilon }_{\bar{{\imath }}}$ zones as

Equation (401)

which is the same for all zones, i. The integral over $\bar{{\imath }}$ is then

Equation (402)

If the left-most zone has a coordinate ${\varepsilon }_{1}=0$, we must modify the above equation to include the whole of the first zone value

Equation (403)

For grid shifts that are larger than the single-zone spacing, $\mathrm{log}(1/\alpha )\gt \mathrm{log}\zeta $, a small modification is required. First, we define a "shift index"

Equation (404)

The β value must now reflect the larger shift, so

Equation (405)

where ${\zeta }^{j}$ refers to ζ to the jth power, which can be used in an extended version of Equation (402):

Equation (406)

The additional zones shifted for ${\varepsilon }_{1}=0$ grids must be added to Equation (403):

Equation (407)

These last two equations reduce to the earlier forms when j = 0.

8.11.3. Interpolation of Two-variable Functions

If we have a function, $f({\varepsilon }_{i},{\varepsilon }_{{\imath }^{\prime} })$, which needs to be remapped to a shifted grid, $f({\varepsilon }_{\bar{{\imath }}},{\varepsilon }_{\bar{{\imath }^{\prime} }})$, we can extend the single-variable method given above. We assume that both the ${\varepsilon }_{i}$ and ${\varepsilon }_{{\imath }^{\prime} }$ grids are the same (i.e., ${\varepsilon }_{i}={\varepsilon }_{{\imath }^{\prime} }$ if $i={\imath }^{\prime} $) and therefore have the same shift and spacing parameters, $(\zeta ,\beta ,1/\alpha ,j)$. We start with a partially shifted version of Equation (402) where i is shifted but ${\imath }^{\prime} $ is not, and then apply the shift to ${\imath }^{\prime} $ to get the shifted double-grid formula:

Equation (408)

For grids shifted such that the shift index $j\gt 0$, we can generalize Equation (408) to

Equation (409)

for cases when $i\gt 1,{\imath }^{\prime} \gt 1$. When either i = 1 or ${\imath }^{\prime} =1$, we need to generalize the equations above:

Equation (410)

Equation (411)

Equation (412)

8.11.4. Interpolation of Chimera Opacities

For the interpolation of opacities in Chimera, we choose to store the raw opacities on a grid, ${\varepsilon }_{}$, that numerically matches the specified E0 grid, which makes $\alpha =a$. Note that this is not a frame transform, just a convenient choice of variables, as all opacities are evaluated in the ${\epsilon }_{0}$ frame of the moving fluid. For single-energy opacities (absorption and emission), f is the inverse mean free path multiplied by ${\epsilon }_{0}^{2}$. For two-energy opacities (scattering and pair processes), f is the Legendre coefficients of the kernels multiplied by ${\epsilon }_{0}^{2}{\epsilon }_{0}{{\prime} }^{2}$. The ${\epsilon }_{0}{{\prime} }^{2}$ term arises from within the collisions integrals, while the ${\epsilon }_{0}^{2}$ arises from integrating the collision integral with the operator $\int {\epsilon }_{0}^{2}d{\epsilon }_{0}$ to conserve the total integral interaction rate.

8.11.5. Local Pool Algorithm

The local pool is implemented for each opacity by first identifying the eight $(\rho ,T,{Y}_{e})$-points needed for each zone using that opacity on each MPI rank by generating an integer hash value that uniquely maps to the potential grid of $(\rho ,T,{Y}_{e})$-points, removing duplicates, and sorting them. The sorted list is then compared to the points in the existing pool to generate a list of points that need to be added to the pool. The list of additions is added to the pool by calling the opacity generation routines for each missing point. The integer hash values for the eight $(\rho ,T,{Y}_{e})$-points for each zone are then regenerated and cross-linked to their pool index numbers. During interpolation, the eight points indicated by the pool index numbers are used to interpolate the kernels for specific cells.

Implementation of this sharing algorithm reduced the opacity generation and management costs to a few percent of the total run time. It had the extra benefit of reducing the opacity interpolation costs by 5%–10%, as the $(\rho ,T,{Y}_{e})$-points needed for interpolating the opacity for one cell frequently fully, or partially, overlap with those needed for the next cell, allowing data stored in processor cache to be reused. This was not possible under the old scheme, as every zone had its own eight points even if they were duplicates of a neighboring cell's points. Implementation was checked with tests that were computed using a 1D reference simulation, with comparisons of critical transport and hydrodynamics variables made at core bounce and at several other times.

9. Comparison with Other Codes

Supernova codes are complex entities, involving the numerical solution of hydrodynamics supplemented by one or more equations of state, neutrino transport with multiple sources of absorption and scattering opacities, nuclear transmutations, and relativistic gravity, with fluid densities ranging over more than ten orders of magnitude. Different numerical techniques have their individual strengths and weaknesses, reinforcing the importance of code validation for a particular choice of techniques and their implementation. In this section, we make a detailed comparison of results of a spherically symmetric simulation performed by Chimera using the code versions utilized in our B-series and C-series simulations with the results of two other codes that have been compared in Liebendörfer et al. (2005): Agile-Boltztran and Prometheus-Vertex. A recent comparison of CCSN codes was presented in O'Connor et al. (2018), using updated physics, particularly the EoS, and a comparison of Chimera with those results will be reported in the future using a later version of Chimera than described here.

9.1. Description of Comparison Simulations and Codes

Agile-Boltztran was developed by an Oak Ridge–Basel collaboration and is a code that solves the general relativistic hydrodynamics equations in a time-implicit fashion, with a dynamical adaptive grid, coupled to a solver for the general relativistic Boltzmann equations for the neutrino distribution functions based on a discrete-ordinates (SN) in angle and finite-difference in space and neutrino energy discretization, all for spherically symmetric spacetimes (Mezzacappa & Bruenn 1993a; Mezzacappa & Messer 1999; Liebendörfer et al. 2002, 2004).

Prometheus-Vertex is a code for multidimensional neutrino radiation hydrodynamics, here used in the context of studies assuming spherical symmetry, developed by the Garching group. It consists of a code to solve the hydrodynamics equations, based on a finite-volume discretization Prometheus (Fryxell et al. 1989), coupled in an operator-split fashion to a code that solves the equations for two-moment neutrino transport, closed using a variable Eddington factor. This factor is derived from a formal solution of a model Boltzmann equation for the counterpart spherically averaged matter configuration corresponding to the actual multidimensional matter configuration for a given time step (Rampp & Janka 2002; Buras et al. 2006b).

From Chimera, we include runs from two versions of the code. The first, Chimera-B, uses the same code as the B-series models reported in Bruenn et al. (2013, 2016) but with microphysics (EoS and opacities) and progenitor (Woosley & Weaver 1995) to match the Liebendörfer et al. (2005) comparison models. The second, dubbed Chimera-C, is the code used for the C-series models (Lentz et al. 2015; E. J. Lentz et al. 2018, in preparation) and differs from Chimera-B in some code refinements, particularly the modification of the transport scheme in the vicinity of large velocity discontinuities (i.e., shocks) to properly account for the large changes in the comoving reference frames described in Section 6.12.

The Agile-Boltztran and Prometheus-Vertex comparisons have been documented in Liebendörfer et al. (2005), as noted above, and in Müller et al. (2010). The progenitor and choice of physics we will use for our comparisons is referred to as model G15 in Liebendörfer et al. (2005), and we also will refer to it as G15. Of the two models used for the comparisons in that work, it is the model that implements the more extensive choice of physics and, thus, more closely approximates the set of physics included in current state-of-the-art CCSN simulations. The progenitor used to initiate the G15 simulations is the 15 ${M}_{\odot }$ progenitor of Woosley & Weaver (1995), widely used in the literature, and representative of the star in the middle- to upper-mass range likely to end its life as a supernova. From the highest densities to a density of $6\times {10}^{7}\,{{\rm{g\; cm}}}^{-3}$, the EoS used in all codes is the compressible liquid drop model of Lattimer & Swesty (1991) with incompressibility modulus K = 180 MeV. It assumes a composition in NSE consisting of neutrons, protons, α-particles, a representative heavy nucleus, electrons, positrons, and photons. At densities below $6\times {10}^{7}\,{{\rm{g\; cm}}}^{-3}$ each code switches to an EoS consisting of a composition of electrons, positrons, photons, nucleons, and nuclei, with the latter two treated as an ideal gas. The detailed treatment of the EoS in the lower-density regime differs among the codes; therefore, our comparisons will be limited to important phenomena occurring above $6\times {10}^{7}\,{{\rm{g\; cm}}}^{-3}$. Agile-Boltztran is fully general relativistic, while both Prometheus-Vertex and Chimera approximate the gravitational potential by including corrections terms due to pressure and energy of the stellar medium and neutrinos, as described in Marek et al. (2006). We plot the results of two Prometheus-Vertex simulations, which we label as Vertex-1 and Vertex-2. The Vertex-1 simulation was performed with gravitational potential "R" of Marek et al. (2006) and is described in Liebendörfer et al. (2005), while the Vertex-2 simulation was performed with gravitational potential "A" and is described in Müller et al. (2010). Both Prometheus-Vertex and Chimera include gravitational redshifting and time dilation in the transport solution but ignore the difference between coordinate and proper radial distances, which has little effect on the transport (Bruenn et al. 2001). The hydrodynamics in both codes is Newtonian.

Neutrinos of all flavors are included in model G15, with the neutrino–matter interactions shown in Table 3. The Agile-Boltztran simulation was performed using 103 adaptive spatial zones from the center to the edge of the included progenitor, which was about 7000 km, and a constant-pressure boundary condition was applied at the outer surface. The neutrino energy grid was resolved with 20 geometrically spaced bins, the first centered at 3 MeV and the last at 300 MeV. The propagation directions were discretized into six angles suitable for Gaussian quadrature. Neutrino flavors were divided into four independently transported species, ${\nu }_{e}$, ${\bar{\nu }}_{e}$, ${\nu }_{\mu \tau }$ (${\nu }_{\mu }$ and ${\nu }_{\tau }$), and ${\bar{\nu }}_{\mu \tau }$ (${\bar{\nu }}_{\mu }$ and ${\bar{\nu }}_{\tau }$). The Prometheus-Vertex simulation was performed on separate radial fluid and transport grids. The fluid grid consisted of 400 zones that moved with the fluid during collapse and were rezoned shortly after core bounce such that inside a radius of 400 km the fluid grid coincides with the transport grid. The latter consisted of an Eulerian radial grid of 235 radial zones spaced logarithmically between 0 and 10,000 km. Forty additional radial zones were added to both grids after 200 ms post-bounce to resolve the steepening density gradient at the surface of the nascent neutron star. The neutrino spectrum was discretized with 19 energy zones between 0 and 380 MeV. Neutrino flavors were divided into three independently transported species, ${\nu }_{e}$, ${\bar{\nu }}_{e}$, and ${\nu }_{x}$ (${\nu }_{\mu \tau }$ and ${\bar{\nu }}_{\mu \tau }$). The Chimera simulations were performed with 512 adaptive radial zones. The zones moved with the fluid during collapse and thereafter adjusted to maintain an approximately constant ${\rm{\Delta }}r/r$ while maintaining a maximum value of ${\rm{\Delta }}\rho /\rho $ to maintain good resolution in the vicinity of the steepening density gradient at the surface of the nascent neutron star as the simulation progresses past bounce. The neutrino spectrum was resolved with 20 energy zones between 0 and 279 MeV with mid-energies of 2.57 and 250 MeV in the first and last zone for the four species: ${\nu }_{e}$, ${\bar{\nu }}_{e}$, ${\nu }_{\mu \tau }$, and ${\bar{\nu }}_{\mu \tau }$. Because the Chimera neutrino energy grid is defined at infinity and the actual grid is radially dependent (the local grid energies are the energies at infinity divided by the local value of the lapse function), the grid expands toward the core center as the latter contracts to higher densities. At bounce, the mid-zone energies at the core center range from 2.92 to 285 MeV, and at 100 ms post-bounce, the mid-zone energies at the core center range from 3.20 to 313 MeV.

Table 3.  Summary of Neutrino Opacities in Model G15

Interaction References
$\nu +{e}^{\pm }\leftrightharpoons \nu +{e}^{\pm }$ Bruenn (1985), Mezzacappa & Bruenn (1993b)
$\nu +A\leftrightharpoons \nu +A$ Bruenn & Mezzacappa (1997), Horowitz (1997)
$\nu +n,p\leftrightharpoons \nu +n,p$ Bruenn (1985), Mezzacappa & Bruenn (1993b)
${\nu }_{e}+n\leftrightharpoons {e}^{-}+p$ Bruenn (1985), Mezzacappa & Bruenn (1993b)
${\bar{\nu }}_{e}+p\leftrightharpoons {e}^{+}+n$ Bruenn (1985), Mezzacappa & Bruenn (1993b)
${\nu }_{e}+A^{\prime} \leftrightharpoons {e}^{-}+A$ Bruenn (1985), Mezzacappa & Bruenn (1993b)
$\nu +\bar{\nu }\leftrightharpoons {e}^{-}+{e}^{+}$ Bruenn (1985), Mezzacappa & Bruenn (1993b)
$N+N\leftrightharpoons N+N+\nu +\bar{\nu }$ Hannestad & Raffelt (1998)
ion–ion correlations Itoh (1975), Horowitz (1997), Bruenn & Mezzacappa (1997)

Download table as:  ASCIITypeset image

9.2. Infall

The evolution of the entropy, electron fraction, and lepton fraction during infall is shown in Figure 38. Prior to the ∼10 ms after shock formation, the Chimera-B and Chimera-C results are essentially identical and are shown simply as "Chimera" in graphs until such times as the differences between them become significant. Before trapping, the entropy evolution is almost identical for all three simulations. Trapping occurs slightly later for Vertex-1 compared with Chimera. Trapping seems to occur for Agile-Boltztran near that of Chimera, but the entropy continues to slowly increase thereafter, likely because of zones moving away from the center to resolve the forming shock. This causes the central zone to encompassing more and more mass of higher entropy, causing the zone-average entropy to rise.

Figure 38.

Figure 38. Evolution of entropy (a), and electron (solid) and lepton (dashed) fractions (b) vs. density in the central zone during infall, for Chimera (red), Agile-Boltztran (black), and Vertex-1 (green) simulations.

Standard image High-resolution image

9.3. Bounce

Figure 39 displays the material and neutrino configurations of model G15 at bounce. The results of all simulations are quite similar at this time. Figure 39(a) plots the neutrino luminosities, with the ${\nu }_{e}$-luminosity for Chimera between 0.8 and 1.3 ${M}_{\odot }$ slightly lower than that for the others—a difference that may be related to the slightly lower core densities (Figure 39(b)) and slightly lower trapped entropy (Figures 39(d) and 38). The ${\bar{\nu }}_{e}$-luminosities are still too small to be shown, and the ${\nu }_{\mu \tau }$-luminosities are just beginning to develop. The ${\nu }_{e}$-rms energies (Figure 39(c)) are almost the same for all models, while the ${\nu }_{\mu \tau }$-rms energies are slightly smaller in the core in the Chimera simulation. The jump in the entropy and electron fraction at an enclosed mass of 1.18 ${M}_{\odot }$ in the Chimera simulation (Figure 39(d)) occurs at about 1.28 ${M}_{\odot }$ in the other two simulations. These jumps in the Chimera simulation appear in the initial model at the same enclosed mass, so it is not a feature that has evolved during infall by Chimera.

Figure 39.

Figure 39. Snapshots at bounce for (a) luminosity, (b) density and velocity, (c) rms-energies, and (d) entropy and ${Y}_{{\rm{e}}}$ of the G15 profile, for the Agile-Boltztran (black), Vertex-1 (green), and Chimera (red) simulations.

Standard image High-resolution image

9.4. Comparisons at 3 ms after Bounce

Figure 40 displays the material and neutrino configurations of Model G15 3 ms after bounce. Again, the results of all simulations are quite similar at this time, but small differences can be observed. The velocity profiles (Figure 40(b)) are almost identical, with the Chimera and Agile-Boltztran shock being slightly farther out in enclosed mass. The Chimera density in the region from 0.5 ${M}_{\odot }$ to the shock is slightly higher than the Vertex-1 density, which in turn is slightly higher than the Agile-Boltztran density. This probably accounts for the similar hierarchy in the ${\nu }_{e}$-rms energy in that region (Figure 40(c)). The Chimera shock at formation is slightly weaker than that of Vertex-1, which is slightly weaker than that of Agile-Boltztran, as inferred by the entropy profile behind the shock (Figure 40(d)), which likely accounts for the differences in the neutrino luminosities at this time (Figure 40(a)).

Figure 40.

Figure 40. As in Figure 39 but for 3 ms after bounce.

Standard image High-resolution image

9.5. Comparisons as a Function of Time

Figures 41(a-d) compare the shock trajectories, neutrino luminosities, and neutrino rms energies computed by the codes as a function of time from bounce to 300 ms post-bounce. Figure 41(a) plots the position of the shocks as a function of time. The outwardly directed "hump" in the shock trajectories at about 170 ms, not seen in the Agile-Boltztran results, is a feature that is due to the passage through the shocks of the silicon layer, with its associated drop in density and consequent reduction in the inwardly directed ram pressure on the shocks. It should be noted that after about 150 ms, the shock exits the region covered by the Lattimer-Swesty EoS and enters a region where the EoS is treated differently by each code, making comparisons problematic. The Chimera-B shock trajectory is very close to that of Agile-Boltztran for the first 60 ms post-bounce, then, after that, rises some 20 km above it. The Chimera-C shock trajectory is initially close to that of Agile-Boltztran, falls below it by up to 5 km from 30 to 80 ms, then stays within 2 km of it from 90 ms to the end of the plot. The Vertex-1 shock trajectory is initially within a few km of that of Agile-Boltztran, then retreats more rapidly after 70 ms post-bounce, falling below it by about 5 km by 100 ms post-bounce. The Vertex-2 shock trajectory follows closely but slightly above that of Agile-Boltztran. Both the Vertex-1 and Vertex-2 shocks exhibit the outwardly directed hump due to the passage through the shock of the silicon layer about 20–30 ms after those of the Chimera shocks.

Figure 41.

Figure 41. Evolution of the (a) model G15 shock radius and the (b) ${\nu }_{e}$, (c) ${\bar{\nu }}_{e}$, and (d) ${\nu }_{x}$ neutrino luminosities and rms energies, for Agile-Boltztran (black), Vertex-1 (green) and Vertex-2 (turquoise), Chimera-B (blue), and Chimera-C (red) simulations.

Standard image High-resolution image

Figures 41(b)–(d) plot, respectively, the ${\nu }_{e}$-, ${\bar{\nu }}_{e}$-, and one of the ${\nu }_{\mu \tau }$-, ${\bar{\nu }}_{\mu \tau }$-luminosities and rms energies as a function of time. The Chimera-B and Chimera-C ${\nu }_{e}$-luminosities exhibit a lower peak at bounce than those of Agile-Boltztran or Vertex-1, 2, but track the Agile-Boltztran luminosities very closely thereafter, except for a slight decline at 170 ms when the shock reaches the large density decrement in the progenitor, which Agile-Boltztran fails to adequately resolve. The Vertex-1 luminosities are above those of the other simulations, reflecting its more rapid shock retraction, while the Vertex-2 ${\nu }_{e}$-luminosities track Agile-Boltztran closely until about 170 ms at which point they, like the Chimera luminosities, fall below. The ${\bar{\nu }}_{e}$-luminosities of Chimera-B and Chimera-C fall between those of Agile-Boltztran and Vertex-1 after bounce, rising about $5\times {10}^{51}$ erg s−1 above those of Agile-Boltztran after bounce and falling below Agile-Boltztran after 150 ms. The Vertex-2 ${\bar{\nu }}_{e}$-luminosities track those of Agile-Boltztran until about 150 ms at which point they too fall below. The ${\nu }_{\mu \tau }$-, ${\bar{\nu }}_{\mu \tau }$-luminosities of Chimera-B and Chimera-C also fall between those of Agile-Boltztran and Vertex-1 after bounce, rising about $5\times {10}^{51}$ erg s−1 above Agile-Boltztran until about 170 ms. Vertex-2 tracks Agile-Boltztran closely until about 170 ms at which point it falls below.

The rms energies exhibit a similar pattern after bounce. The ${\nu }_{e}$- and ${\bar{\nu }}_{e}$-rms energies of both Chimera-B and Chimera-C are between those of Agile-Boltztran and Vertex-1 before bounce, tracking closely those of Agile-Boltztran after bounce and rising slightly above Agile-Boltztran after 190 ms. The ${\nu }_{e}$- and ${\bar{\nu }}_{e}$-rms energies of Vertex-1 veer above those of Agile-Boltztran after bounce while those of Vertex-2 track Agile-Boltztran closely after bounce veering slightly below Agile-Boltztran after 50–100 ms. The ${\nu }_{\mu \tau }$-, ${\bar{\nu }}_{\mu \tau }$-rms energies of Chimera-B, Chimera-C, and Vertex-1 track each other within 1 MeV, all three of which lie about 2 MeV higher than those of Agile-Boltztran, while those of Vertex-2 and Agile-Boltztran lie within 1 MeV of each other.

9.6. Comparisons at 100 ms after Bounce

At 100 ms post-bounce (Figure 42), the Chimera-B and Chimera-C shocks are, respectively, about 20 km and 8 km farther out than that of Agile-Boltztran, the front of the Chimera-C shock being close to that of Agile-Boltztran. The Vertex-2 shock radius is close to that of Agile-Boltztran while the Vertex-1 shock has retreated about 10 km inward. Aside from the positions of the shocks, the densities and velocities as a function of radius, shown in Figure 42(a), agree with each other quite closely. The entropy profiles, shown in Figure 42, also agree with each other modulo the position of the shock. Those of Chimera-B and Chimera-C exhibit entropy wiggles behind the shock. These are not due to the computation of the effective index, ${{\rm{\Gamma }}}_{{\rm{e}}}$ (Buras et al. 2006b), but the use of the Colella & Woodward (1984) suggested parameters for supplying dissipation in the vicinity of strong shocks. These wiggles disappear with a somewhat more aggressive parameter choice for shock dissipation, and this choice is now used in current versions of Chimera. The electron fraction profiles are quite similar up to 30 km and beyond 150 km but are displaced horizontally relative to each other in between these distances. This displacement is also reflected in the entropy profiles but is less obvious, as the plots themselves are more horizontal. The origin of these differences is unclear but may reflect the slight differences in the shock trajectories as a function of time, and the somewhat higher infall velocities in the case of Vertex-1.

Figure 42.

Figure 42. Snapshots of model G15 at 100 ms after bounce for (a) velocity and density, and (b) entropy and ${Y}_{{\rm{e}}}$ for the Agile-Boltztran (black), Vertex-1 (green), Vertex-2 (turquoise), and Chimera-B (blue), Chimera-C (red) simulations.

Standard image High-resolution image

Figure 43 compares the ${\nu }_{e}$, ${\bar{\nu }}_{e}$, and ${\nu }_{x}$ luminosity and rms neutrino energy profiles at 100 ms post-bounce. Compared with the Agile-Boltztran luminosities, Chimera-B and Vertex-1 tend to overestimate these quantities in the region behind the shock, while the agreement between the results from Chimera-C, Vertex-2, and Agile-Boltztran for these quantities is extremely good. Chimera-B and Chimera-C tend to compute larger shock jumps than does Agile-Boltztran, while the shock jumps computed by Vertex-2 and Agile-Boltztran are in good agreement. Both Vertex-1 and Chimera-B compute rms energies larger than those of Agile-Boltztran. The ${\nu }_{e}$ rms energies computed by Chimera-C are quite close to those of Agile-Boltztran, while the ${\nu }_{x}$ rms energies computed by Vertex-2 and Agile-Boltztran are in excellent agreement. In other cases, the rms neutrino energies computed by Chimera-C, Vertex-2, and Agile-Boltztran are typically within an MeV of each other.

Figure 43.

Figure 43. Snapshot of model G15 neutrino properties at 100 ms after bounce, for the Agile-Boltztran (black), Vertex-1 (green), Vertex-2 (turquoise), Chimera-B (blue), and Chimera-C (red) simulations. Neutrino luminosities are plotted on the left in Panels (a), (c), and (e) and rms-energies on the right in Panels (b), (d), and (f), with ${\nu }_{e}$ in Panels (a) and (b), ${\bar{\nu }}_{e}$ in Panels (c) and (d), and ${\nu }_{x}$ (or ${\nu }_{\mu \tau }$) in Panels (e) and (f).

Standard image High-resolution image

10. Summary

This report has documented the development and construction of the Chimera code through its C-series implementations. Chimera has been designed to simulate CCSNe throughout the entire neutrino-driven explosion phase, with outputs that can be used to extract important associated observables, such as element synthesis and dispersal, neutrino signatures, and gravitational radiation. The code couples a multidimensional, PPM-plus-remap, Newtonian hydrodynamics module with radial-ray-plus, multi-group flux-limited diffusion neutrino transport and a multi-species nuclear reaction network. The transport module stems from a general relativistic treatment and currently retains the most important element of general relativity, namely the lapse function, which ensures proper redshifting of neutrinos as they propagate out of the gravitational well. General relativity enters the computation of self gravity through a monopole correction to the Newtonian gravitational potential. Chimera evolves all six neutrino and antineutrino distributions and includes an extensive suite of neutrino weak interactions, including sophisticated treatments of nuclear electron capture and non-isoenergetic scattering on nucleons, as well as industry-standard equations of state.

We have subjected Chimera to a suite of test problems and made comparisons with two other sophisticated neutrino radiation-hydrodynamics codes: (1) Agile-Boltztran, a code based on general relativistic gravity, hydrodynamics, and six-species Boltzmann neutrino kinetics, and (2) Prometheus-Vertex, a code similar in many respects to Chimera but based on a more sophisticated, two-moment neutrino transport scheme closed by a variable Eddington factor computed from approximate Boltzmann kinetics.

Development of Chimera continues to include new features and other enhancements beyond those described herein and will be reported, as appropriately, with the associated results.

The authors would like to thank Bernard Mller for the Vertex-1, 2 data, and Mathias Liebendrfer for posting the Agile-Boltztran data online. We also would like to express our appreciation to an anonymous referee for a careful reading of this paper and for suggestions of some tests that were not included in the original draft of this paper. This work has been supported by a wide range of programs in the 15 years since Chimera's inception. Chimera development has been supported by the NSF PetaApps program (OCI-0749242, OCI-0749248, and OCI-0749204), by the NSF Gravitational Physics Theory program (PHY-0555644, PHY-0652874, PHY-0855315, PHY-1505933, and PHY-1806692), by the NSF Nuclear Theory Program (PHY-0244783 and PHY-1516197), and by the NSF Stellar Astronomy and Astrophysics program (AST-0653376). NASA supported CHIMERA development under the Astrophysics Theory Program (NNH08AH71I, NNH11AQ72I, and 07-ATFP07-0011). Chimera development was also supported by the U.S. Department of Energy Offices of Nuclear Physics and Advanced Scientific Computing Research.

This research used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725, and the National Energy Research Scientific Computing Center, which is supported by the U.S. DOE Office of Science under Contract No. DE-AC02-05CH11231. This research also utilized resources of the NSF TeraGrid provided by the National Institute for Computational Sciences under grant number TG-MCA08X010.

P.M. is supported by the National Science Foundation through its employee IR/D program. The opinions and conclusions expressed herein are those of the authors and do not represent the National Science Foundation.

We wish to thank John C. Hayes for his assistance early in the development of Chimera as well as our many colleagues in the CCSN simulation community, and other of our colleagues, from whom we have benefited greatly through discussions related to supernova simulation.

Footnotes

Please wait… references are loading.
10.3847/1538-4365/ab7aff