Solving the gamma-ray radiative transfer equation for supernovae

TL;DR

This paper introduces a new relativistic radiative-transfer code for gamma-rays in supernova ejecta, enabling accurate and fast modeling of gamma-ray transport and energy deposition, with results comparable to Monte Carlo methods.

Contribution

The paper presents a novel relativistic radiative-transfer code for gamma-ray transport in supernovae, improving accuracy and computational speed over previous methods.

Findings

Code yields results comparable to Monte Carlo simulations.

Differences observed compared to pure absorption models.

Synthetic spectra show asymmetric line profiles at early times.

Abstract

We present a new relativistic radiative-transfer code for $\gamma$-rays of energy less than 5 MeV in supernova (SN) ejecta. This code computes the opacities, the prompt emissivity (i.e. decay), and the scattering emissivity, and solves for the intensity in the co-moving frame. Because of the large expansion velocities of SN ejecta, we ignore redistribution effects associated with thermal motions. The energy deposition is calculated from the energy removed from the radiation field by scattering or photoelectric absorption. This new code yields comparable results to an independent Monte Carlo code. However, both yield non-trivial differences with the results from a pure absorption treatment of $\gamma$-ray transport. A synthetic observer's frame spectrum is also produced from the CMF intensity. At early times when the optical depth to $\gamma$-rays is large, the synthetic spectrum show asymmetric line profiles with redshifted absorption as seen in SN 2014J. This new code is integrated within CMFGEN and allows for an accurate and fast computation of the decay energy deposition in SN ejecta.

Tables

Table 1. Table 1: Example nuclear decay data for the Ni 56 superscript Ni 56 {}^{56}{\rm Ni} → → \rightarrow Co 56 superscript Co 56 {}^{56}{\rm Co} → → \rightarrow Fe 56 superscript Fe 56 {}^{56}{\rm Fe} decay chain. t 1 / 2 subscript 𝑡 1 2 t_{1/2} , Q γ subscript 𝑄 𝛾 Q_{\gamma} , Q th subscript 𝑄 th Q_{\rm th} are the half-life, energy per decay, and thermal energy of the leptons produced. We list the decay line energies E γ subscript 𝐸 𝛾 E_{\gamma} and probabilities for lines with probabilities ≥ \geq 1 per cent. This data and all other nuclear decay data is taken from http://www.nndc.bnl.gov/chart/ .

${}_{}{}^{56}\mathrm{Ni}$ $\to$ ${}_{}{}^{56}\mathrm{Co}$ $\to$ ${}_{}{}^{56}\mathrm{Fe}$
${}_{}{}^{56}\mathrm{Ni}$ $\to$ ${}_{}{}^{56}\mathrm{Co}$		${}_{}{}^{56}\mathrm{Co}$ $\to$ ${}_{}{}^{56}\mathrm{Fe}$
$t_{1/2}$ = 6.075 days		$t_{1/2}$ = 77.233 days
$Q_{\gamma }$ = 1.718 MeV		$Q_{\gamma }$ = 3.633 MeV
$Q_{\mathrm{th}}$ = 0.000 MeV		$Q_{\mathrm{th}}$ = 0.116 MeV
$E_{\gamma }$	Prob.	$E_{\gamma }$	Prob.
(MeV)		(MeV)
0.158	98.8	0.511	38.0
0.270	36.5	0.847	100
0.480	36.5	0.977	1.4
0.750	49.5	1.038	14.0
0.812	86.0	1.175	2.3
1.562	14.0	1.238	67.6
		1.360	4.3
		1.771	15.7
		2.015	3.1
		2.035	7.9
		2.598	17.3
		3.010	1.0
		3.202	3.2
		3.253	7.9
		3.273	1.9

Table 2. Table 2: Listed is the total energy deposition integrated over the whole ejecta ( E dep subscript 𝐸 dep E_{\rm dep} ) and the integrated flux from the synthetic spectrum ( L flux subscript 𝐿 flux L_{\rm flux} ).

	$t_{\exp }$ = 17.4 days			$t_{\exp }$ = 207.0 days
	$E_{\mathrm{dep}}$ (ergs s-1)	$L_{\mathrm{escape}}$ (ergs s-1)	$L_{\mathrm{flux}}$ (ergs s-1)	$E_{\mathrm{dep}}$ (ergs s-1)	$L_{\mathrm{escape}}$ (ergs s-1)	$L_{\mathrm{flux}}$ (ergs s-1)
This work	1.260(43)	9.669(41)	8.466(41)	9.273(40)	1.343(42)	1.359(42)
Hillier & Dessart (2012)	1.279(43)	7.892(41)	7.762(41)	9.291(40)	1.343(42)	1.305(42)

Equations

$$\frac{\gamma (1+\beta \mu )}{c}\frac{\partial I_{\nu }}{\partial t}+\gamma (\mu +\beta )\frac{\partial I_{\nu }}{\partial r}$$

$$+\gamma (1-{\mu }^2)\left[\frac{1+\beta \mu }{r}-\Lambda \right]\frac{\partial I_{\nu }}{\partial \mu }$$

$$-\gamma \nu \left[\frac{\beta (1-{\mu }^2)}{r}+\mu \Lambda \right]\frac{\partial I_{\nu }}{\partial \nu }$$

$$+3\gamma \left[\frac{\beta (1-{\mu }^2)}{r}+\mu \Lambda \right]I_{\nu }$$

$$\Lambda =\frac{{\gamma }^2(1+\beta \mu )}{c}\frac{\partial \beta }{\partial t}+{\gamma }^2(\mu +\beta )\frac{\partial \beta }{\partial r}.$$

