Analytic Closures for M1 Neutrino Transport

TL;DR

This paper evaluates various M1 neutrino transport closures by comparing them with Monte Carlo simulations around protoneutron stars, finding that no single closure is best in all cases, but Minerbo's maximum entropy closure performs relatively well overall.

Contribution

It provides an extensive comparison of M1 closure relations against Monte Carlo benchmarks, highlighting their varying accuracy depending on conditions.

Findings

No closure is universally best across all scenarios.

Minerbo's maximum entropy closure performs relatively well overall.

Closure accuracy depends on matter configuration, neutrino type, and energy.

Abstract

Carefully accounting for neutrino transport is an essential component of many astrophysical studies. Solving the full transport equation is too expensive for most realistic applications, especially those involving multiple spatial dimensions. For such cases, resorting to approximations is often the only viable option for obtaining solutions. One such approximation, which recently became popular, is the M1 method. It utilizes the system of the lowest two moments of the transport equation and closes the system with an ad hoc closure relation. The accuracy of the M1 solution depends on the quality of the closure. Several closures have been proposed in the literature and have been used in various studies. We carry out an extensive study of these closures by comparing the results of M1 calculations with precise Monte Carlo calculations of the radiation field around spherically-symmetric protoneutron star models. We find that no closure performs consistently better or worse than others in all cases. The level of accuracy a given closure yields depends on the matter configuration, neutrino type, and neutrino energy. Given this limitation, the maximum entropy closure by Minerbo (1978) on average yields relatively accurate results in the broadest set of cases considered in this work.

Figures

Tables

Table 1. Table 1: Mean square deviation of the flux and the Eddington factors obtained with closure prescriptions from the analytical solution for the radiative uniform sphere problem. The sum in the formula for the normalized mean square deviation ( 27 ) is taken over radii from r min = 1.0 subscript 𝑟 min 1.0 r_{\rm min}=1.0 to r max = 2.0 subscript 𝑟 max 2.0 r_{\rm max}=2.0 .

Closure prescription	$\delta f_{\nu }\left(r\right)$	$\delta p_{\nu }\left(r\right)$
Kershaw	0.13	0.32
Wilson	0.05	0.14
Levermore	0.10	0.22
ME	0.07	0.17
MEFD	0.07	0.17
Janka_1	0.07	0.13
Janka_2	0.10	0.21
Fit Closure	0.01	0.01

Equations

$$\frac{d{x}^{\alpha }}{d\tau }\frac{\partial \mathcal{F}}{\partial x^{\alpha }}+\frac{d{p}^i}{d\tau }\frac{\partial \mathcal{F}}{\partial p^i}=(-p^{\alpha }{u}_{\alpha })S(p^{\mu },x^{\mu },\mathcal{F}).$$

$$E_{\nu }=\int \epsilon \mathcal{F}(p^{\mu },x^{\mu })\delta (h\nu -\epsilon )d^{3}p,$$

$$F_{\nu }^j=\int p^j\mathcal{F}(p^{\mu },x^{\mu })\delta (h\nu -\epsilon )d^{3}p,$$

$$P_{\nu }^{ij}=\int p^i{p}^j\mathcal{F}(p^{\mu },x^{\mu })\delta (h\nu -\epsilon )\frac{d^{3}p}{\epsilon }.$$

$$\frac{\partial }{\partial t}{E}_{\nu }+\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2{F}_{\nu }^r\right)=S[0{]}_{\nu }$$

$$\frac{\partial }{\partial t}{F}_{\nu }^r+\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2{P}_{\nu }^{rr}\right)=S[1{]}_{\nu }^r,$$

$$\mathrm{Function}(M[0]...M[k+1])=S[k],$$

$$M[k{]}_{{\epsilon }_0}^{{\alpha }_1...{\alpha }_k}=\int {\epsilon }^2\mathcal{F}(p^{\mu },x^{\mu })\delta (h\nu -\epsilon )\frac{p^{{\alpha }_1}}{\epsilon }...\frac{p^{{\alpha }_k}}{\epsilon }\frac{d^{3}p}{\epsilon },$$

$$\displaystyle f(x)=\sum\limits_{N}\frac{1}{N!}\frac{d^{N}f(x)}{dx^{N}}\bigg{|}_{x=x_{0}}(x-x_{0})^{N}.$$

$$P_{\nu }^{ij}=\frac{3p-1}{2}{P}_{\mathrm{thin}}^{ij}+\frac{3(1-p)}{2}{P}_{\mathrm{thick}}^{ij},$$

$$P_{\mathrm{thick}}^{ij}=\frac{1}{3}{E}_{\nu }{\delta }^{ij}$$

$$P_{\mathrm{thin}}^{nn}=E_{\nu }\frac{F_{\nu }^n{F}_{\nu }^n}{|F_{\nu }{|}^2},\mathrm{and}{P}_{\mathrm{thin}}^{ij}=0,\mathrm{if}i\mathrm{or}j\ne n.$$

$$P_{\nu }^{ij}=E{f}^i{f}^j+\frac{E}{3}{\delta }^{ij}\left(1-f^2\right),$$

$$f^i=F_{\nu }^i/E_{\nu }$$

$$p=\frac{1+2f^2}{3}.$$

$$p=\frac{1}{3}-\frac{1}{3}f+f^2.$$

$$p=\frac{3+4f^2}{5+2\sqrt{4-3f^2}}.$$

$$S[\mathcal{F}(\mu )]=(1-\mathcal{F})\log (1-\mathcal{F})+\mathcal{F}\log \mathcal{F},$$

$$e$$

$$f$$

$$\mathcal{F}=\frac{1}{e^{\eta -a\mu }+1},$$

$$p=\frac{1}{3}+\frac{2}{3}(1-e)(1-2e)\chi \left(\frac{f}{1-e}\right),$$

$$\chi (x)=x^2(3-x+3x^2)/5,$$

$$p=\frac{1}{3}\left(1-2f+4f^2\right).$$

$$p=\frac{1}{3}+\frac{2f^2}{15}(3-f+3f^2).$$

$$p=\frac{1}{3}\left[1+a{f}^m+(2-a)f^n\right],$$

