# CASTRO: A New Compressible Astrophysical Solver. III. Multigroup Radiation Hydrodynamics

###### Abstract

We present a formulation for multigroup radiation hydrodynamics that is correct to order using the comoving-frame approach and the flux-limited diffusion approximation. We describe a numerical algorithm for solving the system, implemented in the compressible astrophysics code, CASTRO. CASTRO uses an Eulerian grid with block-structured adaptive mesh refinement based on a nested hierarchy of logically-rectangular variable-sized grids with simultaneous refinement in both space and time. In our multigroup radiation solver, the system is split into three parts, one part that couples the radiation and fluid in a hyperbolic subsystem, another part that advects the radiation in frequency space, and a parabolic part that evolves radiation diffusion and source-sink terms. The hyperbolic subsystem and the frequency space advection are solved explicitly with high-order Godunov schemes, whereas the parabolic part is solved implicitly with a first-order backward Euler method. Our multigroup radiation solver works for both neutrino and photon radiation.

###### Subject headings:

diffusion – hydrodynamics – methods: numerical – radiative transfer## 1. Introduction

Numerical simulations are a useful tool for many radiation hydrodynamic problems in astrophysics. However, numerical modeling of radiation hydrodynamic phenomena is very challenging for a number of reasons. To achieve a highly accurate description of the system usually requires a multidimensional treatment with high spatial resolution. It is also desirable to handle the radiation frequency variable properly, because radiation transport is a frequency-dependent process. Moreover, the numerical algorithms for radiation hydrodynamics must be stable and efficient.

Core-collapse supernovae are such complex phenomena (Bethe, 1990; Kotake et al., 2006; Janka et al., 2007; Janka, 2012) that neutrino radiation hydrodynamics simulations have been indispensable for helping understand the explosion mechanism. Various algorithms have been developed to solve neutrino radiation hydrodynamics equations. Because of the complexities involved in the problem and the limited computer resources available, various tradeoffs have to be made. Many of these algorithms are 1D (e.g., Bowers & Wilson, 1982; Bruenn, 1985; Mezzacappa & Bruenn, 1993; Burrows et al., 2000; Liebendörfer et al., 2004; O’Connor & Ott, 2010), or have used the so-called “ray-by-ray” approach, which does not truly perform multi-dimensional transport (e.g. Burrows et al., 1995; Rampp & Janka, 2002; Buras et al., 2006; Takiwaki et al., 2012). The two-dimensional code of Miller et al. (1993) is a gray flux-limited diffusion (FLD) code that does not include order terms, where is the characteristic velocity of the system. The 2D multigroup multiangle code of Livne et al. (2004) does not include all order terms, and in particular it lacks direct coupling among radiation groups. The two-dimensional algorithm of Hubeny & Burrows (2007) is based on a mixed-frame formulation of a multigroup two-moment system that is correct to order , but it remains to be implemented in a working code. Swesty & Myra (2009) have developed a fully coupled 2D multigroup FLD code based on a comoving frame formulation that is correct to order . Besides the grid-based methods, smoothed particle hydrodynamics (SPH) methods have also been applied to multi-dimensional neutrino radiation hydrodynamics simulations of core-collapse supernovae (e.g., Herant et al., 1994; Fryer & Warren, 2002; Fryer et al., 2006), but the SPH scheme has not yet been extended to multigroup. Recently Abdikamalov et al. (2012) developed a Monte Carlo method for neutrino transport, but hydrodynamics and radiation transport are not yet coupled in the scheme.