$$\frac{dr}{ds}=\gamma (\mu +\beta )\text{and}\frac{d\mu }{ds}=\gamma (1-{\mu }^2)\left[\frac{1+\beta \mu }{r}-\Lambda \right].$$

$$\frac{\partial I_{\nu }}{\partial s}-\nu \Pi \frac{\partial I_{\nu }}{\partial \nu }={\eta }_{\nu }-({\chi }_{\nu }+3\Pi )I_{\nu },$$

$$\Pi =\gamma \left[\frac{\beta (1-{\mu }^2)}{r}+\mu \Lambda \right]$$

$${\chi }_{\nu }^{\mathrm{abs}}={\left(\frac{m_e{c}^2}{h\nu }\right)}^{3.5}{\sigma }_{\mathrm{T}}{\alpha }^{4}8\sqrt{2}\sum _i^{N_{spec.}}{N}_i{Z}_i^5,$$

$${\chi }_{\nu }^{\mathrm{C}}=2\pi r_e^2{N}_e[\left(\frac{1+x}{x^3}\right)\left\{\frac{2x(1+x)}{1+2x}-\log (1+2x)\right\}+$$

$$\frac{\log (1+2x)}{2x}-\frac{1+3x}{(1+2x{)}^2}],$$

$${\eta }_{\nu }^{\mathrm{iso}}=\frac{1}{4\pi }\sum _{i=1}^{N_{isot.}}\sum _{j=1}^{N_{\mathrm{lines}}}\frac{N_i}{{\tau }_i}{E}_{ij}{P}_{ij}\frac{e^{-\Gamma }}{\sqrt{2\pi }{V}_{\mathrm{Dop}}{\nu }_{ij}/c},$$

$${\eta }_{{\nu }^{\prime }}^{\mathrm{s}}(r,{\mathbf{\Omega }}^{\prime })={\int }_0^{\infty }\frac{{\nu }^{\prime }}{\nu }d\nu \oint d\mathbf{\Omega }{\sigma }_s(\nu \to {\nu }^{\prime },\xi )I_{\nu }(r,\mathbf{\Omega }),$$

$$\xi =\mathbf{\Omega }\cdot {\mathbf{\Omega }}^{\prime }=\sqrt{1-{\mu }^2}\sqrt{1-{\mu }^{\prime 2}}\cos (\varphi -{\varphi }^{\prime })+\mu {\mu }^{\prime }.$$

$${\sigma }_s(\nu \to {\nu }^{\prime },\xi )$$

$$\equiv$$

$$\delta \left[\xi -\left(1-\frac{1}{x^{\prime }}+\frac{1}{x}\right)\right]\to$$

$$\frac{\delta \left[\varphi -{\varphi }_1\right]}{\left|\sqrt{(1-{\mu }^2)(1-{\mu }^{\prime 2})-(1-\frac{1}{x^{\prime }}+\frac{1}{x}-\mu {\mu }^{\prime }{)}^2}\right|}+$$

$$\frac{\delta \left[\varphi -{\varphi }_2\right]}{\left|\sqrt{(1-{\mu }^2)(1-{\mu }^{\prime 2})-(1-\frac{1}{x^{\prime }}+\frac{1}{x}-\mu {\mu }^{\prime }{)}^2}\right|},$$

$${\varphi }_1$$

$${\varphi }_2$$

$${\mu }_{1,2}=\frac{a_1\pm \sqrt{a_1^2+4a_2}}{2}$$

$$\text{for}{a}_1=2{\mu }^{\prime }\left[1+\frac{1}{x}-\frac{1}{x^{\prime }}\right]\text{and}$$

$$a_2=\left[\frac{2}{x{x}^{\prime }}+\frac{2}{x^{\prime }}-\frac{1}{x^2}-\frac{2}{x}-\frac{1}{x^{\prime 2}}-{\mu }^{\prime 2}\right].$$

$${\eta }_{{\nu }^{\prime }}^{\mathrm{s}}(r,{\mu }^{\prime })=N_e{r}_o^2{\int }_0^{\infty }\frac{{\nu }^{\prime }}{\nu }d\nu \sigma (\nu ,{\nu }^{\prime }){\int }_{{\mu }_2}^{{\mu }_1}d\mu F(\nu ,{\nu }^{\prime },r,\mu ,{\mu }^{\prime }),$$

$$F(\nu ,{\nu }^{\prime },r,\mu ,{\mu }^{\prime })=\frac{I_{\nu }(r,\mu )}{\sqrt{(1-{\mu }^2)(1-{\mu }^{\prime 2})-(1-\frac{1}{x^{\prime }}+\frac{1}{x}-\mu {\mu }^{\prime }{)}^2}}.$$

$${\int }_{-1}^1\frac{f(x)dx}{\sqrt{1-x^2}}=\sum _{i=1}^n\frac{\pi }{n}f(b_i),$$

$$\left(\begin{array}[]{cc}\mu_{1}&1\\ \mu_{2}&1\end{array}\right)\left(\begin{array}[]{c}c_{1}\\ c_{2}\end{array}\right)=\left(\begin{array}[]{c}1\\ -1\end{array}\right).$$

$$\left(\begin{array}[]{c}c_{1}\\ c_{2}\end{array}\right)=\frac{1}{\mu_{1}-\mu_{2}}\left(\begin{array}[]{c}2\\ -\mu_{1}-\mu_{2}\end{array}\right).$$

$$\frac{d\mu }{\sqrt{(1-{\mu }^2)(1-{\mu }^{\prime 2})-(1-\frac{1}{x^{\prime }}+\frac{1}{x}-\mu {\mu }^{\prime }{)}^2}}\to \frac{dw}{\sqrt{1-w^2}}.$$

