© A. Sarangi 2022

1 Introduction

The rapidly evolving ejecta of core-collapse supernovae (SNe) are known sites for the synthesis of cosmic dust (Bouchet & Danziger 1993; Wooden et al. 1993; Dwek 2006; Szalai & Vinkó 2013; Matsuura 2017; Sarangi et al. 2018b), and are often proposed to be the primary source of dust in high-redshift galaxies (Dwek et al. 2007; Zhukovska et al. 2008; Cherchneff & Dwek 2009; Gall et al. 2011; Gall & Hjorth 2018). The mass and properties of dust present in SN ejecta are quantified by analyzing the spectral energy distribution (SED) at various wavelengths, where the mid-infrared (IR) and far-IR part of the SED is often dominated by emission from dust formed in the ejecta. Observations of several SNe in the mid-IR wavelengths over the past couple of decades found the upper-limit of dust formed in the first 1000 days of the SN explosion to be constrained to within a few times 10^-3 Ms (Sugerman et al. 2006; Ercolano et al. 2007; Kotak et al. 2009; Wesson et al. 2015; Priestley et al. 2020; Szalai & Vinkó 2013; Szalai et al. 2019).

In addition to the emission in the IR, the presence of dust in the ejecta also gives rise to an asymmetry between the red-shifted and blueshifted emission lines in the optical and near-IR wavelengths, since the radiation coming toward us from the far side of the ejecta is subjected to increased absorption and scattering by the dust lying along its path (Bevan & Barlow 2016). The quantity and composition of dust that is required to account for the progressive blueshift of the [OI] doublet in SN 1987A (Lucy et al. 1989) was estimated by Bevan & Barlow (2016). Combining this to the fitting of the mid-IR spectra using a radiative transfer model, it was suggested that the mass of dust formed in the ejecta of SN 1987A was no more than 3 × 10⁻³ M_⊙ at 1000 days, with a predominantly carbon-rich composition (Wesson et al. 2015; Wesson & Bevan 2021; Bevan 2018).

Interestingly, a much larger mass of dust, of about ~0.5 M_⊙, was found to be present in SN 1987A after 23 years of its explosion, derived from the fluxes in the far-IR and submillimeter wavelengths (Matsuura et al. 2011, 2015; Indebetouw et al. 2014). This trend is also evident in the cold ejecta of other SN remnants (SNRs) such as Cas A (Barlow et al. 2010; Arendt et al. 2014; De Looze et al. 2017; Niculescu-Duvaz et al. 2021) and the Crab Nebula (Gomez et al. 2012; Owen & Barlow 2015; De Looze et al. 2019). The apparent dichotomy between the estimated dust masses in the mid-IR and submillimeter wavelengths is perceived as evidence of a gradual increase in dust masses due to continuous dust formation over a couple of decades (Gall et al. 2014; Wesson et al. 2015; Gall & Hjorth 2018).

Dwek & Arendt (2015) and Dwek et al. (2019) propose an alternative scenario where large masses of dust condense within a couple of years after the explosion, when the ejecta are still hot. However, the large masses are not reflected in the SED due to the large opacities of the ejecta in the mid-IR wavelengths. Dwek & Arendt (2015) find a reasonably good fit to the spectra of SN 1987A at days 615 and 775 (Wooden et al. 1993), with a rather large dust mass of ~0.45 M_⊙, dominated by silicates, formed in dense clumps. The early synthesis of O-rich dust, such as silicates, complies well with the characteristic 9.7- and 18-µm feature of silicate dust in the interstellar medium (ISM; van Breemen et al. 2011), which requires high temperatures, typically larger than 1000 K, during its formation (Dwek et al. 2019). In addition, the sharp decrease in intensity of the [Mg I] or [Si I] emission lines (Danziger et al. 1991), concurrent with the increase in the mid-IR emission after day 530, also supports the formation of Mg-silicates in abundance in the gas (Dwek & Arendt 2015). On the contrary, a rapid formation of dust leading to an early saturation is not congruent with the observed trend of a continuous increase in the blueshift of the emission lines (Bevan & Barlow 2016), if a spherical shape is assumed for the ejecta.

Dust forms in the SN ejecta through simultaneous phases of nucleation and condensation (Cherchneff & Lilly 2008; Sarangi et al. 2018b). Dust formation models, in general, support the scenario of efficient dust formation in the first 2–3 yr after the explosion (Nozawa & Kozasa 2013; Sarangi & Cherchneff 2013, 2015; Sluder et al. 2018), since the gas densities in the ejecta are high at the early times, while the warm gas temperatures allow nucleation pathways to overcome activation energy barriers. The final dust masses predicted by the models depend on the assumed density of the gas and on the condensation scheme used.

To account for the extreme nature of the environments typical of SN ejecta, recent models have relied on a chemical kinetic approach to address the chemistry of molecule and dust formation in SNe (Cherchneff & Dwek 2009; Sarangi & Cherchneff 2013, 2015; Sarangi et al. 2018b; Sluder et al. 2018). The most recent models of both Sarangi & Cherchneff (2013, 2015) and Sluder et al. (2018) adopt a similar network of reactions for the gas phase chemistry. For the growth of grains, Sarangi & Cherchneff (2015) assume growth via coagulation and rule out accretion on the surface of the small grains, while Sluder et al. (2018) consider both the processes. Sarangi & Cherchneff (2015) model homogeneous as well as clumpy ejecta, where the geometry is assumed to consist of annular shells with unique abundances. In Sluder et al. (2018), the ejecta is considered as a three-phase system made up of ⁵⁶Ni-rich blobs, the swept-up mass by the blobs in the form of a shell, and the ambient gas. In strong contrast with the observational trends, for SN 1987A, the clumpy model of Sarangi & Cherchneff (2015), and Sluder et al. (2018), both predict a formation of more than 0.01 M_⊙ of dust within a year of the explosion.

When modeling the fluxes from the ejecta dust (Dwek & Arendt 2015; Wesson et al. 2015), the distribution of dust within the He core is always assumed to be uniform or to scale with the densities of the clumps. This therefore ignores the fact that, from the center of the ejecta to the outer edge along any random line of sight, different clumps may have a different efficiency, composition, and timescale for dust production, affecting the overall opacities and fluxes. In addition, the evolution of dust composition, in other words, how the relative masses of each dust type can evolve differently, is also not accounted for. The impact of these factors are further compounded when the geometric shape of the ejecta is not a spherical, isotropic one; for instance, the elongated structure of the evolving ejecta of SN 1987A (Matsuura 2017) is estimated to fit well with an ellipsoidal geometry (Kjær et al. 2010).

The advanced radiative transfer models that address the nebular phase of SNe enable us to better quantify the ejecta conditions, such as the distribution of clumps or the distribution of gas temperatures within the He core (Dessart et al. 2018); (Dessart & Hillier 2020), with a much improved understanding of the molecular cooling, especially by CO molecules (Liljegren et al. 2020). Based on that, in this paper, we adopt a more realistic and higher-resolution model for the evolution of the physical conditions in SN ejecta. Applying the gas phase chemistry developed by Sarangi & Cherchneff (2013, 2015) to the upgraded model of the physical conditions, we derive the distribution and evolution of dust masses in clumpy SN ejecta in the spatial as well as temporal space. Using this generic distribution of dust, we calculate the opacities of the ejecta, and derive the resulting luminosities and fluxes. Having SN 1987A in mind, we also present a case where an ellipsoidal shape is assumed for the ejecta, to estimate the dust distribution, column densities, and fluxes, as a function of the viewing angle. This is the first study to account for the effect of anisotropy on the chemistry of dust formation. Moreover, this is also the first attempt to calculate the fluxes from the dust formed in the SN ejecta using a generic, nonuniform distribution of dust.

The paper is arranged as follows. In Sect. 2 we explain the evolution of physical conditions in the ejecta, in terms of gas densities (Sect. 2.1), gas temperatures (Sect. 2.2), projected dust temperatures (Sect. 2.3), and the velocity-dependent stratification (Sect. 2.4). Following that, Sect. 3 deals with the chemistry of dust formation, and the distribution and evolution of different dust species. In Sect. 4, we discuss the opacities resulting from the dust formed in the ejecta and the expected SEDs from the distribution at various epochs. Section 5 presents the case study of dust formation in anisotropic ejecta of ellipsoidal shape. In Sect. 5.1, we compare the dust distribution and masses between the isotropic and anisotropic cases, while Sect. 5.2 discusses the resulting SEDs as a function of viewing angle. After that, in Sect. 5.3 we present the predicted fluxes at days 615 and 775, in comparison to SN 1987A at the same epochs. Finally, we summarize and discuss our results in Sects. 6 and 7, respectively.

2 Physical conditions of the ejecta

In this section, we describe the physical conditions of SN ejecta with respect to the parameters adopted in this study. After the explosion, typical SN ejecta evolved through a free-expansion phase (homologous expansion), characterized by v ≅ r/t. The ejecta are assumed to be made up of a core of uniform density, and an envelope where the density falls rapidly as a power-law (Truelove & McKee 1999; Dwarkadas 2007). Following the model of Truelove & McKee (1999), we define the velocity of the ejecta core, v_c, as $$v_c={\left[\frac{10\left(n-5\right)}{3\left(n-3\right)}\frac{E}{M_{\text{ej}}}\right]}^{1/2},$$(1)

where E is the explosion energy of the SN, M_ej is the total mass of the ejecta, and n is the power-law exponent that characterizes the density of the ejecta envelope. We model SN ejecta whose progenitor was a 20 M_⊙ main-sequence star, using the post-explosion yields from Rauscher et al. (2002)¹. The mass of the ejecta, M_ej, in our model, is 12 M_⊙ (of which the He core is 4.65 M_⊙), and the explosion energy is taken as 10⁵¹ ergs, with reference to SN 1987A (Ensman & Burrows 1992). The power law for the envelope, n, is taken as 12, typical for red supergiants (Moriya et al. 2013). Using these in Eq. (1), v_c is calculated to be 3300 km s⁻¹. The outer edge of the He core, which is the main site of dust formation in the ejecta, has a velocity of 2650 km s⁻¹. Therefore, the He core is entirely confined within the core of the ejecta and we only focus on that for the purpose of this study.

Fig. 1 Mass and density distribution of SN ejecta are presented in the velocity space. Top left: enclosed mass of the ejecta as a function of the ejecta velocity is shown. The He, O, and Si-cores are marked on the figures, with the respective mass-velocity boundaries (see Sect. 2). Top right: fractional abundances of the elements (H, He, C, O, Ne, Mg, and Si), taken from Rauscher et al. (2002), are shown as a function of ejecta velocity, where the mass-velocity relation can be derived from the figure on the left, as shown in Sect. 2. Bottom left: distribution of gas (number) densities in the ejecta, as the function of velocity at 100 days post-explosion (see Sect. 2.1). The densities for uniform ejecta are shown for reference, but we only model the scenario of clumpy ejecta. Bottom right: volume filling factor for the clumps are shown with reference to Eq. (5).