In this paper, we present a numerical algorithm for solving multi-dimensional multigroup photon and neutrino radiation hydrodynamics equations in the comoving frames that are correct to order . Radiation quantities that are of interest are the radiation energy density, radiation flux, and radiation pressure tensor. It is more accurate to evolve both the radiation energy density and the radiation flux especially when the radiation is highly anisotropic and optically thin. However, the two-moment approach is very expensive in terms of computer time and memory for multi-dimensional multigroup radiation hydrodynamics simulations. Thus, we have adopted the flux-limited diffusion approach (Alme & Wilson, 1973) for computational efficiency. In the FLD approach, only the radiation energy density needs to be evolved over time, and the radiation flux is derived through a diffusion approximation. Furthermore, we use an analytic closure for the radiation pressure tensor. We have adopted a multigroup method, in which the radiation frequency is discretized into multiple groups. The equations we solve are similar to those in Swesty & Myra (2009). However, we have neglected inelastic scattering of neutrinos on matter, which is currently treated as elastic scattering. Since the total cross section for the latter at the neutrino energies near and outside the neutrinospheres are –100 times smaller than the cross sections for scattering off nucleons and nuclei, not including inelastic effects at this stage has a minor effect on the qualitative behavior of collapse, bounce, and shock dynamics. Nevertheless, the approach of Thompson et al. (2003) to handle inelasticity and energy redistribution can be incorporated into our scheme. Another difference between our scheme and the scheme of Swesty & Myra (2009) is that we have developed a Godunov scheme in which radiation is fully coupled into a characteristic-based Riemann solver for the hyperbolic part of the system.

Many multigroup neutrino transport codes use a fixed Eulerian grid. However, it is difficult to achieve sufficient resolution in multi-dimensional simulations with this simple approach. To efficiently achieve high resolution, we use an Eulerian grid with block-structured adaptive mesh refinement (AMR). The AMR algorithm we use has a nested hierarchy of logically-rectangular variable-sized grids and refines simultaneously in both space and time. AMR techniques have been successfully used in multi-dimensional multigroup photon radiation transport (Gittings et al., 2008; van der Holst et al., 2011). However, to the best of our knowledge, this is the first time that AMR has been employed in a multi-dimensional multigroup neutrino solver. Besides AMR, the Arbitrary Lagrange Eulerian (ALE) method is an alternative approach for mesh refinement. We note that ALE has been successfully employed in multi-dimensional multigroup neutrino transport (Livne et al., 2004; Burrows et al., 2007; Ott et al., 2008) and photon transport (Marinak et al., 2001).

We have implemented our algorithm in the compressible astrophysics code, CASTRO. This is the third paper in a series on the CASTRO code and its numerical algorithms. In our previous papers, we describe our treatment of hydrodynamics, including gravity and nuclear reactions (Almgren et al., 2010, henceforth Paper I), and our algorithm for flux-limited gray radiation hydrodynamics based on a mixed-frame formulation (Zhang et al., 2011, henceforth Paper II). Here, we describe our algorithm for flux-limited multigroup radiation hydrodynamics based on a comoving-frame formulation.

In § 2, we present the multigroup radiation hydrodynamics equations that CASTRO solves. We show that the system can be divided into a hyperbolic subsystem (§ 2.3), a subsystem of advection in frequency space (§ 2.4), and a parabolic subsystem for radiation diffusion (§ 2.5). Analytic results for the mathematical characteristics of the hyperbolic subsystem are also presented in § 2.3. The hyperbolic and the frequency space advection steps are solved by an explicit solver, whereas the parabolic subsystems is solved by an iterative implicit solver. We describe the explicit solvers in § 3.1, followed by a description of the implicit solver in § 3.2. We also discuss acceleration schemes that can greatly improve the convergence rate of the iterative implicit solver (§ 3.2.4). In § 4, we present results from a series of test problems. A scaling test is shown in § 5 to demonstrate the scaling behavior of the multigroup group radiation solver. Finally, the results of the paper are summarized in § 6.

## 2. Multigroup Radiation Hydrodynamics

In this section, we present the multigroup radiation hydrodynamics equations that CASTRO solves. Our formulation is based on a comoving-frame approach, unlike our gray radiation hydrodynamics solver presented in Paper II. In the comoving-frame approach, the radiation quantities and the opacities are measured in the comoving frames. The set of equations are correct to order . We then derive the multigroup radiation hydrodynamics equations using frequency space discretization. We also analyze the mathematical characteristics of the system. Here we present the results in terms of neutrino radiation and indicate modifications needed for photon radiation.

### 2.1. Frequency-Dependent Radiation Hydrodynamics Equations

The system of neutrino radiation hydrodynamics that is correct to order can be described by

