Radiative-transfer modeling of supernovae in the nebular-phase. A novel treatment of chemical mixing in spherical symmetry

Supernova (SN) explosions, through the metals they release, play a pivotal role in the chemical evolution of the Universe and the origin of life. Nebular phase spectroscopy constrains such metal yields, for example through forbidden line emission associated with OI, CaII, FeII, or FeIII. Fluid instabilities during the explosion produce a complex 3D ejecta structure, with considerable macroscopic, but no microscopic, mixing of elements. This structure sets a formidable challenge for detailed nonlocal thermodynamic equilibrium radiative transfer modeling, which is generally limited to 1D in grid-based codes. Here, we present a novel and simple method that allows for macroscopic mixing without any microscopic mixing, thereby capturing the essence of mixing in SN explosions. With this new technique, the macroscopically mixed ejecta is built by shuffling in mass space, or equivalently in velocity space, the shells from the unmixed coasting ejecta. The method requires no change to the radiative transfer, but necessitates high spatial resolution to resolve the rapid variation in composition with depth inherent to this shuffled-shell structure. We show results for a few radiative-transfer simulations for a Type II SN explosion from a 15Msun progenitor star. Our simulations capture the strong variations in temperature or ionization between the various shells that are rich in H, He, O, or Si. Because of nonlocal energy deposition, gamma rays permeate through an extended region of the ejecta, making the details of the shell arrangement unimportant. The greater physical consistency of the method delivers spectral properties at nebular times that are more reliable, in particular in terms of individual emission line strengths, which may serve to constrain the SN yields and, for core collapse SNe, the progenitor mass. The method works for all SN types.


Introduction
The complicated 3D structure observed in supernova (SN) remnants from Cas A (Fesen et al. 2006) to the young SN 1987A (Abellán et al. 2017), or the detection of polarized radiation from Type II SNe (Shapiro & Sutherland 1982;Leonard et al. 2006) provide evidence that even standard core-collapse SNe are inherently asymmetric and heterogeneous. The complex 3D distribution of elements likely arises from the neutrino-driven explosion mechanism combined with fluid instabilities triggered by the propagation of the SN shock across the progenitor envelope (Fryxell et al. 1991;Kifonidis et al. 2000;Wongwathanarat et al. 2015;Ono et al. 2020). This structure contrasts with the progenitor chemical distribution, which is thought to exhibit stacked shells of distinct composition (this picture is being challenged by recent simulations; Couch et al. 2015), and with an increasing mean atomic weight towards the denser inner layers of the star, culminating with the Fe core in the innermost layers when the massive star is ripe for core collapse (see, e.g., Arnett 1996).
The complex 3D ejecta structure sets a challenge for the radiative transfer modeling of SN radiation. The common expedient assumes an overall spherical symmetry of the homologously expanding ejecta, but with an enhanced mixing of all species. This mixing is both macroscopic (material from low velocity is advected out to large velocity, and vice versa) and microscopic (the "advected" material is fully mixed with the material at its new location). We refer to this approach as the "old", standard, mixing technique (see, e.g., Dessart et al. 2015). This approximation is adequate for bolometric light-curve calculations. How-ever, because the composition is altered in a nonphysical manner, the opacity and emissivity of the plasma is no longer accurate, which may impact the SN color, the multiband light curves, and the spectral properties.
The shortcomings of a combined macroscopic and microscopic mixing are best seen in nebular-phase spectra. At late times, the plasma cools primarily through line emission, and their cooling power depends sensitively on the local plasma composition. As is well known (Fransson & Chevalier 1989), and recently demonstrated with detailed calculations (Dessart & Hillier 2020), mixing Ca with O-rich material affects [O i] λλ 6300, 6364 emission. The sensitivity of line strengths to microscopic mixing, vividly illustrated by O i, also holds for other emission lines. While some mixing may occur during the ultimate stages of massive star evolution (see, e.g., Collins et al. 2018), it should not be enforced artificially by the need for a simplistic treatment of mixing driven by numerical convenience. Instead, reliable inferences from nebular-phase spectra require an accurate description of chemical mixing in SN ejecta.   Woosley & Heger (2007) at 345 d (the original mass fraction is shown for 56 Ni), arising from a 15 M progenitor star on the main sequence, and having an ejecta kinetic energy of 1.2 × 10 51 erg, a total ejecta mass of 10.96 M , and 0.083 M of 56 Ni initially. To properly identify the eight shells, we also show the profiles for C, N, Mg, and Ca. Right: Composition profile obtained with the new mixing approach presented in this letter, using M sh = 5.36 M , in which only macroscopic mixing is applied by splitting and shuffling ejecta shells identified in the unmixed model shown at left. In both panels (and in Fig. 2), the origin of the x-axis is at 1.83 M , which corresponds to the mass of the compact remnant.  Counterpart to the right panel of Fig. 1, but now using the old, standard, approach for mixing in which the chemical mixing is microscopic and macroscopic and applied to all species (see, e.g., Dessart et al. 2015). The complex composition profiles of Fig. 1  bookkeeping, the method can be set to preserve mass, energy, match the total decay power, and provide a complete nebularspectrum for a given SN ejecta. But by treating each dominant shell independently of the others, any cross-talk between shells is neglected. We have tried various approaches along this route for Type II-Plateau SNe but all failed to deliver a reasonably looking nebular spectrum. The reason is that, even at 200 − 300 d after explosion, there is considerable reprocessing by the whole ejecta of photons injected below 5000 Å (e.g., photons emitted from the hotter He-rich shell material; see also Jerkstrand et al. 2015). The spectrum formation is a global process which requires the modeling of the full ejecta and associated couplings.