2.1 Gas densities in the clumps

We assume that the ejecta are clumpy and the gas densities in any clump (ρ_g) are defined by a homogeneous density (ρ_c) and a volume filling factor (f_V). Following the prescription of Truelove & McKee (1999), the density of the ejecta core is given by $${\rho }_h=\frac{3}{4\pi }\frac{n-3}{n}\frac{M_{ej}}{v_c^3}{t}^{-3},$$(2)

and therefore the density of a clump is simply $${\rho }_g={\rho }_h{f}_V^{-1}.$$(3)

Rauscher et al. (2002) provide the abundances of all the elements (and their isotopes) in closely stratified layers, as a function of the enclosed mass of the ejecta below that layer (assuming spherical geometry). For a freely expanding ejecta, where the velocity of a layer is proportional to r, the enclosed mass can therefore be translated to the velocity of each layer using Eq. (2), as $$v={\left[\frac{3}{4\pi }\frac{m_{\text{enc}}}{t^3{\rho }_h\left(t\right)}\right]}^{1/3}={\left[\frac{n}{n-3}\frac{m_{\text{enc}}}{M_{\text{ej}}}\right]}^{1/3}{v}_c,$$(4)

where m_enc is the enclosed mass below a given layer. Using this, in Fig. 1, the relation between the enclosed mass and the velocity is shown, and the abundances of elements, which are relevant for dust formation, are presented as a function of velocity.

For SN ejecta, it is often estimated that the degree of clumpiness increases toward the inner parts of the He core, which is synonymous with smaller velocities (Wongwathanarat et al. 2015; Dessart et al. 2018). To quantify that, we use the relation between the velocity and the volume filling factor given by Dessart et al. (2018) in our study, so that $$f_V=1+\left(f_{\min }-1\right)e^{-{\left(v-v_{\min }\right)}^2/v_c^2}.$$(5)

The minimum value for the filling factor, f_min, was taken as 0.01, with reference to studies on SN 1987A (Jerkstrand et al. 2011; Ercolano et al. 2007), and v_min corresponds to the velocity of the inner edge of the so-called Si core, which in this case is ~750 km s⁻¹ (see Fig. 1, top right). The volume filling factor is shown as a function of velocity in the bottom right panel of the same figure.

The total number density of the gas in the clumps, calculated using the abundances of elements, the gas densities, and the filling factors, are shown in the bottom left panel of Fig 1. The densities are consistent with those predicted by the 3D SN explosion models (Wongwathanarat et al. 2015; Utrobin et al. 2017), when they are extrapolated to the nebular phase.

Due to the much higher densities in the clumps, compared to the inter-clump medium, we assumed that the formation of dust occurs within the clumps only. Therefore, in the following text, where we deal with the dust formation chemistry and the corresponding dust opacities, the inter-clump region was ignored.

2.2 Gas temperatures

The temperature of the gas in the ejecta evolves through collisional line cooling, initially by atoms and ions, and later by molecules (Kozma & Fransson 1998; Liu et al. 1992; Liu & Dalgarno 1995; Liljegren et al. 2020). Given the rapid evolution of the ejecta, along with the presence of nonthermal processes induced by the radioactivities, it is not straightforward to define an analytical expression for the evolution of gas temperatures.

The cooling of a parcel of gas is often controlled by abundances of the elements in that parcel. Based on the most abundant elements, the He core of the ejecta can be overall subdivided into zones, namely the Si/Fe zone, the O/Si/Mg zone, the O/C zone, and the He/C zone (in order of increasing velocity). The extent of these zones within the He core are shown in Fig. 1 (top left panel). CO molecules form abundantly in the O/C zone (Sarangi & Cherchneff 2013; Liljegren et al. 2020), which is responsible for the rapid cooling of the gas in that region. As Liu & Dalgarno (1995) find, for SN 1987A the temperature of the O/C zones at day 800 is close to 700 K, while at the same time the O/Si zone is still hot with temperatures of ~2200 K.

According to Kozma & Fransson (1998), who calculated the evolution of temperatures in each of these zones between day 200 and day 2000, the zones undergo a rapid cooling phase after a few hundred days of evolution, followed by another phase of leveling-off of the temperature. The inner layers experience the rapid cooling phase earlier, while for the outermost He zone, the rapid cooling phase almost never transpires. This is not completely unexpected, since the heavier elements (inner layers have heavier elements) are characterized by faster radiative cooling rates (Sutherland & Dopita 1995; Sarangi & Slavin 2022).

Dessart & Hillier (2020) calculate the distribution of temperature in the SN ejecta (type II) at a post-explosion epoch of 300 days, as a function of velocity. Their results are also in agreement with Kozma & Fransson (1998), showing that the ejecta at a given epoch are cooler toward the inner layers and the temperature increases with velocity.

For the purpose of this study, we obtained the temperatures as a function of time and abundances from Kozma & Fransson (1998), and as a function of velocity from Dessart & Hillier (2020); combining them, we constructed the matrix of gas temperatures in the ejecta as a function of velocity and time. For the CO forming regions, we consider the temperature to evolve as suggested by Liljegren et al. (2020) (consistent with Liu & Dalgarno 1995) once the CO molecules start to form in the ejecta.

The range of gas temperatures that corresponds to each zone and their evolution from day 100 to day 3000 is shown in Fig. 2. The O/Si zone shows the largest stretch in the range of temperature, as expected, since it is spread over a wide velocity range (see Fig. 1, top left). The He/C zone, despite being larger than 1 M_⊙ in mass, is spread over a relatively small velocity range, given that it is at a larger radius. Therefore, the spread of gas temperature in this zone (at any given epoch) is also limited.

Fig. 2 Range of gas temperatures in different zones and their variation with post-explosion time is shown (see Sect. 2.2 for details). Even though our model does not depend on these coarse zone divisions of the ejecta, this is particularly helpful to understand how the gas temperatures evolve.

2.3 Dust temperatures

Once the dust grains form in the ejecta through nucleation and coagulation (Sarangi & Cherchneff 2015), the temperatures of the grains are not necessarily the same as the kinetic temperature of the ambient gas, since the dust grains are efficient in cooling down via emission, lowering their internal energy. The temperatures of the grains determine their IR fluxes, as well as their growth (in size) via coagulation. Our approach to determine dust temperatures is similar to that of Sluder et al. (2018). However, in this study, we also calculate the evolution and distribution of dust temperatures in the velocity space for each clump in the ejecta, for four different dust species (silicates, amorphous carbon, alumina, and silicon carbide).

A dust grain in the SN ejecta is heated radiatively by the UV photons produced as a result of the decay of radioactive isotopes. In addition, dust is also heated through collisions with the ambient gas, while it cools down via emission. When radiatively heated, the rate of heating is simply the incident flux on a dust grain multiplied to the absorption cross-section, integrated over all the wavelengths (Temim & Dwek 2013). Then again, flux is basically the product of the energy density (at a given radius, at a given time) and the velocity. If we take the number density of photons with energy E at a given time and velocity to be n_γ(E, v, t), then the radiative heating rate (H_r) of a dust grain of species S, of radius a, can be expressed as $$H_{S,r}\left(a,v,t\right)={\displaystyle {\int }_0^{\infty }{n}_{\gamma }\left(E,v,t\right)Ec\pi a^{2}Qs\left(E,a\right)\text{d}E},$$(6)

where Q_S is called the absorption coefficient of the grain for incident energy E and size a. The mass absorption coefficient k(λ, a) (=3Q/4aΩ, Ω being the solid-state density of a grain) of silicates (Ω = 3.3 g cm⁻³), amorphous carbon (Ω = 2.2 g cm⁻³), alumina (Ω = 3.95 g cm⁻³), and silicon carbide (Ω = 3.21 g cm⁻³) are shown in Fig. 3 for a = 0.01 µm (Weingartner & Draine 2001; Koike et al. 1995; Zubko et al. 2004; Pegourie 1988).

As assumed by Sluder et al. (2018), n_γ(E) can be expressed as a Gaussian distribution with a mean energy E₀ and deviation σ₀. Then the energy density (at a given v and t) can be expressed as $$\begin{array}{lll}{u}_{UV}\left(v,t\right)\hfill & =\hfill & D{\displaystyle {\int }_0^{\infty }{e}^{-{\left(E-E_0\right)}^2/2{\sigma }_0^2}E\text{d}E}\hfill \\ \hfill & =\hfill & D{\sigma }_0^2{e}^{-E_0^2/2{\sigma }_0^2}+D{\left(2\pi \right)}^{1/2}{\sigma }_0{E}_0,\hfill \end{array}$$(7)

where D is a normalization coefficient, which is a function of v and t. Taking E₀ = 4.43 eV and σ₀ = 1 eV, with reference to Jerkstrand et al. (2011) and Sluder et al. (2018), the first term of the above expression is much smaller than the second term, and thus can be ignored.

The UV luminosity L_UV, responsible for the heating of the dust grains, is a fraction of the total radioactive energy deposited in the ejecta (L_X) at time t, which is converted to UV photons. The coefficient for the fractional conversion, α, is taken as 0.35 (Sluder et al. 2018; Kozma & Fransson 1992), so that L_UV(t) = αL_X(t). For a thin shell of width ∆r at radius r (=vt) at time t, we can write the energy density as $$u_{\text{UV}}\left(v,t\right)=\frac{\alpha L_X\left(t\right)}{4\pi r^2\text{Δ}r}\frac{\text{Δ}r}{c}=\frac{\alpha L_X\left(t\right)}{4\pi r^{2}c}.$$(8)

Combining Eqs. (7) and 8, we have $$D\left(v,t\right)=\frac{\alpha L_X\left(t\right)}{2^{5/2}{\pi }^{3/2}c{v}^2{t}^2{E}_0{\sigma }_0}.$$(9)

The energy deposited in the ejecta through radioactivities was calculated by Cherchneff & Dwek (2009) following the formalism of Woosley et al. (1989), which was later applied to SN 1987A by Sarangi & Cherchneff (2013) and Sluder et al. (2018). In an identical approach, Seitenzahl et al. (2014) calculated the energy deposited in the ejecta and the masses of all the radioactive species, in accordance with the light curve of SN 1987A until the late epoch of day 4334 post-explosion. For any radioactive species i with total mass m_i produced at the time of explosion, we can express the rate of energy deposition in the ejecta through its radioactive decay (Woosley et al. 1989; Cherchneff & Dwek 2009; Seitenzahl et al. 2014; Sluder et al. 2018) as $$\begin{array}{lll}{L}_{Xi}\left(t\right)\hfill & =\hfill & \frac{m_i{N}_A}{A_i}\frac{e^{-t/{\tau }_{Xi}}}{{\tau }_{Xi}}\left(q_i^l+f\left(t\right)q_i^{\gamma }\right)\hfill \\ f\left(t\right)\hfill & =\hfill & 1-e^{{\tau }_{\gamma }\left(t\right)},{\tau }_{\gamma }\left(t\right)=\frac{3}{4\pi }\frac{{\kappa }_i^{\gamma }{M}_{\text{core}}}{v_c^2{t}^2},L_X={\displaystyle \sum _i{L}_{Xi},}\hfill \end{array}$$(10)