(1) | ||||

(2) | ||||

(3) | ||||

(4) | ||||

(5) |

(Buchler, 1983; Mihalas & Mihalas, 1999; Castor, 2004; Swesty & Myra, 2009). Here , , , and are the mass density, velocity, pressure, and total matter energy per unit mass (internal energy, , plus kinetic energy, ), respectively. is the electron fraction defined to be the net electron number (electrons minus positrons) per baryon. is the external force on the fluid (e.g., gravitational force). , , and are the monochromatic radiation energy density, radiation flux, and radiation pressure tensor at frequency , measured in the comoving fluid frame, respectively. The speed of light is denoted by . The absorption coefficient is denoted by , and the total interaction coefficient including conservative scattering is given by . The emission coefficient is denoted by . For photon radiation is zero, whereas for neutrinos is given by

(6) |

where is the baryon mass, is the Planck constant, and for neutrinos, for neutrinos, and for other flavor neutrinos. Also note that for neutrino radiation with multiple species, the integrals over frequency imply summation over the species. The colon product in indicates contraction over two indices as follows,

(7) |

The system of the radiation hydrodynamics (Equations 1–5) is closed by an equation of state (EOS) for matter, the diffusion approximation for radiation flux, and a closure relation for radiation that computes the radiation pressure tensor from radiation energy density and radiation flux. In the FLD approximation (Alme & Wilson, 1973), the monochromatic radiation flux is written in the form of Fick’s law of diffusion,

(8) |

where is the flux limiter at frequency . We have implemented a variety of choices for the flux limiter (e.g., Minerbo, 1978; Bruenn et al., 1978; Levermore & Pomraning, 1981). For example, one form of the flux limiter in Levermore & Pomraning (1981) is

(9) |

where

(10) |

In the optically thick limit, , , and the flux approaches the classical Eddington approximation . In the optically thin limit, , , and as expected for streaming radiation. For the radiation pressure tensor, we use an approximate closure (Minerbo, 1978; Levermore, 1984)

(11) |

where

(12) |

is the scalar Eddington factor, is the identity tensor of rank 2, and . There are other choices for the scalar Eddington factor. For example, Kershaw (1976) has suggested

(13) |

Note that in the optically thick limit, , and in the optically thin limit, .

Using the FLD approximation and the approximate closure, the equations of the radiation hydrodynamics now become

(14) | ||||

(15) | ||||

(16) | ||||

(17) | ||||

(18) |

### 2.2. Multigroup Radiation Hydrodynamics Equations

To numerically solve Equations (14)–(18), we first discretize the system in frequency space by dividing the radiation into energy groups. We define the radiation energy density in each group as

(19) |

where and are the frequency at the lower and upper boundaries, respectively, for the -th group. For clarity, is denoted by hereafter.

In the multigroup methods, various group mean coefficients are introduced to the system. For example, the Rosseland mean interaction coefficient is defined as

(20) |

where . In this paper, we will simply assume that

(21) | ||||

(22) |

where is the representative frequency for the -th group. Our simple treatment of group mean coefficients is sufficient for the transport of continuum radiation. In particular, it is sufficient for neutrino transport in the core-collapse supernova problem because the monochromatic opacities and emissivities of neutrinos are smooth functions of neutrino frequency. It should, however, be emphasized that our algorithm for radiation hydrodynamics does not rely on the simple treatment of group mean coefficients used in this paper. We also assume that the flux limiter is approximately independent of frequency within each group, and denote the flux limiter and the scalar Eddington factor for the -th group as and , respectively. As for the emissivity, we define

(23) |

where is the group width. However, in the case of photon radiation, we have adopted the polylogarithm function based method of Clark (1987) for computing the group-integrated Planck function. Thus, assuming local thermodynamic equilibrium, we have for photon radiation

(24) |

where .

Using these definitions and approximations, we obtain the multigroup radiation hydrodynamics equations

(25) | ||||

(26) | ||||

(27) | ||||

(28) | ||||

(29) | ||||

where , , denotes .

### 2.3. Mathematical Characteristics of the Hyperbolic Subsystem

