Survival of the Fittest: Numerical Modeling of Supernova 2014C

Initially classified as a supernova (SN) type Ib, $\sim$ 100 days after the explosion SN\,2014C made a transition to a SN type II, presenting a gradual increase in the H${\alpha}$ emission. This has been interpreted as evidence of interaction between the supernova shock wave and a massive shell previously ejected from the progenitor star. In this paper, we present numerical simulations of the propagation of the SN shock through the progenitor star and its wind, as well as the interaction of the SN ejecta with the massive shell. To determine with high precision the structure and location of the shell, we couple a genetic algorithm to a hydrodynamic and a bremsstrahlung radiation transfer code. We iteratively modify the density stratification and location of the shell by minimizing the variance between X-ray observations and synthetic predictions computed from the numerical model. By assuming spherical symmetry, we found that the shell has a mass of 2.6 M$_\odot$, extends from 1.6 $\times 10^{16}$ cm to $1.87 \times 10^{17}$ cm, implying that it was ejected $\sim 60/(v_w/100 {\rm \; km \; s^{-1}})$ yrs before the SN explosion, and has a density stratification decaying as $\sim r^{-3}$. We found that the product of metallicity by the ionization fraction (due to photo-ionization by the post-shock X-ray emission) %and/or the SN UV radiation is $\sim$ 0.5. Finally, we predict that, if the density stratification follows the same power-law behaviour, the SN will break out from the shell by mid 2022, i.e. 8.5 years after explosion.


INTRODUCTION
While mass loss is one of the key mechanisms regulating the evolution of massive stars, a complete understanding of it is still missing, specially during the final phases before the supernova (SN) explosion (e.g., Smith 2014).The mass loss history 100 − 1000 yrs before core-collapse supernova 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 optical, the shock heated gas resulting from the interaction with the environment might be observed in X-ray and radio due to bremsstrahlung and synchrotron radiation from relativistic electrons accelerated at the shock front.
The forward shock of type Ib/c SNe, originated from the collapse of Wolf-Rayet progenitor stars, interacts with the wind of the progenitor star, which has typical mass loss rates of Ṁw ∼ 10 −4 − 10 −6 M yr −1 .The wind is often inhomogenous, as proven by radio emission, showing small flux fluctuations of ∼ a few over timescales of tens-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 ∼ two years before the SN.SN 2001em, initially classified as type Ib SN, presented prominent Hα emission lines at 2.5 yrs.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 16 cm (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;Moriya & Maeda 2014;Chen et al. 2018;Pooley et al. 2019;Dwarkadas et al. 2010;Ben-Ami et al. 2014;Mauerhan et al. 2018;Suzuki et al. 2021).Nevertheless, while in all these cases the intermediate phases of the transition between type Ib and type IIn SN were not observed, this transition has been observed in detail in the SN 2014C.
Discovered by the Lick Observatory Supernova Search (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 about ∼ 100 days after the explosion, showing strong Hα emission (Milisavljevic et al. 2015).The modelling of the optical/UV light curve shows that SN 2014C has a kinetic energy of 1.75 ± 0.25 ×10 51 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 extensivel observed with Xrays Telescope (XRT, (Burrows et al. 2005)) on board the Neil Gehrels Swift Observatory (XRT) (Gehrels et al. 2004) and the Chandra X-ray Observatory (CXO) in the 0.3-10 keV energy band, and by the Nuclear Spectroscopic Telescope Array (NuSTAR) from 3 keV to 79 keV (Margutti et al. 2017;Brethauer et al. 2020Brethauer et al. , 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 raised from 5 × 10 39 erg s −1 to 5 × 10 40 erg s −1 .Then, it peaked ∼ 850 − 1000 days and maintained a nearly constant flux until ∼ 2000 days after the explosion (Brethauer et al. 2020).
Radio observations showed a similar behaviour.Light curves of the SN at 15.7 GHz taken with the Arcminute Microkelvin Imager 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 (Margutti et al. 2017).Furthermore, observations done by using the Very Long Baseline Interferometry found that the shock expansion has already strongly decelerated 384 days after the explosion (Bietenholz et al. 2018(Bietenholz et al. , 2021)).
Altogether, the evolution of optical, radio and X-ray emission have been interpreted by considering the interaction of the SN ejecta with a wind with Ṁ = 5 × 10 −5 M yr −1 around the progenitor star, and a massive shell located at R sh = 5 × 10 16 cm, with an extension of approximately 0.25 R sh and a density of ∼10 6 cm −3 (Milisavljevic et al. 2015;Anderson et al. 2017;Margutti et al. 2017;Bietenholz et al. 2018Bietenholz et al. , 2021)).Assuming a velocity of 10-100 km s −1 for the shell, its position implies that it was ejected ∼ 100-1000 yrs before the SN explosion (Milisavljevic et al. 2015;Anderson et al. 2017;Margutti et al. 2017;Bietenholz et al. 2018Bietenholz et al. , 2021)).Late infrared, optical and UV observations confirm that the SN shock is still interacting with the massive shells five years after the explosion (Tinyanont et al. 2019;Sun et al. 2020).
Radio and X-ray emission from the SN shock are typically described by considering a self-similar behaviour for the dynamics (e.g., Chevalier 1982;Chevalier & Liang 1989;Chevalier 1998).In the case of a supernova interacting with a massive shell, we can not assume a selfsimilar behaviour for the flow as the interaction with the shell leads to the formation of a strong reverse shock (see Figure 1).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) 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 genetic algorithms 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 cou- pled to a genetic algorithm (GA).The GA changes the shell density with time by minimizing the difference between the synthetic X-ray emission (computed by postprocessing the results of the hydrodynamical simulations) and observations of this SN presented in Margutti et al. (2017); Brethauer et al. (2020Brethauer et al. ( , 2021)).In this approach, each density ρ(r i ) (being i = 1, . . ., N ) 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 genetic algorithm 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.

