Learning to exploit z-Spatial Diversity for Coherent Nonlinear Optical
Fiber Communication

Sebastian Jung, Tim Uhlemann, Alexander Span, Maximilian Bauhofer and Stephan ten Brink
Institute of Telecommunications

University of Stuttgart

Pfaffenwaldring 47, 70569 Stuttgart, Germany

Higher-order solitons inherently possess a spatial periodicity along the propagation axis.
The pulse expands and compresses in both, frequency and time domain.
This property is exploited for a bandwidth-limited receiver by sampling the optical signal at two different distances.
Numerical simulations show that when pure solions are transmitted and the second (i.e., further propagated) signal is also processed, a significant gain in terms of required receiver bandwidth is obtained.
Since all pulses propagating in a nonlinear optical fiber exhibit solitonic behavior given sufficient input power and propagation distance, the above concept can also be applied to spectrally efficient Nyquist pulse shaping and higher symbol rates.
Transmitter and receiver are trainable structures as part of an autoencoder, aiming to learn a suitable predistortion and post-equalization using both signals to increase the spectral efficiency.

autoencoder, communication, optical, coherent, nonlinear, chromatic dispersion, soliton

§ INTRODUCTION

The ever-increasing demand for data means that the resources available in today's communications systems must either be used more efficiently or expanded.
Optical fibers in general provide very high bandwidths. For instance, the C-band, commonly used for communication, covers a total bandwidth of around \(4.4~\mathrm{THz}\).
However, the exploitation of this spectrum, or, at least the expansion to higher bandwidths is limited by the required signal processing, e.g., the bandwidths of digital-to-analog (DAC) and analog-to-digital (ADC) conversion, respectively.
These components typically support bandwidths of $100\,\mathrm{GHz}$ to $200\,\mathrm{GHz}$ maximum, which is only a small fraction of the available channel spectrum.
Furthermore, as the optical fiber is a nonlinear medium, the signal broadens as it propagates, occupying additional spectral resources. Summing up, bandwidth limitation is the main driver for the performance degradation, or more positively, one of the highest potentials for higher-power optical communications systems.

One approach to face this bottleneck, is to split the C-band into multiple channels, i.e., in a wavelength-division multiplexing (WDM) setup.
This way, the limited bandwidth of DACs and ADCs is sufficient to cover the entire C-band.
Typically, around $100$ different channels per polarization [1] are used in parallel.
However, there is still a bandwidth limitation, or rather a difference between channel spacing and bandwidth, to reduce the interference of neighboring channels.

In this work, the idea is to exploit additional and available bandwidth of the optical fiber, while still relying on further limited converters. The approach is to take advantage of the propagation properties of signals on the optical fiber, i.e., the solitonic propagation of any waveform given a sufficiently high launch power and distance. Note that the two concepts do not exclude one another.

Higher-order solitons posses a spatial periodicity in bandwidth that can be exploited by sampling the signal at two different points in propagation distance, i.e., the $z$-axis. Thus, information that may be outside the converters' bandwidth at one point may have re-entered the available spectrum at another point in $z$-direction.
This principle is shown in Fig. <ref>.

Intensity of higher-order solitons: Motivation for using \(z\)-spatial diversity

Since there is no known analytical solution to combine the two $z$-spatially different signals obtained in this way, we apply and optimize a neural network (NN) to perform this task.
More precisely, a trainable Autoencoder [2] is used to (i) pre-distort the single transmit signal and (ii) combine and post-equalize the received signals.

It is already well known that Digital-Backpropagation (DBP) is able to compensate for all channel effects if no bandwidth limitation is present and signal-noise interaction can be neglected [3, 4].
Therefore, in this work, we use DBP as a baseline algorithm that does not have the ability to exploit $z$-spatial diversity.
In a previous work, we have already shown that the performance degradation due to the bandwidth limitation of DBP can be partially circumvented by an appropriate Autoencoder-based pre-processing of the transmit signal [5].
Therefore, we use the results of that work to design an Autoencoder that not only compensates for, but also exploits the spectral propagation characteristics of optical signals.
By having designed and trained the spoken Autoencoder, we achieved a gain of over \(+0.4\,\frac{\mathrm{bit}}{\mathrm{s}\cdot\mathrm{Hz}}\) in terms of spectral efficiency.

The remainder of the paper is organized as follows.
Sec. <ref> provides the used channel model and performance metric as well as the basic idea of the concept.
The whole communication setup is provided in Sec. <ref> and the obtained results are presented in Sec. <ref>.
Finally, in Sec. <ref> a summary of the work and an outlook are given.

§ BASICS

Propagation in a single mode fiber is governed by the Nonlinear Schroedinger Equation (NLSE).
Since we want to  exploit the  solitonic  properties  of  propagating  impulses, we assume ideal distributed Raman amplification (IDRA).
Therefore, the attenuation term is dropped and the amplified spontaneous emission (ASE) noise is added instead.

\begin{align}
\frac{\partial q\left(t,z\right)}{\partial z} = \mathrm{j}\frac{\beta_2}{2}\frac{\partial^2 q\left(t,z\right)}{\partial t^2}-\mathrm{j}\gamma\left\vert q\left(t,z\right)\right\vert^2 q\left(t,z\right) + n\left(t,z\right)
\end{align}

