Survival of the Fittest: Numerical Modeling of SN 2014C

Felipe Vargas; Fabio De Colle; Daniel Brethauer; Raffaella Margutti; Cristian G. Bernal

doi:10.3847/1538-4357/ac649d

1. Introduction

While mass loss is one of the key mechanisms regulating the evolution of massive stars, a complete understanding of it is still missing, especially during the final phases before the supernova (SN) explosion (e.g., Smith 2014). The mass-loss history ≲100–1000 yr before core-collapse SN explosions can be inferred from the radio and X-ray emission resulting from the propagation of the SN shock through the circumstellar material (for a review see, e.g., Chevalier & Fransson 2017). While the bulk of the SN ejecta emits in the optical, the shock-heated gas resulting from the interaction with the environment might be observed in X-ray and radio owing to bremsstrahlung and synchrotron radiation from relativistic electrons accelerated at the shock front.

The forward shock (FS) of Type Ib/c SNe, originated from pre-SN stars with high mass loss before collapse, interacts with the wind of the progenitor star, which has typical mass-loss rates of ${\dot{M}}_{w}\sim {10}^{-4}-{10}^{-6}$ M_⊙ yr⁻¹. The wind is often not smooth, as proven by radio emission, showing small flux fluctuations of ∼a few over timescales of tens to hundreds of days after the explosion (e.g., Bietenholz & Bartel 2005; Soderberg et al. 2006; Schinzel et al. 2009; Wellons et al. 2012; Salas et al. 2013; Bietenholz et al. 2014; Corsi et al. 2014; Palliyaguru et al. 2019).

In a few cases, Type Ib/c SNe show signs of much stronger interaction between the ejecta and shells of material ejected before the explosion. For instance, SN 2006jc exploded inside a dense He-rich environment (e.g., Foley et al. 2007), likely produced by an outburst ejected ∼2 yr before the SN event. In addition, SN 2001em, initially classified as a Type Ib SN, presented prominent Hα emission lines at 2.5 yr. Associated with strong radio and X-ray emission, this was interpreted as evidence of the interaction between the SN ejecta and a massive (∼3 M_⊙) hydrogen-rich shell located at ∼7 × 10¹⁶ cm from the progenitor iron core (Chugai & Chevalier 2006; Chandra et al. 2020). Several other SNe show similar signs of early interaction with massive shells (see, e.g., Anupama et al. 2005; Dwarkadas et al. 2010; Ben-Ami et al. 2014; Moriya & Maeda 2014; Chen et al. 2018; Mauerhan et al. 2018; Pooley et al. 2019; Suzuki et al. 2021). Nevertheless, while in all these cases the intermediate phases of the transition between Type Ib and Type IIn SNe were not observed, this transition has been observed in detail in SN 2014C.

Discovered by the Lick Observatory Supernova Search (LOSS; Kim et al. 2014) in the NGC 7331 galaxy at a distance of 14.7 Mpc and initially classified as a Type Ib SN, SN 2014C made a transition to a Type IIn SN ∼100 days after the explosion, showing strong Hα emission (Milisavljevic et al. 2015). The modeling of the optical/UV light curve shows that SN 2014C has a kinetic energy of (1.75 ± 0.25) × 10⁵¹ erg, with an ejecta mass of 1.7 ± 0.2 M_☉ and a nickel mass of 0.15 ± 0.02 M_☉ (Margutti et al. 2017).

SN 2014C has also been extensively observed with the X-ray Telescope (XRT; Burrows et al. 2005) on board the Neil Gehrels Swift Observatory (Gehrels et al. 2004) and the Chandra X-ray Observatory (CXO) in the 0.3–10 keV energy band, as well as by the Nuclear Spectroscopic Telescope Array (NuSTAR) from 3 to 79 keV (Margutti et al. 2017; Brethauer et al. 2020, 2021). Most of the detected X-ray emission is concentrated in the 1–40 keV energy band, while emission below ∼1 keV is strongly absorbed.

The observed X-ray emission increased at 250 days, although the lack of temporal coverage between 100 and 250 days implies that likely the onset of the X-ray increase was at a smaller time. Integrated over the spectral range 0.3–100 keV, it increased from 5 × 10³⁹ erg s⁻¹ to 5 × 10⁴⁰ erg s⁻¹. Then, it peaked at ∼850–1000 days and maintained a nearly constant flux until ∼2000 days after the explosion (Brethauer et al. 2020).

Radio observations showed a similar behavior. Light curves of the SN at 15.7 GHz taken with the Arcminute Microkelvin Imager (AMI) between 16 and 567 days showed that the flux increased rapidly at ∼100–150 days (Anderson et al. 2017). A similar increase (by more than one order of magnitude) was observed by the Very Large Array (VLA; Margutti et al. 2017). Furthermore, observations done by using very long baseline interferometry (VLBI) found that the shock expansion has already strongly decelerated 384 days after the explosion (Bietenholz et al. 2018, 2021).

The evolution of radio, optical, and X-ray emission has been interpreted by considering the interaction of the SN ejecta with a low-density steady wind with ${\dot{M}}_{w}=5\times {10}^{-5}$ M_☉ yr⁻¹ (assuming a wind velocity v_w = 1000 km s⁻¹; see Anderson et al. 2017; Margutti et al. 2017) and a higher-density, extended circumstellar medium (CSM). Although most previous studies agree on the position of the denser shell (e.g., R_shell ∼ (5–6) ×10¹⁶ cm from the progenitor iron core; Milisavljevic et al. 2015; Anderson et al. 2017; Margutti et al. 2017; Tinyanont et al. 2019; Sun et al. 2020), different models have been proposed with respect to the shell thickness, mass, and structure. Consistently with the fact that the interaction lasts several years, Bietenholz et al. (2021), Bietenholz et al. (2018), Margutti et al. (2017), and Milisavljevic et al. (2015) interpreted the data as evidence of the interaction between the SN shock and a thick and massive shell with a size ΔR_shell ∼ R_shell and a mass M_shell ∼ 3.00 M_☉. On the other hand, Harris & Nugent (2020) consider a progenitor with a lower mass-loss rate, which formed a "wall" at R_w ∼ 10¹⁶ cm and resulted in a thin shell with size ΔR_w ≲ 0.25R_w, corresponding to a lower mass, M_w = 0.04–0.31 M_☉. A similar structure, with a thin, high-density shell and an extended outer medium, is also inferred by Tinyanont et al. (2019).