$${\eta }_{{\nu }^{\prime }}^s(r,{\mu }^{\prime })$$

$$\underset{{\mu }^{\prime }\to 1}{\lim }(b_i-c_2)/c_1$$

$$\underset{{\mu }^{\prime }\to -1}{\lim }(b_i-c_2)/c_1$$

$$\delta \left[\xi -\left(1-\frac{1}{x^{\prime }}+\frac{1}{x}\right)\right]\to \delta \left[\pm \mu -\left(1-\frac{1}{x^{\prime }}+\frac{1}{x}\right)\right]$$

Full text

††pagerange: Solving the gamma-ray radiative transfer equation for supernovae–acknowledgements††pubyear: 2016

1 Introduction

Supernovae (SNe) are luminous astrophysical events, and studies of SNe probe stellar evolution of the progenitor and reveal properties of the explosion mechanism. Understanding both spectra and light curves allows us to investigate the physics and retrieve SN ejecta properties. For instance, Type Ia SNe produce large amounts of radioactive material that control the thermal evolution of the ejecta by non-thermal heating.

For decades, the standard paradigm has been that SNe arise through two mechanisms: gravitational core-collapse (CC) and thermonuclear (Type Ia). Historically, SNe have been classified into two spectral types, Type I (no H i lines) and Type II (strong H i lines) (Minkowski, 1941). However, the most successful theory of SNe is that SNe of Type Ib, Ic, and Ibc, and II result from core-collapse – see Colgate & White (1966); Burrows et al. (1995); Janka & Mueller (1996); Mezzacappa et al. (1998). On the other hand, Type Ia SNe are believed to be thermonuclear explosions of carbon-oxygen (CO) white dwarfs (WDs) (Hoyle & Fowler, 1960). These thermonuclear explosions produce large amounts of radioactive material ($\sim$0.6 $M_{\odot}$) (see, e.g., Scalzo et al., 2014), mainly ${}^{56}{\rm Ni}$, which decays into ${}^{56}{\rm Co}$ and then ${}^{56}{\rm Fe}$. Core-collapse SNe (CCSNe) are thought to produce about an order of magnitude less ${}^{56}{\rm Ni}$ than this (\sim$$10^{-2}–$10^{-3}$ $M_{\odot}$), which is what powers the late time light curve (Janka, 2012, see review of CCSNe).

A crucial issue for modelling SNe is the location of radioactive material. If the radioactive material is mixed into the outer ejecta, it will heat it and enhance the ionization. In Type Ia SNe, mixing of ${}^{56}{\rm Ni}$ has been invoked to explain the brightness and colour at very early times (Höflich et al., 1998, 2002; Woosley et al., 2007; Hoeflich et al., 2017). It was also invoked in SN1987A to explain the early detection of X-rays and $\gamma$-rays from 1987A (Pinto & Woosley, 1988b, a; Bussard et al., 1989; The et al., 1990; Dessart et al., 2012, and references therein).

The peculiar Type II SN, SN1987A, is the only CCSN for which we have detected the ${}^{56}{\rm Co}$ decay lines at 847 and 1238 keV (Makino & Moore, 1987; Matz et al., 1988b, c, d, a; Sunyaev et al., 1987; Tanaka, 1988; Cook et al., 1988). After SN1987A was observed, models of expected late time (1–2 years post-explosion due to high initial column densities) $\gamma$-ray and X-ray fluxes and profile shapes calculated from Monte Carlo radiative transfer soon followed (Pinto & Woosley, 1988b, a; Bussard et al., 1989; The et al., 1990, and references therein). To date, SN 2014J is the only Type Ia SN with $\gamma$-ray detections (Churazov et al., 2014, 2015).

Many $\gamma$-ray radiative-transfer codes have utilized Monte Carlo techniques to treat the radiative transfer (Pozdnyakov et al., 1983; Hoeflich et al., 1992; Milne et al., 2004; Sim, 2007; Sim & Mazzali, 2008; Hillier & Dessart, 2012; Summa et al., 2013). Another technique, used by Swartz et al. (1995) and Jeffery (1998), utilizes a grey transfer approach to treat $\gamma$-ray transport. Swartz et al. (1995) finds that the value of $\kappa_{\gamma}=0.06Y_{e}$ cm2 g*-1*, where $Y_{e}$ is the total number of electrons per baryon, best describes the interaction of $\gamma$-rays in the SN ejecta. In contrast, the work presented here is the first of its kind to formally solve the radiative transfer equation for $\gamma$-rays for SN ejecta.

This paper is organized as follows: In Section 2 we outline the method used to calculate the opacity (Compton scattering and X-ray photoelectric absorption) and emissivity (prompt emission and scattering) needed to solve the relativistic radiative transfer equation. The implementation of our method into cmfgen is discussed in Section 2.5. In Section 3 we illustrate our results using a SN Ia ejecta resulting from a delayed-detonation in a Chandrasekhar mass (MCh) WD from Wilk et al. (2018). We also present synthetic $\gamma$-ray/X-ray spectra around bolometric maximum and at nebular times, and compare our results with those from a Monte Carlo calculation and those obtained using the grey approximation. Finally, in Section 4, we summarise our results.

2 Technique

We developed this code for implementation into cmfgen (Hillier & Miller, 1998; Hillier & Dessart, 2012; Dessart et al., 2014), which is a radiative transfer code that solves the spherically symmetric, non-local-thermodynamic-equilibrium (non-LTE), time-dependent, relativistic radiative transfer equation in the co-moving frame (CMF). This work was undertaken as a consistency check of the Monte Carlo (MC) radiative transfer code utilized by cmfgen (Hillier & Dessart, 2012), and to provide an alternative technique to track photons and subsequent Compton scatterings or photon absorption for computation of the energy deposition in SNe. Since the expansion velocities dominate over thermal motions, this work ignores effects of thermal redistribution.