Here \(q\left(t,z\right)\) is the optical baseband signal with time \(t\) and distance \(z\) and \(\beta_2\) and \(\gamma\) are the dispersion and nonlinearity parameters, respectively.
This differential equation is typically solved numerically using the Split-Step Fourier Method (SSFM).
The used implementation can be found in [6].
A standard-single mode fiber (S-SMF) is used as the communication medium in the C-band around \(f_\text{0}=193.55\,\mathrm{THz}\).

The power of the noise term \(n\left(t,z\right)\) is calculated from [7]

\begin{align}
\sigma_{\mathrm{n},\mathrm{ase}}^2\left(\Delta z\right) = n_\text{sp}\left(f_0\right)hf_0\alpha\Delta zB
\end{align}

with spontaneous emission factor \(n_\text{sp}\), Planck's constant \(h\), attenuation \(\alpha\) and bandwidth \(B\).

The transmitters and receivers used in this work are bandwidth-limited due to the DAC and ADC.
Fig. <ref> shows the channel model including the bandwidth limitation.

System model

Since the Kerr-effect  causes  spectral  broadening, information is lost at the receiver.

This bandwidth limitation is also the reason why conventional systems degrade in performance when the Kerr-effect becomes the dominant degradation, i.e., at higher input power.

In our setup, the DAC and ADC have a passband-bandwidth of \(B=20\,\mathrm{GHz}\) and are assumed to be ideal filters with a brickwall spectrum.

The simulations are performed at a frequency of \(f_\text{sim} = 1\,\mathrm{THz}\) with a symbol rate of \(R_\mathrm{S}= 20\,\mathrm{GBd}\).
This results in a resolution of \(f_\text{sim}/R_\mathrm{S}=50\) samples per symbol.

§.§ Performance metric

In contrast to our previous work, we use soft decision spectral efficiency \(\eta\), which is an achievable information rate obtained by evaluating the output of a Gaussian mixture model (GMM)-based demapper with its corresponding likelihoods.
This metric is upper bounded by Shannon's limit [8].

§.§ Spatial diversity

$z$-spatial diversity describes the approach to use the spatial periodicity of higher order solitons to circumvent the bandwidth limitation of DAC and ADC.
Hereby, the signal is sampled at two different distances along the propagation axis $z$.

The signal is sent over an optical fiber that spans the communication distance of \(\ell_1=1000\,\mathrm{km}\).
Afterwards, a \(3\,\mathrm{dB}\)-coupler splits the signal.
One part is passed directly to the digital signal processing (DSP) after the ADC, the other part is optionally amplified using an Erbium-doped fiber amplifier (EDFA), to compensate for the splitter loss, and sent through a second fiber of length \(\ell_2\), ADC, and finally passed to the DSP.
The entire structure is shown in Fig. <ref>.

Realization of spatial diversity

Two hyperparameters remain to be optimized: the fiber length \(\ell_2\) and the use of the EDFA.
These will be optimized later on.

§.§ Spatial diversity with solitons

A proof-of-concept of the $z$-spatial diversity setup was done.
There, second-order solitons were sent through the spatial diversity setup at $5\,\mathrm{dBm}$, with the eigenvalues $\left\{\lambda_1=0.5\mathrm{j}, \lambda_2=\mathrm{j}\right\}$ and a periodic length of $z_\mathrm{p}=483\,\mathrm{km}$.
The first sampling of the signal is performed at \(\ell_1=1000\,\mathrm{km}\).
Compensation of the splitter loss using an EDFA was mandatory since solitons were sent.
In addition, ideal Raman amplification is assumed.

Simulation results for combining both paths using a neural network (NN) are shown in Fig. <ref>.

Simulation results using pure solitons

There, the mutual information versus the length of the second fiber \(\ell_2\) is plotted.
A clear optimum is visible.
It occurs at half of the periodic length, i.e., at the length where the impulse is at the other extreme of its periodic cycle.

§ SETUP

Fig. <ref> shows the TX (top, green box) and RX (bottom, blue box) used to build the Autoencoder (AE) neural network structure.
Tab. <ref> provides the simulation parameters used.

Building blocks of transmitter (above) and receiver (below)

The TX maps the symbol indices onto a geometrically optimized circular constellation as shown in Fig. <ref>.

Used circular constellation

These symbols are then predistorted by a neural network before being upsampled to the simulation frequency.
Pulse shaping is applied and the average power of the signal is normalized.

At the receiver, a chromatic dispersion compensation (CDC) is performed and the signal is downsampled to the symbol rate.
Since the further propagating signal has undergone a different amount of chromatic dispersion, its compensation is performed separately on both paths.
This is indicated by the two parallel paths in the receiver block diagram.
On a symbol-basis, the combined processing of both paths is performed using an NN.
The output is passed on to the trainable demapper, which consists of three fully-connected layers.
Its structure was chosen to give the AE more degrees of freedom during training.
This trainable demapper is used only for training, since a GMM-based and optimized demapper is used for the evaluation.

The architecture for predistortion and combination/compensation is similar to [9] and adapted to the problem at hand.

The pulse shaper and the low-pass filter before downsampling are assumed to be ideal sinc-functions, i.e., with a brickwall spectrum of bandwidth \(B=20\,\mathrm{GHz}\).

System parameters

Property   Symbol   Value

Attenuation                        $\alpha$           \(0.2~\mathrm{dB}\,\mathrm{km}^{-1}\) (IDRA)

Chromatic dispersion               $\beta_2$          \(-21.67~\mathrm{ps}^2\,\mathrm{km}^{-1}\)