Numerical codes and initial conditions
We study the interaction of a 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., González-Casanova et al. 2014;De Colle et al. 2014, 2018a,b).
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 8 cm (the outer edge of the iron core) to ∼ 2 × 10 17 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 8 cm to 10 11 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 16 cm.
In the small scale simulation, we set the density profile of the progenitor star by using the E25 pre-supernova model from Heger et al. (2000).This corresponds to a star which has lost its hydrogen and helium envelope.The resulting Wolf-Rayet star has a mass of 5.25 M and a radius of 3 × 10 10 cm.The computational grid extends radially from 2.2 × 10 8 cm to 6.6 × 10 11 cm.We employ 20 cells at the coarsest level of refinement, with 22 levels of refinement, corresponding to a resolution of 1.6×10 4 cm.The SN energy (E SN ≈ 10 51 erg) is imposed by setting the pressure of the two inner cells of the computational box as p = E SN (Γ ad − 1)/V , being Γ ad = 4/3 the adiabatic index and V the volume of the two cells.Outside the stellar surface, we take ρ = Ṁw /(4πr 2 v w ), being Ṁw = 5 × 10 −6 M yr −1 and v w = 10 8 cm s −1 the mass loss rate and velocity of the wind launched by the Wolf-Rayet star before the collapse.As the velocity of the SN shock front is about two orders of magnitude larger than the Wolf-Rayet wind (i.e., ∼ 10 10 cm s −1 vs. ∼ 10 8 cm s −1 ), 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 to 4.5 × 10 11 cm in 50 seconds.
In the large scale simulations, the computational box goes from 10 10 cm to 5 × 10 17 cm.For r < 4.5 × 10 11 cm we set the density, pressure and velocity by using the values determined from the small scale simulation.For larger radii, we take the density stratification as ρ = Ṁw /(4πr 2 v w ), with Ṁw = 5 × 10 −6 and v w = 10 8 cm s −1 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 9 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.
Trying to reproduce the observed X-ray emission, we have first considered a massive, uniform cold shell located at the 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 M to 2.5 M ), 15 values of ∆R s (ranging from 7.5 × 10 14 up to 2.5 × 10 15 cm), and 11 values of R s (from 1.8 × 10 16 to 7.8 × 10 16 cm).To compare the results of the numerical simulations with observations, we have then employed a ray-tracing code (see Section 2.2).Unfortunately, none of these models give an acceptable fit to the observations, implying that the density of the massive shell is not constant.
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 ten values for the density at 10 differ-ent radii implies running 10 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 other two codes: a radiation transfer code which computes the bremsstrahlung radiation (see Section 2.2), and a genetic algorithm (described in Section 2.3) that automatically and randomly changes the density profile inside the shell by minimazing the variance between the synthetic observations computed from the numerical model and the X-ray observations.