In order to explain the presence of a dense CSM around the progenitor of SN 2014C, two possible scenarios have been suggested. In the first scenario, it is assumed that the progenitor of SN 2014C was a star with a steady mass-loss rate, resulting in a stratified, low-density environment with a density $\rho \propto {\dot{M}}_{w}{r}^{-2}$ (e.g., Anderson et al. 2017; Margutti et al. 2017; Bietenholz et al. 2018, 2021; Harris & Nugent 2020). In this scenario, the mass-loss rate had an abrupt increase roughly tens to hundreds of years before the explosion (with ${\dot{M}}_{w}\sim {10}^{-3}-{10}^{-2}$ M_☉ yr⁻¹ assuming a wind velocity v_w ∼ 100 km s⁻¹), creating a spherically symmetric expanding high-density CSM (with a density ∼10⁶ cm⁻³). Sun et al. (2020) depicted a different scenario, in which the dense CSM results from an 11 M_⊙ star stripped of its outer envelope by binary interaction.

In synthesis, in addition to the progenitor star, the structure, extension, and mass of the environment leading to the X-ray emission have been extensively debated. Clarifying the structure of the dense CSM leading to the X-ray emission, i.e., its density profile and total mass, is the main aim of this paper.

Radio emission and X-ray emission from the SN shock are typically described by considering a self-similar behavior for the dynamics (e.g., Chevalier 1982; Chevalier & Liang 1989; Chevalier 1998). In the case of an SN interacting with a massive shell, we cannot assume a self-similar behavior for the flow, as the interaction with the shell leads to the formation of a strong reverse shock (RS). Then, although an analytical description gives a qualitative picture of the problem, a detailed understanding of the dynamics can be achieved only employing numerical simulations.

Typically, running a large number of models covering a predefined number of parameters provides sufficient detail of the dynamics. When running computational intensive hydrodynamical simulations, however, this approach necessarily limits the parameter space explored. To handle this limitation, in this paper we use an optimization method (a genetic algorithm (GA)) to iteratively determine a relatively large number of parameters by running a limited number of numerical simulations. Although there have been a wide range of applications of GAs in astrophysics (e.g., Charbonneau 1995), it is the first time that this method has been applied by coupling hydrodynamic simulations with radiative transfer, in order to obtain an optimized solution for the CSM density profiles.

Specifically, we run a set of numerical simulations coupled to a GA. The GA iteratively adapts the values of the shell density by minimizing the difference between the synthetic X-ray emission (computed by post-processing the results of the hydrodynamical simulations) and observations of this SN presented in Margutti et al. (2017), Brethauer et al. (2020), and Brethauer et al. (2021). In this approach, the physical space is divided into several parts (according to the epochs of observations), and each part has a density ρ(r_i) (with i = 1, ..., N) that is considered a parameter of the model.

The paper is organized as follows: In Section 2 we describe the hydrodynamics code and the initial conditions of the simulations, the bremsstrahlung radiation transfer code, and the GA employed to solve the optimization problem and to find the density stratification of the shell. In Section 3 we present the results of our numerical calculations. In Section 4 we discuss the limits of the simulations presented and the implications of our findings. Finally, in Section 5 we draw our conclusions.

2. Methods

2.1. Numerical Codes and Initial Conditions

We study the interaction of an SN shock with a massive shell by running a set of one-dimensional (1D), spherically symmetric simulations with the Adaptive Mesh Refinement (AMR) code Mezcal (De Colle et al. 2012). The code solves the special relativistic hydrodynamics equations and has been extensively used to run numerical simulations of astrophysical flows (e.g., De Colle et al. 2014, 2018, 2018; González-Casanova et al. 2014).

We follow the propagation of the SN shock front as it moves through a computational grid covering ∼nine orders of magnitude in space, from ∼2 × 10⁸ cm (the outer edge of the iron core) to ∼2 × 10¹⁷ cm. We do so by running two sets of simulations. First, we follow the propagation of the SN shock front as it moves through the progenitor star, e.g., from ∼10⁸ to 10¹¹ cm. Then, after remapping the results of the small-scale simulation into a much larger computational box, we follow its propagation through the wind of the progenitor star and its interaction with the massive shell located at ≳10¹⁶ cm.

We set the density profile of the progenitor star by using the E25 pre-SN model from Heger et al. (2000). This corresponds to a star that has lost its outer envelope. The resulting star has a mass of 5.45 M_⊙ and a radius of ∼3 × 10¹⁰ cm.

The computational grid extends radially from 2.2 × 10⁸ cm to 6.6 × 10¹¹ cm. We employ 20 cells at the coarsest level of refinement, with 22 levels of refinement, corresponding to a resolution of 1.6 × 10⁴ cm. We set reflecting boundary conditions at the inner boundary and outflow boundary conditions at the outer boundary of the computational box. The SN energy (E_SN ≈ 10⁵¹ erg) is imposed by setting the pressure of the two inner cells of the computational box as $p={E}_{\mathrm{SN}}({{\rm{\Gamma }}}_{\mathrm{ad}}-1)/V$ , where Γ_ad is the adiabatic index and V is the volume of the two cells. We employ a variable adiabatic index Γ_ad, bridging between the relativistic Γ_ad = 4/3 and the Newtonian Γ_ad = 5/3 cases (see Section 2.4 of De Colle et al. 2012).