2.1 Radiative Transfer Equation

We implement the code by solving the relativistic radiative transfer equation along rays as outlined in Olson & Kunasz (1987), Hauschildt (1992), and Hillier & Dessart (2012). The relativistic radiative transfer equation is

[TABLE]

where $\beta=v/c$, $\;\gamma=1/\sqrt{1-\beta^{2}}$, $\;\mu=\cos\theta$, and

[TABLE]

In Eq. 2.1, the specific intensity, emissivity, and opacity (all measured in the CMF) are assumed to be functions of several variables $[I_{\nu}$=$I(t,\,r,\,\mu,\,\nu),\;\eta_{\nu}$=$\eta(t,\,r,\,\mu,\,\nu),\text{ and }\chi_{\nu}$=$\chi(t,\,r,\,\nu)]$. However, if we ignore all time dependence, we can reduce this equation along characteristic rays reducing the partial differential equation with dependent variables $(r,\,\nu,\,\mu)$ to a partial differential equation with dependent variables $(s,\nu)$ (Mihalas, 1980). Time dependence can be neglected as $\gamma$-rays undergo few scatterings before the energy is deposited or the photon escapes to the observer.

From Eq. 2.1, our characteristic equations are

[TABLE]

We can now write the relativistic radiative transfer equation along a characteristic ray as

[TABLE]

where

[TABLE]

In order to solve Eq. 4, we use a backward differencing in frequency (i.e. $\partial\nu=\nu_{i-1}-\nu_{i}$ with $i$ denoting the current frequency). We then solve Eq. 4 by usual means for the formal solution along each ray for each frequency.

2.2 Opacities

Most nuclear decay lines in SNe have energies less than 3.5 MeV (see Table 1). For energies less than this, the dominant opacity source is Compton scattering and photoelectric absorption. Below 100 keV, the dominant opacity is photoelectric absorption and above that it is Compton scattering – (see figure 1 in Milne et al., 2004). This work only incorporates both photoelectric absorption and Compton scattering opacity. We neglect the influence of $e^{-}e^{+}$ pair production opacity because the typical decay $\gamma$-ray energies are less than 3.5 MeV in SNe.

For comparison with the MC method of Hillier & Dessart (2012) which follows that of Kasen et al. (2006), we use a photoelectric absorption opacity given by

[TABLE]

where $m_{e}$ is the electron mass, $\sigma_{\rm T}$ is the Thomson cross section, $\alpha$ is the fine structure constant, $N_{i}$ is the number density of species $i$, and $Z_{i}$ is the atomic number of species $i$. The Compton scattering opacity as given by eq. 7.113 of Pomraning (1973) is

[TABLE]

where $N_{e}$ is the number density of electrons, $r_{e}$ is the classical electron radius, and $x$ is $h\nu/m_{e}c^{2}$.

2.3 Emissivities

The total emissivity in the relativistic radiative transfer equation has two components. The first component is an isotropic prompt emission from nuclear decays, and the second is the scattering emissivity arising from Compton scattering.

2.3.1 Prompt Emission

The simpler of the two, the isotropic emissivity from the prompt decays is given by

[TABLE]

where $\Gamma=\frac{1}{2}[(\nu-\nu_{ij})/(V_{\rm Dop}\nu_{ij}/c)]^{2}$, $N_{i}$ is the number density of the $i$-th species isotope, $\tau_{i}=(t_{1/2})_{i}/\ln(2)$ is the nuclear decay constant for the $i$-th species isotope (see Table 1 for half-lives – $t_{1/2}$ – of ${}^{56}{\rm Ni}$ and ${}^{56}{\rm Co}$), $E_{ij}(\nu_{ij})$ is $j$-th line decay energy (frequency) for the $i$th species isotope, $P_{ij}$ is $j$-th line decay probability for the $i$-th species isotope, and $V_{\rm Dop}$ is the line Doppler velocity width ($\sim$100–200 km s*-1*).

The isotropic emission is the local source of $\gamma$-rays that eventually travel and scatter through the ejecta. Thus, it only needs to be calculated once before the transfer equation is solved.

2.3.2 Scattering Emissivity

Unlike the prompt emission, the scattering emissivity requires more numerical/computational effort and must be calculated concurrently while solving Eq. 4. The difficulty in calculating the scattering emissivity (Eq. 9) is due to the complicated angle and frequency dependence of the anisotropic Klein-Nishina (KN) scattering kernel (Eq. 11). Since we solve the specific intensity along characteristic rays for all impact parameters $p_{i}$, we have a fixed grid of polar angles $\theta_{i}$ (specifically $\mu_{i}=\cos\theta_{i}$) for every radius $r_{i}$ – note azimuthal symmetry is assumed. The scattering emissivity for an outgoing beam of frequency $\nu^{\prime}$ and direction $\bf\Omega^{\prime}$ is generally defined as

[TABLE]

where the prime denotes outgoing, and $\sigma_{s}(\nu\to\nu^{\prime},\xi)$ is the KN scattering kernel for a photon scattering with angle given as

[TABLE]

Following eq. 7.108 of Pomraning (1973), the KN scattering kernel for $x=h\nu/m_{e}c^{2}$ is given by

[TABLE]

Given that we assume $I_{\nu}\neq I_{\nu}(\phi)$, we need to integrate and remove the $\phi$ dependence in Eq. 9. Using the relationship $\delta(f(\phi))=\sum_{i}\delta(\phi-\phi_{i})/|f^{\prime}(\phi_{i})|$ for an arbitrary function $f$ with the zeros $\phi_{i}$, we can transform our $\delta$-function as

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

