An Integral-Based Technique (IBT) to Accelerate the Monte-Carlo Radiative Transfer Computation for Supernovae (2024)

Xingzhuo ChenGeorge P. and Cynthia Woods Mitchell Institute for Fundamental Physics & Astronomy,
Texas A. & M. University, Department of Physics and Astronomy, 4242 TAMU, College Station, TX 77843, USALifan WangGeorge P. and Cynthia Woods Mitchell Institute for Fundamental Physics & Astronomy,
Texas A. & M. University, Department of Physics and Astronomy, 4242 TAMU, College Station, TX 77843, USADaniel KasenAstronomy Department and Theoretical Astrophysics Center, University of California, Berkeley, Berkeley, CA 94720, USA

Abstract

We present an integral-based technique (IBT) algorithm to accelerate supernova (SN) radiative transfer calculations.The algorithm utilizes “integral packets”, which are calculated by the path integral of the Monte-Carlo energy packets, to synthesize the observed spectropolarimetric signal at a given viewing direction in a 3-D time-dependent radiative transfer program.Compared to the event-based technique (EBT) proposed by Bulla etal. (2015), our algorithm significantly reduces the computation time and increases the Monte-Carlo signal-to-noise ratio.Using a 1-D spherical symmetric type Ia supernova (SN Ia) ejecta model DDC10 and its derived 3-D model, the IBT algorithm has successfully passed the verification of: (1) spherical symmetry; (2) mirror symmetry; (3) cross comparison on a 3-D SN model with direct-counting technique (DCT) and EBT.Notably, with our algorithm implemented in the 3-D Monte-Carlo radiative transfer code SEDONA, the computation time is faster than EBT by a factor of $10-30$ , and the signal-to-noise (S/N) ratio is better by a factor of $5-10$ , with the same number of Monte-Carlo quanta.

Supernova

^†^†software: SEDONA(Kasen etal., 2006)

1 Introduction

Type Ia supernovae (SNe Ia) have been widely applied in cosmological studies (e.g. Abbott etal. (2019); Brout etal. (2022)) based on the empirical relation between the light curve properties and the luminosities (e.g. Phillips (1993)).However, the explosion mechanism of SNe Ia is still a mystery; studies of the ejecta structure rely on simulations that include hydrodynamics, nucleosynthesis, and radiative transfer.

Radiative transfer simulations of supernovae calculate the specific intensity and plasma excitation states to obtain spectra that can be compared to optical spectroscopic and photometric observations.Many radiative transfer codes assume that SNe Ia are spherical symmetric and reasonable agreement have been found among the simulation results of different codes (Blondin etal., 2022).

However, imaging of nearby SN Ia remnants (e.g. Ferrazzoli etal. (2023)) and hydrodynamic simulations (e.g. Gamezo etal. (2004); Kromer etal. (2013)) suggest that SNe Ia have an asymmetric 3-D structure.Given the distance to most SNe Ia, this 3-D structure cannot be resolved, but can be probed by spectropolarimetric observations (e.g., Wang etal., 2003; Cikota etal., 2019; Yang etal., 2022).The polarization from SN Ia is primarily due to Thomson scattering of photons, and a non-zero linear polarization signal can arise in asymmetric SN ejecta due to the incomplete cancellation of polarization from photons scattered from different geometrical structures (Wang & Wheeler, 2008).

Monte-Carlo methods, which use energy packets to emulate the propagation of photons in the SN plasma, provide a flexible way to treat complicated physical processes (Lucy, 2002, 2003) and have been applied in 3-D time-dependent radiative transfer codes such as SEDONA (Kasen etal., 2006), and ARTIS (Kromer & Sim, 2009).3-D Monte-Carlo radiative transfer (MCRT) simulations can be computationally expensive, as many photon packets must be propagated to reduce the statistical noise in the synthesized observables.This especially poses a challenge for modeling SN polarization, as the expected signals are typically only at the 1% level or less.To reduce the noise in SN MCRT simulations, Bulla etal. (2015) introduced an “event-based technique” (EBT) for spectropolarimetry synthesis, which has been applied to simulations of SNe Ia (Bulla etal., 2016a, b), superluminous SNe (Inserra etal., 2016), and kilonovae (Bulla etal., 2021).

In this paper, we introduce a new integral-based technique (IBT), which is implemented in the MCRT code SEDONA,to efficiently construct observable spectra from a 3-D MCRT computation.We find that IBT is $10-30$ times faster and $5-10$ times less Monte-Carlo noise than EBT, with the same number of Monte-Carlo quanta.Moreover, IBT can be applied to the ultraviolet and infrared wavelength, and early phase of SNe, which was not calculated with EBT (Bulla etal., 2015).Section 2 reviews the direct-counting technique (DCT) and event-based technique (EBT) that have been used in previous MCRT codes.Section 3 presents the IBT method.Section 4 uses several toy models to verify the spectropolarimetry results from IBT.Section 5 compares the Monte-Carlo simulation error and computation time of the three algorithms: DCT, EBT, and IBT.Section 6 summarizes our main conclusions.

2 Polarized Spectrum Extraction

In the SEDONA code,an optical photon packet (OP) represents a collection of photons with the same frequency, spatial coordinate, and propagation direction.The OPs undergo interaction events due to Thomson scattering, bound-bound transitions, bound-free transitions, and free-free transitions which can change their frequency and direction.

The polarization information of an OP is stored as a dimensionless Stokes vector

\boldsymbol{s}_{OP}=\begin{pmatrix}1\\q\\u\\v\end{pmatrix}\ ,

(1)