The system of multigroup radiation hydrodynamics (Equations 25–29), similar to the system of gray radiation hydrodynamics we have studied in Paper II, has a hyperbolic subsystem,

(30) | ||||

(31) | ||||

(32) | ||||

(33) | ||||

(34) |

where . Since is just an auxiliary variable that does not alter the characteristic wave speeds of the hyperbolic subsystem, we will neglect it in the remainder of this subsection for simplicity. Gravitational force terms are source terms, and thus they are not included in the following analysis. The term in Equation (29) is not included in the analysis of the hyperbolic subsystem here even though it becomes comparable to the term in Equation (34) in the streaming limit where . The first reason for not including the term is that the main purpose of including radiation in the hyperbolic subsystem is to improve the accuracy in the optically thick limit where the radiation can affect the hydrodynamics significantly. The term is not significant in the optically thick limit where . The second reason is that in the streaming limit, where , radiation force has little impact on the motion of matter. Thus, radiation is essentially decoupled from the hydrodynamics of matter in the streaming limit. Moreover, in the streaming limit the diffusion term (Equation 29) will dominate the term and all the terms in Equation (34). Therefore, the term is insignificant in both the diffusion limit and the streaming limit.

In CASTRO, we solve the hyperbolic subsystem with a Godunov method, which utilizes a characteristic-based Riemann solver (see § 3.1 for details). The Godunov method requires that we analyze the mathematical characteristics of the hyperbolic subsystem. For simplicity, let us consider the system in one dimension, which can be written in terms of primitive variables as,

(35) |

where the primitive variables are

(36) |

and the Jacobian matrix is

(37) |

where is the adiabatic index of the matter. This system is hyperbolic because the Jacobian matrix is diagonalizable with real eigenvalues,

(38) |

Here

(39) |

is the radiation modified sound speed, where is the sound speed without radiation and is the contribution to the sound speed from each group. The corresponding right eigenvectors are,

(40) |

and the corresponding left eigenvectors are,

(41) |

where

(42) |

These eigenvectors define the characteristic fields for the one-dimensional system. By computing the product of the right eigenvectors and the gradients of their corresponding eigenvalues (LeVeque, 2002), we find that the first and second fields are genuinely nonlinear corresponding to either a shock wave or a rarefaction wave. The rest of the fields are linearly degenerate corresponding to a contact discontinuity in either density or radiation energy density of one group moving at a speed of . Note that there is also a jump in fluid pressure such that the total pressure, , is constant across the contact discontinuity. Obviously, in three dimensions, there are two additional contact discontinuities for transverse velocities just like the case of pure hydrodynamics.

### 2.4. Advection in Frequency Space

The integrals in Equation (29) form the following system,

(43) |

which can also be written in the conservation law form

(44) |

where and . Here we have used Equation (19). Equation (44) describes the advection of in logarithmic frequency space with an advection speed of . Because is the integrated radiation energy, the system conserves energy.

Various terms in Equation (29) have been split between the hyperbolic subsystem (§ 2.3) and the system for the advection in the frequency space (§ 2.4). There is however another way to look at Equation (29). The evolution of the radiation energy density (Equation 29) without the diffusion term and the source-sink term can also be split as

(45) | ||||

(46) |

Equation (46) also represents an advection in logarithmic frequency space and can be written in the form

(47) |

where . Here the conserved quantity is the number of radiation particles (i.e., ). This new way of splitting will be further discussed later in Section 3.1.

### 2.5. Multigroup Radiation Diffusion

The parabolic part of the system consists of the radiation diffusion and source-sink terms, which were omitted from the discussion of the hyperbolic subsystem (§ 2.3) and advection in frequency space (§ 2.4),

(48) | ||||

(49) | ||||

(50) |

where is the specific internal energy of the fluid, and . The term represents the energy exchange in the comoving frame between the matter and the -th radiation group through absorption and emission of radiation. Note that the terms on the right-hand side (RHS) of Equations (48), (49), and (50) all contain the speed of light . Thus these are order terms in . An implicit treatment for these equations is usually necessary because of their stiffness. For photon radiation, Equation (49) is not needed.