Bremsstrahlung Emission Code
We post-process the results of the numerical simulations by computing synthetic spectra.The specific flux is given by where dΩ = 2π sin θdθ/D 2 , being D = 14.7 Mpc the luminosity distance from the SN 2014C (Freedman et al. 2001), and I ν is the specific intensity.
The observed X-rays 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 where S ν = j ν /α ν is the source function, τ ν = α ν dl is the optical depth, and j ν and α ν are the emissivity and the absorption coefficient respectively.In addition to the bremsstrahlung self-absorption, we also consider photoelectric absorption, which at early times dominates absorption for frequencies 10 18 Hz, so that 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) where Z, n e and n i are the atomic number, 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 solar metallicity and composition for an ideal gas.
Cooling is not included in the numerical simulations as the bremsstrahlung cooling timescale is t ff = 60(T /10 8 K)(n e /10 6 cm −3 ) yrs, so it is much larger than the timescales studied here.-Flowchart showing the genetic algorithm employed to determine the density stratification of the massive shell interacting with the SN 2014C.First, we initialize a population of random densities at different shell radii.We modify the initial population by mutations and cross-over (see the text for details).We run hydrodynamical simulations following the interaction of the SN ejecta with the shell and we compute synthetic X-ray spectra which are compared with the observations.We select the "best" elements of the population using a "fitness function" (the χ 2 test).The process is repeated a thousand times.
commonly used optimization methods4 , 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.
As mentioned before, we made a single small scale simulation for the propagation of the shock-wave in the interior of the star, and we applied the GA 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 shows schematically the GA algorithm implemented.A population formed by 10 elements (the "chromosomes" in the GA terminology) is chosen at the beginning of the iterative process.Each element of the population is defined by setting 11 values of the density at different radii inside the shell (the "genes").Nine of the densities are defined at a radius given by each of the nine observational epochs available.Additionally, we define two densities, one at a smaller and one at a larger radius.At each step, we create 90 new elements of the population.Half of them are defined by randomly choosing two elements of the original population and applying to them cross-over and mutation (see below), while the other half is initialized by directly copying the density values from one random element of the population and modifying it only by mutation.
In the cross-over process, the two "parents" are mixed by choosing randomly a certain number of densities from each element (e.g., the first and second densities from the first element, the third density from the second element 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 element.We do so by setting a Gaussian distribution around the original value of the density ρ 0 , with a width given, in 90% of the cases, by ρ 0 /2, and in 10% of the cases by 10ρ 0 , 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 were then mapped as initial conditions into the HD code.Then, after the simulation was completed, the bremsstrahlung X-ray emission was computed by post-processing the results of the HD calculation.A fitness function (a reduced χ 2 test) was applied to compare the synthetic spectra produced by the model and the observational data at 9 epochs: 308, 396, 477, 606, 857, 1029, 1257, 1574 and 1971 days after explosion.The fittest 10 elements (out of the new 90 and the old 10 elements of the population) were saved and used as initial condition for a new step.This process was repeated for ∼ 100 iterations, for a total number of ∼ 10 4 simulation.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 synchronised the simulations, selected the best elements and managed the cross-over/mutation processes, while the other nodes run (in parallel) each of the 90 hydrodynamics simulation, compute the X-ray spectra and the fitness function (a χ 2 test).

SN shock propagation through the progenitor wind
The propagation of the shock wave through the star has been extensively studied both for the non-relativistic (e.g., Sakurai 1960;Matzner & McKee 1999) and 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 it moves in the stellar envelope, in which the density drops steeper than ρ ∝ r −3 , the shock velocity increases to mildly relativistic speeds (see the blue line in the bottom panel of Figure 3).
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 4 km s −1 , while a small fraction of the mass (corresponding to a kinetic energy ≈ 10 47 ergs) expands with larger velocities (up to v sh ∼ 0.5 c in our simulations).Self-similar solutions describing the propagation of the SN shock through a politropic envelope (with ρ env ∝ (R/r − 1) k ) predict an ejecta density stratification ρ ∝ r −n after the break-out, with n = 7 − 11, depending on the structure of the progenitor star.As we employ a realistic model for the progenitor star structure (which slightly differs from a politrope), we get n ∼ 10 in the inner part of the ejecta, and a transition to a less steep profile as we approach the shock front (with n ∼ 9).
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 (blue lines in Figure 3).As described in section 2.3, 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 3 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 forward shock (FS), which accelerates and heats the WR wind, and the reverse shock (RS), which decelerates and heats the SN ejecta.Following Chevalier (1982), the evolution of the forward and reverse shocks is self-similar, with R ∼ t m , where m = (n − 3)/(n − 2).By taking n ∼ 9, we get m ∼ 8.5, which is consistent with the evolution of the shock wave obtained in our simulation (see Figure 3, bottom panel, and Figure 4).The FS velocity is much larger than the wind velocity.Thus, the postshock temperature achieves values 10 11 K, while the SN bulk temperature quickly drops by adiabatic expansion (heating by Ni 56 and Co 56 decay is not included in the calculation).
In our simulations, the SN shock acceleration stops only when the shock arrives to 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 break-out 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 on particular shock velocity obtained, while it?will depend strongly on the pre-shock structure of the ejecta.The break seen at ∼ 50 days corresponds to the beginning of the interaction with the shell.Bottom panel: shock velocity.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.

SN shock interaction with the massive shell
At ∼ 50 days after the explosion, the shell begins interacting with the massive shell (see Figure 4).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. 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 4 and 5.The shell presents large density fluctuations (see the upper panel of figure 5).Then, the shock propagation leads to the formation of strong reverse shocks which interaction produces the complex shock structure observed in Figure 5. Once interacting with the shell, the shock velocity quickly drops to ∼ 7500 km s −15 , then maintain an approximately constant velocity (see Figure 4).Small fluctuations in the shock velocity are present at ∼ 10 3 days (see the bottom panel of Figure 3), with increases (drops) in velocity by ∼ 3000 km s −1 corresponding to drops (increases) in the density.
Initially, the reverse shock is stronger then the forward shock.Thus, the shocked ejecta is hotter than the shocked wind (Figure 5, blue and green lines in the middle panel).As the ejecta crosses the reverse shock and approaches the reverse shock 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 mostly originated into the shocked ejecta.
A fit to the density profile gives ρ ∝ r −3.00±0.06(although there are large fluctuations), consistently with the constant shock velocity seen in Figure 5 (as E ∼ M v 2 ∼ R 3 ρv 2 implies that the velocity is constant as long as ρ ∝ R −3 and the shock is adiabatic).The shell mass is 2.6 M .This value is consistent with the 3.0 M ± 0.6 M (VLB model) determined by Brethauer et  3, but for larger radii, showing in detail the interaction between the SN ejecta and the shell.The lines corresponds to epoch in which there are X-ray observations available."IC" represents the initial shell density profile.In addition to nine densities corresponding to the epochs with X-ray data available, the amount of unshocked neutral mass is estimated by considering bound-free X-ray absorption.As the absorption does not depend on the density stratification but only on the amount of mass crossed by the X-rays, we show this region as uniform in the figure (dashed line in the top panel).2021) which also assumes spherical symmetry and within the range of typical shell masses observed in type IIn SNe (0.1-10 M , see, e.g., Smith 2017;Branch & Wheeler 2017).