Both $\phi_{1}$ and $\phi_{2}$ exist in $\{0,2\pi\}$ since $\cos(\phi-\phi^{\prime})=\cos(2\pi-[\phi-\phi^{\prime}])$. Each delta contributes an equal value to the integral (see Eq. 9) with respect to $\phi$. Therefore, we have twice the integral of one delta function, giving us a factor of 2.

This transformation changes the $\mu$ integration limits to make sure the Compton relationship holds. To find the new $\mu$ limits, we extremize $\xi$ with respect to $\phi$ evaluated at our roots (i.e. $\partial\xi/\partial\phi|_{\phi=\phi_{i}}=0$). This gives us the constraint that $\phi_{i}-\phi^{\prime}=n\pi$, for an integer $n$. Using this result, we find our new limits on $\mu$ to be

[TABLE]

[TABLE]

[TABLE]

We can then rewrite Eq. 9 after integrating over $\phi$ as

[TABLE]

where

[TABLE]

Note that we have cancelled the one half with the factor of 2 from our $\phi$ integration. If we look at the $\mu$ integrand in Eq. 16 with our new $\mu$ integration limits, we run into a singularity at our limits. However, for this integral we can exploit Gauss-Chebyshev quadrature, which is defined as

[TABLE]

for abscissa $b_{i}=\cos[(2i-1)\pi/2n]$ and integer $n$. In order to get it into the form of Gauss-Chebyshev quadrature, we can make a linear transformation of $\mu$, namely, $w=c_{1}\mu+c_{2}$, for constants $c_{1}$ and $c_{2}$. These constants $c_{1}$ and $c_{2}$ are determined using the integration limits $\mu_{1,2}$ and solving the following linear equation:

[TABLE]

This has the solution

[TABLE]

From the definition of $\mu_{1,2}$ in Eq. 15 and our constants $c_{1,2}$ in Eq. 20, we find that our integrand transforms to

[TABLE]

Finally, Eq. 16 becomes

[TABLE]

This final result for the scattering emissivity is computationally favourable. We avoid having to loop through the large multi-dimensional arrays, thus saving time.

This transformation has certain limiting cases, such as when the $\mu^{\prime}=\pm 1$. In that case, we can look at the problem two ways. First we have that

[TABLE]

Thus, in these cases $I_{\nu}$ is a constant, and the sum in Eq. 22 equals $\pi I_{\nu}(r,\mu=\pm(1-1/x^{\prime}+1/x))$ – note there was a factor of 2 from the $\phi$ integration. The second way to understand these cases goes back to Eqs. 9 and 12. If we look at how the delta function transforms, we have

[TABLE]

In these cases, the $\phi$ integral is $2\pi$, and the $\mu$ integral picks out $\pi I_{\nu}(r,\mu=\pm(1-1/\gamma^{\prime}+1/\gamma))$. Both methods produce the same result.

Our work assumes that there is no contribution from the current frequency to the scattering emissivity (i.e. all photons are down-scattered). This assumption removes coupling between $I_{\nu}(r,\mu)$ and $\eta_{\nu}^{s}(r,\mu)$ at the current frequency $\nu$, so $I_{\nu}(r,\mu)$ can be solved exactly at each frequency when formally integrating Eq. 4.

2.4 Energy Deposition

Once we have solved $I_{\nu}(r,\mu)$ for all depths, we then calculate the energy deposited from scattering at each depth using the following relationship:

[TABLE]

where $E_{\rm lept}$ is the (assumed) local kinetic energy deposition from decay leptons (positrons and electrons). In Eq. 25, the physics we are capturing is the difference between the macroscopic energy lost from the specific intensity and the energy redistributed after scattering. In practice we would only integrate over the range of our frequency grid, which is chosen to cover the physics of the problem.

2.5 Implementation

As previously mentioned, this code is being implemented as part of cmfgen. The code solves Eq. 4 along characteristic rays for a given impact parameter ($p_{i}$) intersecting our radial grid ($r_{i}$). In this set-up we have $N_{D}$ radial grid points and $N_{C}$ core grid points, making $N_{P}=N_{C}+N_{D}$ impact parameters.

Since this code treats $\gamma$-rays from radioactive decays, we begin by reading in nuclear decay data such as nuclear decay energies and their decay probabilities for each unstable isotope included in an ejecta model. Lines with decay probabilities $<$1 per cent are not included in this code. However, like Hillier & Dessart (2012), we scale the decay line probabilities and decay lepton kinetic energies to conserve the total energy released during decay. Table 1 lists the following nuclear decay data: half-life, energy per decay, kinetic energies of leptons produced, and line energies and probabilities for the {}^{56}{\rm Ni}$$\rightarrow$${}^{56}{\rm Co}$$\rightarrow$${}^{56}{\rm Fe} decay chain, which dominates the decay energy for SNe. The annihilation line has a probability of 38 per cent because we assume that each positron produced (with 19 per cent intensity) annihilates without forming ortho-positronium after thermalization. We read in all other supernova data such as the mass fractions of all included species and count either the number of decays since the last time-step (an average) or the instantaneous decay.

After reading in all decay line data, we set up a frequency grid that is equally spaced in a log frequency scale for a given regime such as between lines, across the line, and two regimes for the red Compton tail. Each regime’s spacing is controlled by input parameters to give a desired spectral resolution. A finer frequency grid will produce “narrower” spectral line profiles by reducing numerical diffusion in frequency space as we propagate the photons spatially while solving Eq. 4. For a factor of roughly 3 less frequency points, the Gaussian profiles become broader by $\sim$50 per cent with no more than a few percent difference in the energy deposition.