Figure 1 shows the initial density stratification of the star and of the ejecta at 50 s. The SN energy has been chosen to match that derived by Margutti et al. (2017), while the ejecta mass is larger (3.5 M_⊙ in our simulations vs. 1.7 ± 0.2 M_⊙ inferred by Margutti et al. 2017). Anyway, we notice that (1) the mass inferred depends strongly on the (uncertain) opacity coefficient κ (see Maund 2018, and references therein) and (2) as we set a reflective boundary condition at the inner boundary, our simulations do not include fallback, which could decrease the ejecta mass.

Outside the stellar surface, we take $\rho ={\dot{M}}_{w}/(4\pi {r}^{2}{v}_{w})$ , where ${\dot{M}}_{w}=5\times {10}^{-6}$ M_⊙ yr⁻¹ and v_w = 10⁸ cm s⁻¹ are the mass-loss rate and velocity of the wind launched by the star before the collapse, respectively. As the velocity of the SN shock front is about two orders of magnitude larger than the stellar wind (i.e., ∼10¹⁰ cm s⁻¹ vs. ∼10⁸ cm s⁻¹), we assume that the wind medium is static. The propagation of the SN shock front is followed as it breaks out of the progenitor star and arrives at 4.5 × 10¹¹ cm in 50 s.

In the large-scale simulations, the computational box goes from 10¹⁰ cm to 5 × 10¹⁷ cm. For r < 4.5 × 10¹¹ cm we set the density, pressure, and velocity by using the values determined from the small-scale simulation. We save outputs each ∼20 simulated days and at the times corresponding to the available observational epochs. For larger radii, we take the density stratification as $\rho ={\dot{M}}_{w}/(4\pi {r}^{2}{v}_{w})$ , with ${\dot{M}}_{w}=5\times {10}^{-6}$ and v_w = 10⁸ cm s⁻¹ as in the small-scale simulation. We employ 150 cells, with 20 levels of refinement in the AMR grid, corresponding to a resolution of 2.5 × 10⁹ cm. By running different simulations in which we change the number of levels of refinement between 14 and 22 levels, we verify that 20 levels of refinement are enough to achieve convergence.

2.2. Bremsstrahlung Emission Code

To compare models with observations, we compute the X-ray bremsstrahlung emission coming from the shocked plasma by post-processing the results of the hydrodynamical simulations. We assume that the shocked material is completely ionized⁵ so that it does not contribute to the bound-free opacity and the unshocked shell is neutral. Extending the radio synchrotron emission to X-rays by assuming F_ν ∝ ν^−(p−1)/2 ∝ ν⁻¹ with p ∼ 3, Margutti et al. (2017) showed that the contribution of synchrotron emission to the X-ray flux is negligible (∼2 orders of magnitude smaller than the observed values), suggesting a thermal origin for the X-rays.

The specific flux is given by

$\begin{eqnarray}&&{F}_{\nu }=\int {I}_{\nu }\cos \theta d{\rm{\Omega }},\end{eqnarray} \tag{ 1 }$

where $d{\rm{\Omega }}=2\pi \sin \theta d\theta /{D}^{2}$ , with D = 14.7 Mpc the luminosity distance from SN 2014C (Freedman et al. 2001), and I_ν is the specific intensity.

The observed X-ray radiation is due to thermal bremsstrahlung emission caused by the interaction between the SN shock and a massive shell (Margutti et al. 2017). To determine the specific intensity, we solve the radiation transfer equation

$\begin{eqnarray}&&\displaystyle \frac{{{dI}}_{\nu }}{d{\tau }_{\nu }}={S}_{\nu }-{I}_{\nu },\end{eqnarray} \tag{ 2 }$

where S_ν = j_ν/α_ν is the source function, τ_ν = ∫α_ν dl is the optical depth, and j_ν and α_ν are the emissivity and the absorption coefficient, respectively. The radiation transfer equation is solved by a standard approach (see, e.g., Section 9.2 of Bodenheimer et al. 2007). Shortly, we created a two-dimensional grid in cylindrical coordinates. We mapped the results of our 1D simulations on this grid and integrate along 200 rays. We compute the intensity in the plane of the sky, whose integration leads to the flux (see Equation (1)).

In addition to the bremsstrahlung (free–free) self-absorption, we also consider bound–free absorption, which at early times dominates absorption for frequencies ≲10¹⁸ Hz, so that

$\begin{eqnarray}&&{\alpha }_{\nu }={\alpha }_{\nu ,\mathrm{ff}}+{\alpha }_{\nu ,\mathrm{bf}},\qquad {j}_{\nu }={j}_{\nu ,\mathrm{ff}}.\end{eqnarray} \tag{ 3 }$

To compute the bound–free absorption, we use tabulated cross sections for solar metallicity (Morrison & McCammon 1983). For the bremsstrahlung coefficients we take (e.g., Rybicki & Lightman 1986)

$\begin{eqnarray}&&{j}_{\nu }=\displaystyle \frac{6.8}{4\pi }\times {10}^{-38}{Z}^{2}{n}_{e}{n}_{i}{T}^{-1/2}{e}^{-\tfrac{h\nu }{{kT}}}G\end{eqnarray} \tag{ 4 }$

$\begin{eqnarray}&&{\alpha }_{\nu }=3.7\times {10}^{8}{T}^{-1/2}\displaystyle \frac{{Z}^{2}{n}_{e}{n}_{i}}{{\nu }^{3}}\left(1-{e}^{-\tfrac{h\nu }{{kT}}}\right)G,\end{eqnarray} \tag{ 5 }$

where Z, n_e, and n_i are the atomic number and electron and ion densities, respectively, and G is the Gaunt factor (∼1 for the range of parameters considered here). All other variables have their usual meaning, and everything is in cgs units. The electron/ion densities and the temperature are directly determined from the numerical simulations by assuming that the ejecta is an ideal gas with solar metallicity. We also assume that electrons and ions have the same temperature. Cooling is not included in the numerical simulations, as the bremsstrahlung cooling timescale is t_ff = 60(T/10⁸ K)(n_e/10⁶cm⁻³) yr, so it is much larger than the timescales studied here.