where A_i is the mass number, τ_Xi is called the e-folding decay time, and ${\kappa }_i^{\gamma }$ is the average γ-ray mass-absorption coefficient of the radioactive isotope i. N_A is the Avogadro’s number, τ_γ(t) is the γ-ray optical depth at time t, and $q_i^l$ and $q_i^{\gamma }$ are the energies of photons per decay, produced by the charged leptons and γ-rays, respectively. The mass of the ejecta core, M_core, derived from Eq. (4) (using v = v_c, n =12 and M_ej = 12 M_⊙) is about 9 M_⊙.

We chose the best-fit scenario of Seitenzahl et al. (2014) for the masses of the radioactive species ⁵⁶Co,⁵⁷Co, ⁶⁰Co, and ⁴⁴Ti, listed in Table 1, along with all the parameters from Eq. (11). With that, we have now defined all the quantities required to derive the heating rate of a dust grain, given by Eq. (6), as a function of the velocity of the ejecta where it is located, and the post-explosion time.

Besides radiative heating, dust grains are also subjected to heating by the ambient gas, collisionally. The collisional heating rate of a dust grain was calculated by Dwek (1987); Hollenbach & McKee (1979), and applied to SN environments by Sarangi et al. (2018a). The rate of heating due to collision (Burke & Hollenbach 1983; Sluder et al. 2018), Hc, can be expressed as $$\begin{array}{lll}{H}_c\left(a,v,t\right)\hfill & =\hfill & n_g\left(v,t\right)\times \pi a^2\times {\left[\frac{8k_{\text{B}}{T}_{\text{g}}\left(v,t\right)}{\pi \langle \mu \left(v\right)\rangle }\right]}^{1/2}\hfill \\ \hfill & \hfill & \times f_{\text{ed}}\times \left(2k_{\text{B}}{T}_{\text{g}}\left(v,t\right)-2k_{\text{B}}{T}_{\text{d}}\left(v,t\right)\right)\hfill \\ f_{\text{ed}}\hfill & =\hfill & 0.1+0.35e^{-{\left[\left(T_{\text{d}}+T_{\text{g}}\right)/500\right]}^{1/2}},\hfill \end{array}$$(11)

where k_B is the Boltzmann constant, µ(v) is the mean particle mass of the gas in the ejecta, for the velocity bin v, and T_g and T_d are the gas and dust temperatures, respectively.

Comparing the collisional and radiative heating rates from Eqs. (6) and (12), we find that the rate of radiative heating is far more dominant at all times and at all velocities, and therefore we can safely ignore the impact of collisional heating. In Fig. 4, we show the comparison of the radiative heating rate with the maximum rate of heating due to collision, for all four dust species, at day 800.

The rate of radiative cooling (Temim & Dwek 2013) of a dust grain of species S (with mass absorption coefficient k_S), at temperature T_S, which is analogous to the energy lost by the grain through emission, is given by $$C_S\left(a,T_{\text{d}}\right)=4m_S\left(a\right){\displaystyle \int \pi B_{\lambda }\left(\lambda ,T_S\right)k_S\left(\lambda ,a\right)\text{d}\lambda },$$(12)

where m_S is the mass of the dust grain and B_λ is the Planck function per unit wavelength.

Comparing the radiative heating from Eq. (6) to this cooling rate, the temperature of the dust grains are determined for each velocity bin and time. Figure 5 shows the obtained distribution of dust temperatures for silicates, alumina, amorphous carbon and silicon carbide for a grain size of 0.01 µm. To be clear, this map of dust temperatures in the complete v, t space only reflects the temperatures assuming that dust species were hypothetically present there; in other words, this map does not provide any information on the time and location of dust formation from the chemistry (see Sect. 3 for that).

When accounting for the dust temperatures, we did not include the effect of secondary heating by the IR radiation produced by the dust itself. Given that the coefficient of absorption of IR radiation is much smaller than that of UV (seen in Fig. 3), we do not presume this to be significant; we shall use a more sophisticated radiative transfer model in the future to verify this.

It is important to note that the rapidly evolving SN ejecta in their nebular phase are estimated to not be in chemical or collisional equilibrium (Donn & Nuth 1985; Cherchneff & Dwek 2009; Sarangi et al. 2018b). However, they are still likely to be in radiative equilibrium given the high velocity of the photons; since the dust temperatures are essentially controlled by radiation only, the temperatures derived above are justified.

Fig. 3 Mass absorption coefficients for silicates (Weingartner & Draine 2001), alumina (Koike et al. 1995), amorphous carbon (Zubko et al. 2004), and silicon carbide (Pegourie 1988) are presented, for a grain size of 0.01 µm. For alumina at wavelengths smaller than 0.3 µm and for silicon carbide at wavelengths smaller than 0.2 µm, the values are extrapolated due to lack of data. In addition, for these two dust types, we do not have data for various grain sizes, and therefore these are generalized values for all sizes.

Fig. 4 Radiative heating rate (Eq. (6)) for the four different dust species are shown compared to the (maximum) heating rate of the dust grains via collision with the ambient gas (Eq. (12)), for day 800 and a grain size of 0.01 µm. As it shows, the dust temperatures are controlled by the radiative heating only; that is also the case for other epochs and grain sizes as well.

2.4 Stratification

In Sects. 2.1–2.3, we have explained the physical conditions of the ejecta and their evolution. It is clear that in different parts of the ejecta, characterized by different velocities, there will be variations in the densities of the clumps, as well as in the gas and dust temperatures, all of which are key physical conditions that are critical to the formation and growth of dust grains.

In our model, we stratify the ejecta into 670 thin shells, following the same stratification as that of the post-explosion nucleosynthesis data of a 20 M_⊙ main-sequence star of Rauscher et al. (2002)². Out of those, the He core is confined within the first (inner) 450 shells and we uniquely follow the dust formation chemistry in those shells. In other words, it can be visualized as a set of 450 different groups of clumps, distinguished by their abundances, densities, and gas and dust temperatures. For reference, we identify the clumpy shells by their velocities, since that is unique to each shell throughout the time of their evolution. In the following section, we discuss the chemistry of dust formation and the net chemical yields in light of this model.

3 Dust formation in the ejecta

Dust formation in the SN ejecta is controlled by simultaneous phases of nucleation and condensation (Cherchneff & Dwek 2009, 2010; Sarangi & Cherchneff 2013, 2015). The chemistry of the hot and dense SN ejecta was modeled by Sarangi & Cherchneff (2013), where all the possible chemical reactions and pathways relevant to the atoms, molecules, and molecular clusters present in the gas phase were accounted for. In Sarangi & Cherchneff (2015), the growth of dust grains through stochastic coagulation and coalescence was modeled. We adopted the same formalism in this study, and no changes were made to the chemical reaction network and pathways of dust growth, which were developed by Sarangi & Cherchneff (2013, 2015); for a complete list of all chemical reactions, we refer the reader to Cherchneff & Dwek (2009) and Sarangi & Cherchneff (2015). This means that, by applying the chemical network in this study, we derive the chemical yields in the 450 different groups of clumps, which are identified by their unique abundances of species and unique evolution of gas and dust conditions.

The gas phase nucleation processes lead to the formation of molecules and molecular clusters (Sarangi & Cherchneff 2013). As the He cores of the ejecta are predominantly rich in O, a lot of O-rich diatomic molecules are synthesized in the nebular phase of the SN. Specifically, molecules such as CO, SiO, AlO, O₂, and SO are efficiently produced. Following their formation after about 200 days from the time of explosion, CO molecules retain their large abundances in the ejecta and gradually increase in mass. The total mass of CO molecules formed in the ejecta is 0.37 M_⊙. On the contrary, molecules such as SiO or AlO are increasingly depleted into dust grains with time. Using the formalism developed in Sarangi & Cherchneff (2015), in this study we follow the profiles of the four most abundant dust species in the ejecta, namely, magnesium silicates (of chemical form [Mg₂SiO₄]_n), alumina (of chemical form [Al₂O₃]_n), amorphous carbon (of chemical form [C₂₈]_n) and silicon carbide (of chemical form [SiC]_n; n is an integer generally larger than four). More details on the types of dust species relevant for SN ejecta, and their individual properties, can be found in Sarangi & Cherchneff (2015). For any clump of a given abundance, the required parameters to study the chemistry of dust formation are: the evolution of gas density, the gas temperature, the rates of all possible gas phase reactions, dust temperatures, the solid state density of the dust grains, and the coefficients for the Van der Waals interaction.

The smallest grains of dust, or rather the precursors, are synthesized through gas phase chemical reactions, where the temperatures are expected to be analogous to the gas temperatures. They coagulate and grow larger than 10 Å, and are labeled as dust grains. Coagulation is not well understood in high temperatures, and is generally perceived to be inefficient when the kinetic energies of collisions are large (Hou et al. 2022); we have only assumed coagulation when the gas temperatures are lower than 2000 K.

In Fig. 5, where the dust temperatures are shown, the contours (shown by solid and dashed black lines) refer to the region in velocity space where a particular dust type is present at a given epoch. The left boundary of that contour (shown as a solid black line) is therefore the time when dust formation commences in that velocity bin. This is further understood in Fig. 6, where we show the efficiency of dust formation at various epochs for each of the dust type, in the velocity space.

As a general comment, the inner parts of the ejecta (smaller velocities) tend to be more clumpy as well as cooler in temperature, compared to the outer shells, at any give time. Both these factors enhance the dust formation efficiency and ensure an earlier onset of dust formation. As seen in Fig. 6, the dust first appears in the velocity bin of about 1000 km s⁻¹, while by day 2900, it is spread out over the entire velocity space of the He core (He-core boundary ~2500 km s⁻¹). The majority of the CO molecules form in a narrow band between 2200 and 2350 km s⁻¹. Interestingly, due to faster cooling in those bands, induced by CO molecules, the formation of dust in those regions precedes the formation of dust in the neighboring velocity bins; this can be seen in the panel of day 1000 in Fig. 6. However, a large fraction of the O atoms are locked in CO molecules in that part of the ejecta; hence, the efficiency of O-rich dust production is largely compromised. There is a very narrow (almost) dust-free shell between 2300 and 2400 km s⁻¹, where silicates are formed only in a very small proportion.