Solving Eq. 4 for a given frequency, $k$, introduces difficulty given that the scattering emissivity is an integral over angles at a given depth – see Eqs. 9 and 11. No coupling between $\eta_{k}^{\rm s}$ and $I_{k}$ alleviates some computational difficulty. Full calculation of $\eta_{\nu^{\prime}}^{\rm s}(r,\mu^{\prime})$ requires integration over all angles for a given depth point, but integration of Eq. 4 along rays restricts us to a subset of angles for a given ray – see Section 2.1. It is necessary to map our intensity and emissivity arrays from $(z,p)\,\rightarrow\,(r,\mu)$ in order to perform all scattering calculations. With our assumption that all scattered photons are downgraded in frequency, we solve Eq. 4 exactly from blue to red frequencies. To do this, we calculate the scattering emissivity for all downgraded frequencies from $I_{k}$ and use central differencing quadrature ( $[\nu_{k-1}-\nu_{k+1}]/2$ ) to update $\eta_{j}^{\rm s}$, where $\nu_{j}<\nu_{k}$ – an implicit frequency integration of Eq. 22.

We interpolate $I_{k}$ (using monotonic cubic interpolation) onto a finer equally spaced linear $\mu$ grid in order to use Gauss-Chebyshev abscissa ($b_{i}=\cos\left[(2i-1)\pi/2n\right]$). An equally spaced linear grid allows us to quickly find and select the angle abscissa. After we update $\eta_{\nu^{\prime}}^{\rm s}$ for all possible down-scattered frequencies from $\nu_{k}$, we then map our arrays back into $(z,p)$ and solve for $I_{k+1}$ for all $p$.

In principle, the time required to calculate the scattering emissivity scales as $N_{D}\times N_{\nu}^{2}\times N_{\mu}^{2}$, where $N_{\nu}$ is the number of frequency points and $N_{\mu}$ is the maximum number of angle points equal to $2N_{P}-1$. However, we loop over down-scattered frequencies when calculating $\eta_{\nu^{\prime}}^{\rm s}(r,\mu^{\prime})$, so calculation time will scale less than $N_{D}\times N_{\nu}^{2}\times N_{\mu}^{2}$. Using Gauss-Chebyshev quadrature replaces one loop of length $N_{\mu}$ for a loop of length $N_{\rm Cheb}$. Interpolation onto a monotonic $\mu$ grid circumvents looping to find the abscissa in our arrays and makes its calculation time tractable.

Once we have calculated $I_{\nu}(z,p)$ for all frequencies and impact rays, we map it back into the $(r,\mu)$ space and from Eq. 25 calculate the energy deposited from $\gamma$-rays. This decay energy deposition will then be used and read in by cmfgen as a non-thermal heating source when solving the rate equations coupled to the relativistic radiative transfer equation for lower energy frequencies. We calculate an observer’s frame flux according to Hillier & Dessart (2012), which can be compared to observed $\gamma$-ray spectra of SNe.

3 Results

For this work, we recomputed model CHAN, a Chandrasekhar mass (MCh) WD with 0.62 $M_{\odot}$ of ${}^{56}{\rm Ni}$ initially, from Wilk et al. (2018) at two epochs, 17.4 days after explosion (roughly bolometric maximum) and 207 days after explosion (nebular time – optically thin to $\gamma$-rays). The initial ${}^{56}{\rm Ni}$ mass fraction at 0.75 d is shown in Fig. 1 for model CHAN. We performed the calculations described in this work and compared the results to two other methods cmfgen can use to calculate the energy deposition: (1) MC transport for $\gamma$-rays (Hillier & Dessart, 2012) using 8 000 000 decays and (2) a grey absorption approximation (Swartz et al., 1995) using $\kappa_{\gamma}=0.06Y_{e}$ cm2 g*-1*.

3.1 Runtime

The runtime scaling of our $\gamma$-ray transport code with the number of depth points, number of angles, and the number of frequency points is explained in Section 2.5. To improve efficiency we have made this code parallelizable over depth and have tested it using an Intel(®) Xeon(®) CPU E5-4610 2.40GHz processor and 8 cores. Tests on more modern processors (like Intel(®) Xeon(®) CPU E5-2620 v4 2.10GHz) show an improvement of a factor of two in speed (for the same number of cores). All calculations were performed with $N_{D}=109,\;N_{C}=15,and\;N_{P}=124$. A calculation with a very fine frequency grid resolution ($\sim$26 600 frequency points), needed if very accurate line profiles are to be computed, has a runtime of $\sim$17 hours. However, for most work we are only interested in the energy deposition rate, and we can use a much lower spectral resolution. For example, for a low spectral resolution calculation with $\sim$6 500 frequency grid points the code’s runtime is approximately 45 minutes (the runtime scales roughly as the number of frequency points squared). Even though the number of frequency points has been reduced by a factor of $\sim$4, the energy deposited throughout the ejecta differs by at most 1 per cent from the high resolution calculation.

The runtime on the same machine for a MC calculation on a single processor with 8 000 000 decays (necessary for low statistical noise) is significantly longer than our low spectral resolution calculation. In this case, the MC runtime is roughly 10 hours. Using a factor of 10 less decays per species, the MC calculation runtime is roughly 1 hour. For lower resolution calculations, the runtimes of both codes are somewhat comparable.

3.2 Energy Deposition