2.3. Computing the Shell Parameters from a Grid of Models

As a first attempt to reproduce the observed X-ray emission, we have first considered a massive, uniform cold shell located at radius R_s, with mass M_s and thickness ΔR_s. To determine the best values for these three parameters, we run a grid of 1815 models by using 11 values of M_s (in the range 1.5–2.5 M_⊙), 15 values of ΔR_s (ranging from 7.5 × 10¹⁴ cm up to 2.5 × 10¹⁵ cm), and 11 values of R_s (from 1.8 × 10¹⁶ cm to 7.8 × 10¹⁶ cm). To compare the results of the numerical simulations with observations, we have then employed the ray-tracing code described in Section 2.2.

Unfortunately, none of the models provided a reasonable fit of the data. In particular, shells with larger density (i.e., larger mass and smaller thickness) reproduce well the early observational epochs, while shells with smaller densities (i.e., smaller mass and thicker) reproduced well the late observational epochs. This indicates that the shell density is not constant and has a decreasing density from smaller to larger radii.

We notice that, as the origin of the massive shell is not understood, it could in principle have an arbitrary shape. In this case, to find the density at several radii is a task that is not possible to achieve with a grid of numerical models. For instance, a grid of 10 values for the density at 10 different radii implies running 10¹⁰ simulations. Thus, to determine the density stratification of the shell, we decided to solve the "full" optimization problem. This is done by coupling the Mezcal code with two other codes: a radiation transfer code that computes the bremsstrahlung radiation (see Section 2.2), and a GA (described in Section 2.4) that automatically and randomly changes the density profile inside the shell by minimizing the variance between the synthetic observations computed from the numerical model and the X-ray observations, allowing the density profile to be completely arbitrary.

2.4. Genetic Algorithm

GAs (e.g., Rajpaul 2012) are based on the theory of natural selection. GAs are commonly used optimization methods,⁶ employed also in astrophysics to solve problems with many degrees of freedom, in which finding the optimal solution would be very hard otherwise (e.g., Cantó et al. 2009; De Geyter et al. 2013; Morisset 2016). Nevertheless, this is the first time, as far as we know, that optimization methods are coupled directly to hydrodynamical simulations.

We run a single small-scale simulation following the propagation of the shock wave in the interior of the star. Then, we use the SN ejecta structure given by this simulation to initialize large-scale numerical simulations. The GA (see Figure 2 for a schematic description of the algorithm) is then used to find the best density stratification of the shell by running several simulations of the interaction of the SN ejecta with the massive shell.

**Figure 2.** Flowchart showing the GA employed to determine the density stratification of the massive shell interacting with SN 2014C. First, we initialize a set of shell density profiles. We modify the initial densities by mutations and crossover (see the text for details). We run hydrodynamical simulations following the interaction of the SN ejecta with the shell and compute synthetic X-ray spectra, which are compared with the observations. We select the "best" simulations using a "fitness function" (the χ² test). The process is repeated until convergence.
Download figure:
Standard image High-resolution image

The initial conditions of the large-scale simulations are described in Section 2.1. The shell structure is initialized in the numerical simulations by a power-law interpolation of the densities ρ(r_i) (with i = 1,...,11) defined at 11 radial positions r_i. Nine of the radii r_i are determined directly from the numerical simulations and are given by the location of the shock at the nine available observational epochs. Additionally, we define two densities, one at a smaller radius (r₁) and one at a larger radius (r₁₁). While the value of r₁ = r₂/2 is fixed, r₁₁ is allowed to vary to get the necessary absorption.

Thus, the 11 values of density (and the corresponding 11 radial positions) represent the parameters of the optimization process. We do not impose any constraint on the value of the densities. At the beginning of the iterative process, 10 constant-density shell structures are defined, with equally spaced values (defined in log space) between n = 10² cm⁻³ and n = 10⁸ cm⁻³ and located at a radius r_i = 2 × 10¹⁶ cm to 4 × 10¹⁶ cm. This choice of the initial parameters is completely arbitrary. Indeed, as usual in optimization methods, the final results are independent of the initial choice of the parameters.

Then, we create 90 new shell structures. Half of the densities ρ(r_i) are defined by randomly choosing two shell structures from the initial set and applying to them crossover and mutation (see below), while the other half are initialized by directly copying the density values from one random shell structure and modifying it only by mutation. In GA, the fraction of the parameters that should change as a result of mutation and crossover depends on the particular problem studied. Empirically, we found that, for the optimization problem discussed in this paper, setting equal probability for the crossover and mutation processes was producing reasonably fast results. Anyway, we did not attempt to find the fraction of mutations/crossover that reduces the number of iterations to a minimum. Optimizing it and replacing the GA with more efficient optimization methods could reduce substantially the final computational time.

In the crossover process, two values of density taken from a previous iterative step are mixed by choosing randomly values of densities from each shell structure (e.g., the first and second densities from the first density profile, the third density from the second profile, and so on). This process is inspired by the genetic mixing present in biological evolution. In the mutation process we modify randomly one density in each profile. We do so by setting a Gaussian distribution around the original value of the density ρ₀, with a width given by ρ_i/2 in 90% of the cases and by 10ρ₀ in 10% of the cases, so that in a few cases the system explores density values far away from the initial one (to avoid being trapped by a local minimum).

The shell densities are then mapped as initial conditions into the Mezcal code in each step. Then, after each simulation is completed, the bremsstrahlung X-ray emission is computed by post-processing the results. A fitness function (a reduced χ² test) is applied to compare the synthetic spectra produced by the model and the observational data at nine epochs: 308, 396, 477, 606, 857, 1029, 1257, 1574, and 1971 days after explosion. The fittest 10 simulations (out of the new 90 and the old 10 sets of simulations) are then used as the initial condition for a new step.