Treatment of chemical mixing: standard and new approaches
In this letter, we present a novel treatment of chemical mixing in spherical symmetry, applicable to all SN types, and which requires no adjustment to the radiative transfer, here performed with the nonlocal thermodynamic equilibrium (nonLTE) radiative transfer code CMFGEN (Hillier & Dessart 2012). The adjustment applies instead to the initial conditions for the coasting ejecta. Since the SN ejecta structure cannot be broken laterally in 1D (by having material of distinct composition coexist at the same velocity), the fundamental idea of the method is to break this structure radially. Taking the unmixed SN ejecta, we select the region that covers from the innermost ejecta layers rich in metals to somewhere beyond the base of the H-rich shell, and bounded by the Lagrangian mass M sh . We then split each of the eight dominant shells (concisely tagged H/He, He/N, He/C, O/C, O/Ne/Mg, O/Si, Si/Ca, and Fe/He) into three equal parts and shuffle these 24 sub-shells within the same region limited by M sh . Because the ejecta is coasting, this shuffling in mass space is equivalent to a shuffling in radial or velocity space. In this brute-force approach, the only requirement for CMFGEN is to increase the number of grid points in order to resolve the rapid variation of composition with depth. Rather than using a total of ∼ 80 grid points for a standard nebular-phase model (as in the simulations presented in Dessart & Hillier 2020), we require ∼ 400 grid points, of which 90% are devoted to resolving the ejecta layers within M sh . A detailed presentation of the procedure is given in the Appendix A.
To illustrate the procedure, Fig. 1 shows the undecayed composition versus Lagrangian mass in the s15A ejecta model at 345 d computed with KEPLER (Woosley & Heger 2007). The eight main shells are indicated with an alternate color of grey and white: starting at the remnant mass of 1.83 M , one goes through the shell rich in Fe/He, followed by the Si/Ca shell, the O-rich shell (which we split into three distinct shells, namely O/Si, O/Ne/Mg, and O/C), the He-rich shell (which we split into two distinct shells, namely He/C and He/N), and finally the H/He shell from the progenitor H-rich envelope (the thin He/N/H shell is treated here as being part of the massive H/He shell) all the way to 12.79M , which is the model total mass at collapse.
Splitting each of the eight shells (of distinct composition; the outermost shell involved in the process is limited to the range 4.28 to 5.36 M , while beyond 5.36 M , the H-rich material is left untouched) in three equal-mass sub-shells, we cycle three times through these eight shells and stack them on top of each other, starting from the ejecta base. The result is similar to what is shown in the left panel of Fig. 1 but now the pattern of eight shells is repeated three times with all eight shells being one third of their original mass. This shuffling thus preserves the ejecta mass and kinetic energy but introduces the macroscopic mixing we need. There is no microscopic mixing between shells.  Fig. 3. Comparison of the optical luminosity (no normalization is applied) for models using the standard ("Mix") and the new mixing technique based on the shuffling of shells ("Shuffle", with various choices for M sh ; the fine, resolved, structure in the spectra is real and caused by the rapid variation in composition in these shuffled-shell models). For model with suffix "Rev", the shell stacking order is reversed. While the different models with shuffling appear similar (highlighting some degeneracy with mixing), there is a strong difference with the results obtained with the standard mixing technique, in particular for Hα and [Ca ii] λλ 7291, 7323.
The ejecta composition produced with the new mixing technique is shown for the main species (i.e., H, He, O, Si, and Fe) in the right panel of Fig. 1. For comparison, we show in Fig. 2 the results obtained with the standard mixing technique: the detailed composition segregation is lost in that case, all species being present throughout the ejecta, exhibiting a broad and smooth profile. The new approach avoids these pitfalls and retains the numerical advantages of spherical symmetry. In addition, the 56 Ni is no longer mixed throughout the ejecta so 56 Ni-rich material has a greater density, causing a greater Fe absorption (this property may be mitigated by the 56 Ni bubble effect, which is the expansion caused by decay heating; e.g., Woosley 1988).
The new mixing technique is flexible since one can design many different shell arrangements to mimic macroscopic mixing. Rather than stacking the sub-shells from the ejecta base, one may start from M sh and progress inwards towards the ejecta base. This choice puts some 56 Ni material at larger velocity (the innermost 56 Ni shell would be located at 5.36 M in Fig. 2, while H-rich material would then be present at the lowest ejecta velocities). We may also increase M sh to allow a greater mixing of the H-rich material with the metal-rich regions. Shells could also be split into more than three parts but thinner shells are harder to resolve numerically. On the other hand, splitting the shells in two parts may produce too little mixing. Finally, we could also mix some 56 Ni-rich material to large velocities, beyond M sh , by merely switching a fraction of the 56 Ni-rich shell in the inner ejecta with an H-rich shell of the same mass in the outer ejecta. Results for some of these configurations are presented below.
While the adopted mixing in Fig. 2 may seem weak, the shuffled-shell structure covers the 4π angle, mimicking the presence of numerous "clumps" of the same composition at the same velocity (or mass coordinate). Additional shifts along different directions, which cannot occur in 1D, are not critical since the γrays from radioactive decay fill rather uniformly the entire mixed region -their typical mean free path of ∼ 5 × 10 15 cm in model s15A at 345 d is comparable to the SN radius. Hence, allowing for more shells or even doing a 3D treatment would merely randomize the location of emission and absorption of the ejecta material, but this would merely produce smoother line profiles while not changing much the physics at work. The impact on total line fluxes and such is expected to be small. Although the mixing approach is distinct from that used in the Monte Carlo code SUMO (Jerkstrand et al. 2011), the resulting 3D structure and its influence on the radiative transfer are similar.
Using the coasting ejecta model s15A described above (see also Woosley & Heger 2007), we implement various types of chemical mixing and investigate the resulting radiative properties with CMFGEN. The CMFGEN calculations assume steady state, a SN age of 345 d, and treat the ejecta composition in detail (Fig. 2). We include up to three ionization stages for H, He, C, N, O, Ne, Na, Mg, Al, Sc, Si, S, Ar, K, Ca, Ti, Cr, Fe, Co, and Ni, and only the 56 Ni-decay chain is considered. The ejecta inside M sh , bounded by velocities between 250 and 2500 km s −1 , is resolved uniformly with 320 points, thus with a fixed resolution of 7 km s −1 . Between 2500 km s −1 and the maximum velocity of 7000 km s −1 , we employ 30 points equally spaced on an opticaldepth scale.