and the Stokes vector can be constructed with $\boldsymbol{s}_{OP}$ and the energy of the OP $E_{OP}$ .

Typically, OPs are generated with zero initial polarization as expected for thermal emission.During the propagation of an OP, the linear polarization state changes in a Thomson scattering event according to the Rayleigh scattering phase matrix (Section 1.17, equation 217 of (Chandrasekhar, 1960)).In bound-bound, bound-free, and free-free interactions, the OPs are typically assumed to be fully depolarized.The current version of SEDONA does not include magnetic fields, and the $V$ parameter in the Stokes vector, which represents the circular polarization, is set to zero.

2.1 Direct Counting Technique

A straightforward method of spectral synthesis in MCRT is the direct counting technique (DCT), which sums the energy of OPs that escape the simulation domain within a certain frequency, time, and viewing angle bin.The observed flux is then calculated as

\begin{pmatrix}\mathcal{I}\\\mathcal{Q}\\\mathcal{U}\\\mathcal{V}\end{pmatrix}=\frac{1}{4\pi r^{2}\Delta t\Delta\nu\Delta\Omega}\sum%_{i}E_{OP,i}\boldsymbol{s}_{OP,i}

(2)

where $\boldsymbol{s}_{OP,i}$ is dimensionless Stokes vector of the $i$ -th OP, $E_{OP,i}$ is the lab frame energy of the $i$ -th OP, $\Delta t$ is the size of the arrival time bin, $\Delta\nu$ is the size of the frequency bin, $r$ is the distance to the SN center, $\Delta\Omega$ is the viewing direction bin size.The summation is over all the OPs in the arrival time interval $[t-\Delta t/2,t+\Delta t/2)$ , the frequency interval $[\nu-\Delta\nu/2,\nu+\Delta\nu/2)$ , and the viewing direction interval $[\varphi-\Delta\varphi/2,\varphi+\Delta\varphi/2)$ ; $[\theta-\Delta\theta/2,\theta+\Delta\theta/2)$ .

In the DCT method, the signal-to-noise ratio (S/N) values of the synthetic spectra depend on the number of OPs simulated and the choice of the bins of arrival time, frequency, and viewing direction.Therefore, a higher resolution in viewing angle results in fewer OPs per bin and a lower S/N. The DCT has been the default implementation in the SEDONA code.

2.2 Virtual Packets

Several techniques have been suggested to improve the S/N in MCRT simulations.The idea of a “virtual packet” (VP) in the SN simulation context is discussed by Kerzendorf & Sim (2014) and implemented in the 3-D simulation by Bulla etal. (2015).Generally, VPs are used as follows

•
An observer viewing direction is specified at the start of the MCRT simulation.
•
VPs are created in the SN plasma to represent the emitted and the scattered photons at the corresponding volume-time-frequency bin pointing to the viewing direction.
•
The energy of each VP is reduced by the optical depth along the path to the observer.
•
Finally, the escaped VPs are used to synthesize the spectrum for the specified observer viewing angle.

Event-based techinque (EBT), proposed by Bulla etal. (2015), utilizes VPs to increase the S/N relative to DCT.In the EBT calculation, when an OP has undergone a physical event, a VP is created at the same coordinate and time with the same co-moving frame frequency.If the interaction event is Thomson scattering, then the co-moving frame Stokes vector of the VP is calculated from the Rayleigh scattering rule (Chandrasekhar, 1960):

\bar{\boldsymbol{S}}_{VP}=\bar{\boldsymbol{S}}_{OP}\boldsymbol{P}(\bar{\theta}%_{OP},\bar{\varphi}_{OP},\bar{\theta},\bar{\varphi})\ ,

(3)

where $\bar{\boldsymbol{S}}_{OP}$ is the Stokes vector of the OP in the co-moving frame, $(\bar{\theta}_{OP},\bar{\varphi}_{OP})$ is the OP propagating direction in the co-moving frame, $(\bar{\theta},\bar{\varphi})$ is the viewing direction in the co-moving frame, $\boldsymbol{P}$ is the combination of the rotation matrix and the scattering matrix (See Equation 10 in Bulla etal. (2015)).In the following, co-moving frame quantities are denoted with a bar, whereas non-bar quantities refer to the lab frame.A detailed expression of the $\boldsymbol{P}$ matrix is given in Appendix A.

If the interaction event is a bound-bound transition, bound-free transition, or free-free transition, then an depolarized VP propagating towards the viewing direction is created.The energy of the VP is

\bar{E}_{VP}=\frac{1}{4\pi}\bar{E}_{OP}\ .

(4)

After creation, the VPs propagate along the viewing direction with the speed of light, and are not considered in the subsequent computation of the plasma excitation and ionization states.The total optical depth along the VP trajectory is integrated as

\tau_{tot}=\sum_{j}d_{j}\alpha_{tot,j}\ ,

(5)

where $d_{j}$ is the length of the $j$ -th line segment of the VP trajectory, $\alpha_{tot,j}$ is the total extinction coefficient at the $j$ -th line segment.The lab frame extinction coefficient is calculated from the co-moving frame value

3 The Integral-Based Technique

Our proposed IBT uses integrated packets (IPs), which are an improved version of VPs that can be used to more efficiently construct spectra.Compared to VP, the major changes of IP are:(1) A VP stores the energy and polarization at a single frequency with a four-vector $\boldsymbol{S}_{VP}$ , while an IP stores the energy and polarization at multiple frequency bins with a $4\times N$ tensor of cell radiance $\boldsymbol{R}([I,Q,U,V];[\nu_{1},...,\nu_{N}])$ ;(2) The creation of a VP is based on the physical interactions of OPs, whereas an IP is created in each volume cell for each time step.