This process converged in ∼150 iterations, for a total number of ∼10⁴ simulations, then reducing the number of simulations by a factor of 10⁶ with respect to a grid of simulations covering a similar number of parameters. Each simulation took ∼10 minutes so that the entire process could be completed in less than a day.

The simulations were done on a cluster of CPUs, by using the "Message Passing Interface" library. At each iteration, the master node initialized and synchronized the simulations, selected the best simulations/shell structures, and managed the crossover/mutation processes, while the other nodes ran (in parallel) each of the 90 hydrodynamics simulations and computed the X-ray spectra and the fitness function (a χ² test).

Figure 3 shows the evolution of the fitness function (the normalized χ²) for the best and worst simulation (among the 10 fittest simulations) during each iteration step. The solution reaches a minimum of χ² ≈ 3.03. The evolution of the density profile of the massive shell as a function of radius and iteration is shown in Figure 4. At the beginning of the iteration process, the solution fluctuates substantially between different iterations, to then converge to a final density profile (shown in blue in the figure). The density profile and the value of χ² are discussed in more detail below.

**Figure 4.** Evolution of the number density profiles (in units of cm⁻³) of the outer shell during the iteration process (black lines, with increasing opacity from lower to higher iteration numbers). The final density profile is shown in blue.
Download figure:
Standard image High-resolution image

3. Results

3.1. SN Shock Propagation through the Progenitor Wind

The propagation of the shock wave through the star has been extensively studied for both the nonrelativistic regime (e.g., Sakurai 1960; Matzner & McKee 1999) and the mildly relativistic regime (e.g., Tan et al. 2001; Nakar & Sari 2010; Waxman & Katz 2017). As the shock approaches the surface of the progenitor star and moves into the stellar envelope, in which the density drops steeper than ρ ∝ r⁻³, the shock velocity increases to mildly relativistic speeds.

As a result, once the SN shock breaks out in ∼25 s, most of the mass (and energy) of the SN moves at velocities ∼10⁴ km s⁻¹, while a small fraction of the mass (corresponding to a kinetic energy ≈10⁴⁷ ergs) expands with larger velocities (up to v_sh ∼ 0.5c in our simulations). Self-similar solutions describing the propagation of the SN shock through a polytropic envelope (with ρ_env ∝ (R/r − 1)^k) predict an ejecta density stratification ρ ∝ r⁻ⁿ after the breakout, with n = 7–11, depending on the structure of the progenitor star. We get a nearly constant density profile in the inner ejecta (containing most of the mass; see Figure 5), while the outer ejecta has a steep density profile, i.e., ρ ∼ r^−10.3.

Then, we run a set of simulations by using as input the outcome of the small-scale simulations, i.e., density, pressure, and velocity profiles. As described in Section 2.4, we compute the shell density profile that minimizes the variance between the bremsstrahlung emission computed from the numerical model and the observations. In the following, we discuss the evolution of the SN shock front while it propagates into the progenitor wind and interacts with the shell density profile obtained after the GA algorithm has converged.

Figure 5 shows the evolution of the ejecta before interacting with the outer shell. The expansion of the SN ejecta through the progenitor wind leads to the formation of a double-shock structure, formed by the FS, which accelerates and heats the WR wind, and the RS, which decelerates and heats the SN ejecta. Following Chevalier (1982), the evolution of the FS and RS is self-similar, with R ∼ t^m, where m = (n − 3)/(n − 2). By taking n = 10.3, we get m ∼ 0.88, which is consistent with the evolution of the shock wave obtained in our simulation.

In our simulations, the SN shock acceleration stops only when the shock arrives at the very edge of the progenitor star. This leads to an overestimation of the true shock velocity, as the shock acceleration should stop once the stellar envelope becomes optically thin to the radiation coming from the post-shock region. A detailed calculation of the shock acceleration and breakout is an open problem, which requires a radiation hydrodynamics code. Once the SN ejecta interacts with the shell (see below), the shock velocity quickly drops. Then, the late evolution of the system will be independent of the particular shock velocity obtained, while it will depend on the energy stratification of the SN ejecta, i.e., on the properties of the explosion (energy, ejected mass, fallback, and so on) and of the progenitor star. A detailed study of it could potentially constrain the progenitor star and is left for a future study.

3.2. SN Shock Interaction with the Massive Shell

At ∼50 days after the explosion, the shock begins interacting with the massive shell. Radio and optical observations showed that the interaction started ∼100 days after the explosion (Milisavljevic et al. 2015; Anderson et al. 2017). Then, our simulations overestimate the average shock velocity by a factor of ∼2 with respect, e.g., to radio observation by Bietenholz et al. (2018) or to the shock velocity inferred by Milisavljevic et al. (2015) from absorption lines at 10 days. A lower shock velocity can be due, as mentioned above, to the loss of thermal energy during the shock breakout or, alternatively, to a less steep density profile in the outer layers of the progenitor stars.

The interaction of the ejecta with the shell is shown in Figures 6 and 7. The shell presents large density fluctuations (see the top panel of Figure 7). Then, the shock propagation leads to the formation of strong RSs, whose interaction produces the complex shock structure observed in Figure 7. Once interacting with the shell, the shock velocity⁷ quickly drops to ∼7500 km s⁻¹ and then maintains an approximately constant velocity (see Figure 6). Small fluctuations in the shock velocity are present at ∼10³ days (see the bottom panel of Figure 3), with increases (drops) in velocity by ∼3000 km s⁻¹ corresponding to drops (increases) in the density. The velocities inferred in our model differ from those estimated by VLBI observations (Bietenholz et al. 2018) by ≈50% (slightly larger than 1σ; see the middle panel of Figure 6).