Kerr-effect                        $\gamma$           \(1.27~\mathrm{W}^{-1}\,\mathrm{km}^{-1}\)

Carrier frequency                  \(f_0\)           \(193.55~\mathrm{THz}\)

2*Fiber length      \(\ell_1\)        \(1000~\mathrm{km}\)

\(\ell_2\)   to be optimized

Spontaneous emission factor (Raman)   \(n_{\text{sp},\text{Raman}}\)   \(1.0\)

Spontaneous emission factor (EDFA)   \(n_{\text{sp},\text{EDFA}}\)   \(3.16\)

Simulation frequency               \(f_\text{sim}\)   \(1~\mathrm{THz}\)

Symbol rate                        \(R_\mathrm{S}\)   \(20~\mathrm{GBd}\)

Bandwidth DAC/ADC                  \(B\)   \(20~\mathrm{GHz}\)

Launch power                       \(P_\text{avg}\)   \(-5~\mathrm{dBm}\dots5~\mathrm{dBm}\)

§ RESULTS

Spectral efficiency gain over length of second fiber \(\ell_2\) (SD: spatial diversity w.o. EDFA, SDA: spatial diversity w. EDFA)

To optimize the two hyperparameters mentioned before, the length of the second fiber \(\ell_2\) and the use of an EDFA to compensate for the coupler loss, we tested our setup with different combinations of both.

For the length we considered a range of up to \(500\,\mathrm{km}\) with a spacing of \(20\,\mathrm{km}\) up to \(320\,\mathrm{km}\) plus the lengths \(400\,\mathrm{km}\) and \(500\,\mathrm{km}\).
The denser grid up to \(320\,\mathrm{km}\) was chosen because a second-order soliton propagating with a power of \(5\,\mathrm{dBm}\) has a periodic length of \(483\,\mathrm{km}\) and previous simulations showed a maximum gain at exactly half that length.
Multiplying this by a factor of \(1.5\) should ensure that the maximum is within this range while \(400\,\mathrm{km}\) and \(500\,\mathrm{km}\) are used as trend indicators for further behavior.

Three different power levels were tested for their ideal length: \(0\,\mathrm{dBm}\), \(2.5\,\mathrm{dBm}\) and \(5\,\mathrm{dBm}\).
Fig. <ref> shows the gain in terms of spectral efficiency \(\Delta\eta\) compared to a baseline with exactly the same number of parameters in the neural networks.
This baseline does not exploit spatial diversity but obtains the same information for \(\mathbf{y}_1\) and \(\mathbf{y}_2\).
By comparing the systems with spatial diversity to this baseline, we can see the gain achieved purely by exploiting $z$-spatial diversity .
In the plots, SD (<ref>) and SDA (<ref>) refer to the case with and without splitter loss compensation using an EDFA, respectively.
At first glance, it becomes apparent that SDA outperforms SD at all power levels by up to \(\Delta\eta=+0.15\,\frac{\mathrm{bit}}{\mathrm{s}\cdot\mathrm{Hz}}\).

Both systems, SD and SDA, achieve significant gains of up to \(+0.5\,\frac{\mathrm{bit}}{\mathrm{s}\cdot\mathrm{Hz}}\) over the baseline.
It seems that any small addition of propagation distance in the second fiber produces a large initial gain.
Interestingly, if you compare the curves for SD at \(2.5\,\mathrm{dBm}\) and \(5\,\mathrm{dBm}\), they are almost perfectly aligned and form a sort of plateau after reaching their maximum.
This plateau can also be seen at \(0\,\mathrm{dBm}\), but less pronounced.

Another notable aspect is how the ideal length changes with the input power which is best seen with SDA.
Starting from about \(320\,\mathrm{km}\) at \(0\,\mathrm{dBm}\), the maximum is shifted to a distance of only \(100\,\mathrm{km}\) at \(5\,\mathrm{dBm}\).
This is consistent with soliton theory, where higher power solitons have shorter periodic lengths.

For the \(5\,\mathrm{dBm}\) curve, we observe that the gain achieved with EDFA is partially lower than without EDFA.
There are two possible explanations for this:
Either the maximum gain has already been achieved, indicating that the EDFA is no longer needed at higher input powers.
Thus, the split signal itself still has enough nonlinearity to achieve the same gain.
Another possible explanation is that the optimization did not converge to a good optimum during training.
This theory can be somewhat discounted since we have seen similar maximum gain curves with spatial diversity in simpler setups, i.e., with fewer parameters and simpler constellations.
In addition, the gain for SD and SDA are quite similar up to \(320\,\mathrm{km}\) supporting the first theory.

In summary, in order to get the maximum performance from exploiting $z$-spatial diversity, the coupler loss should be compensated and the length should be optimized.
Although the performance at a non-ideal length is still superior to the case without EDFA.
For all SDA-curves the gain of the SD case can be achieved with the minimum tested length of \(20\,\mathrm{km}\), which is promising and underlines that most of the gain is already achieved with only a small additional propagation distance with the second fiber.

This allows us to compare spatial diversity systems with conventional systems using CDC and DBP.
We use a split DBP [10] where the predistortion for half of the communication distance \(\ell_1\) is applied at the transmitter and the other half is applied at the receiver.
CDC is applied only at the receiver, as in the AE setup, which therefore should provide a lower performance limit.

Fig. <ref> shows the performance in terms of spectral efficiency versus the average launch power for different systems.