The overall efficiency is slightly compromised in the bins where dust formation commences later, which, in general, can be associated to larger velocities; this is also visible in Fig. 6. However, it is important to realize that in the velocity space, the smaller velocities are inner shells, so the total mass of those shells are smaller than the outer ones. Therefore, the inner shells start dust formation earlier with larger efficiencies, but those shells are low in total mass, and therefore end up forming a small mass of dust. When dust formation progresses in the outer shells, the total mass of the shells being larger compromises for the reduced efficiency.

In Fig. 7 we present the complete evolution of each dust type, summed over the all the clumps in the ejecta. We find that the formation of dust commences in the ejecta around 580 days post-explosion. There could be a very small mass of sulfides, or metallic clusters, formed a little earlier than that, but their mass is not expected to exceed 10⁻⁶ M_⊙; we did not follow those in this paper. The earliest grains to form in abundance are the magnesium silicates, followed by alumina, silicon carbide, and amorphous carbon at very late times. The mass of dust grows from 4 × 10⁻⁵ M_⊙ at 600 days to 1.8 × 10⁻³ M_⊙ at day 800, and to 6 × 10⁻³ M_⊙ at day 1000. The dust continues to form until ~day 2800, but after day 1600, the rate of dust formation is considerably reduced. The total dust mass in this scenario is 0.06 M_⊙. Figure 8 shows how the total column density (from the center to the outer edge of the ejecta) of the dust evolve with time, which combines the effect of dust formation with the expansion of the ejecta. Now, we shall briefly discuss the scenario of each dust species individually.

Fig. 5 Resulting dust temperatures for silicates (top left), alumina (top right), carbon (bottom left), and silicon carbide (bottom right) are shown for different days and different velocity bins. The solid black lines refer to loci in terms of time and velocity, when the dust formation starts in a particular bin. The enclosed area inside the contours (dashed black lines) represents the space where that dust species is present (see Sect. 3 for details).

Fig. 6 Efficiency of dust formation, which is the same as the dust-to-gas ratio in the SN ejecta, is shown for days 600, 800, 1000, 1200, 1500, and 2900, as a function of velocity. There is no significant change in the ratios between day 1500 and 2500, and therefore we chose to present the case of day 2900 after day 1500.

Fig. 7 Total mass of grains of silicates, alumina, amorphous carbon, and silicon carbide, as a function of post-explosion time. The sum of all dust species is shown by the solid black line. The results correspond to clumpy, isotropic SN ejecta, whose progenitor was a 20 M⊙ star (at main sequence).

Fig. 8 Column density of dust as a function of time, post-explosion. Despite the nonuniform distribution of dust, this can be correlated to the total opacities (from the center to the outer edge of the ejecta), using the relative abundances of each dust type at a given epoch and their corresponding optical constants (see Sect. 4).

3.1 Silicates

Silicates appear to form in the ejecta about 1.5 yr after the explosion, which is the earliest among the four major dust species. Silicates form in the ejecta over a period that starts from about day 580 and continues until about day 1300. The formation zone spans over a large range in the velocity space (approximately between 900 and 2300 km s⁻¹; see Fig. 6). As seen in Fig. 6, the formation starts in a bin close to v ~ 1000 km s⁻¹; the formation zone generally follows a gradual inside to outside trend, other than for a few exceptions. For instance, the formation of CO molecules in large proportion and the cooling of the gas associated with it, often induce an earlier synthesis of silicates in that region of the ejecta. The mass of silicates grows from 4 × 10⁻⁵ M_⊙ at 600 days to 0.025 M_⊙ at ~1400 days; thereafter its mass remains mostly unchanged (see Fig. 7).

3.2 Alumina

Alumina is the second most abundant dust species to form in the SN ejecta. The final masses are comparable to silicates, but they form much later in time. The synthesis starts at about day 860, and it continues until ~day 1500. Their formation regions in the velocity space overlap with those of silicates, since they are both formed where O atoms are the most abundant in the gas (see Fig. 6). They end up at a final mass of about 0.022 M_⊙ (see Fig. 7).

3.3 Amorphous carbon

The synthesis of amorphous carbon dust provides an interesting scenario. It is found to form in the outer parts of the He core, where the velocities are between 2350 and 2550 km s⁻¹ (Fig. 6). Those parts of the ejecta are He rich in nature, as seen in Fig 1. The helium atoms, when ionized by the fast Compton electrons, produce He⁺, which are extremely detrimental to any stable molecules or clusters present in the gas (Cherchneff & Dwek 2009, 2010; Sarangi & Cherchneff 2013). Eventually, with the decrease in the number of the Compton electrons, ions of an inert gas such as helium tend to quickly recombine back to their neutral state. The timescale for this depends on the density of the gas in the clumps. In this scenario, we find the abundances of He⁺ to be no longer significant after ~900 to 1000 days.

The formation of stable diatomic carbon, C₂, is a crucial bottleneck in the pathways to form carbon chains, which eventually lead to the formation of carbon rings and fullerenes (Sarangi & Cherchneff 2015). In the first ~1000 days, the formation of stable C₂ molecules is hindered by the presence of He⁺. After that, stable C₂ molecules form via the radiative association of two C atoms: $$\begin{array}{lll}C+C\hfill & \to \hfill & C_2+hv\hfill \\ O+C_2\hfill & \to \hfill & C+\text{CO}.\hfill \end{array}$$(13)

In the outer parts of the ejecta, the carbon atoms are more abundant in number compared to oxygen; however, there is a relatively large abundance of ¹⁸O compared to ¹⁶O in those regions (Rauscher et al. 2002). The reaction between neutral O and C₂ (Reaction 13) is found to be a favored chemical path that results in the further breakdown of C₂ into atomic C. In this way, the C atoms end up in a loop with C₂ and CO. Depending on the density, this cycle continues for a rather long time, which does not allow for the formation of larger carbon clusters. Then finally, after a couple of years, all the O atoms end up depleted in CO molecules. Thereafter, rather quickly (within a few days), C chains and clusters are synthesized, leading to the condensation of amorphous carbon dust. Moreover, this series of processes are found to be extremely sensitive to gas densities. Higher gas densities results in an earlier recombination of He ions, and thereby an earlier onset of CO formation that leads to a much faster formation timescale. This also means that a small change in the abundance of oxygen, especially of ¹⁸O, can be significant in altering the timescale of the appearance of carbon dust in the ejecta, and its total mass. There is, of course, a certain degree of uncertainty associated with the ratio of ¹⁸O to ¹⁶O predicted by the SN nucleosynthesis models; on the other hand, observationally it could be possible to identify the O isotope in the CO. This can then be used as a tool to predict the mass and epochs of carbon dust formation in the ejecta.

The formation of carbon dust starts around 2400 days after the explosion and continues until day 2800 (Fig. 7). In this time a mass of about 8 × 10⁻³ M_⊙ is formed, which is strictly controlled by the C-to-O ratio. It found to be always less abundant than O-rich dust species in the gas.

3.4 Silicon carbide

Silicon carbide is the least abundant dust species, which reaches a maximum of about 10⁻⁴ M_⊙ (Fig. 7). SiC dust forms in two phases and in two distinct layers. In the inner layer of the ejecta, where the velocity range is between ~850–950 km s⁻¹, silicon carbide forms around day 900, in a small mass of about 4 × 10⁻⁶ M_⊙ (Fig. 6). Thereafter, post 1900 days, more SiC dust is synthesized in the outer regions of the He core, in the same parts of the ejecta where C dust forms. The production of SiC at times later than 2000 days are potentially in agreement with the isotopic signature of SiC grains (of SN origin) found in meteorites (Liu et al. 2018; Ott et al. 2019).

4 Progressive opacities and infrared spectra

In this section, we discuss the IR radiation and the fluxes generated by the distribution of dust in the SN ejecta. The rate of emission per unit wavelength produced by a single dust grain is shown in Eq. (12), where it is integrated over all wavelengths to find the cooling rate. The emerging radiation from a distribution of dust is, however, not always the total radiation of all dust grains combined, since it also depends on the opacities of the medium as a function of wavelength. Assuming a uniform distribution of dust within a spherical volume, the probability of escape of the IR photons, P_esc(λ), was estimated by Dwek & Arendt (2015) with respect to SN 1987A.

In this study, by modeling the formation of dust in the thinly spaced, clumpy shells in the ejecta, we estimated the generic distribution of each dust species individually, along with their respective distribution of temperatures. This is illustrated in Fig. 9, where we zoom into a small cross-section of the O-rich zone of the ejecta, to show thinly spaced shells characterized by different optical depths (owing to their different dust densities and compositions) and different dust temperatures.

For a thin shell of homogeneous density and thickness ∆r, the optical depth τ_h(λ) and the corresponding escape probability $P_{\text{esc}}^h$ (Inoue et al. 2020) of photon at wavelength λ, can be expressed as $$\begin{array}{lll}{\tau }_h\left(\lambda ,v\right)\hfill & =\hfill & {\displaystyle \sum _S{\displaystyle \sum _a{\displaystyle \int {\rho }_{d,S}\left(v,a\right)k_S\left(\lambda ,a\right)\text{d}r}}}\hfill \\ \hfill & =\hfill & {\displaystyle \sum _S{\displaystyle \sum _a\frac{M_{d,S}\left(v,a\right)}{4\pi r^2\text{Δ}r}{k}_S\left(\lambda ,a\right)\text{Δ}r=\frac{{\displaystyle {\sum }_{S,a}{M}_{d}k}}{4\pi r^2}}}\hfill \\ P_{\text{esc}}^h\hfill & =\hfill & e^{-{\tau }_h},\hfill \end{array}$$(14)

where ρ_d,s and M_d,S are the density and total mass of the dust, respectively, of a given species S and radius a, in that shell. We have summed over all grain sizes a and all dust species S. $P_{\text{esc}}^h$ defined in this way, of course, is a function of the velocity bin (v) of the ejecta; in other words, each shell will have a unique $P_{\text{esc}}^h$.

For a clumpy medium, a suitable approach for deriving the opacities and escape probabilities is to consider a clump as one large dust grain or mega-grain, defined by its individual absorption and emission properties (Városi & Dwek 1999; Inoue et al. 2020). In that scheme, following Inoue et al. (2020), the optical depth of a clump can be written in terms of the corresponding homogeneous optical depth: $$\begin{array}{ccc}{\tau }_{\text{cl}}\left(\lambda \right)={\tau }_h\left(\lambda \right)\frac{r_{\text{cl}}}{r}{f}_V^{-1};& r_{\text{cl}}\simeq \frac{r{f}_V}{10};& {\tau }_{\text{cl}}\simeq 0.1{\tau }_h,\end{array}$$(15)