al. (
To determine the structure of the shell, we left as a free parameter in the GA algorithm the amount of neutral hydrogen still to be crossed by the ejecta (inferred from the bound-free absorption).We get M = 0.38 M solar masses for this component at t = 1971 days, which is represented by a constant density dashed line in the top panel of Figure 5.We notice that the exact structure of this region can not be determined.Nevertheless, it will extend to r ∼ 2.3 × 10 17 cm if the shell density continues dropping as r −3 , in which case (moving at a constant speed as discussed above) the SN shock will break out of it ∼ 8.5 years after the explosion.Late X-ray and radio observations will help to constrain the outer structure of the shell.

Radiation
To compare the model with observations, we compute the X-ray bremsstrahlung emission coming from the shocked material.We assume that the shocked material is completely ionized6 (so that it does not contribute to the bound-free opacity) and that the unshocked shell is neutral.Extending the radio synchrotron emission to X-rays by assuming F ν ∝ ν −(p−1)/2 ∝ ν −1 with p ∼ 3, Margutti et al. (2017) showed that the synchrotron emission gives a negligible contribution by X-ray flux (i.e., ∼ 2 orders of magnitude than the observed values), suggesting a thermal origin for the X-rays.
Figure 6 shows the X-ray spectra at 308,396,477,606,857,1029,1257,1571,1971 days and the best fit obtained by employing the GA (Section 2.3).Dashed lines are computed by assuming solar metallicity in the neutral shell.The best fit, in this case, predicts a larger absorption by a factor 5 at 2 × 10 17 Hz than observed.Solid lines, which fits better the observational data, are computed by assuming that the neutral medium has half of solar metallicity.

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), a density stratification ρ ∝ r −3 and extends from 1.6 × 10 16 cm to 1.87 × 10 17 cm.While uncommon in type Ib SNe, SN IIn show evidence of strong interaction, with shell masses of 0.1 − 10 M (see, e.g.Smith 2017;Branch & Wheeler 2017, and references therein).Also, many type IIn SNe show an X-ray emission inconsistent with a density profile ρ ∝ r −2 , implying a steeper density stratification (Dwarkadas & Gruszko 2012).While sharing many characteristics with SN IIn, in the case of the SN 2014C the shell is located at a much larger distance, implying that it was ejected ∼ 60/(v w /100 km s −1 ) yrs before the SN explosion.
Harris & Nugent (2020) presented a series of numerical simulations of a SN ejecta interacting with a wall located at r 0 ∼ 1.6 × 10 16 cm, and with a CSM wind at larger radii.They deduce a small mass for the wall (∼ 0.04 − 0.31 M ), in apparent contradiction with previous estimations (e.g., Margutti et al. 2017).Actually, the density of the CSM wind is assumed to be a factor 4-7 smaller than the density of the wall at r 0 .The CSM wind, then, corresponds to a mass loss Ṁw ∼ 10 −2 M yr −1 (assuming v w = 10 8 cm s −1 ).Our results (which, we stress, have been obtained by a "blind" fit, i.e. without assumptions on the final shell structure) clarify this inconsistency, showing that the wall and the CSM wind are both part of the same extended structure (see figure 5), with the densest shell material located at ∼ 2 × 10 16 (see Figure 5) in agreement with Harris & Nugent 2020, 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 the SN 2014C are remarkably similar to those inferred for the SN 2001em.The X-ray, radio and Hα emission from SN 2001em have 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 16 We do not include the shaded grey area in the fit, as it is dominated by Fe emission lines.Curves at different times are rescaled so they do not overlap with each other.cm, and expanding with a velocity of 5800 ± 10 4 km s −1 (Bietenholz & Bartel 2007), which is consistent with the 7500 km s −1 inferred here.
The best fit to the data is achieved by considering a bound-free cross-section, corresponding to half solar metallicity.This low 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).
An equivalent effect is obtained considering that the X-ray emission and/or UV radiation coming from the shocked shell and SN ejecta can partially ionize the unshocked material, producing an effect similar to a drop in the metallicity, as the optical depth is ∝ n H 0 σ ν , being n H 0 the density of neutral hydrogen and σ ν the boundfree cross-section (which depends on metallicity).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 as the bremsstrahlung emission is ∝ n 2 e .A detailed calculation of the ionization of the shell is left for a future work.
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 −2 M yr −1 (v w /100 km s −1 ).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 ejections can be achieved only by detailed theoretical models coupled to a larger sample of observed interacting supernovae.
The density fluctuations have a periodicity of ∼ 4 × 10 16 cm.A similar periodicity has been observed in the radio emission from SN 1979C (Weiler et al. 1991) and have 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 16 /v w ∼ 10 yrs (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).Small scales inhomogeneities are expected.The density profile shown in Figure 5 show large scale density fluctuations.Furthermore, it is likely that the medium is, at some scale, clumpy.If the shell is not perfectly homogeneous before interacting with the SN ejecta, the interaction will amplify the inhomogeneities, leading to a multi-phase medium with denser/colder regions in pressure equilibrium with more tenuous/hotter regions, which is consistent with the strong Fe emission line observed at ∼ 6.5 keV (Margutti et al. 2017).Also if the shell is initially nearly perfectly homogeneous, the contact discontinuity which 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 dishomogeneities in the post-shock region which can not be captured by our numerical simulations (but see Harris & Nugent 2020 for an approximated treatment of RT instabilities in one dimensional simulations).As the bremsstrahlung emission is ∝ Z 2 n 2 , 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.