Comparison of spectral efficiency curves versus launch power for different system configurations (AEC: only compensation, AEP: only predistortion, AEPC: predistortion and compensation, SD: spatial diversity w.o. EDFA, SDA: spatial diversity w. EDFA)

We only consider the interval from \(-5\,\mathrm{dBm}\) to \(5\,\mathrm{dBm}\), since the (split) DBP already reaches the maximum achievable performance up to about \(-3\,\mathrm{dBm}\).
The plot shows the performance of three different AE systems: predistortion only at the transmitter (AEP, <ref>), compensation only  at the receiver (AEC, <ref>) and the full system with both predistortion and compensation (AEPC, <ref>).
Spatial diversity was enabled where applicable, i.e., in systems with compensation.

The comparison of the spatial diversity systems (SD/SDA) with their respective baselines shows that by exploiting the solitonic propagation properties, significant gains can be obtained over the whole power range if we are not limited by saturation.
The gain achieved for AEPC-SDA (<ref>) compared to the baseline (<ref>) corresponds to a shift of about \(1\,\mathrm{dBm}\) in terms of average launch power.
Accordingly, the split DBP can already be beaten at \(2.5\,\mathrm{dBm}\) and \(5\,\mathrm{dBm}\) for the DBP operating at \(40\,\mathrm{GHz}\) and \(f_\text{sim}\), respectively.

§ SUMMARY AND CONCLUSION

In this paper, the use of $z$-spatial diversity for state-of-the-art communication systems has been investigated.
We have shown that significant gains in terms of spectral efficiency (up to almost \(+0.5\,\frac{\mathrm{bit}}{\mathrm{s}\cdot\mathrm{Hz}}\)) can be achieved by combining the two signals.
Simulations show that even higher gains can be achieved by compensating the coupler loss with an EDFA.

Currently, the main issue with this approach is the required length of the additional fiber.
It has already been shown that a length of \(20\,\mathrm{km}\) is sufficient to achieve gains over the baseline.
Fibers designed specifically for this use-case, e.g. highly nonlinear fibers (HNLF) could improve performance for this or even shorter lengths, making the system more practical.
Combining this approach with WDM may also be interesting, since such systems require hard bandwidth limitations and DBP is suboptimal there, leaving more potential for improvement.

§ ACKNOWLEDGMENT

The authors would like to thank D. J. Faul for his preliminary investigation of $z$-spatial diversity with pure solitons.

[1] R. Dar, P. J. Winzer: ”On the Limits of Digital Back-Propagation in Fully Loaded WDM Systems” , IEEE Photonics Technology Letters 28(11), 1253–1256, 2016
[2] T. O’Shea and J. Hoydis, "An Introduction to Deep Learning for the Physical Layer," in IEEE Transactions on Cognitive Communications and Networking, vol. 3, no. 4, pp. 563-575, Dec. 2017
[3] C. Paré, A. Villeneuve, P.-A. Bélanger, and N. J. Doran, "Compensating for dispersion and the nonlinear Kerr effect without phase conjugation,"  Opt. Lett. 21, 459-461, 1996
[4] R.-J. Essiambre and P. J. Winzer, "Fibre Nonlinearities in Electronically Pre-Distorted Transmission", 31st European Conference on Optical Communication, Glasgow, 2005
[5] T. Uhlemann, A. Span, S. Dörner, an S. ten Brink: ”Introducing $\gamma$-lifting for Learning Nonlinear Pulse Shaping in Coherent Optical Communication”, Photonic Networks, 23th ITG-Symposium, 2022
[6] J. Hoydis, S. Cammerer, F. Ait Aoudia, A. Vem, N. Binder, G. Marcus, and A. Keller: ”Sionna: An Open-Source Library for Next-Generation Physical Layer Research”, arXiv preprint , Mar. 2022
[7] N. J. Muga, M. C. Fugihara, M. F. S. Ferreira, and A. N. Pinto: ”ASE Noise Simulation in Raman Amplification Systems”, 4, 2009
[8] G. Kramer, M. I. Yousefi, and F. R. Kschischang: ”Upper bound on the capacity of a cascade of nonlinear and noisy channels”, 2015 IEEE Information Theory Workshop (ITW) , 1–4, 2015
[9] P. J. Freire, V. Neskornuik, A. Napoli, B. Spinnler, N. Costa, G. Khanna, E. Riccardi, J. E. Prilepsky, and S. K. Turitsyn: ”Complex-Valued Neural Network Design for Mitigation of Signal Distortions in Optical Links” , Journal of Lightwave Technology 39(6), 1696–1705, 2021
[10] D. Lavery, D. Ives, G. Liga, A. Alvarado, S. J. Savory, and
P. Bayvel, “The Benefit of Split Nonlinearity Compensation for
Single-Channel Optical Fiber Communications,” IEEE Photonics
Technology Letters, vol. 28, no. 17, pp. 1803–1806, Sep. 2016


# Survival of the fittest: Numerical modeling of Supernova 2014c

Felipe Vargas11affiliation: Instituto de Ciencias Nucleares, Universidad
Nacional Autónoma de México, A. P. 70-543 04510 D. F. Mexico 22affiliation:
Universidade Federal do Rio Grande, Av. Italia km 8 Bairro Carreiros, Rio
Grande, RS, Brazil , Fabio De Colle11affiliation: Instituto de Ciencias
Nucleares, Universidad Nacional Autónoma de México, A. P. 70-543 04510 D. F.
Mexico , Daniel Brethauer33affiliation: Center for Interdisciplinary
Exploration and Research in Astrophysics (CIERA) and Department of Physics and
Astronomy, Northwestern University, Evanston, IL 60208 , Raffaella
Margutti33affiliation: Center for Interdisciplinary Exploration and Research
in Astrophysics (CIERA) and Department of Physics and Astronomy, Northwestern
University, Evanston, IL 60208 & Cristian G. Bernal22affiliation: Universidade
Federal do Rio Grande, Av. Italia km 8 Bairro Carreiros, Rio Grande, RS,
Brazil