Results
Figure 3 compares the optical spectra (shown as a luminosity L λ , and thus without any scaling) for the CMFGEN calculations based on the standard technique for mixing (i.e., which is both microscopic and macroscopic), and various incarnations of the new mixing technique that shuffles ejecta shells in mass (or in velocity) space. Although all four models have the same 56 Ni mass, the different mixing introduces variations in decay-power absorbed (and thus bolometric luminosity) at the 10% level.
Two striking results emerge from Fig. 3, whose interpretation can be facilitated by inspection of the ejecta and radiation properties for one of the shuffled-shell models shown in Fig. 4. First, the main difference between models is limited to Hα and [Ca ii] λλ 7291, 7323, with little impact on [O i] λλ 6300, 6364. It is complicated to identify the core reason for this because the physics is very nonlinear and complex. With the standard mixing technique, the decay power is absorbed by a nearly homoge- neous metal-rich inner ejecta, from which [Ca ii] forbidden lines perform most of the cooling. In contrast, in the new shuffledshell approach, the [Ca ii] forbidden lines form mostly within the Si-rich layer (Fig. 4). Being of low mass (see Fig. 1), this Si-rich layer absorbs a small fraction of the total decay power, which drastically limits the maximum flux in [Ca ii] λλ 7291, 7323. In the "Mix" model, this maximum flux can be huge, and potentially equal to the total decay power available, because the Ca abundance is allowed to be large throughout the ejecta (Fig. 2).
The second striking result is that the various choices for the shuffling of shells do not lead to much diversity. The spectral properties are therefore weakly sensitive to the details of mixing (as long as the mixing of 56 Ni is comparable; see Dessart & Hillier 2020), which is poorly constrained in SN ejecta. This degeneracy arises because the emergent radiation is primarily dependent on the decay power absorbed by each shell. By allowing for the presence of 56 Ni-rich shells at various locations out to about 2000 km s −1 , nonlocal energy deposition allows decay power to permeate nearly uniformly throughout this extended region (see Fig. 4) -the specific location of the various other shells becomes irrelevant. Furthermore, the key physics operating in each of these shells (opacity, emissivity) is primarily dependent on the shell mass and composition, which is the same in all three incarnations shown in Fig. 3 since it is fixed for a given progenitor (here the unmixed s15A model).
Although one may argue that the standard mixing technique and our new technique yield similar-looking spectra, they are quantitatively different, and the new technique has physical consistency that the standard technique was severely lacking. With the new technique, we can rely on our results, which are no longer biased by an unphysical mixing of species. For example, our new simulations yield very complex temperature and ionization profiles, reflecting the different composition and coolants in the various shells. We also have confidence on the origin of the emission from forbidden lines, while the treatment of the full ejecta takes full consideration of any cross-talk between individual shells. In the shuffled-shell model with M sh = 5.36 M , the Hα line forms throughout the H-rich ejecta layers (down to the innermost layer which is placed around 1400 km s −1 ), the [O i] λλ 6300, 6364 line forms primarily in the O-rich shell, while [Ca ii] λλ 7291, 7323 forms in the Si-rich shell, where Ca is neutral -over-ionization of Ca is in part responsible for the lack of Ca ii emission from the outer H-rich ejecta layers. Figure 5 shows a comparison between our shuffled-shell s15A model with M sh = 5.36 M and the observations of SN 2004et. Model and observation are only offset by 0.1 mag in V-band magnitude after correction for distance and reddening. The agreement in the spectral properties is also satisfactory, suggesting that model s15A is a suitable representation of the ejecta properties of SN 2004et. One discrepancy is the lack of Hα emission at low velocity. When we reverse the shell stacking order, more H is present in the inner ejecta but this causes only a modest increase in the Hα strength (see Fig. 3). Allowance for clumping might help resolve this small discrepancy.