## 3. Numerical Methods

CASTRO uses a nested hierarchy of logically-rectangular, variable-sized grids with simultaneous refinement in both space and time. The AMR algorithms for hyperbolic and radiation diffusion equations in CASTRO have been described in detail in Papers I and II, respectively. The multigroup solver does not require any changes in the AMR algorithms. Thus, we will not describe them again in this paper.

In this section, we describe the single-level integration scheme. For each step at a single level of refinement, the state is first evolved using an explicit method for the hyperbolic subsystem (Equations 30–34) and advection in frequency space (§ 3.1). Then an implicit update for radiation diffusion and source-sink terms is performed (§ 3.2). In § 3.2.4, we describe two acceleration schemes for speeding up the convergence rate of the implicit solver for multigroup radiation diffusion equations.

### 3.1. Explicit Solver for Hyperbolic subsystem and Advection in Frequency Space

The hyperbolic subsystem (§ 2.3) is treated explicitly. Our Godunov algorithm for the hyperbolic subsystem of the multigroup radiation hydrodynamics is an extension of the hyperbolic solver for the gray radiation hydrodynamics presented in Paper II , which in turn is based on the hydrodynamics algorithm presented in Paper I of this series. We refer the reader to Papers I and II for a detailed description of the integration scheme, which supports a general equation of state, self-gravity, nuclear reactions, the advection of nuclear species and electron fraction. Here we will only present a brief overview of the Godunov scheme.

The time evolution of the hyperbolic subsystem can be written in the form

(51) |

where denotes the conserved variables, and is their flux. The conserved variables are defined at cell centers. We predict the primitive variables, including , , , , , where , from cell centers at time to edges at time and use an approximate Riemann solver to construct fluxes, , on cell faces. This algorithm is formally second order in both space and time. The time step is computed using the standard Courant-Friedrichs-Lewy (CFL) condition for explicit methods, and the sound speed used in the computation is now the radiation modified sound speed (Equation 39). The time step is also limited by maximally allowed changes in temperature and electron fraction in the implicit solver. Additional constraints for the time step are applied if additional physics, such as burning, is included. CASTRO solves the hyperbolic subsystem of radiation hydrodynamics with an unsplit piecewise parabolic method (PPM) with characteristic tracing and full corner coupling (Miller & Colella, 2002). An approximate Riemann solver based on the Riemann solver of Bell et al. (1989) and Colella et al. (1997), that utilizes the eigenvectors of the system, is used to compute the Godunov states and fluxes at the cell interface.

Advection in frequency space (§ 2.4) can also be treated explicitly because its time scale is comparable to that of the hyperbolic subsystem. One approach is to adopt an operator-splitting method and to evolve Equation (44) with a Godunov method. This is the approach taken by van der Holst et al. (2011) in the CRASH code for photon radiation. The term can be treated as a source term of Equation (44). The success of this approach for a number of test problems has been demonstrated by van der Holst et al. (2011). We have also found that this approach can produce satisfactory results for all the photon test problems in Section 4. However, we have found that erroneous results were obtained for core-collapse supernova simulations (§ 4.2), probably due to operator-splitting errors. One difficulty in core-collapse supernova simulations is that the divergence of velocity can be significant because there is a strong shock and the typical velocity in front of the shock is –. In the case of core-collapse supernovae, the negative divergence of velocity will shift neutrinos towards high frequency (Equation 44). The stronger conservative scattering at higher frequency can further amplify the neutrino energy density due to the work term (Equation 34). These two processes can combine to cause unphysical results due to splitting errors. However, the origin of the various terms involved in these two processes is the term in Equation (5). To avoid operator-splitting errors, it is therefore desirable to treat the term without splitting it. On the other hand, the term can also make significant contribution to the hyperbolic subsystem in the optically thick limit, and we do not wish to sacrifice accuracy of the hyperbolic subsystem by taking the traditional approach of leaving out radiation from the Riemann solver for the hydrodynamics.