###### Abstract

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⊙, 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) 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.

###### Subject headings:

supernovae: individual: SN2014C - stars: mass loss - X-rays: individual: SN
2014C - circumstellar matter - methods: numerical

##  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, specially
during the final phases before the supernova (SN) explosion (e.g., Smith,
2014). The mass loss history $\lesssim 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 $\dot{M}_{w}\sim 10^{-4}-10^{-6}$ M⊙ yr-1. The wind
is often inhomogenous, as proven by radio emission, showing small flux
fluctuations of $\sim$ 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 $\sim$ two years before the SN.
SN 2001em, initially classified as type Ib SN, presented prominent H$\alpha$
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 ($\sim 3$ M⊙) hydrogen-rich shell located at $\sim 7\times 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 $\sim 100$ days after
the explosion, showing strong H$\alpha$ emission (Milisavljevic et al., 2015).
The modelling of the optical/UV light curve shows that SN 2014C has a kinetic
energy of 1.75 $\pm$ 0.25 $\times 10^{51}$ erg, with an ejecta mass of 1.7
$\pm$ 0.2 M☉ and a Nickel mass of 0.15 $\pm$ 0.02 M☉ (Margutti et al., 2017).

SN 2014C has also been extensivel observed with X-rays 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., 2020, 2021).
Most of the detected X-ray emission is concentrated in the 1-40 keV energy
band, while emission below $\sim 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\times 10^{39}$ erg s-1 to $5\times 10^{40}$ erg
s-1. Then, it peaked $\sim 850-1000$ days and maintained a nearly constant
flux until $\sim 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 $\sim$ 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, 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
$\dot{M}=5\times 10^{-5}$ M☉ yr-1 around the progenitor star, and a massive
shell located at R${}_{\rm sh}=5\times 10^{16}$ cm, with an extension of
approximately 0.25 Rsh and a density of $\sim$106 cm-3 (Milisavljevic et al.,
2015; Anderson et al., 2017; Margutti et al., 2017; Bietenholz et al., 2018,
2021). Assuming a velocity of 10-100 km s-1 for the shell, its position
implies that it was ejected $\sim 100$-$1000$ yrs before the SN explosion
(Milisavljevic et al., 2015; Anderson et al., 2017; Margutti et al., 2017;
Bietenholz et al., 2018, 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 self-similar 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 coupled 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 post-processing
the results of the hydrodynamical simulations) and observations of this SN
presented in Margutti et al. (2017); Brethauer et al. (2020, 2021). In this
approach, each density $\rho(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.

Figure 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.

##  2\. Methods

###  2.1. 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, 2018b).

We follow the propagation of the SN shock front as it moves through a
computational grid covering $\sim$ nine orders of magnitude in space, from
$\sim 2\times 10^{8}$ cm (the outer edge of the iron core) to $\sim 2\times
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 $\sim 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 $\gtrsim 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\times
10^{10}$ cm. The computational grid extends radially from $2.2\times 10^{8}$
cm to $6.6\times 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\times 10^{4}$ cm. The SN energy ($E_{\rm SN}\approx 10^{51}$ erg) is
imposed by setting the pressure of the two inner cells of the computational
box as $p=E_{\rm SN}(\Gamma_{\rm ad}-1)/V$, being $\Gamma_{\rm ad}=4/3$ the
adiabatic index and $V$ the volume of the two cells. Outside the stellar
surface, we take $\rho=\dot{M}_{w}/(4\pi r^{2}v_{w})$, being
$\dot{M}_{w}=5\times 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., $\sim 10^{10}$ cm s-1 vs.
$\sim 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\times 10^{11}$ cm in 50 seconds.

In the large scale simulations, the computational box goes from $10^{10}$ cm
to $5\times 10^{17}$ cm. For $r<4.5\times 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
$\rho=\dot{M}_{w}/(4\pi r^{2}v_{w})$, with $\dot{M}_{w}=5\times 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\times 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 $\Delta 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_{\odot}$ to 2.5 $M_{\odot}$), 15 values of $\Delta R_{s}$
(ranging from $7.5\times 10^{14}$ up to $2.5\times 10^{15}$ cm), and 11 values
of $R_{s}$ (from $1.8\times 10^{16}$ to $7.8\times 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 different 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.

###  2.2. Bremsstrahlung Emission Code

We post-process the results of the numerical simulations by computing
synthetic spectra. The specific flux is given by

$F_{\nu}=\int I_{\nu}\cos\theta d\Omega\;,$ (1)

where $d\Omega=2\pi\sin\theta d\theta/D^{2}$, being $D=14.7$ Mpc the
luminosity distance from the SN 2014C (Freedman et al., 2001), and $I_{\nu}$
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

$\frac{dI_{\nu}}{d\tau_{\nu}}=S_{\nu}-I_{\nu}\;,$ (2)

where $S_{\nu}=j_{\nu}/\alpha_{\nu}$ is the source function,
$\tau_{\nu}=\int\alpha_{\nu}dl$ is the optical depth, and $j_{\nu}$ and
$\alpha_{\nu}$ are the emissivity and the absorption coefficient respectively.
In addition to the bremsstrahlung self-absorption, we also consider photo-
electric absorption, which at early times dominates absorption for frequencies
$\lesssim 10^{18}$ Hz, so that

$\alpha_{\nu}=\alpha_{\nu,{\rm ff}}+\alpha_{\nu,{\rm bf}}\;,\qquad
j_{\nu}=j_{\nu,{\rm ff}}\;.$ (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)

$\displaystyle j_{\nu}=\frac{6.8}{4\pi}\times
10^{-38}Z^{2}n_{e}n_{i}T^{-1/2}e^{-\frac{h\nu}{kT}}G$ (4)
$\displaystyle\alpha_{\nu}=3.7\times
10^{8}T^{-1/2}\frac{Z^{2}n_{e}n_{i}}{\nu^{3}}\left(1-e^{-\frac{h\nu}{kT}}\right)G\;,$
(5)

where $Z$, $n_{e}$ and $n_{i}$ are the atomic number, electron and ion
densities respectively, and $G$ is the Gaunt Factor ($\sim 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_{\rm ff}=60(T/10^{8}{\rm\,K})(n_{e}/10^{6}{{\rm
cm}^{-3}})$ yrs, so it is much larger than the timescales studied here.

###  2.3. Genetic Algorithm

Figure 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 $\chi^{2}$ test). The process is
repeated a thousand times.

Genetic algorithms (GA hereafter, e.g., Rajpaul, 2012) are based on the theory
of natural selection. GA are commonly used optimization methods111We have
employed GA in this paper but the results obtained are independent on the
particular optimization method chosen., 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 $\rho_{0}$, with a width given, in 90% of the
cases, by $\rho_{0}/2$, and in 10% of the cases by $10\rho_{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 $\chi^{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 $\sim 100$ iterations, for a total number of $\sim$
104 simulation. Each simulation took $\sim 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 $\chi^{2}$ test).

##  3\. Results

###  3.1. 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 $\rho\propto 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 $\sim 25$ s, most of the mass
(and energy) of the SN moves at velocities $\sim 10^{4}$ km s-1, while a small
fraction of the mass (corresponding to a kinetic energy $\approx 10^{47}$
ergs) expands with larger velocities (up to $v_{\rm sh}\sim 0.5$ c in our
simulations). Self-similar solutions describing the propagation of the SN
shock through a politropic envelope (with $\rho_{\rm env}\propto(R/r-1)^{k}$)
predict an ejecta density stratification $\rho\propto 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\sim 10$ in the inner part of the ejecta,
and a transition to a less steep profile as we approach the shock front (with
$n\sim 9$).

Figure 3.— Time evolution 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\lesssim 1.6\times 10^{16}$ cm) and the outer shell (at $\gtrsim
1.6\times 10^{16}$ cm).

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\sim t^{m}$, where
$m=(n-3)/(n-2)$. By taking $n\sim 9$, we get $m\sim 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 post-shock temperature achieves values $\gtrsim 10^{11}$
K, while the SN bulk temperature quickly drops by adiabatic expansion (heating
by Ni56 and Co56 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.

###  3.2. SN shock interaction with the massive shell

At $\sim 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 $\sim 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 $\sim 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 $\sim
7500$ km s-1222From the Hydrogen line, a velocity of $\sim$ half of this value
is inferred (Milisavljevic et al., 2015). This is consistent with the hydrogen
lines being produced in the post shock clumpy medium or by recombination of
the upstream medium., then maintain an approximately constant velocity (see
Figure 4). Small fluctuations in the shock velocity are present at $\sim
10^{3}$ days (see the bottom panel of Figure 3), with increases (drops) in
velocity by $\sim 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.

Figure 4.— _Top panel_ : position of the SN shock as function of time. The
break seen at $\sim 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 $\sim 0.35$c to $\sim 0.15$c. After
interacting with the shell, it drops to $\sim 0.025$c.

A fit to the density profile gives $\rho\propto r^{-3.00\pm 0.06}$ (although
there are large fluctuations), consistently with the constant shock velocity
seen in Figure 5 (as $E\sim Mv^{2}\sim R^{3}\rho v^{2}$ implies that the
velocity is constant as long as $\rho\propto R^{-3}$ and the shock is
adiabatic). The shell mass is 2.6 M⊙. This value is consistent with the 3.0 M
$\pm$ 0.6 M⊙ (VLB model) determined by Brethauer et al. (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).

Figure 5.— The same as Figure 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).

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\sim 2.3\times 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 $\sim 8.5$ years after the explosion. Late X-ray
and radio observations will help to constrain the outer structure of the
shell.

###  3.3. 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 ionized333This assumption is justified by the large post-shock
temperature and the presence of strong photo-ionizing X-ray and UV emission.
(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_{\nu}\propto\nu^{-(p-1)/2}\propto\nu^{-1}$ with $p\sim 3$,
Margutti et al. (2017) showed that the synchrotron emission gives a negligible
contribution by X-ray flux (i.e., $\sim 2$ orders of magnitude smaller 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 $\lesssim 5$
at $2\times 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.

Figure 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.

##  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), a density stratification $\rho\propto r^{-3}$ and extends from
$1.6\times 10^{16}$ cm to $1.87\times 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 $\rho\propto 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 $\sim 60/(v_{w}/100{\rm\;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}\sim 1.6\times 10^{16}$ cm,
and with a CSM wind at larger radii. They deduce a small mass for the wall
($\sim 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 $\dot{M}_{w}\sim 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 $\sim 2\times 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$\alpha$ 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\times 10^{16}$ cm, and expanding with a velocity of $5800\pm
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 $\sim 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 $\propto n_{H^{0}}\sigma_{\nu}$, being $n_{H^{0}}$ the
density of neutral hydrogen and $\sigma_{\nu}$ the bound-free 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 $\propto n_{e}^{2}$. 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 $\sim 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 $\sim 4\times 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 $\gg 4\times 10^{16}/v_{w}\sim 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
$\sim$ 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 $\propto 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.

##  5\. 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_{\rm 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 $\rho\propto r^{-3}$. We also found that the shell is very
extended, with a size $\gtrsim 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 time-scales 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 environment 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$\alpha$ 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.

We thank Luc Binette, Cesar Fernández Ramírez, Leonardo Ferreira and Claudio
Toledo Roy for useful discussions. FV and FDC 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 No. AST-1909796 and AST-1944985.
Support for this work was provided by the National Aeronautics and Space
Administration through Chandra Award Number 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.

## References

  * Anderson et al. (2017) Anderson, G. E., Horesh, A., Mooley, K. P., et al. 2017, MNRAS, 466, 3648
  * Anupama et al. (2005) Anupama, G. C., Sahu, D. K., Deng, J., et al. 2005, ApJ, 631, L125
  * Ben-Ami et al. (2014) Ben-Ami, S., Gal-Yam, A., Mazzali, P. A., et al. 2014, ApJ, 785, 37. doi:10.1088/0004-637X/785/1/37
  * Bietenholz & Bartel (2005) Bietenholz, M. F. & Bartel, N. 2005, ApJ, 625, L99
  * Bietenholz & Bartel (2007) Bietenholz, M. F. & Bartel, N. 2007, ApJ, 665, L47
  * Bietenholz et al. (2014) Bietenholz, M. F., De Colle, F., Granot, J., et al. 2014, MNRAS, 440, 821
  * Bietenholz et al. (2018) Bietenholz, M. F., Kamble, A., Margutti, R., et al. 2018, MNRAS, 475, 1756
  * Bietenholz et al. (2021) Bietenholz, M. F., Bartel, N., Kamble, A., et al. 2021, MNRAS, 502, 1694. doi:10.1093/mnras/staa4003
  * Branch & Wheeler (2017) Branch, D. & Wheeler, J., C. Supernova Explosions (New York: Springer-Verlag Berlin Heidelberg)
  * Brethauer et al. (2020) Brethauer, D., Margutti, R., Milisavljevic, D., Bietenholz, M., arXiv:2020.2020
  * Brethauer et al. (2021) Brethauer, D., Margutti, R., Vargas, F., et al., in preparation (2021)
  * Burrows et al. (2005) Burrows, D. N., Hill, J. E., Nousek, J. A., et al. 2005, Space Sci. Rev., 120, 165. doi:10.1007/s11214-005-5097-2
  * Cantó et al. (2009) Cantó, J., Curiel, S., & Martínez-Gómez, E. 2009, A&A, 501, 1259
  * Chandra et al. (2020) Chandra, P., Chevalier, R. A., Chugai, N., et al. 2020, arXiv:2008.13724
  * Charbonneau (1995) Charbonneau, P. 1995, ApJS, 101, 309. doi:10.1086/192242
  * Chen et al. (2018) Chen, T.-W., Inserra, C., Fraser, M., et al. 2018, ApJ, 867, L31
  * Chevalier (1982) Chevalier, R. A. 1982, ApJ, 259, 302
  * Chevalier & Liang (1989) Chevalier, R. A. & Liang, E. P. 1989, ApJ, 344, 332. doi:10.1086/167802
  * Chevalier (1998) Chevalier, R. A. 1998, ApJ, 499, 810
  * Chevalier & Fransson (2017) Chevalier, R. A. & Fransson, C. 2017, Handbook of Supernovae, 875
  * Chugai & Chevalier (2006) Chugai, N. N., & Chevalier, R. A. 2006, ApJ, 641, 1051
  * Corsi et al. (2014) Corsi, A., Ofek, E. O., Gal-Yam, A., et al. 2014, ApJ, 782, 42
  * De Colle et al. (2012) De Colle, F., Granot, J., López-Cámara, D., et al. 2012, ApJ, 746, 122
  * De Colle et al. (2014) De Colle, F., Raga, A. C., Contreras-Torres, F. F., et al. 2014, ApJ, 789, L33
  * De Colle et al. (2018a) De Colle, F., Kumar, P., & Aguilera-Dena, D. R. 2018, ApJ, 863, 32
  * De Colle et al. (2018b) De Colle, F., Lu, W., Kumar, P., et al. 2018, MNRAS, 478, 4553
  * De Geyter et al. (2013) De Geyter, G., Baes, M., Fritz, J., et al. 2013, A&A, 550, A74
  * de Jager et al. (1988) de Jager, C., Nieuwenhuijzen, H., & van der Hucht, K. A. 1988, A&AS, 72, 259
  * Dwarkadas et al. (2010) Dwarkadas, V. V., Dewey, D., & Bauer, F. 2010, MNRAS, 407, 812. doi:10.1111/j.1365-2966.2010.16966.x
  * Foley et al. (2007) Foley, R. J., Smith, N., Ganeshalingam, M., et al. 2007, ApJ, 657, L105
  * Dwarkadas & Gruszko (2012) Dwarkadas, V. V. & Gruszko, J. 2012, MNRAS, 419, 1515
  * Freedman et al. (2001) Freedman, W. L., Madore, B. F., Gibson, B. K., et al. 2001, ApJ, 553, 47
  * Gehrels et al. (2004) Gehrels, N., Chincarini, G., Giommi, P., et al. 2004, ApJ, 611, 1005. doi:10.1086/422091
  * González-Casanova et al. (2014) González-Casanova, D. F., De Colle, F., Ramirez-Ruiz, E., et al. 2014, ApJ, 781, L26
  * Harris & Nugent (2020) Harris, C. E. & Nugent, P. E. 2020, ApJ, 894, 122. doi:10.3847/1538-4357/ab879e
  * Heger et al. (2000) Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368
  * Kim et al. (2014) Kim, M., Zheng, W., Li, W., et al. 2014, Central Bureau Electronic Telegrams, 3777, 1
  * Kuriyama & Shigeyama (2020) Kuriyama, N. & Shigeyama, T. 2020, A&A, 635, A127. doi:10.1051/0004-6361/201937226
  * Leitherer (2010) Leitherer, C. 2010, Hot and Cool: Bridging Gaps in Massive Star Evolution, 425, 171
  * Margutti et al. (2017) Margutti, R., Kamble, A., Milisavljevic, D., et al. 2017, ApJ, 835, 140
  * Matzner & McKee (1999) Matzner, C. D. & McKee, C. F. 1999, ApJ, 510, 379
  * Mauerhan et al. (2018) Mauerhan, J. C., Filippenko, A. V., Zheng, W., et al. 2018, MNRAS, 478, 5050. doi:10.1093/mnras/sty1307
  * Milisavljevic et al. (2015) Milisavljevic, D., Margutti, R., Kamble, A., et al. 2015, ApJ, 815, 120
  * Morisset (2016) Morisset, C. 2016, arXiv:1612.04250
  * Moriya & Maeda (2014) Moriya, T. J. & Maeda, K. 2014, ApJ, 790, L16. doi:10.1088/2041-8205/790/2/L16
  * Morrison & McCammon (1983) Morrison, R., & McCammon, D. 1983, ApJ, 270, 119
  * Nakar & Sari (2010) Nakar, E. & Sari, R. 2010, ApJ, 725, 904. doi:10.1088/0004-637X/725/1/904
  * Palliyaguru et al. (2019) Palliyaguru, N. T., Corsi, A., Frail, D. A., et al. 2019, ApJ, 872, 201
  * Pooley et al. (2019) Pooley, D., Wheeler, J. C., Vinkó, J., et al. 2019, ApJ, 883, 120
  * Rajpaul (2012) Rajpaul, V. 2012, arXiv:1202.1643
  * Rybicki & Lightman (1986) Rybicki, G. B., & Lightman, A. P. 1986, Radiative Processes in Astrophysics
  * Sakurai (1960) Sakurai, A. 1960, Commun. Pure Appl. Math., 13, 353
  * Salas et al. (2013) Salas, P., Bauer, F. E., Stockdale, C., & Prieto, J. L. 2013, MNRAS, 428, 1207
  * Schinzel et al. (2009) Schinzel, F. K., Taylor, G. B., Stockdale, C. J., et al. 2009, ApJ, 691, 1380
  * Smith & Arnett (2014) Smith, N. & Arnett, W. D. 2014, ApJ, 785, 82. doi:10.1088/0004-637X/785/2/82
  * Smith (2014) Smith, N. 2014, ARA&A, 52, 487
  * Smith (2017) Smith, N. 2017, Handbook of Supernovae, 403
  * Soderberg et al. (2006) Soderberg, A. M., Chevalier, R. A., Kulkarni, S. R., & Frail, D. A. 2006, ApJ, 651, 1005
  * Sun et al. (2020) Sun, N.-C., Maund, J. R., & Crowther, P. A. 2020, MNRAS, 497, 5118. doi:10.1093/mnras/staa2277
  * Sun et al. (2020) Sun, N.-C., Maund, J. R., & Crowther, P. A. 2020, MNRAS, doi:10.1093/mnras/staa2277
  * Suzuki et al. (2021) Suzuki, A., Nicholl, M., Moriya, T. J., et al. 2021, ApJ, 908, 99. doi:10.3847/1538-4357/abd6ce
  * Tan et al. (2001) Tan, J. C., Matzner, C. D., & McKee, C. F. 2001, ApJ, 551, 946
  * Tinyanont et al. (2019) Tinyanont, S., Lau, R. M., Kasliwal, M. M., et al. 2019, ApJ, 887, 75
  * Waxman & Katz (2017) Waxman, E. & Katz, B. 2017, Handbook of Supernovae, 967. doi:10.1007/978-3-319-21846-5_33
  * Weiler et al. (1991) Weiler, K. W., van Dyk, S. D., Panagia, N., et al. 1991, ApJ, 380, 161. doi:10.1086/170571
  * Wellons et al. (2012) Wellons, S., Soderberg, A. M., & Chevalier, R. A. 2012, ApJ, 752, 17
  * Yalinewich & Portegies Zwart (2019) Yalinewich, A. & Portegies Zwart, S. 2019, ApJ, 872, L3. doi:10.3847/2041-8213/ab01e4