Discussion of pros and cons of the method
The complex 3D structure of core-collapse SN ejecta, which typically exhibit asymmetry from small to large scales, would require high spatial and angular resolution together with 3D nonLTE radiative transfer. This is beyond current capabilities, which warrants the development of simplified approaches. Our approach has several caveats, but these are offset by the benefits.
At intermediate nebular times of a few hundred days, γ-rays have a large mean free path, and hence the arrangement of 56 Nirich material in shells, blobs, or fingers is not crucial for determining the energy deposition. With time passing, γ-ray escape increases until the dominant power source arises from positrons. In between these regimes, the adopted shell arrangement may influence the results, and the limitations of the approach may become more apparent. The transition from full γ-ray trapping to positron escape occurs on a time scale that varies considerably between SN types. For modeling SNe II at ∼ 300 d or SNe Ibc at ∼ 150 d, this is irrelevant so our technique should be robust in these cases. For the faster expanding SN Ia ejecta, γ-ray escape occurs much sooner. However, their chemical segregation is less severe than in core-collapse SNe, with the dominance of 56 Ni-rich material and Si-rich material (unburnt material plays little role at nebular times). Tests are needed to investigate the influence of the adopted shell arrangement for each SN type.
The radiation from one shell is modified as it crosses other shells. Blanketing below 5000 Å will lead to some frequency redistribution, in particular at earlier times. The arrangement in shells may overestimate this reprocessing since emission from inner shells is forced to cross all overlying shells before escape, while a more realistic configuration with fingers or blobs could facilitate this photon escape. But if SN ejecta are made of numerous blobs and fingers covering all sight lines, the difference with our simple radial shell arrangement might be small.
Our approach is fundamentally 1D and the emerging radiation is the same for any distant observer. This assumption will break down in the presence of large-scale asymmetries. Our technique, for example, will not capture the essential physics if the 56 Ni distribution is offset to one side. This could influence where the γ-ray power is absorbed, perhaps favoring some shells over other shells (for example, the Si-rich material might also be preferentially present on the same side as the 56 Ni, while the O-rich material could lie on the opposite side). We note that asymmetry will not affect the total flux in optically-thin lines (for the same decay power absorbed), but only their line profile.
Although presently ignored, we need to allow for compression (clumping) and depression (rarefaction) of the various shells in order to mimic the bubble effect caused by 56 Ni heating and the associated compression of the surrounding material. This will influence the ionization of the gas, for example enhancing the recombination in compressed regions (e.g., Ca in the H-rich material).