where we have used the limiting relations for a radius of a clump, r_cl, from Dessart et al. (2018). Assuming most of the dust is located inside the clumps, we can derive the escape probability of the individual shells, $P_{\text{esc}}^{\mathrm{sh}}$ (Inoue et al. 2020), in the following form: $$\begin{array}{lll}{\tau }_{\text{eff}}\left(\lambda \right)\hfill & =\hfill & {\tau }_h\left(\lambda \right)P_{\text{esc}}^h\left({\tau }_{\text{cl}}\right)={\tau }_h\left(\lambda \right)e^{-0.1{\tau }_h\left(\lambda \right)}\hfill \\ P_{\text{esc}}^{sh}\left(\lambda ,v\right)\hfill & =\hfill & P_{\text{esc}}^h\left({\tau }_{\text{eff}}\right)=e^{-{\tau }_{\text{eff}}\left(\lambda ,v\right)}.\hfill \end{array}$$(16)

Finally, using Eq. (16), we calculate the probability of escape of an IR photon from the entire ejecta, to reach the observer. If the photon originates from a clump located in a given shell identified by velocity v, then P_esc can be expressed as the cumulative product of the probability of escape of all the shells lying along its path, starting from that shell to the outer edge of the He core. Therefore, we can write $$P_{\text{esc}}\left(\lambda ,v\right)={\displaystyle \prod _v^{v_c}{P}_{\text{esc}}^{sh}\left(\lambda ,v\right)=\exp \left[-{\displaystyle \sum _v^{v_c}{\tau }_{\text{eff}}\left(\lambda ,v\right)}\right].}$$(17)

We have termed P_esc(λ, v) as the progressive escape probability, since it is expressed in the velocity space (along with λ and t). Considering the generic distribution of dust is nonuniform, this formulation is the way to account for the dust density distribution at any give time, when estimating the chances of a photon originating from a given clump in a shell, to reach the observer.

We find that the sizes (radii) of all dust species range between 0.001 and 0.1 µm. In the infrared, the wavelengths are much larger compared to the average grain radii (λ_IR ≫ a). Within this Rayleigh limit, the luminosities of the dust grains are not dependent on the grain sizes. The same can be said for the optical depths in the IR regime of the individual clumps. For simplicity, a single grain of radius 0.01 µm was assumed.

In Fig. 10 we show P_esc for days 800, 1000, 1500, and 2000 for the Spitzer Space Telescope Infrared Array Camera (IRAC) bands (λ = 3.6, 4.5, 5.8, 8.0 µm). As the figures show, at the earlier epochs, the opacities are large (P_esc small) only in the inner parts of the ejecta. With time, the opacities are found to increase more and more in the outer velocity bands, as the formation zones of silicates and alumina shift outward. Overall, the escape probabilities are shown to decrease consistently until about day 1500, after which there is a decline in the opacity, and therefore an increase in P_esc. This results from a combined effect of the slowing down of the dust production and the expansion of the ejecta. Figure 8 shows the dust column densities, which also present similar trends and can easily be translated to the total opacities from the center of the ejecta to the outer edge of the He core. The 8-µm band clearly shows larger opacities than the other three wavelengths, owing to the larger values of the absorption coefficients of O-rich dust grains close to that wavelength (see Fig. 3).

The IR flux from SN ejecta is the IR luminosity (amount of energy that escapes the ejecta per unit time), divided by 4πL², with L being the distance of the SN from the Earth. The observable IR luminosity is expressed as $$\begin{array}{lll}{L}_{\lambda }^{\text{obs}}\left(t\right)\hfill & =\hfill & {\displaystyle \sum _v{\displaystyle \sum _{S,a}4M_{d,S}\left(v,t,a\right)k_S\left(\lambda ,a\right)}}\hfill \\ \hfill & \hfill & \times \pi B_{\lambda }\left(\lambda ,T_S\left(v,t,a\right)\right)P_{\text{esc}}\left(\lambda ,v,t\right),\hfill \end{array}$$(18)

where we have summed over all dust species S, and all grain sizes and temperatures, over the complete velocity space. The total energy emitted by dust per unit time, which can also be called the unattenuated luminosity, is given by the same expression, without the P_esc term.

In Fig. 11, the unattenuated as well as the observable luminosities produced by the dust grains in the SN ejecta are shown for days 700, 800, 900, 1000, 1200, 1500, 2000, 2500, and 2900. As expected, the ratio of the two luminosities changes over time, and the contrast is largest between 1200 and 1500 days. At longer wavelengths, typically larger than ~100 µm, the unattenuated and the observable luminosities almost always concur, since the dust opacities are much lower in the far-IR wavelengths compared to the mid-IR bands.

Observationally, the detection or non-detection of the strong silicate feature at 9.7 µm (Simpson 1991) often motivates the estimation of the possible dust composition (Ercolano et al. 2007; Wesson et al. 2015; Wesson & Bevan 2021). Our results suggest that the composition of dust in the ejecta is heavily dominated by O-rich dust at all times. However, we want to stress that, when there is a local maximum in the absorption coefficient at any wavelength band (for example around 9.7 µm for silicates), while the IR luminosity in that band increases, the optical depths at that wavelength also increase simultaneously. Therefore, this leads to a strong attenuation at those wavelengths. The combination of the two (increased luminosity, increased attenuation), often leads to a complete suppression of the strong features in the spectra. This was also pointed out by Dwek & Arendt (2015). Figure 11 confirms this as well, where we can clearly see the strong features originating from the O-rich dust components in the unattenuated luminosities until day 1000, while the observable luminosity never shows any noticeable feature. After day 1000, dust formation is accelerated in the outer parts of the O/Si zone, and the mass-weighted dust temperatures are already too cool to radiate strongly in the 10–20 µm range.

As we see in Fig. 11, post day 2500, the ejecta becomes increasingly thin to the self-emitted IR radiation amorphous carbon dust appears in the outermost layers of the He core at about the same time. However, they do not form an outer dusty shell that is thick enough to block the radiation coming from the inner O-rich dust grains.

Fig. 9 Schematics of the geometry is shown, where the ejecta is characterized by different zones according to the abundances of elements, while each zone is stratified into several layers, defined by unique gas densities in the clumps, as well as unique dust and gas temperatures. As a result, the efficiency of dust production varies widely among the layers, and so do the opacities. The zoomed-in view of the small random rectangular cross-section of the gas indicates the unique dust densities, dust temperatures, and opacities in its different layers (despite being entirely in the O core).

Fig. 10 Wavelength-dependent progressive escape probabilities (Eq. (17)), which are a function of the local velocity, are presented for days 800, 1000, 1500, and 2000, for four different wavelengths (3.6, 4.5, 5.8, and 8.0 µm) that correspond to the Spitzer Space Telescope IRAC bands.

Fig. 11 Radiated energy per unit time (in solar luminosities) emitted by dust in the SN ejecta as a function of wavelengths, for nine different epochs (see Eq. (18)). The unattenuated luminosity (dashed line) corresponds to the total energy per unit time (from the dust), while the observable luminosity (solid line) is the luminosity that escapes the system, after taking the progressive escape probabilities into account, and will be the one that can be observed (see Sect. 4 and Eqs. (17) and (18) for details).

Fig. 12 Geometry of anisotropic ejecta and the corresponding zones, in the form of an ellipsoid. To emulate SN 1987A, we chose the velocity ratios to be 2:3:1 (with reference to Kjær et al. 2010), where the velocity distribution along the mean axes (taken as the x axis in the figure) remain the same as the isotropic model.

5 Elliptical geometry

In this section, we discuss the formation and distribution of dust in the ejecta, and the resulting IR emission, when the geometry of the ejecta is assumed to be anisotropic.

By analyzing the spectral information obtained from the 3D spectroscopy of SN 1987A, Kjær et al. (2010) found that the 2D projection of the ejecta is elliptical in shape with an axial ratio of approximately 1:1.8. The observations can be fit well using an ellipsoidal geometry where the axes are in ratios of 1:2:3.

We adopted the same geometry for the ejecta, where the ellipsoid is defined by the ratio of velocities along the x, y, and z axes, which are taken to be 2:3:1, as shown in Fig. 12. The velocities along the x direction are kept identical to the velocities of the spherical ejecta that were derived in Sect. 2. With respect to that, the velocities along the y axis are then given by v_y = 1.5v_x, and along z axis by υ_z = 0.5v_x. Assuming a conservation of column densities, the gas density therefore scales along the y and z directions as n_y = (2/3)n_x and n_z = 2n_x, respectively, while n_x (along the x direction) remains the same as the gas densities in the spherical case.

In Sect. 2.4, we explained the stratification of the ejecta in the velocity space, where the He core is divided into 450 thinly spaced spherical, clumpy shells, identified by their velocities. Adopting the same approach, each shell in this case is converted into a thin ellipsoidal volume, where the velocities scale in the chosen ratio of 2:3:1, along the three axes. To be specific, the clumps at velocity v of the spherical model are now distributed in a range of velocities from 0.5v to 1.5v; all these clumps have identical abundances, but the gas densities, and the gas and dust temperatures are now distributed over a range as a function of their radial velocities. The spread of gas temperatures in each zone of the ellipsoidal ejecta are shown in Fig. 13, in comparison to Fig. 2 of the isotropic model. The O/C zone and the He/C zone, which are toward the outer part of the ejecta, clearly shows a larger range in gas temperatures compared to Fig. 2. In the isotropic model, the O/C layer was a thin shell, so the range of temperatures was also narrow. However, that thin shell in an ellipsoidal geometry now covers a range of velocities along different axes, which explains the larger spread in the temperature distribution. The same is true also for the He/C zone as well.

Within each ellipsoidal shell, the individual clumps can be identified by their polar coordinates. We associate the clumps with their respective velocities within each shell, as shown for a hypothetical clump at point P in the figure. A bin size of 80 was chosen for each shell, which means the ejecta were stratified into 450 × 80 different groups of clumps (450 ellipsoidal shells, each divided into 80 bins). Fortunately, for each of these shells, all the clumps within a single shell have the same mass fraction of elements, and only vary in physical conditions. Hence, the computation is relatively less strenuous, compared to an actual 3D hydrodynamic model, where any random parcel of gas has unique abundances. In the following sections, we discuss the chemical yields, estimated assuming this anisotropic geometry.

Fig. 13 Distribution and evolution of temperature with respect to the various zones are shown for the ellipsoidal ejecta. Comparing to the same distribution and evolution (Fig. 2) for the spherical model, it is evident that the temperatures are spread to wider ranges in this case, as expected.

5.1 Comparing dust masses

Here we describe the chemistry and dust masses formed in the ejecta. The same formalism, as described in Sect. 3 is used, where we trace the distribution and evolution of dust grains of silicates, amorphous carbon, alumina, and silicon carbide.

Along the smallest axis (which is the z axis in our model) of the ellipsoidal ejecta, the clumps are denser and the gas is found to cool faster. Both these factors considerably enhance the rate and efficiency of dust production along that axis, which has the smallest velocities. Conversely, along the large axis, the gas is found to be hotter (compared to the other two directions), as well as less dense, leading to a delayed dust formation.