3.1 The Integrated Packets

The IP computation is based on a universal logarithmically spaced frequency grid $[\nu_{0},\nu_{1}...\nu_{N}]$ .The frequency ratio between the adjacent pixels is a constant $C$ , satisfying

\nu_{n+1}/\nu_{n}=C.

(8)

Using this frequency grid, the Doppler shift of the variables could be simplified to the shifting of the index of the frequency bin, which saves computation time.

We define the cell radiance $\boldsymbol{R}(\vec{x},t,\nu)$ in the lab frame for a cell located at coordinate $\vec{x}$ and time $t$ as the energy radiated per frequency bin at frequency $\nu$ toward the viewing direction, the units of $\boldsymbol{R}(\vec{x},t,\nu)$ is $\rm erg\ Hz^{-1}\ sr^{-1}$ .The formal solution of the cell radiance is:

	$\displaystyle\boldsymbol{R}(\vec{x},t,\nu)=\Delta v\Delta t_{s}\left(\frac{%\bar{\nu}}{\nu}\right)^{-2}\left(\bar{\boldsymbol{j}}_{emi}+\bar{\boldsymbol{j%}}_{sc}\right)$		(9)
	$\displaystyle\bar{\boldsymbol{j}}_{emi}=\bar{j}_{emi}(\vec{x},t,\bar{\nu})%\begin{pmatrix}1\\0\\0\\0\end{pmatrix}$		(10)
	$\displaystyle\bar{\boldsymbol{j}}_{sc}=\alpha_{T}\int_{\Omega}\bar{\boldsymbol%{I}}(\vec{x},t,\bar{\theta}_{in},\bar{\varphi}_{in},\bar{\nu})\boldsymbol{P}(%\bar{\theta}_{in},\bar{\varphi}_{in},\bar{\theta},\bar{\varphi})d\Omega,$		(11)

where $\bar{\boldsymbol{I}}$ is the specific intensity as a four-vector of the Stokes parameters; $\Delta t_{s}$ is the size of the simulation time step; $\bar{\boldsymbol{j}}_{emi}$ and $\bar{\boldsymbol{j}}_{sc}$ are the four-vectors of the emission terms due to depolarized isotropic emission and scattered light from Thomson scattering, respectively;the $\bar{j}_{emi}$ function is the depolarized isotropic emission term, which includes all the emission from bound-bound, bound-free, and free-free transitions;the $\boldsymbol{P}(\bar{\theta}_{in},\bar{\varphi}_{in},\bar{\theta},\bar{\varphi})$ function is a combination of the Rayleigh scattering phase matrix and rotation matrix discussed in Appendix A; $\alpha_{T}$ is the extinction coefficient for Thomson scattering;the integral is over all the incident directions ( $\bar{\theta}_{in},\bar{\varphi}_{in}$ ); the $(\bar{\nu}/\nu)^{-2}$ term is the correction of the Doppler effect (Castor e.g., 1972, eq.(1–3); see also Mihalas 1978, p.31,33,495–496).

For each IP, the cell radiance $\boldsymbol{R}$ is calculated on the universal frequency grid and stored as a $4\times N$ vector.During the propagation of an IP, $\boldsymbol{R}$ is reduced by the extinction coefficient on the trajectory.The spectral synthesis is the summation of $\boldsymbol{R}$ over the IPs at a specific arrival time.Section 3.2 and Section 3.3 illustrate the computation of $\boldsymbol{R}$ in a MCRT simulation during the creation of an IP.Section 3.4 illustrates the update of $\boldsymbol{R}$ during the propagation of IP and spectral synthesis.

3.2 Isotropic Emission

The isotropic emission term, including bound-bound, bound-free, and free-free interactions, can be calculated by the summation of the OPs from the same physical events

\Delta v\Delta t_{s}\left(\frac{\bar{\nu}}{\nu}\right)^{-2}\bar{j}_{emi}(\vec{%x},t,\bar{\nu})=\left(\frac{\bar{\nu}}{\nu}\right)^{-2}\sum_{i}\frac{\bar{E}_{%OP,i}}{4\pi\Delta\bar{\nu}}\ ,

(12)

where the summation is over all the OPs within the frequency bin $[\bar{\nu}-\Delta\bar{\nu}/2,\bar{\nu}+\Delta\bar{\nu}/2)$ , the volume cell $\Delta v$ , and the time step bin $[t_{s},t_{s}+\Delta t_{s})$ .Only OPs that are newly created or undergone the same physical events are taken into account.

In the present SEDONA calculation, we include bound-free and free-free transitions and adopt the line expansion opacity approximation (Kasen etal. 2006, eq.(8)) for bound-bound transitions which are treated as purely absorptive.Therefore, the isotropic emission term is calculated from the current plasma state as $\bar{j}_{emi}(\vec{x},t,\bar{\nu})=B(\bar{\nu})\bar{\alpha}_{abs}(\vec{x},t,%\bar{\nu})$ , where $B(\bar{\nu})$ is the blackbody spectrum and $\bar{\alpha}_{abs}$ is the total extinction coefficient including bound-bound, bound-free, and free-free transitions (Thomson scattering is not included in this term).

3.3 The Path Integral

The cell radiance from Thomson scattering is calculated from the path integral of the OPs