CONCLUSIONS
In this paper, we presented hydrodynamical simulations of the strongly interacting SN 2014C.First, we follow the propagation of a 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.3, we initialize the shell with a uniform density n shell = 10 7 cm −3 .We follow 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 compute the bremsstrahlung emission using the algorithm described in Section 2.2, and compare 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 and metallicity.In particular, we get a mass of 2.6 M for the shell and a density profile ρ ∝ r −3 .We also found that the shell is very extended, with a size 10 17 cm.If the shell stratification continues with the same slope, the SN shock will break out of it nearly 8 yrs after the explosion, i.e. during 2022.
Radio and X-ray emission allows us to understand the mass loss history of core-collapse SNe progenitor on timescales which 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 envi-ronment once data at several epochs are available, as in the case of SN 2014C.The X-ray emission tracks the forward and reverse shock emission, depending on the density of the environment and the ejecta velocity.The Hα emission tracks the shocked shell and the unshocked medium fotoionized by the X-ray and UV radiation.All together, 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.

Fig. 1 .
Fig. 1.-Schematic evolution of a supernova shock interacting with an external shell (not-to-scale).(A) the supernova explosion produces a outwardly propagating shock; (B) the reverse shock becomes stronger as the shock front collides with the external shell, moving back towards (in mass coordinates) the supernova center.
Fig.2.-Flowchart showing the genetic algorithm employed to determine the density stratification of the massive shell interacting with the SN 2014C.First, we initialize a population of random densities at different shell radii.We modify the initial population by mutations and cross-over (see the text for details).We run hydrodynamical simulations following the interaction of the SN ejecta with the shell and we compute synthetic X-ray spectra which are compared with the observations.We select the "best" elements of the population using a "fitness function" (the χ 2 test).The process is repeated a thousand times.

Fig. 3 .
Fig.3.-Timeevolution of the SN shock as it moves through the wind of the progenitor star.From top to bottom: density, temperature and velocity profiles of the SN shock as it interacts with the wind of the progenitor star (at r 1.6 × 10 16 cm) and the outer shell (at 1.6 × 10 16 cm).
Fig. 4.-Top panel: position of the SN shock as function of time.The break seen at ∼ 50 days corresponds to the beginning of the interaction with the shell.Bottom panel: shock velocity.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.
Fig.5.-The same as Figure3, but for larger radii, showing in detail the interaction between the SN ejecta and the shell.The lines corresponds to epoch in which there are X-ray observations available."IC" represents the initial shell density profile.In addition to nine densities corresponding to the epochs with X-ray data available, the amount of unshocked neutral mass is estimated by considering bound-free X-ray absorption.As the absorption does not depend on the density stratification but only on the amount of mass crossed by the X-rays, we show this region as uniform in the figure (dashed line in the top panel).
Fig.6.-Comparison between X-ray observations and synthetic observations computed for the best model.Full lines represents solar metallicity while dashed lines are for half solar metallicity.We do not include the shaded grey area in the fit, as it is dominated by Fe emission lines.Curves at different times are rescaled so they do not overlap with each other.