To avoid the operator-splitting without sacrificing the accuracy of the hyperbolic subsystem, we have adopted a new approach based on splitting the radiation energy density equation into Equations (45) and (46). Note that the radiation diffusion and source-sink terms will be treated separately by an implicit solver (§ 3.2). Our new approach still utilizes the Godunov scheme for the hyperbolic subsystem governed by Equations (30)–(34). After the Godunov states and fluxes at the cell interfaces are obtained, all variables in the hyperbolic system except the radiation energy density are updated as usual. For the radiation energy density, two steps are taken. The radiation energy density is first updated according to Equation (45) with its RHS computed using the Godunov states. Then another Godunov solver is used to evolve Equation (47), which is equivalent to Equation (46). The conversion between and is performed according to . The Godunov states at time of the first Godunov solver are used in computing and to obtain the wave speed in the frequency space advection. The second Godunov scheme for the frequency space advection is a standard approach based on the method of lines. For time integration we use the third-order total variation diminish Runge-Kutta scheme of Shu & Osher (1988). We use the piecewise linear method with a monotonized central slope limiter (van Leer, 1977) to reconstruct and in order to achieve high order in (logarithmic frequency) space. The reconstructed variables are used to compute the HLLE flux (Harten et al., 1983; Eisenstat, 1981). The frequency space advection has its own CFL condition for stability. Since CASTRO uses the hydrodynamic CFL condition with other additional constraints for the time steps, if necessary, subcycling in time is employed in order to satisfy the frequency space CFL condition.

### 3.2. Implicit Solver for Radiation Diffusion and Source-Sink Terms

The implicit solver evolves the radiation and matter according to equations (48), (49), and (50) with the results of the explicit update (§ 3.1) as the initial conditions. Our scheme is based on a first-order backward Euler method.

#### 3.2.1 Outer Newton Iteration

We solve the parabolic subsystem (Equations 48, 49, and 50) iteratively via Newton’s method. We define

(52) | ||||

(53) | ||||

(54) |

where the superscript denotes the state following the explicit update, and . The desired solution is for , and to all be zero. We approach the solution by Newton iteration, and in each iteration step we solve the following system

(55) |

Here , and , where the superscript denotes the stage of the Newton iteration. The exact Jacobian matrix in equation (55) is very complicated. For simplicity, we neglect the derivatives of the diffusion coefficient . This is a common practice in radiation transport, and does not appear to significantly degrade the convergence rate. With the approximate Jacobian matrix, we can eliminate and from equation (55) and obtain the following equation for ,

(56) |

To reduce clutter, we have dropped the superscript for , , , and . The newly introduced variables in equation (56) are given by

(57) | ||||

(58) |

where

(59) | ||||

(60) | ||||

(61) | ||||

(62) | ||||

(63) |

Here all the derivatives are computed from the results of the previous Newton iteration (i.e., the state).

#### 3.2.2 Inner Iteration

We now present our algorithm for solving equation (56), which is actually a set of equations coupled through the summation terms on the RHS of equation (56). The system is usually too large to be solved directly. Instead, we solve it by decoupling the groups and utilizing an iterative method. Note that the iteration here for solving equation (56) is different from the Newton iteration for solving the entire parabolic subsystem. The inner iteration here is embedded inside the outer Newton iteration step. In our algorithm, we avoid the coupling by replacing the radiation energy density in the coupling terms with its state from the previous inner iteration. Thus, each step of the inner iteration solves

(64) |

where is the inner iteration index. Here, we have dropped the superscript for and to reduce clutter. For the first inner iteration, the state of the previous inner iteration is set to the result of the previous outer iteration. The inner iteration for equation (64) is stopped when the maximum of on the computational domain is below a preset tolerance (e.g., ).

CASTRO can solve equations in the following canonical form (Paper II)

(65) |

where are cell-centered coefficients and , , and are centered at cell faces. We use the same approach for spatial discretization at the AMR coarse-fine interface as in Almgren et al. (1998). For equation (64) the and coefficients are not used. The same canonical form works for Cartesian, cylindrical, and spherical coordinates so long as appropriate metric factors are included in the coefficients and the RHS. CASTRO uses the hypre library (Falgout & Yang, 2002; hypre Code Project, 2012) to solve the linear system.