\Delta v\Delta t_{s}\left(\frac{\bar{\nu}}{\nu}\right)^{-2}\bar{\boldsymbol{j}%}_{sc}=\frac{\sum_{i,j}\bar{E}_{i}\bar{d}_{i,j}\bar{\alpha}_{T}\bar{%\boldsymbol{s}}_{i,j}\boldsymbol{P}(\bar{\theta}_{i,j},\bar{\varphi}_{i,j},%\bar{\theta},\bar{\varphi})}{4\pi\Delta\bar{\nu}}

(13)

where $\bar{d}_{i,j}$ is the co-moving frame length of the $j$ -th line segment of the $i$ -th OP trajectory in the volume cell and in the time step, the summation is over all the trajectories within the frequency bin $[\bar{\nu}-\Delta\bar{\nu}/2,\bar{\nu}+\Delta\bar{\nu}/2)$ , the volume cell $\Delta v$ , the time step bin $[t_{s},t_{s}+\Delta t_{s})$ .

For each volume cell and at each time step, an IP is created.The initial coordinate of the IP is randomly chosen in the volume cell, and the initial time of the IP is randomly chosen in the time step.The cell radiance $\boldsymbol{R}$ is calculated on the universal frequency grid and stored in IP as a $4\times N$ tensor.

3.4 The Propagation of Integral Packets

After creation, the IPs propagate toward the viewer with the speed of light, and the cell radiance $\boldsymbol{R}$ is reduced during propagation.Similar to the propagation of the VPs, the trajectory of an IP is split into line segments by the edges of time steps and volume cells.Moreover, the trajectory is further split to satisfy the criterion:

\left|\frac{\bar{\nu}_{begin}-\bar{\nu}_{end}}{\bar{\nu}_{end}}\right|\leq C-1,\

(14)

where $\bar{\nu}_{begin}$ is the co-moving frame frequency at the beginning of the line segment, $\bar{\nu}_{end}$ is the co-moving frame frequency of the same lab frame frequency $\nu$ at the ending of the line segment, the constant $C$ is defined in Equation 8.

In the calculation of VP in Equation 5 and 7, the total opacity is accumulated over the trajectory, which computationally requires one more double precision floating point number per VP to store the opacity.This process reduces the computation time on exponentials at the cost of memory space.While in the propagation of IPs, the cell radiance is updated at every line segment on the trajectory without accumulating the total opacity over the trajectory

\boldsymbol{R}_{i,j+1}(\nu)=\boldsymbol{R}_{i,j}(\nu)e^{-d_{i,j}\alpha_{tot,i,%j}(\nu)}\ ,

(15)

where $d_{i,j}$ is the length of the $i$ -th IP at the $j$ -th line segment over the trajectory, $\alpha_{tot,i,j}(\nu)$ is the total extinction coefficient of the $i$ -th IP at the $j$ -th line segment and at lab frame frequency $\nu$ .The final spectrum is the sum of the IPs:

\begin{pmatrix}\mathcal{I}\\\mathcal{Q}\\\mathcal{U}\\\mathcal{V}\end{pmatrix}=\frac{1}{4\pi r^{2}\Delta t}\sum_{i}\boldsymbol{R}_{i%,final}\ ,

(16)

where $\boldsymbol{R}_{i,final}$ is the final cell radiance when the $i$ -th IP escapes the simulation domain.The summation is over all IPs in the arrival time interval $[t-\Delta t/2,t+\Delta t/2)$ .Note that the grid size of arrival time $\Delta t$ should be much larger than the time step size $\Delta t_{s}$ to avoid grid mismatch issue, while this issue is not significant in DCT or EBT calculations.

4 Model Validation

An Integral-Based Technique (IBT) to Accelerate the Monte-Carlo Radiative Transfer Computation for Supernovae (1)

In this section, we prepare two toy models to validate the IBT computation results.The first model is the SNIa delayed-detonation model DDC10.DDC10 is calculated from hydrodynamic simulations in Blondin etal. (2013) and adopted as a benchmark in Blondin etal. (2022) to evaluate the differences between radiative transfer codes.We remap the 1-D model into a $60\times 60\times 60$ 3-D Cartesian grid.retaining the spherical symmetry and with the outer boundary velocity limited to 25,000 km/s.Figure 1 shows the abundance of several elements and the density profile of the DDC10 model.The second model is an ellipsoidal model based on the DDC10 model; the y and z coordinates are reprojected with $y_{new}=1.3y$ and $z_{new}=z/1.3$ , in order to introduce asphericity.This model is also realized on a $60\times 60\times 60$ Cartesian grid, with a boundary velocity of 32,500 km/s.The geometry of the model is shown in the right panel of Figure 1.

In the following text, the spherical symmetric DDC10 model is denoted as D1D, and the modified ellipsoidal 3-D model is denoted as D3D.Both the models are treated with 3-D time-dependent radiative transfer calculations in SEDONA starting from 0.976 days after the explosion.Before 6.6 days, the simulation time step is a logarithmic with $\Delta t_{s}/t_{s}=0.03$ .After 6.6 days, the simulation time step is $\Delta t_{s}=0.2$ days.The arrival time bin size is $\Delta t=1$ day, which satisfies $\Delta t\geq 5\Delta t_{s}$ .The frequency grid ranges from $1\times 10^{14}$ Hz to $5\times 10^{15}$ Hz (30,000 $\AA$ to 600 $\AA$ ) with $\Delta\nu/\nu=0.002$ .

4.1 Spherical Symmetry