Fig. 2 compares the ratio of the energy deposition at 17.4 days calculated using our new radiative transfer code to that computed with the MC code, and to that obtained using the grey absorption approximation, as a function of velocity for a MCh WD. Fig. 2 shows that our work is in agreement with the MC method within 3 per cent at $<$20 000 km s*-1*. Below 3000 km s*-1*, the MC method is subject to statistical noise and has discrepancies with this work due to a “${}^{56}{\rm Ni}$ hole” where little radioactive material is mixed in. Beyond 20 000 km s*-1*, MC statistical noise is the source of the discrepancy between the two codes. The error between 5 000 and 11 000 km s*-1* is partially numerical since doubling the value of $N_{D}$ caused the error in this region to decrease. However, doubling the value of $N_{\nu}$ gave minimal improvement. The error is also likely sensitive to the interpolation techniques. However, despite our best efforts, we were unable to reduce the discrepancy below 1 per cent.

Table 2 lists the total integrated energy deposition over the whole ejecta at this epoch and shows that the two methods agree within $\sim$2.5 per cent. Fig. 2 also shows the ratio of the non-thermal energy deposited to that of local energy released from nuclear decays. We see that beyond 12 000 km s*-1*, the energy deposition comes from the inner ejecta as the $\gamma$-ray photons scatter. In the region between 12 000–20 000 km s*-1* where many optical and diagnostic lines are formed, this work is consistent within 2.5–3 per cent to that of Hillier & Dessart (2012). In the same region, the energy deposition computed using the grey approximation diverges from the other methods.

Fig. 3 is the same as Fig. 2, except now at 207 days and without the ratio of the energy deposition to the local energy emitted being plotted. Despite the fundamental differences in the approach each code uses to calculate the energy deposition, the MC code and our work agree to within 1 per cent (highlighted in Table 2). At late times and for ejecta velocities less than 10 000 km s*-1*, the grey approximation is inconsistent with the two other methods by more than $\sim$5 per cent. At this epoch, important strong cooling lines form at velocities $\leq$10 000 km s*-1*, so wrongly estimating the energy deposited may affect the ionization structure and/or flux in strong cooling lines.

3.3 Synthetic Spectra

From the CMF at the outer boundary, we can transform the specific intensity into the observer’s frame to produce a synthetic $\gamma$-ray spectrum – see section 11 of Hillier & Dessart (2012). Fig. 4 shows our resulting synthetic spectra calculated at two epochs, 17.4 and 207 days post-explosion. At 17.4 days, the dominant decay luminosity begins to switch from ${}^{56}{\rm Ni}$ to ${}^{56}{\rm Co}$, and the spectrum shows strong lines from both ${}^{56}{\rm Ni}$ and ${}^{56}{\rm Co}$ – see Table 1. However, at 207 days, all the ${}^{56}{\rm Ni}$ has decayed and the synthetic spectra is dominated by ${}^{56}{\rm Co}$ decay lines. At both epochs, synthetic spectra from our work and the MC method are in good agreement. The total integrated flux listed in Table 2 shows that the two methods are within $\sim$9 per cent at 17.4 days and $\sim$4 per cent at 207 days.

Both the MC method and our radiative transfer code produce synthetic spectra with predicted asymmetric profiles with absorption on the red side of the emission line. These asymmetric profiles are not uncommon. They are predicted and seen in X-ray line profiles for massive stars (Macfarlane et al., 1991; Owocki & Cohen, 2001; Cohen et al., 2010, 2014). They are also a product of dust scattering (Romanik & Leung, 1981), and have been modelled for dust in the ejecta of SN1987A (Bevan & Barlow, 2016). Electron scattering opacity also produces blue shifted asymmetric profiles for some optical lines in Type II SNe (Dessart & Hillier, 2005). Many previous theoretical studies have predicted the anticipated asymmetric $\gamma$-ray line profiles, notably Burrows & The (1990); Mueller et al. (1991); Hoeflich et al. (1992, 1993, 1994); Maeda (2006). Since Compton scattering is a continuum opacity, the optical depth of the red side of the line is higher because the path length is larger to the far side of the ejecta. We expect our profiles to exhibit the same effect. In Fig. 5 we highlight two ${}^{56}{\rm Co}$ decay lines at 1038 and 1238 keV. Fig. 5 shows that at 17.4 days our profiles are asymmetric as the optical depth to $\gamma$-rays is large, causing most of the emission to be in the blue side of the line profile; whereas, at 207 days its optical depth is low, and the profile is symmetric.

As can be seen from Figs. 4 and 5, the profiles produced by the MC calculation are somewhat narrower than those produced by our $\gamma$-ray transfer approach. This arises from numerical diffusion as we propagate photons from the inner regions to the outer boundary of the model (as the calculation is done in the CMF, the photons are propagated in both frequency and space). Numerical diffusion can be reduced by increasing the grid resolution, or by reducing the extent of the outer boundary (especially relevant at 207 days). The best approach would be to utilise the computed scattering emissivities in an observer’s frame calculation, but given the lack of high quality observed data we have not implemented such a procedure.

3.4 Comparison to SN2014J

To compare our work to the observations from SN 2014J, we computed $\gamma$-ray synthetic spectra using our ejecta model at 75 days. Fig. 6 shows the 847 and 1238 keV line profiles as a function of velocity, comparing synthetic spectra computed from this work and the MC code at 75 days. The 847 keV line centroid velocities are $-$1836 and $-$2158 km s*-1* for this work and the MC method respectively. Similarly, the 1238 keV line centroid velocities are $-$960 and $-$1348 km s*-1*. These 847 keV line results are consistent with the values measured for the $\gamma$-ray spectrum obtained for SN2014J. Churazov et al. (2014, in fig. 4 and table 1) shows that the 847 keV cobalt line is slightly blueshifted with a velocity of $-$1900$\pm$1600 km s*-1*. However, our work disagrees with the 1238 keV cobalt line. Table 1 of Churazov et al. (2014) shows this line to have a peak velocity shift of $-$4300$\pm$1600 km s*-1*. This line should have a smaller blueshifted velocity relative to the 847 keV line since the optical depth is lower due to the smaller cross section at 1238 keV. Given the very large errors on the mean shifts, the disagreement may simply be statistical. The fiducial model plotted in fig 4 of Churazov et al. (2014) shows a less blueshifted 1238 keV line profile.