Conclusions
We have presented a new technique for the treatment of chemical mixing in SN ejecta for the purpose of radiative-transfer modeling. It relies on the shuffling in mass space (or equivalently velocity space) of the pristine ejecta shells produced in a 1D explosion model having reached homologous expansion. Each shell retains its original composition, thus no microscopic mixing is introduced, while the shuffling of shells mimics the process of macroscopic mixing in velocity space. Because the method is 1D, it is amenable to detailed steady-state nonLTE radiative transfer calculations in grid-based codes. The need for several hundred grid points increases the memory requirements and the computing times by a factor of a few tens. To circumvent this hurdle, the calculation for a given explosion model is done in steps. First, we obtain a converged model using the standard mixing approach. The solution is then used for the initial guess on level populations, temperature, electron density etc for a shuffled-shell model using a Doppler line width set to 50 km s −1 (for discussion on the role of the line Doppler width, see Dessart & Hillier 2020). Once converged, this model is then rerun with a line Doppler width fixed at 10 km s −1 . The whole process takes several days, which is affordable for a steady-state calculation : only one model (rather than a full time-sequence of models) is needed to produce a spectrum at a given SN age. The method requires some care but is straightforward to implement and can be applied to any SN type, either a core collapse or a thermonuclear explosion.
With this more physically-consistent mixing technique, it becomes possible to perform reliable nebular-phase calculations for all SN types with CMFGEN. It is also of interest to test the impact of this different mixing technique during the photospheric phase, for example to assess the role of He i line excitation in SNe Ib, or the behavior of H recombination during the plateau phase of SNe II-Plateau.
To prepare the initial model for CMFGEN, we used the unmixed s15A ejecta model at 345 d computed with KEPLER (Woosley & Heger 2007), whose undecayed composition is shown in the left panel of Fig. 1. So, unlike numerous studies in which we used a combination of our own stellar evolution and explosion models (as in, for example, , the initial conditions for the present CMFGEN calculations are based exclusively on that KEPLER model s15A at 345 d. Obviously, the new mixing technique can be applied to any unmixed coasting ejecta model, whatever its origin, and could also be used for calculations at other times during the coasting phase (i.e., when dynamical effects have abated).
At such a late time, the ejecta is essentially in homologous expansion (the current radius R of each mass shell is huge relative to its initial radius), that is the velocity V of each mass shell is simply V = R/t, where t is the time since explosion. For the CMFGEN input ejecta model, we enforced homologous expansion exactly, as usual (this requirement arises from the numerical method; see details in Hillier & Dessart 2012).
Before applying our mixing technique to the ejecta, we also suppress the density jump caused by the reverse shock formed at the He-core edge by forcing the density to be constant throughout the region we intend to mix (i.e., this density is just the total mass divided by the total volume of the ejecta between the innermost ejecta layer and M sh ). Changing the inner density profile of 1D explosion models is routinely done by other radiative transfer modelers (see, for example, Fig. 1 of Jerkstrand et al. 2012, where the density is similarly forced to be constant throughout the inner ejecta). This density jump is an artifact of the imposed spherical symmetry. Indeed, the Rayleigh-Taylor instability that grows at the He core edge after shock passage completely erases this jump (see, for example, Fig. 8 of Utrobin et al. 2017 for a comparison between 1D and 3D predictions). Obviously, after changing the density in the inner ejecta regions, we have to recompute the Lagrangian mass, which is trivial. This change in density structure at low velocity affects the velocity versus mass, but only modestly since the inner ejecta contains a small fraction of the total kinetic energy. Indeed, the total ejecta kinetic energy increases by 3%, while the kinetic energy within M sh , which is originally only 4% of the total, is raised by 80%. The velocity and density before and after this process are shown, versus Lagrangian mass, in Fig. A.1. This change in density profile (i.e., both versus radius and Lagrangian mass), is independent of the mixing step that we describe next.
We then selected the region of the ejecta that we wished to mix. In the example shown in the right panel of Fig 1, we selected the region between the innermost layers at a Lagrangian mass of 1.83 M (i.e. the remnant mass) and the Lagrangian mass M sh , here chosen to be 5.36 M . The heart of the method is to swap a chunk δM at some depth M 1 with the same δM at some other depth M 2 , both within M sh . Doing this obviously conserves mass (and the yields of all isotopes) since we are merely swapping chunks of identical mass between two different locations. It requires no change to the density profile (one does not even need to know the density structure to do this operation). The velocity structure is obviously unchanged in the process since it is uniquely defined by homology. Because the velocity increases with Lagrangian mass, this shuffling of shells leads to a shuffling in velocity space, which is precisely what we aim to achieve so that, for example, 56 Ni becomes present at higher velocities (Wongwathanarat et al. 2015).
There are many ways to design the shell swapping. One could, for example, choose to start at the ejecta base and swap a given δM with the same δM located at M sh , and then proceed until we reach M sh /2 (i.e., halfway through the regions to mix). But this is not desirable because some shells are thin and others thick. So, for example, with an adopted ∆M of a few 0.1 M , all the 56 Ni would end up just below M sh . This situation would capture very poorly the predictions for 56 Ni mixing, which suggest the ubiquitous presence of 56 Ni from low to high velocities (Wongwathanarat et al. 2015). Similarly, these simulations suggest that blobs and fingers of distinct composition (i.e., not just the 56 Ni material from the 56 Ni-rich shell) are thoroughly mixed within an extended region of the ejecta.
Hence, our choice of shell shuffling was driven by the need to shuffle all mass shells, so that material from the thin 56 Nirich shell and the thicker O-rich shell get mixed by comparable amounts. This way, fractions of the corresponding shells are made present in multiple regions of the ejecta. Our approach was therefore to split each of the eight dominant shells (concisely tagged H/He, He/N, He/C, O/C, O/Ne/Mg, O/Si, Si/Ca, and Fe/He) into three equal parts and shuffle these 24 sub-shells within the same region limited by M sh . Starting from the ejecta base, we cycle three times through these eight shells, placing each sub-shell on top of the previous one. When we are done, we have reached the Lagrangian mass M sh , beyond which the ejecta is unchanged from the original KEPLER model. The result is similar to what is shown in the left panel of Fig. 1 but now the pattern of eight shells is repeated three times with all eight shells being one third of their original mass but appearing at three different locations (this is most easily seen for the massive O-rich shell, colored in olive in Fig. 1). This shuffling (which is equivalent to swapping chunks of mass between two different locations) thus preserves the ejecta mass and kinetic energy but introduces the macroscopic mixing we need. There is no microscopic mixing between shells.
To facilitate the procedure and the CMFGEN modeling, each of the eight shells are made homogeneous initially to avoid having additional composition gradients within shells. A very small mixing is also applied to all species to reduce the composition gradient at each shell interface in the original 1D model (see left panel of Fig. 1).