Because the linear polarization components cancel out, a spherical symmetric model should produce zero polarization.Using IBT, we calculate the spectropolarimetry of the D1D model to test if the algorithm recovers zero polarization.Figure 2 shows the spectra and the polarization percentage at different times after the explosion.At each time step, $6\times 10^{6}$ OPs are generated to represent the energy release from radioactive decay and gamma-ray scattering.We find that the polarization signal at all times and wavelengths is consistent with the expected zero polarization.Moreover, the simulation error between 7 days and 31 days and between 2000 and 10000 $\AA$ wavelength range (within which most SNe Ia spectropolarimetry observations are made (Cikota etal., 2019)) is as low as $\sim 0.2\%$ .Notably, the spectra in the ultraviolet wavelength between 7 days and 21 days also show high S/N, despite the emergent flux being low in this band.In the early phase of SNe Ia, most of the OPs are generated below the photosphere, therefore the cell radiance $\boldsymbol{R}$ above the photosphere has low S/N.As a result, the spectrum shows low S/N at 3-4 days after the explosion.Moreover, the S/N in the infrared wavelength is lower than that in the optical wavelength, because most of the OPs are in the optical wavelength.Due to the limited number of OPs in the ultraviolet wavelength at late phase, the ultraviolet spectrum at 30-31 days and after is dominated by Monte-Carlo noise.In the IBT simulation of this paper, this phenomenon is only observed below 2000 $\rm\AA$ and can be alleviated by increasing the number of OPs.

An Integral-Based Technique (IBT) to Accelerate the Monte-Carlo Radiative Transfer Computation for Supernovae (2)

An Integral-Based Technique (IBT) to Accelerate the Monte-Carlo Radiative Transfer Computation for Supernovae (3)

4.2 Mirror Symmetry

In this section, we use IBT to calculate the spectropolarimetry at two mirror symmetric viewing directions of the D3D model.At each time step we generate $1\times 10^{7}$ OPs.Based on the mirror symmetry of the D3D model in the X-Z, Y-Z, and X-Y planes, spectropolarimetry in the viewing directions $(cos(\theta)=0.51122,\ \varphi=\pi/6)$ and $(cos(\theta)=0.51122,\ \varphi=11\pi/6)$ should have the same Stokes $\mathcal{Q}$ terms and inverted $\mathcal{U}$ terms.Figure 3 shows the spectropolarimetry at the two symmetric viewing directions.We notice the simulation results are consistent with the theoretical expectations of the symmetry at above 2000 $\rm\AA$ wavelength from early phase (7-8 days after the explosion) to late phase (50-51 days after explosion).Similar to the spherical symmetric validation results in Figure 2, the ultraviolet spectropolarimetry below 2000 $\rm\AA$ after $\sim$ 25 days are dominated by Monte-Carlo noise.Moreover, ultraviolet spectropolarimetry at 7-8 days and 15-16 days shows an exceptional high S/N.

An Integral-Based Technique (IBT) to Accelerate the Monte-Carlo Radiative Transfer Computation for Supernovae (8)

An Integral-Based Technique (IBT) to Accelerate the Monte-Carlo Radiative Transfer Computation for Supernovae (9)

4.3 Cross Comparison

In this section, we compare the simulation results from the EBT, IBT, and DCT methods using the D3D model.Figure 4 shows the spectra and linear polarization time sequence calculated by the three methods.In the DCT calculation, the viewing direction is split into a $10\times 10$ bin and $4.6\times 10^{8}$ OPs are generated per time step to achieve a reasonable S/N.For the EBT and IBT calculations, we only generate $4\times 10^{6}$ OPs per time step.Figure 4 shows that the IBT calculation produced spectra and polarization with a S/N comparable to a DCT calculation that used $\sim 115$ times more OPs.This noise reduction holds from 6 days to 33 days after the explosion and from ultraviolet to infrared wavelengths.

Due to the high opacity at the early phase of SN, the OPs undergo multiple scattering events near the core of SN, which results in a rapidly growing number of VPs in the EBT calculation and an increasing memory usage expense.Therefore, in this example EBT calculation we started the calculation at 25 days after the explosion due to the limit of hardware.

An Integral-Based Technique (IBT) to Accelerate the Monte-Carlo Radiative Transfer Computation for Supernovae (12)

An Integral-Based Technique (IBT) to Accelerate the Monte-Carlo Radiative Transfer Computation for Supernovae (13)

5 Computational Performance

5.1 Computational Speed Comparison

The computational speed of the three methods is measured on one 48-core compute node of Grace supercomputer in TAMU HPRC ¹¹1https://hprc.tamu.edu/.Each node has 2 sockets of Intel Xeon 6248R CPUs, in total 48 cores, and has 384 GBs of RAM memory.The majority of the computation time in a 3-D time-dependent radiative transfer computation consists of two components: (1) the computation of the plasma state, level population, and opacity of the grid; (2) the propagation of energy packets including OPs, VPs, or IPs.The choice of spectral synthesis methods only affects the propagation time of the energy packets.

An Integral-Based Technique (IBT) to Accelerate the Monte-Carlo Radiative Transfer Computation for Supernovae (15)

Figure 5 shows the computation time of particle propagation at each time step as a function of the time after explosion using the EBT and IBT methods as described in Section 4.3, and the corresponding DCT computation time with $4\times 10^{6}$ OPs per time step, the same number of OPs as the EBT and IBT results.The time needed to compute the plasma state and opacity does not change drastically at different time steps.In contrast, the computation time for packet propagation increases in the first $\sim$ 3 days, because the optical depth is high and the OPs cannot easily escape the SN ejecta and undergo repeated interaction events.With the expansion of the SN ejecta, the optical depth decreases, and the accumulated OPs escape the SN ejecta, therefore the computation time decreases.In particular, the computation time for EBT is a factor of $\sim$ 30 higher than DCT around 25 days after the explosion, and is a factor of $\sim$ 10 higher than DCT around 60 days.In contrast, the introduction of IBT only increases the computation time by a factor of 0.3-0.5 of the DCT computation time, throughout the time after the SN explosion, to reach a comparable S/N as EBT.