**Figure 6.** Top panel: position of the SN shock front as a function of time. The break seen at ∼50 days corresponds to the beginning of the interaction with the shell. The thin vertical lines corresponds to the epochs observed in X-rays. The radii inferred by VLBI observations (Bietenholz et al. 2018) at 8.4 and 22 GHz are also included in the plot. Middle panel: shock velocity as a function of time. During the self-similar phase, the shock velocity drops from ∼0.35c to ∼0.15c. After interacting with the shell, it drops to ∼0.025c. The inset shows a zoom-in corresponding to t ≲ 2000 days. Velocities inferred from VLBI observations (Bietenholz et al. 2018) are also included in the plot. Bottom panel: evolution of the FS and RS (in mass coordinates).
Download figure:
Standard image High-resolution image

The bottom panel of Figure 7 shows the evolution of the FS and RS in mass coordinates. Before interacting with the dense CSM, the low-mass post-shock structure remains thin, with the FS and RS very close to each other. While the ejecta interact with the dense CSM, more and more mass is accumulated in the post-shock region, separating the two shocks when seen in mass coordinates. By the end of the simulation (∼2000 days after the explosion), most of the CSM mass and about half of the ejecta mass have been shocked.

**Figure 8.** X-ray light curve at 0.1–12.5 keV. The observational data (stars) are compared with the models computed by the GA. The three models shown correspond to a cross section computed assuming solar metallicity (indicated by "Model 2" in the figure), a model computed assuming half of that cross section ("Model 1"), and a simplified r⁻³ density profile model, where the formation of the bumps was not allowed, showing that the bumps are necessary to fit observational data.
Download figure:
Standard image High-resolution image

The velocity plot (third panel from the top in Figure 7) illustrates the evolution of the different components, i.e., the unshocked CSM (the stationary medium at large radii), the unshocked ejecta (the region with increasing velocities at small radii), and the shocked CSM and ejecta (the region with fluctuating velocities between the unshocked ejecta and CSM).

Initially, the RS is stronger than the FS. Thus, the shocked ejecta is hotter than the shocked wind, i.e., the progenitor wind accelerated and heated by the SN shock front (Figure 7, blue and green lines in the middle panel). As the ejecta crosses the RS and approaches the RS velocity, the RS temperature drops, becoming smaller than the FS temperature. Thus, the bremsstrahlung specific emission (larger at smaller temperatures) is initially dominated by the shocked wind and later the bremsstrahlung emission is dominated by the shocked ejecta. We also notice that the unshocked ejecta mainly cools by adiabatic expansion.

A fit to the density profile (see Figure 5, bottom panel) gives ρ ∝ r^{−3.00±0.06}, consistently with the constant shock velocity seen in Figure 7 (as E ∼ Mv² ∼ R³ ρ v² implies that the velocity is constant as long as ρ ∝ R⁻³ and the shock is adiabatic), although the density profile is clearly more complex than what is predicted by a power-law fit, as it presents fluctuations larger than ∼1 order of magnitude, apparently being formed by the joining of three different shells that are centered at 4 × 10¹⁶ cm, 7.5 × 10¹⁶ cm, and 13.75 × 10¹⁶ cm. The shell mass is 2.6 M_⊙. This value is consistent with the 3.0 ± 0.6 M_⊙ (VLB model) determined by Brethauer et al. (2021), which also assumes spherical symmetry and is within the range of typical shell masses observed in Type IIn SNe (0.1–10 M_⊙; see, e.g., Branch & Wheeler 2017; Smith 2017).

One of the parameters of the GA was the 11th density specified in the shell structure (corresponding to the dashed region of Figure 7). This was the average density of shell material not yet reached by the SN shock front. This density allows us to determine the amount of neutral hydrogen still to be crossed by the ejecta (inferred from the bound-free absorption). We get M = 0.38 M_⊙ for this component at t = 1971 days, which is represented by a constant-density dashed line in the top panel of Figure 7. We notice that the exact structure of this region cannot be determined. Nevertheless, it will extend to r ∼ 2.3 × 10¹⁷ cm if the shell density continues dropping as r⁻³, in which case (moving at a constant speed as discussed above) the SN shock will break out of it ∼8.5 yr after the explosion. Late X-ray and radio observations will help to constrain the outer structure of the shell. We also notice that Milisavljevic et al. (2015) constrained the density of the unshocked CSM to be <10⁷ cm⁻³, which is consistent with the results shown here, and the temperature to be ≈(2–8) × 10⁴ K.

3.3. Radiation

Figure 8 shows the total integrated flux over time and Figure 9 shows the X-ray spectra at 308, 396, 477, 606, 857, 1029, 1257, 1571, and 1971 days and the radiation emitted as a result from the models obtained by employing the GA (Section 2.4). Dashed lines are computed by assuming solar metallicity in the neutral shell, whereas the solid lines consider half of the cross section of the material used in the solar-metallicity model. We include in the figure also the X-ray emission computed by considering a simple r⁻³ density profile obtained by using as parameters of the GA the normalization factor (i.e., the parameter A in the equation ρ = A r⁻³) and the position and width of the shell. The fit is clearly poor in this case, implying that the "bumpy" structure inferred from the GA algorithm is not due to an overfitting of the X-ray data. We also notice that the presence of small-scale fluctuations in the X-ray observations is not generated in our simple 1D model and can be due to the presence of a clumpy medium, which would require a multidimensional study, or to the presence of unresolved lines (this explains the normalized χ² larger than 1 obtained in our optimizing procedure).

The best fit, in this case, predicts a larger absorption by a factor ≲5 at 2 × 10¹⁷ Hz than observed. Solid lines, which fit better the observational data, are computed by reducing the cross section by a factor of two and then running a new set of simulations using the GA. We notice that the shell density structure obtained from this calculation is nearly identical to the first one, showing the robustness of the result. The implications of this result will be discussed in the next section.

Finally, Figure 8 shows the X-ray flux integrated over the range of frequencies used in the GA, i.e., 0.1–12.5 keV, excluding the region between 5.3 and 7.8 keV, where the emission from Fe lines dominates the X-ray flux. The "bump" models are in general agreement with the observations, and also it is possible to see that observations cannot be explained by a simple r⁻³ density profile.

4. Discussion

