We present an integralbased technique IBT algorithm to accelerate supernova SN radiative transfer calculations. The algorithm utilizes integral packets which are calculated by the path integral of the MonteCarlo energy packets to synthesize the observed spectropolarimetric signal at a given viewing direction in a 3D timedependent radiative transfer program. Compared to the eventbased technique EBT proposed by Bulla et al. 2015 our algorithm significantly reduces the computation time and increases the MonteCarlo signaltonoise ratio. Using a 1D spherical symmetric type Ia supernova SN Ia ejecta model DDC10 and its derived 3D model the IBT algorithm has successfully passed the verification of 1 spherical symmetry 2 mirror symmetry 3 cross comparison on a 3D SN model with directcounting technique DCT and EBT. Notably with our algorithm implemented in the 3D MonteCarlo radiative transfer code SEDONA the computation time is faster than EBT by a factor of 1030 and the signaltonoise SN ratio is better by a factor of 510 with the same number of MonteCarlo quanta.
['Astrophysics - High Energy Astrophysical Phenomena']
An IntegralBased Technique IBT to Accelerate the MonteCarlo Radiative Transfer Computation for Supernovae
{'An Integral-Based Technique (IBT) to Accelerate the Monte-Carlo Radiative Transfer Computation for Supernovae': 'Xingzhuo Chen , 1 Lifan Wang, 1 and Daniel Kasen 2  \n1 George P. and Cynthia Woods Mitchell Institute for Fundamental Physics & Astronomy,  \nTexas A. & M. University, Department of Physics and Astronomy, 4242 TAMU, College Station, TX 77843, USA 2 Astronomy 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 et al. (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.  \nKeywords: Supernova", '1. INTRODUCTION': "Type Ia supernovae (SNe Ia) have been widely applied in cosmological studies (e.g. Abbott et al. (2019); Brout et al. (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.  \nRadiative 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 et al. 2022).  \nCorresponding author: Lifan Wang  \nlifan@tamu.edu  \nHowever, imaging of nearby SN Ia remnants (e.g. Ferrazzoli et al. (2023)) and hydrodynamic simulations (e.g. Gamezo et al. (2004); Kromer et al. (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 et al. 2003; Cikota et al. 2019; Yang et al. 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).  \nMonte-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 3D time-dependent radiative transfer codes such as SEDONA (Kasen et al. 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 espe-  \ncially 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 et al. (2015) introduced an 'event-based technique' (EBT) for spectropolarimetry synthesis, which has been applied to simulations of SNe Ia (Bulla et al. 2016a,b), superluminous SNe (Inserra et al. 2016), and kilonovae (Bulla et al. 2021).  \nIn 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 et al. 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.  \nThe polarization information of an OP is stored as a dimensionless Stokes vector  \ns OP =      1 q u v      , (1)  \nand the Stokes vector can be constructed with s OP and the energy of the OP E OP .  \nTypically, 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  \nthe 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  \n     I Q U V      = 1 4 πr 2 ∆ t ∆ ν ∆Ω ∑ i E OP,i s OP,i (2)  \nwhere 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, ∆ t is the size of the arrival time bin, ∆ ν is the size of the frequency bin, r is the distance to the SN center, ∆Ω is the viewing direction bin size. The summation is over all the OPs in the arrival time interval [ t -∆ t/ 2 , t +∆ t/ 2), the frequency interval [ ν -∆ ν/ 2 , ν + ∆ ν/ 2), and the viewing direction interval [ φ -∆ φ/ 2 , φ + ∆ φ/ 2); [ θ -∆ θ/ 2 , θ +∆ θ/ 2).  \nIn 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 3D simulation by Bulla et al. (2015). Generally, VPs are used as follows  \n- · An observer viewing direction is specified at the start of the MCRT simulation.\n- · 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.\n- · The energy of each VP is reduced by the optical depth along the path to the observer.\n- · Finally, the escaped VPs are used to synthesize the spectrum for the specified observer viewing angle.  \nEvent-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):  \n¯ S V P = ¯ S OP P ( ¯ θ OP , ¯ φ OP , ¯ θ, ¯ φ ) , (3)  \nwhere ¯ S OP is the Stokes vector of the OP in the comoving frame, ( ¯ θ OP , ¯ φ OP ) is the OP propagating direction in the co-moving frame, ( ¯ θ, ¯ φ ) is the viewing direction in the co-moving frame, P is the combination of the rotation matrix and the scattering matrix (See Equation 10 in Bulla et al. (2015)). In the following, comoving frame quantities are denoted with a bar, whereas non-bar quantities refer to the lab frame. A detailed expression of the P matrix is given in Appendix A.  \nIf 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  \n¯ E V P = 1 4 π ¯ E OP . (4)  \nAfter 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  \nτ tot = ∑ j d j α tot,j , (5)  \nwhere d j is the length of the j -th line segment of the VP trajectory, α 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  \nα tot,j ( ν ) = (¯ ν/ν )¯ α tot,j (¯ ν ) , (6)  \nwhere ¯ ν and ν are the co-moving frame frequency and the rest-frame frequency respectively.  \nThe VP is moved in small line segments, the length of which are determined by finding the minimum distance among the following  \n- · The VP reaches the boundary of a volume cell.\n- · The VP reaches the end of a time step.\n- · The co-moving frame frequency of the VP reaches the boundary of the frequency grid of α tot .  \nThe final spectrum is constructed by summing the VPs  \n     I Q U V      = 1 4 πr 2 ∆ t ∆ ν ∑ i E V P,i s V P,i e -τ tot,i , (7)  \nthe summation is over all the VPs in the frequency bin [ ν -∆ ν/ 2 , ν + ∆ ν/ 2), and arrival time interval [ t -∆ t/ 2 , t +∆ t/ 2).  \nThe EBT can be accelerated by only emitting VPs with frequencies within a desired spectral range (e.g., 3500 -10000 ˚ A ) for the synthetic observables, and by removing the VP once the integrated optical depth reaches a critical value τ 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 ∼ 4 while not affecting accuracy.", '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 S V P , while an IP stores the energy and polarization at multiple frequency bins with a 4 × N tensor of cell radiance R ([ I, Q, U, V ]; [ ν 1 , ..., ν 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 [ ν 0 , ν 1 ...ν N ]. The frequency ratio between the adjacent pixels is a constant C , satisfying  \nν n +1 /ν n = C. (8)  \nUsing 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.  \nWe define the cell radiance R ( ⃗x, t, ν ) in the lab frame for a cell located at coordinate ⃗x and time t as the energy radiated per frequency bin at frequency ν toward the viewing direction, the units of R ( ⃗x, t, ν ) is erg Hz -1 sr -1 . The formal solution of the cell radiance is:  \nR ( ⃗x, t, ν ) = ∆ v ∆ t s ( ¯ ν ν ) -2 ( ¯ j emi + ¯ j sc ) (9)  \n¯ j emi = ¯ j emi ( ⃗x, t, ¯ ν )      1 0 0 0      (10)  \n¯ j sc = α T ∫ Ω ¯ I ( ⃗x, t, ¯ θ in , ¯ φ in , ¯ ν ) P ( ¯ θ in , ¯ φ in , ¯ θ, ¯ φ ) d Ω , (11)  \nwhere ¯ I is the specific intensity as a four-vector of the Stokes parameters; ∆ t s is the size of the simulation time step; ¯ j emi and ¯ j sc are the four-vectors of the emission terms due to depolarized isotropic emission and scattered light from Thomson scattering, respectively; the ¯ j emi function is the depolarized isotropic emission term, which includes all the emission from bound-bound, bound-free, and free-free transitions; the P ( ¯ θ in , ¯ φ in , ¯ θ, ¯ φ ) function is a combination of the Rayleigh scattering phase matrix and rotation matrix discussed in Appendix A; α T is the extinction coefficient for Thomson scattering; the integral is over all the incident directions ( ¯ θ in , ¯ φ in ); the (¯ ν/ν ) -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).  \nFor each IP, the cell radiance R is calculated on the universal frequency grid and stored as a 4 × N vector. During the propagation of an IP, R is reduced by the extinction coefficient on the trajectory. The spectral synthesis is the summation of R over the IPs at a specific arrival time. Section 3.2 and Section 3.3 illustrate the computation of R in a MCRT simulation during the creation of an IP. Section 3.4 illustrates the update of 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  \n∆ v ∆ t s ( ¯ ν ν ) -2 ¯ j emi ( ⃗x, t, ¯ ν ) = ( ¯ ν ν ) -2 ∑ i ¯ E OP,i 4 π ∆¯ ν , (12)  \nwhere the summation is over all the OPs within the frequency bin [¯ ν -∆¯ ν/ 2 , ¯ ν +∆¯ ν/ 2), the volume cell ∆ v , and the time step bin [ t s , t s +∆ t s ). Only OPs that are newly created or undergone the same physical events are taken into account.  \nIn the present SEDONA calculation, we include bound-free and free-free transitions and adopt the line  \nexpansion 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 ¯ j emi ( ⃗x, t, ¯ ν ) = B (¯ ν )¯ α abs ( ⃗x, t, ¯ ν ), where B (¯ ν ) is the blackbody spectrum and ¯ α abs is the total extinction coefficient including bound-bound, bound-free, and freefree 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  \n∆ v ∆ t s ( ¯ ν ν ) -2 ¯ j sc = ∑ i,j ¯ E i ¯ d i,j ¯ α T ¯ s i,j P ( ¯ θ i,j , ¯ φ i,j , ¯ θ, ¯ φ ) 4 π ∆¯ ν (13)  \nwhere ¯ 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 [¯ ν -∆¯ ν/ 2 , ¯ ν +∆¯ ν/ 2), the volume cell ∆ v , the time step bin [ t s , t s +∆ t s ).  \nFor 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 R is calculated on the universal frequency grid and stored in IP as a 4 × 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 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:  \n∣ ∣ ∣ ∣ ¯ ν begin -¯ ν end ¯ ν end ∣ ∣ ∣ ∣ ≤ C -1 , (14)  \nwhere ¯ ν begin is the co-moving frame frequency at the beginning of the line segment, ¯ ν end is the co-moving frame frequency of the same lab frame frequency ν at the ending of the line segment, the constant C is defined in Equation 8.  \nIn 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  \ntrajectory without accumulating the total opacity over the trajectory  \nR i,j +1 ( ν ) = R i,j ( ν ) e -d i,j α tot,i,j ( ν ) , (15)  \nwhere d i,j is the length of the i -th IP at the j -th line segment over the trajectory, α tot,i,j ( ν ) is the total extinction coefficient of the i -th IP at the j -th line segment and at lab frame frequency ν . The final spectrum is the sum of the IPs:  \n     I Q U V      = 1 4 πr 2 ∆ t ∑ i R i,final , (16)  \nwhere 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 -∆ t/ 2 , t +∆ t/ 2). Note that the grid size of arrival time ∆ t should be much larger than the time step size ∆ t s to avoid grid mismatch issue, while this issue is not significant in DCT or EBT calculations.', '4. MODEL VALIDATION': 'In this section, we prepare two toy models to validate the IBT computation results. The first model is the SN Ia delayed-detonation model DDC10. DDC10 is calculated from hydrodynamic simulations in Blondin et al. (2013) and adopted as a benchmark in Blondin et al. (2022) to evaluate the differences between radiative transfer codes. We remap the 1-D model into a 60 × 60 × 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 . 3 y and z new = z/ 1 . 3, in order to introduce asphericity. This model is also realized on a 60 × 60 × 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.  \nIn the following text, the spherical symmetric DDC10 model is denoted as D1D, and the modified ellipsoidal 3Dmodel 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 ∆ t s /t s = 0 . 03. After 6.6 days, the simulation time step is ∆ t s = 0 . 2 days. The arrival time bin size is ∆ t = 1 day, which satisfies ∆ t ≥ 5∆ t s . The frequency grid ranges from 1 × 10 14 Hz to 5 × 10 15 Hz (30,000 ˚ A to 600 ˚ A ) with ∆ ν/ν = 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 × 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 ˚ A wavelength range (within which most SNe Ia spectropolarimetry observations are made (Cikota et al. 2019)) is as low as ∼ 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 R above the photosphere has low S/N. As a result, the spectrum shows low S/N at 34 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 ˚ A and can be alleviated by increasing the number of OPs.', '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 × 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 ( θ ) = 0 . 51122 , φ = π/ 6) and ( cos ( θ ) = 0 . 51122 , φ = 11 π/ 6) should have the same Stokes Q terms and inverted 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 ˚ A 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 ˚ A after ∼ 25 days are dominated by MonteCarlo noise. Moreover, ultraviolet spectropolarimetry at 7-8 days and 15-16 days shows an exceptional high S/N.  \nFigure 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 × 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 ∼ 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 ∼ 30 higher than DCT around 25 days after the explosion, and is a factor of ∼ 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.  \n<!-- image -->  \nFigure 1. Left panel: the mass fractions of several isotopes in the DDC10 model. Middle panel: the density structure of the DDC10 model. Homologous expansion is assumed, so the radius is equal to velocity times time. Right panel: the geometry of the modified DDC10 model with 3-D structures (D3D model).  \n<!-- image -->', '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 × 10 bin and 4 . 6 × 10 8 OPs are generated per time step to achieve a reasonable S/N. For the EBT and IBT calculations, we only generate 4 × 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 ∼ 115 times more OPs. This noise reduction holds from 6 days to 33 days after the explosion and from ultraviolet to infrared wavelengths.  \nDue 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.', '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 1 . 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  \nmethods only affects the propagation time of the energy packets.', '5.2. Signal to Noise Ratio': 'In an MCRT simulation, the output spectral flux error scales roughly as σ ∝ √ 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 (10002000 ˚ A ), optical wavelength range (2000-10000 ˚ A ), and infrared wavelength range (10000-25000 ˚ A ). Second, we measure the simulation root mean squared error (RMSE) of the Stokes parameters Q and U using the following equation:  \nFigure 6 shows the relationship between the number of OPs per time step and the resulting RMSE at 1011 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 ∝ N -0 . 5 relation in most of the time and  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 2. The spectral flux ( I ) and the polarization percentage ( Q / I , U / I ) of the spherical symmetric D1D model at different times. The spectra are calculated by IBT with 6 × 10 6 OPs per simulation time step. The arrival time bin relative to the explosion time is labeled in the title of each subplots. The polarization percentages of Q / I and U / I are averaged using a 7-pixel bin.  \n<!-- image -->  \nRMSE = √ √ √ √ √ ∑ N ν ( Q ( ν ) I ( ν ) -0 ) 2 + ( U ( ν ) I ( ν ) -0 ) 2 N ν , (17)  \nwhere N ν 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 ∝ N -0 . 5 .  \nwavelength ranges. The DCT result also reproduces the relation RMSE ∝ 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 ∼ 4 × 10 5 , the inverse square root relation is not observed due to the large simulation error. The result from EBT also reproduces the RMSE ∝ 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 ∝ N -0 . 5 relation.  \nA 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  \n<!-- image -->  \nFigure 3. The spectropolarimetry of the D3D model at the two symmetric viewing directions ( cos ( θ ) = 0 . 51122 , φ = π/ 6) and ( cos ( θ ) = 0 . 51122 , φ = 11 π/ 6). The arrival time bin at each panel is labeled in the title of figures. In the sub-figures of Stokes U / I term, the flipped values at one of the viewer points is shown as dashed line, in order to make direct comparisons. The polarization percentages of Q / I and U / I are averaged using a 7-pixel bin.  \n<!-- image -->  \npackets 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).', '6. CONCLUSION': 'We present IBT, an algorithm to efficiently synthesize the spectropolarimetry signal from 3-D MCRT simulations. Comparing to EBT (Bulla et al. 2015), the present IBT method is upgraded in the following aspects:  \n- · 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 ∼ 3 days after the explosion). However, IBT uses the cell radiance R ( ⃗x, t, ν ) 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.\n- · The frequency grid in IBT is logarithmically spaced (Equation 8), which simplifies the Doppler shift of cell radiance R ( ⃗x, t, ν ) 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.  \nIn the tests on the D3D model, both IBT and EBT have successfully reproduced the spectral time sequence with spectropolarimetry from 1000 ˚ A to 25000 ˚ A, and the S/N is comparable to DCT with ∼ 115 more MonteCarlo OPs. The computation time for IBT is ∼ 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 et al. 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.  \n<!-- image -->  \n<!-- image -->  \nFigure 4. The spectra and linear polarization of the D3D model at the viewing direction cos ( θ ) = 0 . 5, φ = 3 . 4455, using DCT (black), EBT (green), and IBT (red) methods. The DCT calculation is performed with 4 . 6 × 10 8 OPs per time step (which is 115 times more than EBT or IBT, marked in the legend), while EBT and IBT calculation generates 4 × 10 6 OPs per time step. The arrival time bin is shown on top of each panel. Due to the memory limit, the EBT calculation is only started at 25 days after the explosion, and the EBT results are not shown at 6-7 days and 18-19 days. The Monte-Carlo error bar of the DCT result is estimated by the number of escaped OPs at each bin. Due to the lack of OPs, the spectrum between 1000 ˚ A and 1500 ˚ A at 32-33 days from DCT is missing. The polarization percentage is binned with 7 pixels.  \n<!-- image -->  \nFigure 5. Comparison of particle propagation computation time at different time steps using the DCT (black), EBT (green), and IBT (red) methods. Each time step generates 4 × 10 6 OPs in the three simulations. The computation time of the plasma state and opacity is also shown in the figure with ultramarine color. Note that the EBT calculation is not activated until 25 days after the explosion, and the green curve before day 25 is overlapped to the black curve.  \n<!-- image -->  \nIn 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 ∼ 30 times higher than the DCT spectrum S/N and 5 -10 times higher than the EBT spectrum.  \nChen et al. (2020, 2024) developed the Artificial Intelligence Assisted Inversion (AIAI) method which enables theoretical models of SNe Ia 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 SNe Ia with detailed radiative transfer code to improve the use of SNe Ia 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 SNe Ia into studies of the physics of thermonuclear explosions and the use of SNe Ia as cosmological distance candles.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 6. The RMSE of DCT, EBT, and IBT simulation results on the D1D model, as a function of the number of OPs generated per time step. The polarization spectra are binned with 7 pixels. Upper panel shows the RMSE in the time bin of 10-11 days after the explosion, lower panel shows the RMSE in the time bin of 37-38 days after the explosion. From left to right, the wavelength range is ultraviolet, optical, and infrared. The dashed line shows the theoretical inverse square root relation between RMSE and the number of packets ( RMSE ∝ N -0 . 5 ) for comparison.  \n<!-- image -->  \nSoftware: SEDONA(Kasen et al. 2006)  \n- X.C. is supported by the grant NSF-AST 1817099. 1\n- The authors would like to thank Prof. J. C. Wheeler 2\n- from University of Texas at Austin and Prof. David 3\n- Jeffery from University of Nevada at Las Vegas for sup4\n- portive discussions. Portions of this research were con5\n- ducted with the advanced computing resources provided 6\n- by Texas A&M High Performance Research Comput7\n- ing. This work used FASTER super computer at TAMU 8\n- HPRC through allocation PHY240215 from the Ad9\n- vanced Cyberinfrastructure Coordination Ecosystem: 10\n- Services & Support (ACCESS) program, which is sup11\n- ported by National Science Foundation grants 2138259, 12\n- 2138286, 2138307, 2137603, and 2138296(Boerner et al. 13\n- 2023). 14', 'A. RAYLEIGH SCATTERING': 'The Thomson scattering between free electron and photon is calculated by the Rayleigh scattering phase matrix:  \nP R (Θ) = 3 4      cos 2 Θ+1 cos 2 Θ -1 0 0 cos 2 Θ -1 cos 2 Θ+1 0 0 0 0 2 cos Θ 0 0 0 0 2 cos Θ      , (A1)', 'Chen et al.': "where Θ 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:  \nL ( ϕ in ) =      1 0 0 0 0 cos 2 ϕ in sin 2 ϕ in 0 0 -sin 2 ϕ in cos 2 ϕ in 0 0 0 0 1      , (A2)  \nwhere ϕ in is the rotation angle between the observer's reference axis and the scattering plane reference axis. A similar rotation matrix of L ( π -ϕ out ) is also applied to the emergent light, in order to convert the polarization state to the observer's reference axis.  \nWe define ¯ S as the multiplication of the packet energy and the dimensionless Stokes vector:  \n¯ S =      ¯ I ¯ Q ¯ U ¯ V      = ¯ E ¯ s . (A3)  \nTherefore, the full Rayleigh scattering formula in the SN radiative transfer simulation is written as:  \n¯ S out = L ( π -¯ ϕ out ) P R ( ¯ Θ) L ( ¯ ϕ in ) ¯ S in . (A4)  \nNote 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 P ( ¯ θ in , ¯ ϕ in , ¯ θ out , ¯ ϕ out ) = L ( π -¯ ϕ out ) P R ( ¯ Θ) L ( ¯ ϕ in ).  \nReplacing the axis rotation angles and scatter angle ( ¯ ϕ in , ¯ ϕ out , ¯ Θ) with the propagating direction of the incident OP ( ¯ θ in , ¯ φ in ) and the propagating direction of the emergent VP/IP ( ¯ θ , ¯ φ ), the expression of the operator P ( ¯ θ in , ¯ ϕ in , ¯ θ, ¯ ϕ ) becomes (Chandrasekhar 1960):  \n¯ I out = 3 8 π (( ( l, l ) 2 +( l, r ) 2 ) ¯ I in + ( ( r, l ) 2 +( r, r ) 2 ) ¯ Q in +(( l, l )( r, l ) + ( l, r )( r, r )) ¯ U in ) (A5)  \n¯ Q out = 3 8 π (( ( l, l ) 2 -( l, r ) 2 ) ¯ I in + ( ( r, l ) 2 -( r, r ) 2 ) ¯ Q in +(( l, l )( r, l ) -( l, r )( r, r )) ¯ U in ) (A6)  \n¯ U out = 3 8 π ( 2( l, l )( l, r ) ¯ I in +2( r, r )( r, l ) ¯ Q in +(( l, l )( r, r ) + ( r, l )( l, r )) ¯ U in ) (A7)  \n¯ V out = 3 8 π (( l, l )( r, r ) -( r, l )( l, r )) ¯ V in , (A8)  \nwhere ( l, l ), ( r, l ), ( l, r ), ( r, r ) are:  \n( l, l ) = sin ( ¯ θ ) sin ( ¯ θ in ) + cos ( ¯ θ ) cos ( ¯ θ in ) cos ( ¯ φ in -¯ φ ) (A9)  \n( r, l ) = cos ( ¯ θ ) sin ( ¯ φ in -¯ φ ) (A10)  \n( l, r ) = -cos ( ¯ θ in ) sin ( ¯ φ in -¯ φ ) (A11)  \n( r, r ) = cos ( ¯ φ in -¯ φ ) . (A12)", 'REFERENCES': 'Abbott, T. M. C., Allam, S., Andersen, P., et al. 2019,  \nApJL, 872, L30, doi: 10.3847/2041-8213/ab04fa  \nBlondin, S., Dessart, L., Hillier, D. J., & Khokhlov, A. M.  \n2013, MNRAS, 429, 2127, doi: 10.1093/mnras/sts484'}