5.2 Signal to Noise Ratio

In an MCRT simulation, the output spectral flux error scales roughly as $\sigma\propto\sqrt{N}$ , where $N$ is the number of photon packets.We therefore use the simulation results on the D1D model to measure the simulation error of the three spectral retrieval methods.First, we split the output spectra into ultraviolet wavelength range (1000-2000 $\AA$ ), optical wavelength range (2000-10000 $\AA$ ), and infrared wavelength range (10000-25000 $\AA$ ).Second, we measure the simulation root mean squared error (RMSE) of the Stokes parameters $\mathcal{Q}$ and $\mathcal{U}$ using the following equation:

RMSE=\sqrt{\sum_{N_{\nu}}\frac{\left(\frac{\mathcal{Q}(\nu)}{\mathcal{I}(\nu)}%-0\right)^{2}+\left(\frac{\mathcal{U}(\nu)}{\mathcal{I}(\nu)}-0\right)^{2}}{N_%{\nu}}}\ ,

(17)

where $N_{\nu}$ is the number of frequency indicies in the selected wavelength range, the summation is over all frequency indicies in the selected wavelength range.Note that the theoretical solution of the D1D model is no polarization which results in the zero term in the equation.The RMSE observes the inverse square root relation with the number of packets, $RMSE\propto N^{-0.5}$ .

An Integral-Based Technique (IBT) to Accelerate the Monte-Carlo Radiative Transfer Computation for Supernovae (16)

An Integral-Based Technique (IBT) to Accelerate the Monte-Carlo Radiative Transfer Computation for Supernovae (17)

Figure 6 shows the relationship between the number of OPs per time step and the resulting RMSE at 10-11 days after the explosion using DCT or IBT, and at 37-38 days after the explosion using DCT, EBT, or IBT.We notice the result from IBT could reproduce the $RMSE\propto N^{-0.5}$ relation in most of the time and wavelength ranges.The DCT result also reproduces the relation $RMSE\propto N^{-0.5}$ in the optical and infrared wavelength ranges when the number of packets is large enough.However, in the ultraviolet wavelength, and in the optical wavelength with the number of packets per step smaller than $\sim 4\times 10^{5}$ , the inverse square root relation is not observed due to the large simulation error.The result from EBT also reproduces the $RMSE\propto N^{-0.5}$ relation in optical ind infrared wavelength ranges at 37-38 days after the explosion.In the ultraviolet wavelength range at 37-38 days after the explosion, none of the methods could simulate good S/N spectra to reproduce the $RMSE\propto N^{-0.5}$ relation.

A direct comparison of the computation efficiency of the three methods could be made using the RMSE values measured in the optical and infrared wavelength range.The RMSE of IBT is smaller than that of DCT by a factor of $20-30$ , leading to a factor of $400-900$ less packets needed to calculate to reach similar S/N.The RMSE of EBT is smaller than that of DCT by a factor of $5-10$ , resulting in a factor of $25-100$ less packets to calculate to achieve a comparable S/N.The accuracy comparison between EBT and IBT is consistent with the reports in Table 1 of Bulla etal. (2015).

6 Conclusion

We present IBT, an algorithm to efficiently synthesize the spectropolarimetry signal from 3-D MCRT simulations.Comparing to EBT (Bulla etal., 2015), the present IBT method is upgraded in the following aspects:

•
In EBT, a VP is created at each scattering event of an OP, which could lead to an uncontrolled number of VPs and limit the application of EBT on high-opacity plasma models (i.e. SNe Ia at $\sim$ 3 days after the explosion). However, IBT uses the cell radiance $\mathbf{R}(\vec{x},t,\nu)$ as a proxy of all VPs in a volume cell of the plasma model, and limits the number of IPs to save memory and computation time in tracing the particles.
•
The frequency grid in IBT is logarithmically spaced (Equation 8), which simplifies the Doppler shift of cell radiance $\mathbf{R}(\vec{x},t,\nu)$ to re-indexing the pixel. Therefore, the computation of Equation 15 is simplified to vector operations. This upgrade further accelerates the computation of IP propagation.

In the tests on the D3D model, both IBT and EBT have successfully reproduced the spectral time sequence with spectropolarimetry from 1000Åto 25000Å, and the S/N is comparable to DCT with $\sim$ 115 more Monte-Carlo OPs.The computation time for IBT is $\sim$ 0.3 of DCT with the same number of OPs per viewing direction, while the computation time for EBT is a factor of $10-30$ of DCT per viewing direction.With an opacity limit and wavelength limit to accelerate the computation, the EBT computation time is still about a factor of two of DCT per viewing direction(Bulla etal., 2015).Moreover, IBT has successfully calculated the spectral time sequence from 1 day to 60 days after the explosion, while EBT failed before day 25 due to the limited memory space of the computing facility.

In the tests on the D1D model, IBT has successfully reproduced the zero polarization signal as the theoretical predictions from the spherical symmetry.Using the same number of Monte-Carlo OPs, the S/N of the IBT spectrum is $\sim$ 30 times higher than the DCT spectrum S/N and $5-10$ times higher than the EBT spectrum.

