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

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

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.

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 et al. (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 et al. (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):

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 et al. (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

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

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

where $\bar{\nu}$ and $\nu$ are the co-moving frame frequency and the rest-frame frequency respectively.

The VP is moved in small line segments, the length of which are determined by finding the minimum distance among the following

• •

  The VP reaches the boundary of a volume cell.

• •

  The VP reaches the end of a time step.

• •

  The co-moving frame frequency of the VP reaches the boundary of the frequency grid of $\alpha_{tot}$.

The final spectrum is constructed by summing the VPs

the summation is over all the VPs in the frequency bin $[\nu-\Delta\nu/2,\nu+\Delta\nu/2)$, and arrival time interval $[t-\Delta t/2,t+\Delta t/2)$.

The EBT can be accelerated by only emitting VPs with frequencies within a desired spectral range (e.g., $3500-10000\ \AA$) for the synthetic observables, and by removing the VP once the integrated optical depth reaches a critical value $\tau_{tot}=10$. We do not apply these optimizations in the present calculations, but Bulla et al. (2015) find that they can decrease the computational time by a factor of $\sim 4$ while not affecting accuracy.

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

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:

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.

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

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 et al. 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).

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

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.

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:

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

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:

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.

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 et al., 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.

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.

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.

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.

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.

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:

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}$.

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 et al. (2015).