The optimization method employed in this work allows us to determine the detailed structure of the shell. As described in the previous section, the shell has a mass of 2.6 M_⊙ (2.2 M_⊙ of shocked and 0.4 M_⊙ of still unshocked gas) and a density stratification that at first order can be approximated as ρ ∝ r⁻³ and extends from 2 × 10¹⁶ cm to 2 × 10¹⁷ cm. While uncommon in Type Ib SNe, Type IIn SNe show evidence of strong interaction, with shell masses of 0.1–10 M_⊙ (see, e.g., Branch & Wheeler 2017; Smith 2017, and references therein). Also, many Type IIn SNe show an X-ray emission inconsistent with a density profile ρ ∝ r⁻², implying a steeper density stratification (Dwarkadas & Gruszko 2012). While sharing many characteristics with Type IIn SNe, in the case of SN 2014C the shell is located at a much larger distance, implying that it was ejected ∼ 60/(v_w/100 km s⁻¹) yr before the SN explosion.

Several papers have inferred the properties of the massive CSM from the observations. Based on the models of Chevalier & Liang (1989), Margutti et al. (2017) modeled the behavior of the X-ray light curve by considering three components: a low-density wind, a high-density shell, and an outer wind with an unknown density profile. Margutti et al. (2017) modeled the high-density shell as a ≈1–1.5 M_⊙, constant-density medium (which is excluded by our numerical simulations). Tinyanont et al. (2019) and Harris & Nugent (2020) presented a similar three components model. These works focused on the high-density shell and outlying material. Tinyanont et al. concluded that ∼1 M_⊙ of material has been shocked by the SN 2014c ejecta by 400 days. Harris & Nugent (2020) presented numerical simulations of an SN ejecta interacting with a wall located at r₀ ∼ 2 × 10¹⁶ cm and with a CSM wind at larger radii. They deduce a small mass for the wall (∼0.04–0.31 M_⊙), while the CSM wind density is assumed to be a factor 4–7 smaller than the density of the wall at r₀. The CSM wind, then, corresponds to a mass loss ${\dot{M}}_{w}\sim {10}^{-2}$ M_⊙ yr⁻¹ (assuming v_w = 10⁸ cm s⁻¹).

Our results (which, we stress, have been obtained by a "blind" fit, i.e., without assumptions on the final shell structure) clarify the structure of the massive CSM. In particular, our results imply that the wall and the CSM wind considered by Harris & Nugent (2020) are both part of the same extended structure (see Figure 7), with the densest shell material located at ∼2 × 10¹⁶ cm and the outer shell density dropping quickly with radius. Thus, we conclude that the same event is responsible for the ejection of the full massive shell.

The parameters determined for SN 2014C are remarkably similar to those inferred for SN 2001em. The X-ray, radio, and Hα emission from SN 2001em has been interpreted as evidence of interaction with a 3 M_⊙ hydrogen-rich shell (Chugai & Chevalier 2006). VLBI observations showed that the shell is located at 7 × 10¹⁶ cm and expanding with a velocity of 5800 ± 10⁴ km s⁻¹ (Bietenholz & Bartel 2007), which is consistent with the 7500 km s⁻¹ inferred here.

The X-ray emission is strongly absorbed at low frequencies. To compute it, we have used a tabulated cross section that assumes solar metallicity and have assumed that the upstream medium is neutral. We have also shown that employing a smaller cross section leads to a better fit at low frequencies. This has interesting implications for the properties of the upstream shell material.

There are two physical main mechanisms that can modify the absorption cross section. In the first one the metallicity can be lower than solar (which, by the way, would also lead to a lower X-ray flux). We notice that a lower metallicity is in contradiction with observations of the Fe line. The prominent Fe emission line at 6.7–6.9 keV is consistent with a metallicity larger by a factor of ∼5 with respect to solar metallicity. This can be reconciled with our results by assuming that the medium is clumpy, with high-density (and lower-temperature) regions responsible for most of the Fe emission (Margutti et al. 2017). In the second one the upstream medium can be partially ionized owing to the strong ionizing X-ray flux coming from the post-shock region, as the optical depth is $\propto {n}_{{{\rm{H}}}^{0}}{\sigma }_{\nu }$ , where ${n}_{{{\rm{H}}}^{0}}$ is the density of neutral hydrogen and σ_ν is the bound-free cross section. This is indeed very likely considering that the X-ray emission and/or UV radiation coming from the shocked shell and SN ejecta can partially ionize the unshocked material. In this case, we expect the mass of the shell to be larger than the value obtained in this paper, although by a small factor since the bremsstrahlung emission is $\propto {n}_{e}^{2}$ . A detailed calculation of the ionization of the shell, left for a future work, will help to clarify which among these two effects is dominant and will potentially help to determine the ionization degree and metallicity of the dense CSM.

Different origins for the shell have been considered. The possibility of the shell being due to a massive wind ejection (e.g., de Jager et al. 1988; Leitherer 2010; Kuriyama & Shigeyama 2020) is unlikely, as it would correspond to an extremely large mass-loss rate of ∼10⁻² M_⊙ yr⁻¹ (v_w/100 km s⁻¹). Other possibilities include a sudden outburst some time before collapse, which would remove the most external layer of the star, where almost all hydrogen is found (Smith & Arnett 2014), or binary system interactions in which the envelope of the most massive star has been stripped away (Sun et al. 2020). A better understanding of the origin of this ejection can be achieved only by detailed theoretical models coupled to a larger sample of observed interacting SNe. Anyway, future models exploring the origin of this massive shell will need to account for the density profile, size, and mass obtained in this work.