Chen etal. (2020, 2024) developed theArtificial Intelligence Assisted Inversion (AIAI) method which enables theoretical models of SNeIa to be derived using the observational data as input constraints.The published AIAI is based on 1-D models.The studies demonstrate that it is possible to combine empirical models of SNeIa with detailed radiative transfer code to improve the use of SNeIa as cosmological probes.In the future, we will extend the 1-D models to 3-D time-dependent models.The proposed IBT is the ideal algorithm to efficiently generate a time sequence of spectrophotometric and spectropolarimetric data bases which can be used to train AI models, this can be an important step in assimilating the extensive observational data of SNeIa into studies of the physics of thermonuclear explosions and the use of SNeIa as cosmological distance candles.

X.C. is supported by the grant NSF-AST 1817099.The authors would like to thank Prof. J. C. Wheeler from University of Texas at Austin and Prof. David Jeffery from University of Nevada at Las Vegas for supportive discussions.Portions of this research were conducted with the advanced computing resources provided by Texas A&M High Performance Research Computing.This work used FASTER super computer at TAMU HPRC through allocation PHY240215 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by National Science Foundation grants 2138259, 2138286, 2138307, 2137603, and 2138296(Boerner etal., 2023).

Appendix A Rayleigh Scattering

The Thomson scattering between free electron and photon is calculated by the Rayleigh scattering phase matrix:

\boldsymbol{P}_{R}(\Theta)=\frac{3}{4}\begin{pmatrix}cos^{2}\Theta+1&cos^{2}%\Theta-1&0&0\\cos^{2}\Theta-1&cos^{2}\Theta+1&0&0\\0&0&2cos\Theta&0\\0&0&0&2cos\Theta\\\end{pmatrix}\ ,

(A1)

where $\Theta$ is the angle between the incident light and the emergent light.Before the calculation of Rayleigh scattering phase matrix, the reference axis of the incident light should be rotated into the scattering plane, and the polarization state is changed by the rotation matrix:

\boldsymbol{L}(\phi_{in})=\begin{pmatrix}1&0&0&0\\0&cos2\phi_{in}&sin2\phi_{in}&0\\0&-sin2\phi_{in}&cos2\phi_{in}&0\\0&0&0&1\\\end{pmatrix}\ ,

(A2)

where $\phi_{in}$ is the rotation angle between the observer’s reference axis and the scattering plane reference axis.A similar rotation matrix of $\boldsymbol{L}(\pi-\phi_{out})$ is also applied to the emergent light, in order to convert the polarization state to the observer’s reference axis.

We define $\bar{\boldsymbol{S}}$ as the multiplication of the packet energy and the dimensionless Stokes vector:

\bar{\boldsymbol{S}}=\begin{pmatrix}\bar{I}\\\bar{Q}\\\bar{U}\\\bar{V}\end{pmatrix}\ =\bar{E}\bar{\boldsymbol{s}}\ .

(A3)

Therefore, the full Rayleigh scattering formula in the SN radiative transfer simulation is written as:

\bar{\boldsymbol{S}}_{out}=\boldsymbol{L}(\pi-\bar{\phi}_{out})\boldsymbol{P}_%{R}(\bar{\Theta})\boldsymbol{L}(\bar{\phi}_{in})\bar{\boldsymbol{S}}_{in}\ .

(A4)

Note that the Rayleigh scattering happens in the co-moving frame, and all the variables in the above equation are measured in the co-moving frame.In Equation 3, Equation 9, and Equation 13, the operator is defined as $\boldsymbol{P}(\bar{\theta}_{in},\bar{\phi}_{in},\bar{\theta}_{out},\bar{\phi}%_{out})=\boldsymbol{L}(\pi-\bar{\phi}_{out})\boldsymbol{P}_{R}(\bar{\Theta})%\boldsymbol{L}(\bar{\phi}_{in})$ .

Replacing the axis rotation angles and scatter angle ( $\bar{\phi}_{in}$ , $\bar{\phi}_{out}$ , $\bar{\Theta}$ ) with the propagating direction of the incident OP ( $\bar{\theta}_{in}$ , $\bar{\varphi}_{in}$ ) and the propagating direction of the emergent VP/IP ( $\bar{\theta}$ , $\bar{\varphi}$ ), the expression of the operator $\boldsymbol{P}(\bar{\theta}_{in},\bar{\phi}_{in},\bar{\theta},\bar{\phi})$ becomes (Chandrasekhar, 1960):

$\displaystyle\bar{I}_{out}$	$\displaystyle=$	$\displaystyle\frac{3}{8\pi}\left(\left((l,l)^{2}+(l,r)^{2}\right)\bar{I}_{in}+%\left((r,l)^{2}+(r,r)^{2}\right)\bar{Q}_{in}+\left((l,l)(r,l)+(l,r)(r,r)\right%)\bar{U}_{in}\right)$	(A5)
$\displaystyle\bar{Q}_{out}$	$\displaystyle=$	$\displaystyle\frac{3}{8\pi}\left(\left((l,l)^{2}-(l,r)^{2}\right)\bar{I}_{in}+%\left((r,l)^{2}-(r,r)^{2}\right)\bar{Q}_{in}+\left((l,l)(r,l)-(l,r)(r,r)\right%)\bar{U}_{in}\right)$	(A6)
$\displaystyle\bar{U}_{out}$	$\displaystyle=$	$\displaystyle\frac{3}{8\pi}\left(2(l,l)(l,r)\bar{I}_{in}+2(r,r)(r,l)\bar{Q}_{%in}+\left((l,l)(r,r)+(r,l)(l,r)\right)\bar{U}_{in}\right)$	(A7)
$\displaystyle\bar{V}_{out}$	$\displaystyle=$	$\displaystyle\frac{3}{8\pi}\left((l,l)(r,r)-(r,l)(l,r)\right)\bar{V}_{in}\ ,$	(A8)