$$\delta Y(X)=\sqrt{\frac{1}{N_X}\sum _{X_{\min }}^{X_{\max }}{\left[1-\frac{Y(X_i)}{Y_0(X_i)}\right]}^2}.$$

$$\bar{\delta }Y=\frac{\sum w_i\delta Y_i}{\sum w_i},w_i=S_i/S_{\max },$$

$$\mathcal{F}(r,\mu )=B\left[1-e^{-\kappa Rs(r,\mu )}\right],$$

$$\displaystyle s(r,\mu)=\left\{\begin{array}[]{lll}\displaystyle\frac{r}{R}\mu+g(r,\mu)&&\mathrm{if}\ \ r<R,\quad-1\leq\mu\leq 1\\ &\\[-5.0pt] \displaystyle 2g(r,\mu)&&\mathrm{if}\ \ r\geq R,\quad\sqrt{1-\left(\frac{R}{r}\right)^{2}}\leq\mu\leq 1\\ &\\ \displaystyle 0&&\mathrm{otherwise}\end{array}\right.$$

$$g(r,\mu )=\sqrt{1-{\left(\frac{r}{R}\right)}^2(1-{\mu }^2)}.$$

Full text

††pagerange: Analytic Closures for M1 Neutrino Transport–References

1 Introduction

Neutrinos play an important role in core-collapse supernovae (CCSNe), coalescence of binary neutron stars and many other astrophysical phenomena. Their collective behavior is described by the distribution function that obeys the Boltzmann transport equation. The state of the radiation field is characterized by spatial coordinates, the propagation direction (two angles), energy, and time, making the problem seven-dimensional in the most general case. Many astrophysical systems have dense and opaque central regions surrounded by transparent low-density envelopes. Radiation moves within the dense central regions via diffusion and, when it leaks into the outer regions, it streams freely. The transport equation has a parabolic character in the former region, while in the latter region, it has a hyperbolic character (e.g., Mihalas & Mihalas, 1984; Pomraning, 1983). In order to model such systems accurately, the solution techniques must handle not only the two different regimes, but also the transition between the two. In the presence of scattering, the collision terms on the right-hand side of the Boltzmann equation contains the angular moments of the specific intensity, which makes the Boltzmann equation an integro-differential equation. These aspects make solving the transport equation a challenging computational task. For this reason, often one has to resort to approximations and simplifications to make the problem tractable.

One way of simplifying the problem is to assume spherical or axial symmetries to reduce the number of spatial dimensions. There are many situations where such assumptions have been employed. Spherically-symmetric calculations have been performed by, e.g., Mezzacappa & Bruenn (1993a, b, c); Yamada (1997); Rampp & Janka (2000); Liebendörfer et al. (2005); Sumiyoshi et al. (2005). 2D axisymmetric simulations have been performed by, e.g., Livne et al. (2004); Ott et al. (2008); Brandt et al. (2011); Skinner et al. (2016); Burrows et al. (2016). In addition, ray-by-ray approaches, in which multi-dimensional transport problem is approximated as a set of one-dimensional problems along radial rays, has been used widely (e.g., Burrows et al., 1995; Buras et al., 2006; Marek & Janka, 2009; Mueller et al., 2011; Müller et al., 2012; Bruenn et al., 2013, 2016; Lentz et al., 2015; Melson et al., 2015a; Melson et al., 2015b; Müller & Janka, 2015; Müller, 2015). Most realistic problems, however, do not possess spatial symmetries. For these problems, the transport equation has to be solved in full three dimensions. The pioneering attempts to solve three-dimensional Boltzmann equation have been already taken (e.g., Radice et al., 2013; Sumiyoshi & Yamada, 2012; Sumiyoshi et al., 2015), yet the computational cost remains too high for solving it in realistic scenarios.

In order to further reduce the cost, one can approximate the Boltzmann equation either by neglecting the time dependence (steady-state solution) and/or energy dependence (gray schemes). The simplest treatment of the transport problem is the "light bulb" approach, in which simple parametrized neutrino heating and cooling rates are imposed (e.g., Ohnishi et al., 2006; Murphy & Burrows, 2008; Radice et al., 2016). The less cruder approximation, the so-called leakage/heating scheme, has been used extensively in the literature (e.g., Ruffert et al., 1996; Rosswog & Liebendörfer, 2003; O’Connor & Ott, 2011; Ott et al., 2012, 2013a, 2013b; Mösta et al., 2014; Abdikamalov et al., 2015; Deaton et al., 2013). The leakage/heating scheme evaluates the local neutrino energy and number emission rates, which are then locally subtracted from the matter. A fraction of the emitted energy is deposited back as neutrino heating in the “gain” region outside the protoneutron star (PNS)(e.g., O’Connor & Ott, 2010; Ott et al., 2013b).

In this paper we focus on an alternative approach, employing the reduction of the angular degrees of freedom of the problem, called the moment scheme. The simplest version of the moment scheme is the diffusion approximation. One takes the zeroth moment of the transport equation, which yields an equation that contains the zeroth and first moments of the distribution function. For non-static moving media, the third moment is also present (e.g., Just et al., 2015). These three moments represent the energy density, flux, and pressure tensor of radiation, respectively. In the optically thick limit, the first moment can be approximated using the gradient of the zeroth moment via the Fick’s law, while the second moment can be approximated as one-third of the zeroth moment (e.g., Pomraning, 1983). These relations allows us to “close” the zeroth moment of the transport equation. While the resulting equation is far simpler than the original Boltzmann equation, the diffusion approximation is not valid in the free-streaming regime and yields inaccurate results such as acausal flux. This can be mitigated by using the flux limiter (e.g., Smit et al., 2000; Burrows et al., 2000) or by using other advanced prescriptions (e.g., Dgani & Janka, 1992; Scheck et al., 2006; Müller & Janka, 2015). The way to obtain a more accurate solution is to incorporate higher-order moments.

The first moment of the transport equation results in an equation containing up to the second moments of the distribution function. In general relativity and for non-static media, the third moment is also present (Shibata et al., 2011; Cardall et al., 2013; Just et al., 2015). Together with the zeroth moment of the transport equation the system has four sets of unknowns: the zeroth, first, second, and third moments. There are two commonly used approaches for closing the system. In the first method, using the first and second moments as given, one can express the source terms of the Boltzmann equation due to interaction with matter as functions of only space, time, and momentum coordinates (Burrows et al., 2000; Rampp & Janka, 2002). This transforms the Boltzmann equation from a non-linear integro-differential equation into a linear differential equation. The solution of this simpler equations yields higher moments to close the original system of the lowest two moment equations, from which we can obtain updated zeroth and first moments. This procedure is iterated until convergence. Depending on the method for obtaining the closure, the approach can yield the full solution of the Boltzmann equation. This method is usually called the variable Eddington tensor method (Burrows et al., 2000; Rampp & Janka, 2002; Cardall et al., 2013).

Another approach for closing the system is by expressing the second and the third moments in terms of the lower-order moments using approximate analytical relations. This results in a closed system of two equations for the zeroth and first moments. Originally proposed by Pomraning (1969); Kershaw (1976); Levermore (1984), such methods are often called the M1 methods (Smit et al., 1997; Pons et al., 2000; Audit et al., 2002) or the “algebraic Eddington factor” methods (Just et al., 2015).

A common way to derive a closure relation by interpolating between optically thick and optically thin limits. In these limits, the second and third moments can be expressed precisely in terms of the zeroth and first moments111More specifically, in the optically thin limit, one can derive an expression for the second and third moments for a freely propagating radiation beam (e.g., Shibata et al., 2011).. The Eddington factor, which ranges from $1/3$ to $1$ between these limits, serves as an interpolation parameter. The functional form of the Eddington factor in terms of the local energy density and flux of radiation is a relation often called as the M1 closure in the literature. Once this relation is established, the system is closed (e.g., Shibata et al., 2011)222This is not the only approach for closing the system that exists in the literature. In principle, if the second and third moments are assumed to be functions of the local values of the lowest two moments, then the second moment can be expressed in terms of the latter two via the Eddington factor (e.g., Just et al., 2015). Similarly to the second moment, the third moment can be expressed in terms of the lowest two via the ”third-order counterpart” of the Eddington factor. To close the system, this factor must be expressed in terms of the radiation energy density and flux, which has been achieved by, e.g., Vaytet et al. (2011) and Just et al. (2015) for two different closures..

The M1 approach is particularly suitable for problems with not too complex geometries such as CCSNe and remnants neutron star mergers. In these problems, the radiation field is often arranged in such a way that there exists some approximate relationship between higher-order and lower-order moments. That said, not all problems possess such properties. A prominent example is a collision of radiation beams coming from multiple sources. For this problem, the closure relation depends not only on the local values of the first and second moments, but also on the spatial distribution of the radiation sources. In general, if the second and third moments are assumed to be functions of the zeroth and first moments, then the former two must be symmetric with respect to rotation around the direction of the radiation flux (e.g., Cardall et al., 2013). For problems with such symmetries, the M1 approach offers excellent compromise between computational cost and accuracy.

Moreover, the moment equations constitute a hyperbolic system, which allows us to utilize a wide variety numerical methods developed for solving hyperbolic system of conservations laws (e.g., Godunov-type methods) to calculate the solution of the transport problem (Pons et al., 2000). For these reasons and because of their relatively modest computational cost, such methods recently gained significant popularity in astrophysics.

The M1 method has been applied to a wide range of problems such as core-collapse supernovae (e.g., O’Connor & Couch, 2015; Roberts et al., 2016; Kuroda et al., 2016), evolution of protoneutron stars (Pons et al., 2000), accretion disks (e.g., Shibata & Sekiguchi, 2012; Foucart et al., 2015; Just et al., 2015; Just et al., 2016), and many more in AGN accretion literature. A number of analytical closures have been suggested in the literature. The accuracy of the M1 solution depends on the quality of the closure used and it is a priori unclear which closure yields the best results for a given problem. While the quality of individual closures has been examined in different contexts (e.g., Janka, 1991; Smit et al., 1997; O’Connor, 2015), a systematic analysis for neutrino transport has been performed only by Janka (1991), Smit et al. (2000), and Just et al. (2015)333Note that there is a relation between the M1 scheme and the flux-limited diffusion approximation and each flux limiter is associated with an M1 closure relation (Levermore, 1984; Smit et al., 2000). The quality of some of the flux limiters has been explored by, e.g., Burrows et al. (2000); Just et al. (2015) using the flux-limited diffusion framework for neutrino transport in the context of core-collapse supernovae.. The aim of this work is to extend these two works, to consider a wider selection of M1 closures, to verify them using a wider range of test problems that are relevant to neutrino transport, and to present a quantitative assessment of their quality.

In this work, we evaluate the quality of various closures proposed in the literature by comparing the radiation field distribution in and around radiating objects obtained with the M1 method with the one obtained analytically or with Monte Carlo. We consider two types of radiating objects: a uniform sphere with a sharp surface and a protoneutron star with a hot envelope obtained from core-collapse simulations. These two objects possess the opaque central radiating source surrounded by a transparent envelope, i.e. the important characteristics common to many astrophysical systems. Since our goal is to study the quality of the analytical closures and in order not to contaminate our results with errors emanating from other sources such as hydrodynamics and non-linear radiation-matter coupling, we consider only static matter configurations in our tests. For simplicity, we limit ourselves to spherical symmetry and ignore spacetime curvature around PNSs. Implications of these assumptions will be discussed in Section 5.

The uniform sphere problem consists of an opaque radiating sphere with a sharp surface surrounded by vacuum, and it has an analytical solution (cf. Section 4.2). In the PNS problem, we take three different post-bounce configurations (obtained from simulations of Ott et al. 2008) of a $20M_{\odot}$ progenitor star at 160, 260, and 360 ms after bounce. We obtain precise solution of this problem by performing Monte Carlo radiation transport calculations using the code of Abdikamalov et al. (2012) (cf. Section 4.3). These solutions are compared to M1 solutions obtained using the GR1D code (O’Connor & Ott, 2011, 2013; O’Connor, 2015) available at http://www.GR1Dcode.org.

We consider seven different closures. These are the closures by Kershaw (1976), Wilson et al. (1975), Levermore (1984), the classical maximum entropy closure of Minerbo (1978), and the maximum entropy closure with the Fermi-Dirac distribution by Cernohorsky & Bludman (1994). In addition, we consider two closures by Janka (1991) that are constructed by fitting closure relations to exact Monte Carlo solutions of the radiation field around PNSs (Janka, 1991).

This paper is organized as follows. In Section 2, we give a theoretical overview of the neutrino transport problem and the M1 scheme. In Section 3, we give a brief description of the seven closures we study in this work. Section 4 presents the details of the test problems. We also describe our tools for systematic quantitative assessment of the quality of the closure relations and present the results of our analysis. Our conclusions are provided in Section 5. In Appendix A, we describe the two codes that we use in our analysis: the GR1D code for M1 transport and the Monte Carlo code of Abdikamalov et al. (2012).

2 Boltzmann Equation and M1 method

Neutrinos are described by the distribution function ${\cal F}$, which characterizes the number of neutrinos in a phase-space volume element and which obeys the relativistic Boltzmann equation (e.g., Lindquist, 1966; Mezzacappa & Matzner, 1989):

[TABLE]

Here, $\tau$ is an affine parameter of the neutrino trajectory, $u^{\mu}$ is the four-velocity of the medium, and $p^{\mu}$ is the four-momentum of radiation, from which one can obtain the neutrino energy in the rest frame of the medium via relation $\varepsilon=-p^{\alpha}u_{\alpha}$. The Greek indices $\mu,\alpha=0,1,2,3$ run over space-time components and the Latin indices $i=1,2,3$ runs over the spatial components. $S(p^{\mu},x^{\mu},{\cal F})$ is the collision term that describes the interaction of radiation with matter via absorption, emission and scattering. The evaluation of $S(p^{\mu},x^{\mu},f)$ is a domain of a separate field of study and it is beyond the scope of this work (e.g., Bruenn, 1985; Burrows et al., 2006). In this work, we treat neutrinos as massless particles and fix units using $\hbar=c=1$.

The zeroth, first, and second moments of the distribution function represent the energy density,

[TABLE]

the radiation flux,

[TABLE]

and the radiation pressure,

[TABLE]

Here, $E_{\nu}$, $F^{j}_{\nu}$, and $P^{ij}_{\nu}$ are the functions of neutrino energy $\varepsilon=p^{0}=|\vec{p}|$. In order to obtain the total energy density, flux, and pressure, one has to integrate (2)-(4) over energy, as discussed in, e.g., Thorne (1981); Novikov & Thorne (1973).

The zeroth and the first moments of the Boltzmann equation constitute the system of equations for $E_{\nu}$ and $F^{j}_{\nu}.$ In Minkowski space, spherical symmetry, and neglecting the velocity of the medium, these two evolution equations are

[TABLE]

where $S[0]_{\nu}$ and $S[1]_{\nu}$ are accordingly the zeroth and first moments of the collision term $S(p^{\alpha},x^{\beta},{\cal F})$. Note that, since we consider flat spacetime and static matter, the third moment does not appear in this equation. As we pointed out above, this system is not closed. There are two equations (5) and (6), but three unknowns $E_{\nu}$, $F^{r}_{\nu}$, and $P^{rr}_{\nu}$. This is a simple reflection of the fact that, although the system (5)-(6) is obtained from the Boltzmann equation, it does not capture all the information contained in the Boltzmann equation. To capture the complete information, one has to solve the complete system, which can be expressed as

[TABLE]

where

[TABLE]

is the moment of order $k$. Note that subscripts and superscrips are omitted in this equation to avoid clutter. This is an infinite system of an infinite number of unknowns $M[k]$, which is not feasible to solve in practice.

This situation is somewhat analogous to the Taylor expansion. Function $f(x)$ can be represented through the infinite sum

[TABLE]

This allows one to express the value of $f(x)$ at an arbitrary point $x$ through its properties at a given point $x_{0}$. In order to calculate $f(x)$ exactly, one has to incorporate all the terms in the infinite series. The low-order terms yield accurate results only in the vicinity of the point $x_{0}$. Similar is true for the moments of the distribution function. When we constrain ourselves to the first few moments, we sacrifice the accuracy of our description. To capture all the information contained in the distribution function one needs to employ the whole infinite set of moments.

The M1 approach used in the literature is based on an interpolation of the radiation pressure $P^{ij}$ between optically thick and thin limits (e.g., Shibata et al., 2011)

[TABLE]

where $P^{ij}_{\mathrm{thick}}$ and $P^{ij}_{\mathrm{thin}}$ are the radiation pressure in these limits. In the former limit, radiation is in thermal equilibrium with matter and is isotropic. This results in $F_{\nu}^{i}=0$ and

[TABLE]

for the gas of ultrarelativistic particles such as photons and neutrinos (Mihalas & Weibel-Mihalas, 1999). In the free-streaming limit, radiation propagates like a beam along a certain direction $n$ and exerts pressure only along this direction, giving us $F^{n}_{\nu}=E_{\nu}$ and $F^{i\neq n}_{\nu}=0$ and

[TABLE]

The parameter $p$ in equation (10) is known as the variable Eddington factor and it plays the role of the interpolation factor between the two regimes. The functional form of $p$ in terms of the lower moments is referred to as the M1 closure.

Equation (10) is derived based on the assumption that the radiation is symmetric around the direction parallel to the flux. While the assumption is valid for the spherically symmetric matter and radiation distributions, it does not always hold in more general cases. Colliding radiation beams emanating from different sources is a prominent example. Therefore, equation (10), even before we fix the form of the function $p$, already contains an approximation.

Note that equation (10) in its modern form is often cited as derived by Levermore (1984) in the literature. While it is true, Kershaw (1976) also proposed the interpolation between thick and thin regimes like (10). He then suggested using the simplest among such interpolation

[TABLE]

where

[TABLE]

and $f^{2}=f^{i}f_{i}$. This relation is equivalent to equation (10) with a closure relation $p=(1+2f^{2})/3$, which is known as the Kershaw closure. Even earlier, a formulation similar to M1 was discussed by Pomraning (1969).

3 Closures

In this section, we present a list of seven different closures most commonly used in the literature and describe their main properties.

3.1 Kershaw Closure

The Kershaw (1976) closure is a simple interpolation between the optically thick ($f\to 0$) and the optically thin ($f\to 1$) limits. In the spherically symmetric, case the Kershaw closure reads

[TABLE]

This closure is shown with the solid red line in Fig. 1. In the following, we refer to this closure as the Kershaw closure.

3.2 Wilson Closure

Wilson et al. (1975) and LeBlanc & Wilson (1970) presented a flux-limiter for neutrino diffusion, which corresponds to the closure

[TABLE]

Physically, this expression is equivalent to an interpolation of the diffusive and free-streaming fluxes via harmonic averaging (Smit et al., 2000). This ensures correct diffusive and free-streaming limits, but may yield imprecise results in the intermediate regime. This closure is shown with the solid yellow line in Fig. 1. Hereafter, we refer to this closure as the Wilson closure.

3.3 Levermore Closure

The Levermore closure can be derived assuming that the radiation is isotropic and satisfies the Eddington closure ($P^{ij}_{\nu}=P^{ij}_{\rm thick}$ or $p=1/3$ everywhere) in the “rest frame” of radiation, i.e., in the frame in which the radiation flux is zero (Levermore, 1984; Sa̧dowski et al., 2013):

[TABLE]

This closure is shown with the solid green line in Fig. 1. We refer to this closure as the Levermore closure.

3.4 MEFD: Maximum Entropy Closure for Fermionic Radiation

The idea to use the maximum entropy principle to construct the closure relation was first suggested by Minerbo (1978), who applied it to photons assuming a classical distribution. Later, Cernohorsky & Bludman (1994) applied it to fermions using the Fermi-Dirac distribution.

By maximizing the entropy functional

[TABLE]

under the constraints that the dimensionless zeroth and first moments,

[TABLE]

are given, one can obtain a distribution of radiation in terms of Lagrange multipliers $\eta$ and $a$ (e.g., Smit et al., 2000):

[TABLE]

where $\mu=\cos\theta$. The second moment of equation (21) yields $p$ as a function of $\eta$ and $a$. The closure relation is obtained by expressing $\eta$ and $a$ in terms of $e$ and $f$ through inversion of $e(\eta,a)$ and $f(\eta,a)$:

[TABLE]

where $\chi(x)=1-3/q(x)$ and $q(x)$ is the inverse Langevin function $L(q)\equiv\coth q-1/q$. The lowest-order polynomial approximation to function $\chi(x)$ that has the correct free-streaming and diffusive limits is

[TABLE]

which is $\sim 2\%$ accurate (Cernohorsky & Bludman, 1994; Smit et al., 2000). The substitution of this approximation into equation (22) yields an analytical closure that is a function of both $f$ and $e$. We refer to this closure as the MEFD closure. It is shown in Fig. 1 as a a series of curves for $e=0.1,\ 0.3,\ 0.5,\ 0.7$ and $0.9$ with dashed lines. Note that, in the limit of maximum packing, the MEFD closure reduces to (e.g., Smit et al., 2000)

[TABLE]

This closure, shown with the bottom curve in Fig. 1, represents one boundary of the MEFD closure. The other boundary is the classical limit of this closure, the Maximum Entropy (ME) closure, discussed below.

3.5 ME: Maximum Entropy Closure in the Classical Limit

The classical limit of the MEFD closure is the closure by Minerbo (1978). It can be obtained from equations (22)-(23) by formally taking the $e\rightarrow 0$ limit:

[TABLE]

This closure is shown with the solid sky blue line in Fig. 1. We refer to this closure as the ME closure.

3.6 Janka Closures

Based on extensive Monte Carlo neutrino transport calculations in PNS envelopes, Janka (1991, 1992) constructed several analytic fits to energy-averaged radiation fields, which were parametrized as

[TABLE]

where $a$, $n$, and $m$ are the fitting parameters. We consider two closures corresponding to sets $\{a=0.5,b=1.3064,n=4.1342\}$ Janka_1 and $\{a=1,b=1.345,n=5.1717\}$ Janka_2. The former is obtained by combining the MC outputs for electron neutrinos from two matter distribution models corresponding to extended hot shocked mantle and compact post bounce configuration. The latter closure is obtained from the $\nu_{\mu}$ radiation field of the matter configuration at $300\,\mathrm{ms}$ after bounce. These two closures are shown in Fig. 1 with dark and bright purple colors, respectively.

4 Results

In order to asses the quality of M1 results, we consider the radiation field in and around the uniform sphere (Section 4.2) and a set of protoneutron star models (Section 4.3). The former case has an analytical solution, while the latter is calculated with the MC method using the code of Abdikamalov et al. (2012). Both of these problems have the central opaque region and outer transparent envelope common to many astrophysical sources.

4.1 Quantitative Estimate of Accuracy

In order to estimate the accuracy of the M1 results, we use the normalized mean square deviation and the spectrum-weighted mean square deviation. The former is defined as

[TABLE]

Here, $Y$ stands for any quantity we want to compare (e.g., energy density, flux factor, and etc.), while $Y_{0}$ is the “exact” value of this quantity obtained from the analytical solution or a Monte Carlo calculation. $X$ is a variable on which both $Y$ and $Y_{0}$ depend (e.g., the radial coordinate) and $X_{i}$ are its discrete values ranging from $X_{\rm min}$ to $X_{\rm max}$. Thus, $\delta{Y}$ provides an estimate of how well the closure solution approximates the exact solution in the entire range from $X_{\rm min}$ to $X_{\rm max}$.

The spectrum-weighted mean square deviation is defined as

[TABLE]

where $i$ is the index of the neutrino energy group and $\delta Y_{i}$ is defined by equation (27) for each energy group independently. The spectral weights $w_{i}$ are obtained using the spectral energy density $S_{i}$ at energy $\varepsilon_{i}$ and the peak value of spectral energy density $S_{\rm max}$. In our analysis of the spectral weighted quantities, we restrict ourselves to the energies lying near the spectral peak. More specifically, we consider only the energy groups with spectral energy densities greater than $0.3S_{\rm max}$ in order to cut out low statistics energy groups.

4.2 Uniform Sphere

The uniform sphere problem consists of a static homogeneous and isothermal sphere of radius $R$ surrounded by vacuum. Matter inside the sphere can absorb and emit radiation. This problem has an analytical solution and possesses important physical and numerical characteristics. The central opaque source with transparent outer regions are characteristics of many astrophysical systems, while the sharp surface represents a serious challenge for many numerical techniques. For this reason, this problem is often used as a test problem for radiation transport codes (Schinder & Bludman, 1989; Smit et al., 1997; Rampp & Janka, 2002; O’Connor, 2015)

The analytical solution for the distribution function is given as

[TABLE]

where $r$ is the radial coordinate, $R$ is the radius of the sphere, $\mu=\cos\theta$,

[TABLE]

and

[TABLE]

Inside the sphere the absorption coefficient $\kappa$ and emissivity $B$ are constants. Outside the sphere, there is no emission and absorption. For our test, we use $\kappa R=7500$ and $B=1$, which ensures that radiation is fully isotropic inside the sphere and a tiny region $\sim 1/\kappa\ll R$ separates it from the free-streaming regime outside the sphere.

The flux factor $f$ (top panel) and the Eddington factor $p$ (bottom panel) as a function of the radial coordinate are shown on Fig. 3. The dim gray line represents the analytical solution, while the rest of the lines represent the solutions obtained with M1 approximations. The values of normalized mean square deviations of these solutions from the analytical result are listed in Table 1.

As we see, all closures perform poorly for this problem. The Kershaw closure yields significantly worse results than the rest of the closures. The normalized mean square deviation of the flux and Eddington factors obtained with closure prescriptions from analytical result are $0.13$ and $0.32$, respectively (cf. Table 1). The next worse performers are the Levermore and Janka_2 closures. The normalized deviations of the Eddington factor are $0.22$ and $0.21$ for these two closures, respectively. The best performers are the Wilson and the Janka_1 closures, for which the normalized deviations of the Eddington factor are $0.14$ and $0.13$. The rest of the closures yield intermediate results.

Such a poor result of the analytical closures has a simple explanation. The true closure for this problem, which can be reconstructed directly from the solution (29), is shown with a solid black line in Fig. 3. This closure is significantly different from all the aforementioned closures, as we can glean by comparing Fig. 3 with Fig. 3. Therefore, the reason why our closures yield inaccurate results is simply because these closures are different from the true closure for this problem.

Interestingly, the true closure in Fig. 3 can be fit well with the following simple analytical expression

[TABLE]

which is shown with the orange dashed line in Fig. 3. If we perform M1 calculations with this closure, it reproduces the analytical result (29) with excellent $\sim 1\%$ accuracy (cf. Table 1). Note that this fit closure is not expected to perform well for any other matter distributions except the uniform sphere. It is specific to this particular model problem.

4.3 Protoneutron Star

In this section, we evaluate the ability of the M1 closures to model the neutrino radiation field around spherically symmetric models of protoneutron stars (PNSs) formed in core-collapse supernovae. We take three different post-bounce configurations (obtained from 2D radiation-hydrodynamics simulations of Ott et al. 2008) of a $20M_{\odot}$ progenitor star at 160 ms, 260 ms, and 360 ms after bounce. We average the 2D profiles of Ott et al. (2008) over angle to obtain spherically symmetric configurations of PNSs. The radial profiles of density, temperature, and electron fraction are shown in upper, center, and bottom panels of Fig. 4. The spectra of neutrino luminosity obtained from MC code at the radius of $100\,\mathrm{km}$ are shown in Fig. 5. The top, center, and bottom panels represent the PNS models at 160, 260, and 360 ms after bounce, respectively.

The “exact” solution of the problem is obtained from MC calculations using the code of Abdikamalov et al. (2012), while the M1 results are obtained using the GR1D code (O’Connor, 2015). We evolve our time dependent MC code until we reach steady-state neutrino radiation field for each model of PNS. This field is then averaged over many timesteps until we get rid of the stochasticity in the MC solution. Since the MC method does not use any approximations in the solution procedure (Abdikamalov et al., 2012), the solution obtained this way is exact for a given configuration of matter (i.e., a given configuration of opacities and emissivities). In order to ensure the consistency of the results, the two codes use identical microphysical inputs. Both use the Shen et al. (1998) equation of state (EOS) table and a NuLib opacity table (O’Connor, 2015) that was generated with the same EOS table. In both codes, we use 48 logarithmic energy groups ranging from $0.5\,\mathrm{MeV}$ to $200\,\mathrm{MeV}$. In the MC code, we use $150$ radial logarithmically spaced zones with the central resolution of $0.2\,\mathrm{km}$. We have performed extensive resolution tests in order to ensure that our results are convergent.

We examine the quality of the closures in the free streaming, semi-transparent, and opaque regimes. We separate these regimes based on the value of the flux factor $f$. We choose the transparent regime as the one where $0.9<f<1$, the semi-transparent as $0.5\leq f\leq 0.9$ and the opaque as $0.2\leq f\leq 0.5$. We neglect the region of low $f$ because the MC results suffer from noise in the highly diffusive region. It is more appropriate to separate the different regimes based on $f$ rather than, e.g., the value of the radial coordinate, because at a given radius, neutrinos of different energies behave differently. For example, low-energy neutrinos have lower opacity and hence they propagate more freely compared to higher-energy neutrinos at the same radius.

The radial profiles of the flux factor $f$ obtained using the seven closures and from the Monte Carlo code (dotted line) for three different types of neutrinos and three PNS configurations at 160, 260, and 360 ms after bounce are shown in Fig. 6. All flux factors are measured at the neutrino energy groups corresponding to the peak luminosity at $100\,\mathrm{km}$ for each neutrino type (cf. Fig. 5). As we can see, all M1 closures yield qualitatively correct results.

For a more precise quantitative estimate, we utilize the spectrum-weighed deviations of the flux factor and energy density from the MC results, equation 28. The spectrum-weighted deviation $\bar{\delta}e$ of the energy density $e$ in the transparent regime for the seven different closures, for the three neutrino types, and the PNS models at $160$, $260$, and $360$ ms are shown on the top left panel of Fig. 7. The deviations $\bar{\delta}e$ are calculated using formula (28), in which the spectrum is taken at $100\,\mathrm{km}$. In order to verify that our results are not too dependent on spectral information at different locations, we have calculated $\bar{\delta}e$ using spectra at several different radii and obtain results very similar to $\bar{\delta}e$ that use spectra at $100\,\mathrm{km}$. This suggests that the values of $\bar{\delta}e$ presented in Fig. 7 are robust measures of errors of the M1 closures for the radiation field around PNSs.

As we can see, different closures yield different levels of accuracy depending on the neutrino type and PNS model. No single closure performs consistently better (or worse) than other closures in all cases. That said, the Wilson and the Levermore closures perform better than others in most cases in the transparent regime, followed by the ME and MEFD closures. The Janka_1 and the Janka_2 closures exhibit $\bar{\delta}e\gtrsim 0.04$ in most cases, which is worse than $\bar{\delta}e$ for the rest of the closures. This is a remarkable result because the Janka_1 and Janka_2 closures were constructed from fitting to the exact solution for the neutrino radiation field around PNS. This demonstrates that a closure constructed for one PNS model (with a given EOS and opacity table) does not necessarily yield a good result for all other PNS models.

The behavior of the deviation of the flux factor $\bar{\delta}f$ is shown on the top right panel of Fig. 7. In this case, the Levermore closure performs better than other closures in almost all cases and yields small $\bar{\delta}f\sim 0.004$. The Janka_2 closure performs only slightly worse than this closure, yielding $0.004\lesssim\bar{\delta}f\lesssim 0.009$. However, this is $\sim 2-3$ times smaller than what the rest of the closures yield, which is a surprising result because the Janka_2 closure yields relatively poor result for the energy density compared to most of the closures, as we discussed in the previous paragraph.

The Kershaw closure produces the least accurate $f$ compared to the other closures in all models except the PNS models at $160\,\mathrm{ms}$ after bounce, yielding deviations of $\bar{\delta}f\sim 0.015$ in all cases. The Wilson closure yields intermediate results in most situations except for the PNS model at $160$ ms after bounce, for which it yields the largest deviation of $\bar{\delta}f\sim 0.015$. This again shows that a closure that yields the best result for the energy density does not necessarily yield the best result for the flux factor. As we see below, the reverse of this statement is also true.

The spectrum-weighed deviations $\bar{\delta}e$ and $\bar{\delta}f$ for the semi-transparent regime are shown in the left and right center panels of Fig. 7, respectively. In this regime, both the ME and MEFD closures often – but not always – yield the smallest deviations $\bar{\delta}e$. The Janka_1 and Janka_2 closures yield the largest $\bar{\delta}f$ of $\sim 0.04-0.06$ in all cases. The Wilson closure yields the smallest $\bar{\delta}e$ in most cases, but produces $\bar{\delta}f\sim 0.03$, which is roughly the mean of the values of $\bar{\delta}f$ produced by all of the closures. On the other hand, the Kershaw closure yields the smallest $\bar{\delta}f$ of $\sim 0.01-0.02$ in most cases, but yield relatively large $\bar{\delta}e$ of $\sim 0.035-0.05$.

In the opaque region, the situation is significantly different. The Kershaw closure, which often yields the largest $\bar{\delta}e$ and $\bar{\delta}f$ in the transparent and $\bar{\delta}e$ in the semi-transparent regimes, produces the smallest $\bar{\delta}e$ and $\bar{\delta}f$ of $\lesssim 0.05$ in the opaque regime. The Wilson closure yields slightly worse results than the Kershaw closure with $\bar{\delta}e\sim\bar{\delta}f\sim 0.05$. The Janka_1 and Janka_2 closures yield the largest errors ($\bar{\delta}e\sim 0.1$ and $\bar{\delta}f\sim 0.1-0.15$). This again shows that a closure that yields good results for one model of PNS does not necessarily produce good results for other PNS models.

In all cases, the ME and MEFD maximum entropy closures yield almost identical results. This is direct consequence of the fact the former closure is the classical limit of the latter and the radiation field, in the regions we considered, is predominately in the classical regime.

Overall, these two closures yield relatively good results in all cases and never result in the largest deviations compared to the other closures. This, in combination with the fact that the ME closure is simpler and requires fewer operations to compute than the MEFD closure, makes the former a more attractive option for neutrino transport applications involving protoneutron stars.

These results demonstrate that no single closure performs the best or the worst in all cases. Whether a specific closure is “correct” for a given problem depends on the parameters of the problem, such as the matter configuration (i.e., profiles of density, temperature, and composition) and the neutrino type. We find a similar level of differences between deviations of $e$ and $f$ for different neutrino energy groups (not shown here). Fig. 8 shows the functional form of the closures $p(f)$ extracted from the MC simulations for different neutrino types, different PNS models and different energy groups. We can clearly see variations between closures $p(f)$ corresponding to different problems.

5 Conclusion

We conducted a systematic, quantitative study of the accuracy of analytical closure relations for two-moment neutrino radiation transport schemes commonly used in the literature. We considered the neutrino field around two sets of radiating objects: the uniform radiative sphere and PNS models at $160$, $260$, and $360$ ms after core bounce obtained from simulations of Ott et al. (2008). In all cases, the matter configuration is assumed to be static. This restriction allows us to focus on the quality of the closure relations and exclude other sources of errors such as those stemming from non-linear radiation-matter coupling.

We considered seven different closures. These are the closures by Kershaw (1976), Wilson et al. (1975), Levermore (1984), and the maximum entropy closures of Cernohorsky & Bludman (1994) and Minerbo (1978). In addition, we considered two closures that are constructed by fitting to exact Monte Carlo solutions of the radiation field around PNSs by Janka (1991).

We find that no single closure, among those studied here, is consistently better or worse than any other. A closure that yields accurate results in one case may not yield as good results in other situations. The level of accuracy that a given closure yields varies for different quantities.

Given this limitation of the closures, the maximum entropy closure by Minerbo (1978) and Cernohorsky & Bludman (1994), which yield almost identical results, often yield better results among all the closures studied. These two closures never yielded the worst results compared to all other closures. In this sense, these two are a safe choice, as they are less likely to yield extremely erroneous result over a wide range of problems and variables. Since the Minerbo (1978) closure is simpler to compute than the one by Cernohorsky & Bludman (1994), we conclude that the former closure is the most attractive choice for problems involving neutrino transport around PNSs.

In this work, we assumed spherical symmetry, static matter, and flat spacetime. These assumptions limit the scope of the implications of our results. In particular, the hydrodynamic and general relativistic effects introduce shifts in the energy spectrum of neutrinos, which alter the moment equations (e.g, Cardall et al., 2013; Just et al., 2015). Also, in non-spherically-symmetric cases, we can have radiation sources at different spatial locations. Notable examples are the accretion disks formed in neutron star mergers and hot spots near the PNS surfaces in the context of core-collapse supernovae. The interaction of radiation beams from such sources cannot be modeled within our spherically-symmetric setup. It is a priori unclear the extent to which our results are valid in such cases. This will the subject of a future study.

Acknowledgements

We thank Evan P. O’Connor, Christian D. Ott and John Wendell for their input. We also thank Adam Burrows, David Radice, Sherwood Richers, and Luke Roberts for helpful discussions. This work was supported by IGPPS program at Los Alamos National Laboratory (2012-2015) and partially by the Sherman Fairchild Foundation. EM is grateful to David and Barbara Groce for their kindness and support. EA acknowledges support from NU ORAU and Social Policy grants. LANL report number LA-UR-17-20390.

Appendix A Methods

A.1 The GR1D Radiation-Hydrodynamics Code

We employ the M1 radiation transport solver that is part of the GR1D radiation-hydrodynamics Code (O’Connor & Ott, 2011, 2013; O’Connor, 2015) available at http://www.GR1Dcode.org. The matter in our calculations is static and the metric is Minkowski. There is no coupling between different energy groups. For these conditions, GR1D implements M1 closures by setting the pressure tensor according to equation (10). Under these conditions, the transport equations read

[TABLE]

where $\eta_{\nu}$, $\kappa_{\nu,a}$, and $\kappa_{\nu,s}$ are the neutrino emissivity, absorption opacity, and scattering opacity, respectively. To numerically solve these equations we discretize the neutrino spectrum into 48 energy groups, logarithmically spaced between 0.5 MeV and 180 MeV. For each energy group and species, we compute the closure ($p=P_{\nu}^{rr}/E$), the spatial flux terms ($\partial_{r}r^{2}F_{\nu}^{r}$ and $\partial_{r}r^{2}P_{\nu}^{rr}$), and the values of the neutrino interaction coefficients ($\eta_{\nu}$, $\kappa_{\nu,a}$, and $\kappa_{\nu,s}$) explicitly at the beginning of the time step (denoted via the index $(n)$). We then use a first order, implicit/explicit time integration method to solve for the values of the energy and momentum density at time $t+\Delta t$ (or the $(n+1)$ time step),

[TABLE]

The explicit calculation of the spatial flux remains valid in the diffusion limit due to corrections applied to the Riemann solution in high optical depth regions (O’Connor, 2015). The neutrino interaction coefficients are computed using NuLib, an open-source neutrino interaction library available at http://www.nulib.org (O’Connor, 2015). For PNS calculations, we include elastic scattering of neutrinos on nucleons, and coherent elastic scattering on alpha particles and heavy nuclei as contributions to $\kappa_{\nu,s}$. Charged current absorption of electron neutrinos on neutrons and heavy nuclei and electron antineutrinos on protons is included in $\kappa_{\nu,a}$ and in $\eta_{\nu}$ via Kirchhoff’s law. For heavy lepton neutrinos we determine the emissivity $\eta_{\nu}$ from pair processes (electron-positron annihilation and nucleon-nucleon Bremsstrahlung) and an effective absorption opacity via an approximation that works well for supernova conditions (O’Connor, 2015).

A.2 The Monte-Carlo Neutrino Transport Code

In order to asses the quality of the M1 closures, we compare the M1 results to Monte Carlo radiation transport calculations using the code of Abdikamalov et al. (2012). Here, we outline some salient aspects of such methods, while more in-depth discussion can be found in Abdikamalov et al. (2012).

Monte Carlo methods have been used for many applications (e.g., Janka & Hillebrandt, 1989; Janka, 1991, 1992; Burrows & van Riper, 1995; Wolf et al., 1999; Keil et al., 2003; Lucy, 2005; Kasen et al., 2006; Densmore et al., 2007; Dolence et al., 2009; Wollaeger et al., 2013; Richers et al., 2015). Such methods use sequences of pseudo-random numbers to simulate the transport of radiation using the concept of Monte Carlo particles. Each MC particle represents a group of physical particles with a given location, direction, and energy (or frequency). The number of physical radiation particles (i.e., photons or neutrinos) represented by a MC particle is called the weight of the MC particle. The smaller the weight, the larger the number of MC particles that are needed to model a given problem, which means higher precision at the price of higher computational cost.

The spatial problem domain is divided into a number of interconnected cells, and matter in each cell has its own temperature, density, and composition. Using this information, one calculates the number of MC particles that has to be emitted in each cell within a timestep. These particles are then randomly placed in each cell by randomly sampling their coordinates. The propagation direction is sampled randomly with isotropic distribution. The frequencies/energies of particles are chosen using the form of the energy-dependence of the emissivity function.

Once the MC particles are placed in each cell, they are then transported within a timestep. This can be achieved by calculating three distances for each particle: the distance to collision (absorption or scattering), the distance to the boundary of its cell, and the distance the particle would travel until the end of timestep if it were to travel freely (i.e., assuming that no collision happens). What happens to a particle depends on which of these distances is the smallest. If the distance to collision is the smallest, then the particle is either absorbed or moved to its new location and scattered. The probabilities of absorption and scattering are proportional to their relative opacities. If the distance to a cell boundary is the smallest, then the particle is moved to the boundary of the next cell and transported further in the new cell by calculating a new set of three distances. If the particle crosses the outer boundary of the computational domain, it is removed from the system. Finally, if the third distance is the smallest, then the particle moves by that distance within its cell. Once this step is accomplished for all of the MC particles, the remaining particles, i.e., the ones that have not been absorbed or left the domain, are stored in memory as a preparation for the next timestep. At the next timestep, the cells are populated with newly emitted MC particles, which together with the MC particles remaining from the previous step are then transported further within the second timestep. Subsequent timesteps are performed in the same manner until the end of the simulation.