The top panel of Fig. 14 shows the evolution of dust masses in this scenario and its comparison with the isotropic model. The onset of dust formation is about 440 days, leading to the formation of a small amount of dust, typically less than 10⁻⁵ M_⊙. By 600 days the total mass of dust increases to 9 × 10⁻⁵ M_⊙, and then to ~ 10⁻³ M_⊙ by the end of 800 days post-explosion. After 1000 days, the total mass of dust in the ejecta is 5.5 × 10⁻³ M_⊙, and finally it ends up forming about 0.059 M_⊙ of dust.

The dust composition is dominated by silicates in the first ~1400 days, after which alumina and silicates are found in almost equal proportion. The final mass of both these species are between 0.02 and 0.03 M_⊙. Carbon dust appears as soon as 1200 days after the explosion, and continues to form until 2800 days. For silicon carbide, the formation starts after about 900 days and continues slowly over a long time, finally ending up at a mass close to 2 × 10⁻⁴ M_⊙.

In comparison to the isotropic case, the evolution of dust masses differs in the following way: (a) the formation of dust starts at about 450 days post-explosion in the ellipsoidal ejecta, compared to 580 days in the spherical scenario; (b) the trend of an earlier synthesis of dust is valid for both the O-rich dust species in the ellipsoidal case; (c) the timescale for the onset of carbon dust production changes significantly from 2400 days in the isotropic model to 1200 days in the anisotropic scenario; (d) the final mass of carbon dust produced in the anisotropic ejecta is 5 × 10⁻³ M_⊙, compared to 8 × 10⁻³ M_⊙ in case of the isotropic model, and so even though carbon dust forms earlier in the ellipsoidal scenario, the final mass is smaller; (e) the mass of silicon carbide in the ellipsoidal ejecta is larger than the mass in the spherical case by a factor of 1.5.

The bottom panel of Fig. 14 compares the evolution of each dust species for the two cases. We find that the final dust masses in the two cases are almost identical. The epochs for the onset of formation and the initial rate of production vary to some extent for all the dust species. This is mainly induced by the geometry, where within the same ellipsoidal shell, some of the clumps see much earlier dust production than others. When this is aggregated over all the shells, the earliest and the latest epochs for dust formation in the ejecta get spread out over a larger span of time compared to an isotropic ejecta model. In Sect. 3.3, we stressed the fact that the formation of amorphous carbon is highly sensitive to the physical conditions of the ejecta. In the same context, here we see that in the denser clumps, which are along the slowest axis (z axis), carbon dust forms as soon as 1200 days post-explosion, while in the clumps along the fastest axis (y axis), the formation of C dust does not take place even at the end of day 3000, when we terminate our simulations. This results in an overall smaller mass of C dust at 3000 days in the ellipsoidal scenario.

Even though the masses of dust are not significantly different in the two cases, the distribution of dust inside the ejecta shows strong anisotropy. In Fig. 15 the evolution of dust column densities (from the center to the edge of the He core) are presented for days 700, 1000, 1500, 2000, and 2900 along all the directions in each of the three axial planes. The figures clearly indicate that along the xz or the yz plane, the column densities increase rapidly with a decrease in θ (θ = 0 is the z direction). The velocity vz being the slowest velocity, the dust is more concentrated in that direction, as previously explained. The contrast is more visible in the earlier epochs, because the onset of dust formation also varies widely along the planes. In the yz plane, the variation in column density with θ is most prominent, since the velocities along the y and z axes are the largest and smallest, respectively. The xy plane has relatively lower columns of dust and also later epochs of formation, compared to the other two.

Fig. 14 Mass of dust for the ellipsoidal case is presented. Top: evolution of total dust masses, in the case of anisotropic ejecta, are shown for each dust species, and the sum is shown by the solid-line. Bottom: mass of dust produced in the spherical and ellipsoidal cases are compared for each dust species.

Fig. 15 Dust column densities are shown for the ellipsoidal model for five different epochs. In the top left panel, the chosen polar coordinates are shown in terms of θ and φ for an ellipsoidal shell. The velocity of any clump at point P inside the shell is calculated. In the top right panel and the bottom panel, we show the evolution of dust column densities (from the center to the outer edge of the ejecta) along various directions as a function of the viewing angle, in the three axial planes. See Sect. 5.1 for details.

5.2 Anisotropic IR spectra

Based on the resulting distribution of dust, we calculated the opacities and progressive escape probabilities of the medium along various directions. Using Eq. (17), the P_esc(λ, v_x, θ, ϕ) can be estimated, which tells how many of the IR photons of wavelength λ, originating in a thin ellipsoidal shell (defined by velocity along the x axis, v_x), escape the ejecta to reach the observer who is in the (θ, ϕ) direction. When calculating the optical depth τ_h(λ, v_x, θ, ϕ) using Eq. (14), the total mass of dust cannot be used in the expression, since the distribution of dust is anisotropic. Instead, the dust densities in a given direction should be integrated over all the shells along its path, as it is also shown in Eq. (14).

In Fig. 16, the progressive escape probability along the three axes of the ellipsoidal ejecta are shown for days 1000 and 1500, for the Spitzer Space Telescope IRAC bands λ = 3.6, 8.0 µm. To summarize the trend, initially dust forms in the clumps centered around the z axis. However, the emission from those dust grains cannot easily escape in the z direction since the opacities become large very quickly. On the contrary, the same radiation can easily pass through the He core in the x or y direction, since the clumps along the path in those directions are not very dusty in the earlier epochs. With time and dust formation, the opacities along other directions also start to increase, although this trend remains true for later epochs as well, since the dust is always more concentrated toward the smaller-velocity axis. Moreover, along the faster axes, even if dust forms with equal efficiency, as the densities are lower and the clumps are more spread out, the opacities cannot get much larger. All these factors collectively result in the variation in P_esc with velocities along different directions, which is seen in Fig. 16. For the spherical ejecta, the outer-edge of the He core is at 2650 km s⁻¹, which is the same for the x direction in this scenario, while the velocities are half of that and 1.5 times that along the z and y directions, respectively. In Fig. 16 we also show the escape probabilities as a function of normalized He-core velocities, which actually reflect the efficiency of dust production along the three axes of the ellipsoidal ejecta, normalized to the same scale.

The observable luminosity, given by Eq. (18), is again a function of the direction (θ, ϕ), in correlation with the directional escape probabilities. However, the unattenuated luminosity has no direction dependence since it is just the total amount of energy radiated by all the dust grains together, per unit time. In Fig. 17, we show the observable luminosities along the three axes, compared to the unattenuated luminosities, for nine epochs. In qualitative terms, the variation in luminosity with wavelength, which is explained in Sect. 4, is also applicable to the anisotropic ejecta.

Looking at Fig. 17, it is evident that the observable luminosity in the z direction is always much lower than in the other two axes, given the large opacities in that direction. Moreover, there are no particular spectral features visible along that axis. For days 700 and 800, the observable luminosity along the y direction coincides with the unattenuated luminosity, which implies that the ejecta is optically thin along the y axis, since the clumps aligned along that axis are yet to form dust. The 9.7- and 18-µm spectral features of silicate dust (Simpson 1991) are therefore visible in the spectra along the xy plane until about day 900, when the distribution of dust is very skewed toward the smallest velocity axis only.

With time, dust forms in all the clumps, with varied efficiency. The expanding ejecta also become optically thinner. Typically, after 2400 days, we find the observable luminosity in all the directions to be within comparable range. After 2800 days, the ratio between the unattenuated and the observable luminosities in all the directions approaches unity.

Fig. 16 Progressive escape probabilities, corresponding to the anisotropic case, are shown for two epochs, days 1000 and 1500, for the 3.6- and 8.0-µm Spitzer IRAC bands. The other two IRAC bands, 4.5 and 5.8 µm, also show a similar behavior. The escape probabilities are presented in the top panel as a function of ejecta velocities, while in the bottom panel, we show the probabilities as functions of the normalized He-core velocities (normalized with respect to the velocity of the outer edge of the He core in that direction). The escape probability in this scenario is dependent on the viewing angle; here we present the three axial directions. See Sect. 5.2 for more details.

Fig. 17 Unattenuated luminosities, as well the observable luminosities, along the three axial planes are shown for nine different epochs. For reference, in our chosen geometry, the z axis is the smallest axis and the y axis is the largest axis of the ellipsoid. See Sect. 5.2 for the description.

5.3 SN 1987A at days 615 and 775

The spectra of SN 1987A at post-explosion epochs 615 and 775 are presented by Wooden et al. (1993). Using these data, the masses and the composition of dust formed in the ejecta are constrained (Ercolano et al. 2007; Wesson et al. 2015; Dwek & Arendt 2015; Wesson & Bevan 2021).

In this paper, we study the formation and distribution of dust, and their temperatures, in a bottom-up approach. In Fig. 18, we plot the IR fluxes from our model alongside the observed fluxes of SN 1987A taken from Wesson et al. (2015). For days 615 and 775, the predicted fluxes for the isotropic model and the three axial directions of the anisotropic model are presented. It is important to note that these figures are not parameter fits, but instead just the co-plotting of the data and the predictions of our model. Using our results, we only show the contribution of dust to the fluxes in the IR in Fig. 18. The UV-optical fluxes needs to be added to it in order to obtain the complete SED. For instance, at day 615 the observed fluxes at wavelengths smaller than 2 µm, shown in the figure, can only be accounted when we consider the contribution of the optical luminosities.

The total amount of dust formed in our calculations at days 615 and 775 for the isotropic model are 1.7 × 10⁻⁴ M_⊙ and 1.5 × 10⁻³ M_⊙, respectively; for the anisotropic case, these are are 2.1 × 10⁻⁴ M_⊙ and 9.7 × 10⁻⁴ M_⊙. A distance of 51.4 kpc was assumed for SN 1987A from the Earth.

When the ejecta are assumed to be anisotropic, we find that, between days 600 and 800, newly formed dust is mainly confined along the clumps that are aligned close to the z direction. Therefore, the opacities are much enhanced along the z direction at these times, which results in a reduced flux escaping along that direction. On the contrary, when the viewing angle is along the x or y direction, the dusty clumps are aligned in a perpendicular cross-section to the observer, and thus the fluxes from all those clumps reach the observer almost unattenuated. This is visible in the top panel of Fig. 18, where along the x or y direction, the predicted fluxes are larger compared to the observations, while the fluxes along the z direction, and also the flux produced by the isotropic ejecta, are smaller than the observed fluxes.

An ellipsoidal ejecta with a tilt of 25° to the plane of the sky is best fit to the observations of SN 1987A (Kjær et al. 2010). Using the polar coordinates given in Fig. 15, we find that for day 615, the fluxes obtained in our calculations are closest to the observations at a viewing angle of θ = 16°, ϕ = 65°. For day 775, these are at θ = 21°, ϕ = 69°. Both the viewing angles being close to each other, in a way, justifies their credibility.