The dense CSM presents density fluctuations. It is roughly formed by three consecutive shells, the first one from (2–5) × 10¹⁶ cm, the second from (5–10) × 10¹⁶ cm, and the third from (10–16) × 10¹⁶ cm, with average shell density declining like r⁻³. The density fluctuations have a characteristic length scale of ∼4 × 10¹⁶ cm. A similar length scale has been observed in the radio emission from SN 1979C (Weiler et al. 1991) and has been interpreted as evidence of a binary system in which the orbital motion modulates the wind density (see Yalinewich & Portegies Zwart 2019), which interacts with the stellar outburst. If this is also the case of SN 2014C, it would imply that the binary system is very detached (as the binary period is ≫ 4 × 10¹⁶/v_w ∼ 10 yr ${({v}_{w}/{10}^{8})}^{-1}$ ). The companion star would then be not responsible for the loss of the envelope of the primary star. An alternative explanation is that the density fluctuations seen in the GA fit are the direct result of a modulation in the outburst from the progenitor star.

Finally, we notice that the simulations presented here assume that the shell is spherically symmetric. This is consistent with VLBI observations (Bietenholz et al. 2018). On the other hand, the velocities estimated by the GA are smaller than those inferred from VLBI observations (and, puzzling, larger than those inferred by the Hα emission line, corresponding to roughly a few thousand kilometers per second). The density profile shown in Figure 7 shows large-scale (i.e., much larger than the size of our computational cells) density fluctuations that, in spherical symmetry, correspond to fluctuations among contiguous shells of material. Actually, it is likely that the medium is clumpy, with the clumps necessarily breaking the spherical symmetry. In this case, the interaction with the SN ejecta will amplify the inhomogeneities (see, e.g., Guo et al. 2012; Velázquez et al. 2017), leading to a multiphase medium with denser/colder regions in pressure equilibrium with more tenuous/hotter regions, which could explain the strong Fe emission line observed at ∼6.5 keV (Margutti et al. 2017).

A more realistic scenario, with density stratification in the shell not only in the radial but also along the polar direction, could also explain the apparent discrepancies in the velocities inferred by the different emission processes, with faster material, moving along lower-density regions of the shell, emitting in radio and denser regions, moving at lower speeds, emitting in X-rays and radio. Such a complex (but possibly more realistic) scenario requires multidimensional simulations and is left for a future study.

In addition, if the shell is initially nearly perfectly homogeneous, the contact discontinuity that separates the shocked SN ejecta from the shocked shell is prone to Rayleigh–Taylor (RT) instabilities, so that we can expect the formation of plasma filaments and inhomogeneities in the post-shock region, which cannot be captured by our numerical simulations (but see Harris & Nugent 2020 for an approximated treatment of RT instabilities in 1D simulations). As the bremsstrahlung emission is ∝ Z² n², inhomogeneities in the shell and mixing with the higher-metallicity ejecta lead to a larger emissivity, implying that the mass of the shell should be taken as an upper limit.

5. Conclusions

In this paper, we presented hydrodynamical simulations of the strongly interacting SN 2014C. First, we follow the propagation of an SN shock through the progenitor star. Then, by using as input the outcome of the small-scale simulation (i.e., density, pressure, and velocity profiles), we run a large set of simulations. As described in Section 2.4, we initialize the shell with a uniform density n_shell = 10⁷ cm⁻³.

We followed the propagation of the SN shock as it interacts with the wind launched by the progenitor Wolf-Rayet star and with the massive shell. We computed the bremsstrahlung emission using the algorithm described in Section 2.2 and compared the results with observations. At each step, we run a large number of simulations changing the shell density profile. As a result, we determine the shell structure. In particular, we get a mass of 2.6 M_⊙ for the shell and a density profile roughly dropping as ρ ∝ r⁻³, although presenting large density fluctuations. We also found that the shell is very extended, with a size ≳10¹⁷ cm. If the shell stratification continues with the same average slope, the SN shock will break out of it nearly 8 yr after the explosion, i.e., during 2022.

Although this paper represents one of the most detailed models of the interaction between the SN ejecta and the surrounding massive CSM, several ingredients are still missing from our model. In particular, we have explored only one pre-SN structure; we have modeled only the X-ray emission, while radio emission (especially images obtained by VLBI; see, e.g., Bietenholz et al. 2021) contains detailed information about the evolution of the shock and its interaction with the medium; we have not included a detailed treatment of the photoionization of the upwind medium, which can affect the bound-free absorption. Finally, and most importantly, we assumed that the shell structure is spherically symmetric.

Radio and X-ray emission allows us to understand the mass-loss history of core-collapse SN progenitors on timescales that are impossible to study by direct observations. As we have shown in this paper, optimization methods can be used, coupled with hydrodynamical simulations, to model the density stratification of the environment once data at several epochs are available, as in the case of SN 2014C. The X-ray emission tracks the FS and RS emission, depending on the density of the environment and the ejecta velocity. The Hα emission tracks the shocked shell and the unshocked medium photoionized by the X-ray and UV radiation. Altogether, a detailed fit of the different components can help us to get a better understanding of this system. Then, coupled with detailed modeling of the radio emission, this analysis can allow us to determine the microphysical parameters as a function of time (which are usually degenerate with the density of the environment and ejecta velocity), giving us direct information on the particle acceleration process. In this paper, we describe this technique by analyzing the X-ray bremsstrahlung emission. The extension to radio and optical emission will be considered in a future study.

We thank the anonymous referee for the helpful comments and constructive remarks on this manuscript. We thank Luc Binette, Cesar Fernández Ramírez, Leonardo Ferreira, and Claudio Toledo Roy for useful discussions. F.V. and F.D.C. acknowledge support from the UNAM-PAPIIT grant AG100820. We acknowledge the computing time granted by DGTIC UNAM (project LANCAD-UNAM-DGTIC-281). R.M. acknowledges support by the National Science Foundation under award Nos. AST-1909796 and AST-1944985. Support for this work was provided by the National Aeronautics and Space Administration through Chandra award No. GO9-20060A issued by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under contract NAS8-03060.

Software: Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), Mezcal (De Colle et al. 2012).

Survival of the Fittest: Numerical Modeling of SN 2014C

Article metrics

Author e-mails

Author affiliations

Author notes

ORCID iDs

Dates

Abstract

1. Introduction

2. Methods