where $(l,l)$ , $(r,l)$ , $(l,r)$ , $(r,r)$ are:

$\displaystyle(l,l)$	$\displaystyle=$	$\displaystyle sin(\bar{\theta})sin(\bar{\theta_{in}})+cos(\bar{\theta})cos(%\bar{\theta_{in}})cos(\bar{\varphi_{in}}-\bar{\varphi})$	(A9)
$\displaystyle(r,l)$	$\displaystyle=$	$\displaystyle cos(\bar{\theta})sin(\bar{\varphi_{in}}-\bar{\varphi})$	(A10)
$\displaystyle(l,r)$	$\displaystyle=$	$\displaystyle-cos(\bar{\theta_{in}})sin(\bar{\varphi_{in}}-\bar{\varphi})$	(A11)
$\displaystyle(r,r)$	$\displaystyle=$	$\displaystyle cos(\bar{\varphi_{in}}-\bar{\varphi})\ .$	(A12)

References

Abbott etal. (2019)Abbott, T.M.C., Allam, S., Andersen, P., etal. 2019, ApJ, 872,L30, doi:10.3847/2041-8213/ab04fa
Blondin etal. (2013)Blondin, S., Dessart, L., Hillier, D.J., & Khokhlov, A.M. 2013,MNRAS, 429, 2127, doi:10.1093/mnras/sts484
Blondin etal. (2022)Blondin, S., Blinnikov, S., Callan, F.P., etal. 2022, A&A, 668,A163, doi:10.1051/0004-6361/202244134
Boerner etal. (2023)Boerner, T.J., Deems, S., Furlani, T.R., Knuth, S.L., & Towns, J. 2023, inPractice and Experience in Advanced Research Computing, 173–176
Brout etal. (2022)Brout, D., Scolnic, D., Popovic, B., etal. 2022, ApJ, 938, 110,doi:10.3847/1538-4357/ac8e04
Bulla etal. (2015)Bulla, M., Sim, S.A., & Kromer, M. 2015, MNRAS, 450, 967,doi:10.1093/mnras/stv657
Bulla etal. (2016a)Bulla, M., Sim, S.A., Pakmor, R., etal. 2016a, MNRAS,455, 1060, doi:10.1093/mnras/stv2402
Bulla etal. (2016b)Bulla, M., Sim, S.A., Kromer, M., etal. 2016b, MNRAS,462, 1039, doi:10.1093/mnras/stw1733
Bulla etal. (2021)Bulla, M., Kyutoku, K., Tanaka, M., etal. 2021, MNRAS, 501, 1891,doi:10.1093/mnras/staa3796
Castor (1972)Castor, J.I. 1972, ApJ, 178, 779, doi:10.1086/151834
Chandrasekhar (1960)Chandrasekhar, S. 1960, Radiative transfer
Chen etal. (2020)Chen, X., Hu, L., & Wang, L. 2020, ApJS, 250, 12,doi:10.3847/1538-4365/ab9a3b
Chen etal. (2024)Chen, X., Wang, L., Hu, L., & Brown, P.J. 2024, ApJ, 962, 125,doi:10.3847/1538-4357/ad0a33
Cikota etal. (2019)Cikota, A., Patat, F., Wang, L., etal. 2019, MNRAS, 490, 578,doi:10.1093/mnras/stz2322
Ferrazzoli etal. (2023)Ferrazzoli, R., Slane, P., Prokhorov, D., etal. 2023, ApJ, 945, 52,doi:10.3847/1538-4357/acb496
Gamezo etal. (2004)Gamezo, V.N., Khokhlov, A.M., & Oran, E.S. 2004, Phys.Rev.Lett., 92, 211102,doi:10.1103/PhysRevLett.92.211102
Inserra etal. (2016)Inserra, C., Bulla, M., Sim, S.A., & Smartt, S.J. 2016, ApJ, 831,79, doi:10.3847/0004-637X/831/1/79
Kasen etal. (2006)Kasen, D., Thomas, R.C., & Nugent, P. 2006, ApJ, 651, 366,doi:10.1086/506190
Kerzendorf & Sim (2014)Kerzendorf, W.E., & Sim, S.A. 2014, MNRAS, 440, 387,doi:10.1093/mnras/stu055
Kromer & Sim (2009)Kromer, M., & Sim, S.A. 2009, MNRAS, 398, 1809,doi:10.1111/j.1365-2966.2009.15256.x
Kromer etal. (2013)Kromer, M., Fink, M., Stanishev, V., etal. 2013, MNRAS, 429, 2287,doi:10.1093/mnras/sts498
Lucy (2002)Lucy, L.B. 2002, A&A, 384, 725, doi:10.1051/0004-6361:20011756
Lucy (2003)—. 2003, A&A, 403, 261, doi:10.1051/0004-6361:20030357
Mihalas (1978)Mihalas, D. 1978, Stellar atmospheres (San Francisco: W.H.Freeman andCompany)
Phillips (1993)Phillips, M.M. 1993, ApJ, 413, L105, doi:10.1086/186970
Wang & Wheeler (2008)Wang, L., & Wheeler, J.C. 2008, ARA&A, 46, 433,doi:10.1146/annurev.astro.46.060407.145139
Wang etal. (2003)Wang, L., Baade, D., Höflich, P., etal. 2003, The Astrophysical Journal,591, 1110, doi:10.1086/375444
Yang etal. (2022)Yang, Y., Yan, H., Wang, L., etal. 2022, The Astrophysical Journal, 939, 18,doi:10.3847/1538-4357/ac8d5f