Fig. 18 Fluxes obtained from our calculations (assuming the SN to be at the distance of SN 1987A) for days 615 and 775, are co-plotted with the near- and mid-IR data of SN 1987A (shaded region in yellow) from Wesson et al. (2015). In the top panel, the fluxes for the isotropic model, as well as the three axial directions of the anisotropic model, are shown for the two epochs. In the bottom panel, we show the viewing angle with respect to the coordinates shown in Fig. 15, where we find the fluxes emerging from the anisotropic ejecta match well with the observations. For the details, see Sect. 5.3.

6 Summary

In this study, we have revisited the scenario of dust production in SN ejecta, based on an improved understanding of the ejecta temperature, and their cooling, along with the nature of clumpiness within the ejecta cores. The evolution and distribution of dust in the He core was derived, along with the respective opacities. The emerging IR emission from the dust was modeled, based on the generic, nonuniform distribution of the dust grains that we estimated together with their respective temperature distributions and progressive opacities. For the first time, we show the impact of anisotropy on the chemistry of dust formation by choosing an ellipsoidal geometry for the SN ejecta that imitates SN 1987A.

In a few points below, we summarize the results of this study:

(1) The composition of dust in SN ejecta is dominated by O-rich components (silicates and alumina) at all times. For this study, we considered the progenitor star to be 20 M_⊙ at main sequence. Sarangi & Cherchneff (2013) found that the relative mass of O-rich dust increases with an increase in progenitor mass.

(2) The formation of dust commences at around 580 days post-explosion, if the ejecta are assumed to be isotropic, while it starts at day 450 when an anisotropic, ellipsoidal geometry is adopted for the ejecta.

(3) The rate of dust production is relatively rapid in the period from year 2 to 5, after the explosion of the SN, typically averaging about 6 × 10⁻⁵ M_⊙ day⁻¹; following that, the rates decline significantly.

(4) The formation of carbon dust is especially sensitive to the physical conditions of the ejecta. In our model, we find that in the isotropic model, carbon dust forms almost 6.5 yr after the explosion. The clumps that form carbon dust are also rich in C¹⁸ O molecules, whose formation precedes the synthesis of the dust grains.

(5) In the spherical ejecta, the mass of dust evolves from less than 10⁻⁵ M_⊙ at 600 days to about 10⁻³ M_⊙ at 800 days, to ~0.06 at 1600 days. In the ellipsoidal scenario, the formation of dust starts earlier. However, post day 800, the total mass of dust present in the ejecta at any given time are comparable to the spherical model, with a variation of no more than 5%.

(6) In this study, we have derived the distribution of dust in the velocity space. For the isotropic scenario, we find that the dust is initially confined in the inner regions only, where the velocities are about 1000 km s⁻¹. With time, the distribution spreads toward the outer layers. One exception to this trend is the region where CO molecules form in abundance and induce extra cooling, resulting in early synthesis of dust. In the outermost layer, where the velocities are larger than 2300 km s^-1, dust forms as late as 2400 days post-explosion.

(7) Even though the total mass of dust in the ellipsoidal ejecta is comparable to the spherical model, the distribution of dust is highly anisotropic in nature. We find that the clumps that are aligned close to the smallest axis of the ellipsoid are more suitable for dust production, both in terms of timescales and efficiency. This results in a larger concentration of dust in the slower clumps at a given time, if the clumps with identical abundances are compared.

(8) To account for the nonuniform distribution of dust and its temperature, we derived the opacities as a function of velocity, thereby defining the progressive escape probability of a photon originating from a clump at velocity v to the outer-edge of the ejecta. The opacities were found to increase due to dust formation until about day 1600, followed by a gradual decline.

(9) Even though the dust is always O-rich, we show that for the isotropic ejecta, the spectral features of silicates or alumina, which are between 10–18 µm, will not be present in the observed SED after 700 days. This is simply because when there is a local maxima in the emission, the opacities are also large in that wavelength. Moreover, the emission from the inner part of the ejecta is increasingly blocked, as the dust formation progresses outward with time. The effective temperatures of the observable fluxes, estimated from our calculations, are already too low after day 1000 to have strong emission around the 10-µm band.

(10) In the anisotropic ellipsoidal geometry, the polar distribution of dust results in relatively much lower opacities along the faster velocity-axes; in other words, dust forms in clumps, which are in a plane perpendicular to the line of sight, when observed along the large axis of the ellipsoid. Therefore, the predicted observable luminosities are found to be more than an order of magnitude larger than the luminosity along the small axis.

(11) By comparing the fluxes from SN 1987A at days 615 and 775 with our predicted SEDs, we estimate the most likely alignment of the SN ejecta with our line of sight.

7 Discussion

In Sarangi & Cherchneff (2013, 2015), the gas densities derived from Nozawa et al. (2010) were used to construct a homogeneous density model for the ejecta, and a clumpy model was based on a density enhancement to that, using the filling factors given by Jerkstrand et al. (2011). In this paper, we adopted a revised density structure for the ejecta, which is consistent with the dynamics of SN ejecta (Truelove & McKee 1999), the nature of clumpiness used by the radiative transfer models of the nebular phase (Dessart et al. 2018), and which are in good agreement with the multidimensional explosion models of SNe (Wongwathanarat et al. 2015; Utrobin et al. 2017; Gabler et al. 2021). We find that the densities derived from Nozawa et al. (2010) are actually more consistent with clumpy ejecta. Hence, a density enhancement over that leads to an artificial increase in the densities, thereby suggesting a very rapid formation of dust in the first couple of years in the clumpy model of Sarangi & Cherchneff (2015), which saturates at a total budget of 0.13 M₀. In this paper, our results show that only a small mass of dust is formed in the first couple of years in the clumpy ejecta. Dust formation continues until day 2800, though the first 1600 days mark the most active phase. In addition, compared to Sarangi & Cherchneff (2013, 2015), in this paper we have assumed a more realistic evolution of gas and dust temperatures, which are characterized by limited cooling in the first couple of years. The cooling only escalates between year two and year three from the time of explosion, depending on the zone. This results in a delay in the synthesis of dust, thereby lowering densities at the epochs of dust formation and reducing efficiencies compared to Sarangi & Cherchneff (2015). However, to be clear, the revised estimations of the gas densities in the clumps are primarily responsible for the lower dust yields in this study (0.06 M₀), compared to the larger, rapidly growing dust masses predicted by the clumpy model of Sarangi & Cherchneff (2015).

Sluder et al. (2018) predict the synthesis of both O-rich and C-rich dust grains in abundance within the first year of explosion; the dust masses are as large as ~0.4 M_⊙ within a 1000 days. In this paper, we predict much lower dust masses in the first 1000 days, which does not exceed 0.01 M_⊙ and is almost entirely made up of O-rich dust. We find that dust formation starts in the SN ejecta between days 450 and 550, which is in reasonable agreement with the rapid decrease in the line intensities of Mg or Si about the same time (Danziger et al. 1991), implying the depletion of atoms in the gas phase due to dust formation. The model of Sluder et al. (2018) suggests that the rapid increase in dust masses in the first two years from the time of explosion is strongly aided by the accretion of metals (Hirashita & Kuo 2011) on the surface of the small, newly formed grains. In line with Sarangi & Cherchneff (2013, 2015), we continue to argue that the surface reactions are likely to be important only at much longer timescales, and when the gas is much cooler. Our knowledge of accretion on the grain surfaces is largely based on ISM-like environments (Hirashita & Kuo 2011; Hirashita 2012) where chemical equilibrium can be assumed, and thus its application to the much hotter and rapidly evolving SN ejecta is questionable. The rate of sticking on the surface of a grain is expected to considerably decline at high temperatures (Blum 2004). Moreover, following the treatment of ISM dust by Guhathakurta & Draine (1989), it can be suggested that the smallest grains are much more prone to lose atoms from its surface at higher temperatures, which are analogous to the time of their genesis from the gas phase. Most importantly, the detection of molecules (Lepp et al. 1990; Kotak et al. 2009; Matsuura 2017) and atoms (Laming & Temim 2020) present in the gas also speaks against rapid growth of dust via accretion of all the residual gas-phase species on the surface of the grains.

We find that the dust masses in the SN ejecta do not exceed 0.06 M_⊙. This is almost an order of magnitude lower than the masses reported to be present in the SNRs when probed in the far-IR and submillimeter wavelengths (Gomez et al. 2012; Indebetouw et al. 2014; Matsuura et al. 2015). We do not rule out the possibility of accretion when the gas is cool enough after about five years from explosion. However, the grain growth via cold accretion is expected to have loosely bound surfaces that are easily susceptible to destruction by the reverse shock (Dwek & Arendt 2015). The dust destruction efficiency by the reverse shock in Cas A is estimated to be between 70-90% (Slavin et al. 2020; Priestley et al. 2021), which could be correlated to the amount of dust formed through nucleation compared to the dust that grows through accretion at late times. Another important factor could be a significant cooling of the gas due to dust formation, which might have a snowball effect on the dust production. The likelihood of these scenarios could only be determined through further analysis.

In the last few years, most papers dealing with dust in SNe have taken it as a general statement that dust formation models always predict a lot of dust to form very early after explosion, which contradicts the observations. We would like to clarify here that the actual contradiction between theoretical and observational predictions is that theoretical models do not support a scenario of continuous dust production over decades. In the first 1000 days, the most relevant phase for mid-IR observations of any SN ejecta, the mass of dust produced is not necessarily very large, as shown in this paper. Given the compact nature of the ejecta at these early times, even a small mass of dust is sufficient to make the opacities large, as shown in this paper, thereby suppressing the intensity of the emerging flux as explained by Dwek et al. (2019). The generic distribution of dust, nonuniform in nature, is also expected to impact the asymmetries in the emission lines (Bevan & Barlow 2016) differently, compared to a uniform scenario. In the future, we plan to test that for the isotropic and anisotropic cases.

In line with Dwek & Arendt (2015), we show in this paper that even if the ejecta dust is largely O-rich (silicates, alumina), the relatively large opacities efficiently suppress the strong silicate features (between 10 and 20 µm) of the mid-IR wavelengths. In addition, after about 1200 days, the ejecta dust cools down further so that the emission in that wavelength range is not significant. Therefore, the SED itself cannot be conclusively used to assign the composition of dust in the ejecta. In this study, we stress the importance of the generic, nonuniform distribution of dust, and also the dependence on the viewing angle when estimating the emerging fluxes.