Not only do our profiles agree with those measured by Churazov et al. (2014), but our flux levels also agree. Adjusting our flux at 75 d in Fig. 6 for a distance of 3.5 Mpc to M82, our flux levels are roughly $\sim 9\times 10^{-6}$ photons cm*-2* s*-1* keV*-1* for our 847 keV line, consistent with the flux levels shown in fig 1 of Churazov et al. (2014).

3.5 Grey transfer

Since SN radiative transfer calculations are already time-intensive, it is beneficial to use a simple and fast prescription to calculate the energy deposited by $\gamma$-rays in the ejecta. The grey absorption method of Swartz et al. (1995) (hereafter S95) is one such fast procedure to calculate the energy deposition. Comparing the results of both the MC radiative transfer and this work from Figs. 2 and 3, we see that the grey approximation of S95 would require a time varying grey opacity factor as well as one that varies spatially. Simply using the mass absorption coefficient, $\kappa_{\gamma}=0.06Y_{e}$ cm2 g*-1*, does not reproduce the energy deposition the other methods produce.

Fig. 7 shows the ratio of the calculated energy deposition of this work to that calculated using the grey transfer from S95. However, we show the energy deposition ratio for varying coefficients for $\kappa_{\gamma}$ at 17.4 and 207 days in the grey approximation. At 17.4 days, we see that $\kappa_{\gamma}=0.07Y_{e}$ cm2 g*-1* matches to our work below 10 000 km s*-1*, while $\kappa_{\gamma}=0.09Y_{e}$ cm2 g*-1* more accurately reproduces the energy deposition beyond 10 000 km s*-1*. For low values of the grey absorption (i.e. 0.05$Y_{e}$ cm2 g*-1*) too little energy is deposited in the inner ejecta, which is instead deposited in the outer region. However, increasing the constant in the grey absorption coefficient still produces too much absorption in the outer ejecta.

At nebular times, we see that a lower value of $\kappa_{\gamma}=0.05Y_{e}$ cm2 g*-1* reproduces the energy deposition of the other methods when the ejecta is optically thin to $\gamma$-rays. However, Fig. 7 shows too much energy being deposited into the outer ejecta beyond 12 000 km s*-1* at nebular times. More sophisticated approaches like that of Jeffery (1998) may be needed to model different parts of the ejecta as a function of time.

These values of $\kappa_{\gamma}$ required to reproduce the energy deposition are roughly aligned with those of Maeda (2006). Maeda (2006) argue that $\kappa_{\gamma}=0.027$ cm2 g*-1* best reproduces a light curve of their spherically symmetric F model. With $Y_{e}\approx 0.5$, the result of Maeda (2006) is consistent with the value $\kappa_{\gamma}=0.05Y_{e}$ cm2 g*-1* that we claim agrees with our nebular energy deposition. However, our work demonstrates that values of much higher $\kappa_{\gamma}$ are required to reproduce the energy deposition at early times (see Fig. 7). Woosley et al. (1986) claim a higher value (because of the centrally located distribution of ${}^{56}{\rm Ni}$ and sensitivity to angle averaging effects along density gradients) of $\kappa_{\gamma}=0.07$ cm2 g*-1* reproduces the energy deposition function.

4 Conclusion

We have presented a new code that solves the relativistic radiative transfer equation for $\gamma$-rays, taking into account opacity, prompt radioactive decay emissivity, and scattering emissivity. In computing the scattering emissivity, we assume that all photons are downgraded in energy and ignore any thermal redistribution effects since the expansion velocities dominate the transfer. From the specific intensity, we are able to produce an observer’s frame spectrum as well as the energy deposition consistent with that of the MC code of Hillier & Dessart (2012).

For a low spectral resolution ($\sim$6 500 frequency grid points) calculation, our new code has the advantage of running in approximately 45 minutes using parallel processing with 8 CPUs. Low resolution calculations result in at most 1 per cent differences in calculated energy deposition within the ejecta compared to the higher spectral resolution. In comparison to the MC code, with 8 000 000 decays needed to achieve low statistical noise, the code runs in approximately 10 hours on the same machine using one CPU. In terms of the integrated energy deposition, the two codes agree within 3 per cent at early times and within 1 per cent at late times.

We have shown that this code produces the expected line profiles. When the optical depth to $\gamma$-rays is large, the red side of the line has a higher optical depth than the blue side, and thus most of the emission comes out in the blue side of the line profile – see Figs. 4, 5, 6, and Churazov et al. (2014, fig. 4 and table 1).

This code will be publicly available and serves (along with all other MC $\gamma$-ray radiative transfer codes) to improve the astrophysics community’s constraints on nucleosynthetic yields as well as the stratification of nuclear material in SN ejecta. We are currently limited by observations of $\gamma$-rays from SNe, so future observations of $\gamma$-rays will uncover a previously untapped opportunity to understand more about the nature of SNe and their progenitors.

acknowledgements

DJH acknowledges partial support from STScI theory grant HST-AR-12640.001-A, and DJH and KDW thank NASA for partial support through theory grant NNX14AB41G.