The carbon-rich dust components, silicon carbide and amorphous carbon, are found to form must later than the O-rich ones. amorphous carbon dust appears in the ejecta after about 6–7 yr from the time of explosion in the isotropic model. In the anisotropic case, it starts forming after 4 yr from explosion, but only in small amounts. The synthesis of SiC grains takes place in two distinctly separate regions of the ejecta: the outermost and the innermost layers of the He core. In the outermost layer, the formation timescale is comparable to that of amorphous carbon, and in the innermost part of the core, a small mass of SiC forms after about day 1000. This result is especially relevant to account for the estimated delay in the formation of the SiC grains of SN origin found in the presolar grains in meteorites (Liu et al. 2018; Ott et al. 2019), and also to explain the isotopic ratios of elements found in the presolar grains through mixing between the layers of the SN ejecta (Xu et al. 2015).

The chemical kinetic models of dust formation in SNe argue that the production of dust grains starting from gas phase reactions is unlikely to take place at very late times, when the gas temperatures in the SN ejecta are low, since several gas phase reactions require sufficient energy to overcome activation barriers (Sarangi et al. 2018b). In the network leading to the formation of amorphous carbon dust, the formation of C₂ holds the bottleneck; we find that the formation of C₂ molecules takes place after about 1000 days post-explosion, once the He⁺ recombines and importantly when the gas is still warm. However, as we explain in Sect. 3.3, due to a very slow rate of CO production, the condensation of the C₂ molecules into C chains does not transpire until later than day 2400. However, since C₂ was already formed when the gas was still warm, the C chains could condense even at low temperatures to form clusters and eventually carbon dust at such late times. As an additional note, we also want to point out that the conversion of C chains to amorphous carbon dust is still quite swift, once it starts; therefore, this late formation of C dust itself cannot be taken as evidence of the gradual growth of dust over a long span of time. Also, given that the more abundant O-rich dust grains are formed much earlier in time, the late formation of C dust does not provide any significant boost to the overall dust production rate at late times. However, as we mentioned earlier in the text, in the future, we wish to better understand the possibility of accretion of molecules and residual atoms on the surface of these dust grains in cold ejecta, after about a decade from the explosion, to verify if that can amplify the dust formation rates at late times.

In this study, we find the grain sizes to vary between 0.001 and 0.1 µm. When estimating the fluxes, we did not use a distribution of grain sizes, but instead used a fixed value of 0.01 µm, since these small grains are within the Rayleigh limit (typically when a is smaller than 0.1 µm), and we expect only minor differences for the change in grain sizes. In the future, we shall reconsider the possibility of cold accretion and derive the modified grain size distribution taking that into account.

In this paper, we did not test the dependence of net dust production on metallicities, while we acknowledge its importance. The post-explosion elemental and isotopic yields are correlated to the explosion energy and amount of ⁵⁶Ni produced at the explosion. Our current understanding of this correlation (⁵⁶Ni mass to explosion energy to element abundances) is largely incomplete. In the future, when we study the effect of metallicity, we want to simultaneously address this issue.

By analyzing a large variety of isotropic and anisotropic clumps in the velocity space, this study is set as a foundation stone to address the chemistry of 3D ejecta. In this approach, we did not account for any mixing between layers. This model in the future can be coupled to the shuffled-shell technique to imitate macroscopic mixing (Dessart et al. 2021). From our understanding of ejecta chemistry, we can predict that even though the mixing may alter the relative timescale for the formation of certain dust species, we do not anticipate a large change in the final dust masses or the timespan for dust production in the ejecta.

Acknowledgements

We are deeply indebted to Dr. Isabelle Cherchneff, for the development of the SN ejecta chemistry together with Arka Sarangi for previous studies, which is now applied in this paper. We also acknowledge the very useful inputs from Dr. Anders Jerkstrand and Dr. Eli Dwek during the course of this research. Finally, we like to thank the anonymous referee and the editor of this article for his/her precious time. This work was supported by a grant from VILLUM FONDEN, PI Jens Hjorth (project number 16599).

References

All Tables

All Figures

Fig. 1 Mass and density distribution of SN ejecta are presented in the velocity space. Top left: enclosed mass of the ejecta as a function of the ejecta velocity is shown. The He, O, and Si-cores are marked on the figures, with the respective mass-velocity boundaries (see Sect. 2). Top right: fractional abundances of the elements (H, He, C, O, Ne, Mg, and Si), taken from Rauscher et al. (2002), are shown as a function of ejecta velocity, where the mass-velocity relation can be derived from the figure on the left, as shown in Sect. 2. Bottom left: distribution of gas (number) densities in the ejecta, as the function of velocity at 100 days post-explosion (see Sect. 2.1). The densities for uniform ejecta are shown for reference, but we only model the scenario of clumpy ejecta. Bottom right: volume filling factor for the clumps are shown with reference to Eq. (5).

In the text

	Fig. 2 Range of gas temperatures in different zones and their variation with post-explosion time is shown (see Sect. 2.2 for details). Even though our model does not depend on these coarse zone divisions of the ejecta, this is particularly helpful to understand how the gas temperatures evolve.
In the text

Fig. 3 Mass absorption coefficients for silicates (Weingartner & Draine 2001), alumina (Koike et al. 1995), amorphous carbon (Zubko et al. 2004), and silicon carbide (Pegourie 1988) are presented, for a grain size of 0.01 µm. For alumina at wavelengths smaller than 0.3 µm and for silicon carbide at wavelengths smaller than 0.2 µm, the values are extrapolated due to lack of data. In addition, for these two dust types, we do not have data for various grain sizes, and therefore these are generalized values for all sizes.

In the text

	Fig. 4 Radiative heating rate (Eq. (6)) for the four different dust species are shown compared to the (maximum) heating rate of the dust grains via collision with the ambient gas (Eq. (12)), for day 800 and a grain size of 0.01 µm. As it shows, the dust temperatures are controlled by the radiative heating only; that is also the case for other epochs and grain sizes as well.
In the text

Fig. 5 Resulting dust temperatures for silicates (top left), alumina (top right), carbon (bottom left), and silicon carbide (bottom right) are shown for different days and different velocity bins. The solid black lines refer to loci in terms of time and velocity, when the dust formation starts in a particular bin. The enclosed area inside the contours (dashed black lines) represents the space where that dust species is present (see Sect. 3 for details).

In the text

	Fig. 6 Efficiency of dust formation, which is the same as the dust-to-gas ratio in the SN ejecta, is shown for days 600, 800, 1000, 1200, 1500, and 2900, as a function of velocity. There is no significant change in the ratios between day 1500 and 2500, and therefore we chose to present the case of day 2900 after day 1500.
In the text

	Fig. 7 Total mass of grains of silicates, alumina, amorphous carbon, and silicon carbide, as a function of post-explosion time. The sum of all dust species is shown by the solid black line. The results correspond to clumpy, isotropic SN ejecta, whose progenitor was a 20 M⊙ star (at main sequence).
In the text

	Fig. 8 Column density of dust as a function of time, post-explosion. Despite the nonuniform distribution of dust, this can be correlated to the total opacities (from the center to the outer edge of the ejecta), using the relative abundances of each dust type at a given epoch and their corresponding optical constants (see Sect. 4).
In the text

Fig. 9 Schematics of the geometry is shown, where the ejecta is characterized by different zones according to the abundances of elements, while each zone is stratified into several layers, defined by unique gas densities in the clumps, as well as unique dust and gas temperatures. As a result, the efficiency of dust production varies widely among the layers, and so do the opacities. The zoomed-in view of the small random rectangular cross-section of the gas indicates the unique dust densities, dust temperatures, and opacities in its different layers (despite being entirely in the O core).

In the text

	Fig. 10 Wavelength-dependent progressive escape probabilities (Eq. (17)), which are a function of the local velocity, are presented for days 800, 1000, 1500, and 2000, for four different wavelengths (3.6, 4.5, 5.8, and 8.0 µm) that correspond to the Spitzer Space Telescope IRAC bands.
In the text

Fig. 11 Radiated energy per unit time (in solar luminosities) emitted by dust in the SN ejecta as a function of wavelengths, for nine different epochs (see Eq. (18)). The unattenuated luminosity (dashed line) corresponds to the total energy per unit time (from the dust), while the observable luminosity (solid line) is the luminosity that escapes the system, after taking the progressive escape probabilities into account, and will be the one that can be observed (see Sect. 4 and Eqs. (17) and (18) for details).

In the text

	Fig. 12 Geometry of anisotropic ejecta and the corresponding zones, in the form of an ellipsoid. To emulate SN 1987A, we chose the velocity ratios to be 2:3:1 (with reference to Kjær et al. 2010), where the velocity distribution along the mean axes (taken as the x axis in the figure) remain the same as the isotropic model.
In the text

	Fig. 13 Distribution and evolution of temperature with respect to the various zones are shown for the ellipsoidal ejecta. Comparing to the same distribution and evolution (Fig. 2) for the spherical model, it is evident that the temperatures are spread to wider ranges in this case, as expected.
In the text

	Fig. 14 Mass of dust for the ellipsoidal case is presented. Top: evolution of total dust masses, in the case of anisotropic ejecta, are shown for each dust species, and the sum is shown by the solid-line. Bottom: mass of dust produced in the spherical and ellipsoidal cases are compared for each dust species.
In the text

Fig. 15 Dust column densities are shown for the ellipsoidal model for five different epochs. In the top left panel, the chosen polar coordinates are shown in terms of θ and φ for an ellipsoidal shell. The velocity of any clump at point P inside the shell is calculated. In the top right panel and the bottom panel, we show the evolution of dust column densities (from the center to the outer edge of the ejecta) along various directions as a function of the viewing angle, in the three axial planes. See Sect. 5.1 for details.

In the text

Fig. 16 Progressive escape probabilities, corresponding to the anisotropic case, are shown for two epochs, days 1000 and 1500, for the 3.6- and 8.0-µm Spitzer IRAC bands. The other two IRAC bands, 4.5 and 5.8 µm, also show a similar behavior. The escape probabilities are presented in the top panel as a function of ejecta velocities, while in the bottom panel, we show the probabilities as functions of the normalized He-core velocities (normalized with respect to the velocity of the outer edge of the He core in that direction). The escape probability in this scenario is dependent on the viewing angle; here we present the three axial directions. See Sect. 5.2 for more details.

In the text

	Fig. 17 Unattenuated luminosities, as well the observable luminosities, along the three axial planes are shown for nine different epochs. For reference, in our chosen geometry, the z axis is the smallest axis and the y axis is the largest axis of the ellipsoid. See Sect. 5.2 for the description.
In the text

Fig. 18 Fluxes obtained from our calculations (assuming the SN to be at the distance of SN 1987A) for days 615 and 775, are co-plotted with the near- and mid-IR data of SN 1987A (shaded region in yellow) from Wesson et al. (2015). In the top panel, the fluxes for the isotropic model, as well as the three axial directions of the anisotropic model, are shown for the two epochs. In the bottom panel, we show the viewing angle with respect to the coordinates shown in Fig. 15, where we find the fluxes emerging from the anisotropic ejecta match well with the observations. For the details, see Sect. 5.3.

In the text