Cloud Formation by Supernova Implosion

The deposition of energy and momentum by supernova explosions has been subject to numerous studies in the past few decades. However, while there has been some work that focused on the transition from the adiabatic to the radiative stage of a supernova remnant (SNR), the late radiative stage and merging with the interstellar medium (ISM) have received little attention. Here, we use three-dimensional, hydrodynamic simulations, focusing on the evolution of SNRs during the radiative phase, considering a wide range of physical explosion parameters ( nH,ISM∈0.1,100cm−3 and ESN∈1,14×1051erg ). We find that the radiative phase can be subdivided in four stages: A pressure-driven snowplow phase, during which the hot overpressurized bubble gas is evacuated and pushed into the cold shell; a momentum-conserving snowplow phase that is accompanied by a broadening of the shell; an implosion phase where cold material from the back of the shell is flooding the central vacuum; and a final cloud phase, during which the imploding gas is settling as a central, compact overdensity. The launching timescale for the implosion ranges from a few 100 kyr to a few Myr, while the cloud formation timescale ranges from a few to about 10 Myr. The highly chemically enriched clouds can become massive (M cl ∼ 103–104 M ⊙) and self-gravitating within a few Myr after their formation, providing an attractive, novel pathway for supernova-induced star and planet formation in the ISM.


INTRODUCTION
It has long been recognized that supernovae (SNe) play an important role in maintaining the balance and structure of the interstellar medium (ISM).Even though there is only about one supernova per 100 M ⊙ of formed stars, due to their enormous energy output, these destructive events can have an enormous impact on their surroundings.
Of particular interest is the concept of positive SN feedback or triggered star formation, where a strong shock wave compresses the gas, leading to further collapse, fragmentation and eventually the formation of new stars.This processes has been predicted in numerous theoretical works (e.g.Dwarkadas et al. 2017;Krause et al. 2018;Herrington et al. 2023) and has recently been confirmed with observations from the Gaia mission (e.g.Zucker et al. 2022;Miret-Roig et al. 2022;Ratzenböck et al. 2023).
SNR evolution in a uniform medium has been studied at great length using analytical models (e.g.Woltjer 1972;Gaffet 1978;Ostriker & McKee 1988), and numerical simulations in one (e.g.Chevalier 1974;Cioffi et al. 1988;Fierlinger et al. 2016), two (e.g.Blondin et al. 1998;Ntormousi et al. 2011;Meyer et al. 2023) and three dimensions (e.g.Kim & Ostriker 2015;Makarenko et al. 2023), which have lead to a comprehensive picture comprised of a series of different stages, characterized by different deceleration parameters q = −d 2 R/dt 2 and conserved quantities.In the first stage, known as the free expansion phase the SNR expands with constant velocity until the reverse shock has fully thermalized the ejecta (Truelove & McKee 1999).In the next so-called Sedov-Taylor (ST) phase (Sedov 1959;Taylor 1950), the SNR expands adiabatically as cooling losses are still negligible and thus energy is conserved.The ST phase ends, when radiative cooling losses become important and a thin, cold shell forms at the shock front.
After the shell has formed, the SNR keeps expanding in what is known as the pressure-driven snowplow (PDS) phase (Cox & Anderson 1982;Ostriker & McKee 1988).During the PDS, the hot bubble is rapidly evacuated as hot material is pushed into the shell (Gaffet 1983;Cioffi et al. 1988;Kim & Ostriker 2015).Once the bubble pressure has dropped below that of the shell, the PDS ends and transitions into a momentum conserving snowplow (MCS) phase (Cioffi et al. 1988;Thornton et al. 1998).It has been claimed that the SNR evolution ends, when the shock velocity becomes comparable to the typical velocity dispersion of the ambient medium (Cioffi et al. 1988;Draine 2011;Faucher-Giguère et al. 2013;Krumholz et al. 2018) and the shock merges with the ISM.However, the details of the merging have only received little attention and it remains unclear how the evacuated bubble carved out by the blast wave is refilled.
In order to address this gap, in this work, we utilize three-dimensional, hydrodynamical simulations with cooling to study the late radiative stage of SNR evolution and the onset of fade out.We thus provide a more complete picture for the later stages of SNR evolution.We show that as the SNR shell pressure approaches the ISM pressure, the SNR implodes, filling the central cavity.This implosion leads to the formation of a dense compact, cloud in its center, which has the potential to form new stars and thus provides a novel pathway for triggered star formation.
The remainder of this paper is organized as follows.In section 2 we describe the numerical scheme and the setup of our simulation suite.In section 3 we present the results of our numerical simulations.We discuss the limitations and implications of our results in section 4. Finally, we summarize our findings and conclude in section 5.In the appendix we present a number of tests, related to the question of numerical convergence and the adopted treatment for radiative cooling.

Numerical Methods
We utilize the adaptive mesh refinement (AMR) code ramses (Teyssier 2002) to simulate the hydrodynamic evolution of blast waves in a uniform density medium, including radiative cooling.ramses is solving the system of hydrodynamic equations utilizing a second-order unsplit Godunov method (MUSCL scheme) on a finite volume, cartesian grid.Variables at the cell interfaces are reconstructed from the cell-centered values using the HLLC Riemann solver (Toro et al. 1994) with MinMod total variation diminishing scheme.Cooling is solved for the default courty cooling function implemented in ramses, which provides a basic treatment of primordial chemistry, metal line cooling and heating due to ultraviolet background (UVB) radiation.
In our fiducial simulation suite we consider a cubic computational domain with a side length of L = 1024 pc and periodic boundaries, which we refine with l min = 7 to l max = 11 refinement levels, corresponding to a spatial resolution of ∆x max = 8 pc and ∆x min = 0.5 pc, respectively.However, we ensure that the resolution criterion of Kim & Ostriker (2015) is fulfilled and accordingly increase the resolution in runs, where the expected radius at shell formation would not be resolved with at least 10 grid cells.
Initially, the simulation domain is filled with uniform density gas with log n H cm −3 ∈ {−1, 0, 1, 2} at solar metallicity and an initial temperature set to be close to cooling equilibrium.In the domain center, we initialize the explosive ejecta uniformly within a spherical region of radius R inj ∼ 5 ∆x min .We inject E SN = 10 51 erg and M ej = 5 M ⊙ per SN, corresponding to an initial ejecta temperature of T SN ∼ 10 9 K.We ensure that the injection region is maximally refined, by statically refining the central R ref = 50 ∆x min with the maximum resolution.
In order to ensure that the shock and the bubble are maximally refined, while only as little as possible of the surrounding medium is refined, we advect a passive scalar variable Z ej with the injected mass and maximally refine all cells where Z ej > 10 −15 .The criterion might fail, if numerical errors in the advection of pristine cells trigger the criterion or if cells just behind the shock are not polluted enough to trigger refinement.However, we have checked that both of these cases do not occur frequently enough to cause any serious problems.In the Appendix A we discuss the role of the AMR in more detail.Ploeckinger & Schaye (2020).
1 Equation 3Besides different densities we also consider different explosion strengths, mimicking the feedback from a single stellar population with N SN ∈ {1, 5, 14} massive stars exploding all at once, by simply injecting N SN times as much mass and energy.We thus label a model with N SN = x, log n H = y and l max = z as Nx ny Lz.While this simplistic approach can capture some aspects of clustered feedback, it is worth noting that studies that take into account the time delay between explosions find some qualitative differences, such as an increased momentum per SN (Walch & Naab 2015;Gentry et al. 2019) and a longer lived hot bubble (Kim et al. 2017).
All models are run until t = 14 Myr at which point the largest bubbles are reaching the domain limits.We reran some of the models in a smaller domain for a shorter time span, but with a much higher frequency of snapshots.For these models we add the suffix HC to the name and the part of the name referring to the resolution refers to the equivalent refinement level for the fiducial domain, i.e.L11 refers to ∆x min = 0.5 pc in both the fiducial and the HC runs.
In appendix C we will discuss the effect that different cooling functions could have on our results.To this end, we rerun a few of the models with different cooling tables taken from Ploeckinger & Schaye (2020) and label the models with suffices corresponding to the respective alternate cooling model.
Finally, in section B we are discussing the results obtained for the N1 n2 model at different resolutions.
A list of all the different models and their properties is given in table 1.

Data Analysis
In order to quantify the global evolution of the gas, we distinguish between bubble and shell gas and further within these components, differentiate between gas which is moving radially outward (v r > 0) and inward (v r < 0).We distinguish between the SNR and the ISM using the passive scalar, i.e. gas with Z ej > 10 −10 is considered part of the SNR.The bubble is defined as SNR gas that is either hot (T > 2 × 10 4 K) or diffuse (n H < 10 −2 cm −3 ), while the shell is all SNR gas that is not part of the bubble.A summary of the classification is given in Figure 1.Kim & Ostriker (2015) use a similar criterion for the bubble gas, considering only the temperature of the gas.
Figure 1.Schematic overview of the different gas components, described in section 2.2.The classification differentiates between the ISM (gray) and the SNR (color), consisting of a bubble (shades of red) and a shell (shades of blue), which itself might be out-or inflowing.
The addition of the density criterion only becomes important at late times, when the bubble has cooled below 10 4 K at which point the temperature criterion alone would fail to differentiate between the bubble and the shell.
The partition of the SNR into a bubble and shell, allows us to measure the shell formation timescale t sf , which denotes the time when the cold shell at the shock front forms.We follow Kim & Ostriker (2015), who define the numerically measured shell formation timescale t n sf as the time at which the mass of the hot bubble reaches its maximum.

RESULTS
In this section we describe a mechanism for the formation of a cloud through the implosion of a radiative SNR.In section 3.1 we give a brief overview of the physical mechanism.We provide a detailed description of our simulation results for a single SN in a high density ISM in section 3.3.In section 3.4 we extend our analysis to the whole simulation suite and investigate the dependence of the relevant timescales on the explosion parameters in section 3.4.1.We describe the properties of the implosion clouds in section 3.5.Finally, in section 3.2 we describe a model for the launching of the implosion.

Schematic Overview
Figure 2 gives a schematic description of the evolution of the SNR after shell formation, which is separated into four stages.
In the first stage, shortly after shell formation the interior of the SNR is still hot and overpressurized relative to the isothermal shell, which is kept at about T ∼ 10 4 K.During this stage the rarefied bubble material tries to expand and as a result is pressed into the shell.This phase corresponds to the modified PDS phase described by Cioffi et al. (1988) and Kim & Ostriker (2015).After about 3 − 5 shell formation timescales the bubble has been evacuated and its pressure has dropped below that of the shell.
At this point, the expansion of the SNR is entirely inertia driven, corresponding to the MCS.The mass of the evacuated bubble is negligibly small and its density is many orders of magnitude below the ambient density.The shell which is now overpressurized relative to both the ISM and the bubble, begins to broaden, leading to a gradual reduction in the shell pressure and an accelerated weakening of the shock.Meanwhile, the shell begins to get deformed and fragmented due to thin-shell overstability and nonlinear thin-shell instability (Vishniac 1983;Blondin et al. 1998).
Once the shell pressure reaches pressure equilibrium with the ISM a reflected version of the outgoing shock wave is launched driving cooling material from the shell back into the center, refilling the warm, evacuated cavity.We refer to this reflected wave as Implosion or Backflow.The implosion is very similar to the socalled "negative phase" in the context of terrestrial blast waves (see e.g.Glasstone & Dolan 1977) and appears to be a purely hydrodynamic realization of the hydromagnetic Rayleigh-Taylor instability (RTI) described by Breitschwerdt et al. (2000).
After a few Myr, the backflow reaches the center of the SNR and collides with the backflowing gas from all directions.The colliding gas piles up in the center and is reflected, forming a slowly expanding cloud.The cloud keeps accreting material from the backflowing gas reaching a mass of 10 3 − 10 4 M ⊙ within ∼ 10 Myr.As the cloud is directly formed from the SN ejecta, it is highly chemically enriched.

Model for the Launching of the Backflow
In the previous subsection we have given an overview of the different phases of the radiative stage.Here we describe a model for estimating the relevant timescales.
Right after shell formation, we assume that the bubble is following the modified PDS described by Kim & Ostriker (2015).In their description the thermal energy Opacity alpha is scaled with the logarithm of the density and is set to zero for densities in the range nH ∈ (95, 105) cm −3 in order to remove the foreground.The octant facing the observer has been made transparent.
of the bubble evolves as where t sf and R sf are given by (Kim & Ostriker 2015) Here E 51 = E SN / 10 51 erg .The average bubble pressure is then given by Meanwhile, the temperature of the shell remains roughly constant at T shell ∼ 10 4 K and the compression ratio of the shell is χ ∼ 10, leading to a shell pressure of P shell, PDS ∼ 10 5 n 0 k B K cm −3 .The PDS phase ends, when the pressure in the shell and bubble is equal, at The radius of the SNR at this time is After the PDS phase has ended, the momentum of the shell remains constant and the MCS phase begins.During this phase the radius of the SNR evolves ∝ t 1/4 and correspondingly the shock velocity evolves ∝ t −3/4 .If one assumes a constant compression ratio and that the temperature of the shell is proportional to the square of the velocity, as one would expect for a strong shock, one finds for the pressure during the MCS phase The SNR implodes, when the pressure of the shell approaches the pressure of the ISM.In the standard RAM-SES cooling prescription at solar metallicity, which assumes collisional ionization equilibrium the pressure on the cooling-equilibrium curve for a given density is approximately (see e.g. Figure 17 of Kim et al. (2023)) Thus, the launching timescale can be inferred as and the radius of the SNR at this time is Finally, cloud formation happens once the imploding shell has reached the center.The cloud formation timescale is thus the combination of the launching and the crossing timescale where V in is a characteristic inflow velocity.The results in section 3.4.1 indicate that this velocity is independent of the explosion energy, but depends on the ambient density in a complicated way.

Single Supernova in a High Density Medium
In this subsection we describe the time evolution of the SNR formed by a single SN in a stationary, uniform, high density medium in cooling equilibrium.
We first give a qualitative overview of our ultra high resolution model N1 n2 L14 with a maximum grid resolution of ∆x = 0.0625 pc.In Figure 3 the density (top panels) and radial mass flux (bottom panels) is shown in slices through the xy-plane at four different points in time corresponding to the four different stages of SNR evolution after shell formation.In order to provide a quantitative reference, radial profiles of various physical quantities at the same points in time are shown in Figure 4.During the first stage, right after shell formation, the shell reaches a maximum compression ratio of χ ∼ 10.Its width is comparable to the resolution limit ∆R ∼ ∆x = 0.0625 pc.The pressure of the hot interior is about an order of magnitude higher than that of the rapidly cooling shell, and as a result bubble material is condensing onto the shell.
After 0.1 Myr, the bubble has been mostly evacuated (n H, Bubble /n H, ISM < 10 −4 ) and the shell has thickened considerably (∆R ∼ 3 pc), leading to an overall reduction in the compression ratio (χ ∼ 2.5) of the shell.The pressure in the bubble has dropped significantly to about 10 % of the ISM pressure, while the pressure of the shell is still overpressurized with respect to the ISM.During this stage the mass flux is concentrated within the shell, with only little outward mass flux from inside the bubble.The shell is subject to thin-shell overstability and nonlinear thin-shell instability (Vishniac 1983(Vishniac , 1994;;Blondin et al. 1998), resulting in ripples on the shell's surface.The differently shaded curves correspond to t = 0.006, 0.1, 1 and 10 Myr, respectively.In panel (c) solid (dashed) lines correspond to outward (inward) mass flux.
We note that we do not explicitly seed perturbations that would drive these instabilities.Instead, they arise from grid scale perturbations due to the mapping of the sphere onto a Cartesian grid and numerical instabilities such as the carbuncle instability (see e.g.appendix C of Stone et al. 2008).We refer the interested reader to appendix B, where we discuss in more detail the dependence of these artifacts on the resolution and how they might affect our results.
After 1 Myr the flow just behind the shell has reversed and is now flowing inward, starting to slowly fill up and cool the bubble with material from the backside of the shell.Meanwhile, the pressure of the shell has dropped to a level comparable to the ISM pressure.The shell continues to broaden (∆R ∼ 7 pc) and the compression ratio (χ ∼ 1.5) continues to drop, approaching unity.During these first three stages, most of the mass of the bubble gas is composed of ejecta material.
Finally, after 10 Myr the inward flow has reached the center and is compressed into a compact, dense, expanding cloud with a constant density of about 3 times the ISM density and a size of about 5 pc.The mass fraction of the ejecta in the cloud is up to an order of magnitude larger than in the rest of the SNR.The outer radius of the cloud is bounded by more inflowing material, which is slightly underdense relative to the ISM.Meanwhile, the shell has broadened to about ∆R ∼ 20 pc and the compression ratio is only slightly above unity.The shell instabilities have lead to complex substructure within the shell.The shell is composed of many ∼ pc size blisters, which are bounded by outflowing, overdense shells and filled with inflowing, underdense gas.
In order to describe the launching mechanism of the implosion in more detail, in Figure 5 we show the time evolution of various globally computed quantities for the model N1 n2 L13 HC with a maximum grid resolution of ∆x = 0.125 pc.Here we use the L13 model, because due to storage limitations it was not feasible to run an L14 model with high output cadence.The different panels show the time evolution of mass, momentum, kinetic & thermal energy, pressure, volume and ejecta mass, calculated for the different gas components described in Figure 1.
Extensive quantities are computed by summing up the contributions from each cell belonging to the respective gas components.The volume averaged pressure of each gas component i is calculated as where γ = 5/3 is the adiabatic index, E th,i is the thermal energy and V i is the volume.
After an initial relaxation period the solution approaches the energy-conserving ST phase, during which the entire SNR is hot and ∼ 70 % of the energy is thermal, in agreement with analytical calculations (Taylor 1950;Sedov 1959).
Shortly before shell formation, after t ∼ 10 −3 Myr, a reverse pulse emerges, is reflected in the center and merges with the shock again.
At shell formation t sf ∼ 3 × 10 −3 Myr the bubble mass reaches a maximum and the shell mass begins to increase.Similarly the momentum of the bubble gas reaches a maximum as the momentum of the shell starts to increase.The total thermal energy begins to drop steeply, while the kinetic energy remains constant.The thermal energy is dominated by the bubble while the shell carries negligible amounts of thermal energy.On the contrary the kinetic energy of the bubble drops rapidly and is taken over by that of the shell.The pressure of the shell after its formation is initially roughly constant and much lower than that of the bubble, which however drops rapidly.The volume of the SNR is dominated by the hot bubble, which after shell formation initially decreases until it reaches a fixed volume filling factor of about ∼ 2/3, with a corresponding volume filling factor of the shell of about ∼ 1/3.Most of the ejecta stay within the bubble and only slowly get incorporated into the shell.This behavior is in line with the modified version of the PDS phase described in 1D by Cioffi et al. (1988) and in 3D by Kim & Ostriker (2015).2).In panel (c) the thermal (solid) and kinetic (dashed) energy is shown.The horizontal green dashed line in panel (d) marks the equilibrium pressure of the ISM.The bubble pressure is initially very high (≳ 10 9 K cm −3 ) and thus outside of the axis limits.
After about 10 −2 Myr most of the SNRs mass, momentum and kinetic energy is carried by the still entirely outward moving shell.Until this point the momentum has still been increasing, but at this point the increase stops, marking the beginning of MCS phase.The kinetic energy starts to decrease as t ∝ t −0.75 , consistent with the analytical expectation.The pressure of the shell begins to decrease due to the effect of radiative cooling, while the pressure of the bubble keeps decreasing rapidly.About 10 % of the ejecta are now in the shell.
After 2 × 10 −2 Myr the pressure of the bubble falls below that of the shell.At this point the volume filling factor of the bubble begins to decrease and that of the shell correspondingly has to increase, corresponding to a relative broadening of the shell.At this point the majority of the ejecta mass is in the shell.
The backflow emerges after t launch ∼ 0.3 Myr at the same time as the pressure of the shell approaches the ambient pressure.This inward flow is quite different from the series of reflected sound waves described by Cioffi et al. (1988).While the sound waves are associated with the bubble and carry negligible amounts of mass, the implosion is associated with the shell and carries relatively large amounts of mass, which are growing.The pressure of the outward moving shell levels off at the ISM pressure, while that of the inward moving com-ponent increases.The emergence of the backflow leads to a slight decrease in total radial momentum.
After 1.5 Myr, about 1 % of the shell mass is moving inward, carrying about 1 % of the ejecta back to the center.At this point, the bubble has disappeared entirely.

Universality of the Mechanism
In Figure 6 we show a rescaled version of Figure 5 for the models N1 n0 L11 HC, N1 n1 L11 HC and N1 n2 L11 HC.Time is measured in units of the shell formation timescale t n sf .Mass, momentum and volume are normalized to their value at t n sf and pressure is normalized to the ISM value in cooling equilibrium.
As expected from the self-similarity during the ST phase (Sedov 1959), all models exhibit a very similar time evolution before shell formation.In the models N1 n0 L11 HC and N1 n1 L11 HC there is already some small amount of cold gas during this phase.This is due to the method used to extract SNR gas, which might include a negligible number of unshocked cells.
After shell formation, the time evolution in all models is qualitatively the same as for N1 n2 L13 HC, though timescales in units of the shell formation timescale may differ.In the models N1 n0 L11 HC and N1 n1 L11 HC the pressure of the outward moving shell starts to decrease after the bubble pressure falls below it after 3 -4 t n sf , while in N1 n2 L13 HC it already starts to decrease already slightly before that.We interpret this difference as follows: The pressure can either decrease by radiative cooling or adiabatic expansion.When the pressure of the bubble drops below the shell pressure before radiative cooling becomes important, the shell will cool adiabatically.On the other hand, if the cooling timescale is shorter than the timescale to reach pressure equilibrium between the shell and the bubble, the shell will start to cool radiatively before the bubble pressure has dropped.In both cases the shell pressure will start to decrease, albeit with a slightly different scaling.It is therefore no surprise that in the run with a higher density and therefore a shorter cooling timescale, the pressure starts to drop slightly earlier.
In all models, once the pressure of the shell approaches that of the ISM, a steadily growing backflow is launched.

Timescales
Having established the universal emergence of backflows, once the shell pressure approaches that of the ambient medium, we may investigate next, how the timescales for launching the backflow t launch and subsequently forming a central overdensity t cf depend on the explosion parameters.
We define the launching timescale as the earliest time, when the inflowing shell mass exceeds 0.1 M ⊙ .We choose this threshold, because especially in the runs, where the shell is resolved with many cells, instability of the shell itself can lead to eddies within the shell, which lead to a small amount of inflowing shell gas that is not associated to the implosion.The cloud formation timescale is defined as the earliest time after launching, when the density in the innermost radial bin exceeds the ambient density.
Due to the limited temporal resolution of the snapshots, the events actually occur somewhere between the first snapshot when the above conditions are met and the previous snapshot.We therefore report the arithmetic mean of the two time points and indicate the time interval between the two snapshots with error bars.
In Figure 7 the timescales are shown as functions of the ambient density (left panels) and the explosion energy (right panels).The markers are colored by the respective other variable.
We find launching timescales between a few hundred kyr and few Myr.There is a negative trend with density and a positive trend with energy.The scaling and normalization are in rough agreement with Equation 10, with factor of ∼ 2 differences.
At high densities differences might arise, because in the derivation of Equation 10we have neglected the role of radiative cooling below T ∼ 10 4 K for the shell gas (see also discussion of Figure 6), which would lead to a shallower scaling of t PDS (Eq. 6) with density, which in turn manifests as a steeper scaling of t launch .
At the low density end, the differences might arise, because here P shell, PDS is already quite similar to P ISM, eq and therefore the assumed scaling might not apply since the shock is already quite weak at t PDS .
The cloud formation timescale is on the order of several Myr to over 10 Myr for our simulations with the highest explosion energies and lowest densities.At n H = 0.1 cm −3 no overdense clouds have formed by the end of our simulations.There is no clear trend with density.The cloud formation timescale becomes shorter for SNRs in higher density environments, though for high explosion energies the timescale appears to level off and it even slightly increases for the highest densities.The timescale generally increases with energy, roughly scaling like t cf ∝ E 0.3 51 .The scaling with the explosion energy appears to be the same as the scaling of the radius of the SNR at launching R launch (E 51 , n H ), as predicted by Equation 11.Equation 12 then implies that the implosion velocity V in is independent of explosion energy and depends only on density.

Cloud Properties
In the previous section we have shown that for a wide range of explosion parameters, the SNRs implode and form a cloud in their center.Here we summarize the properties and evolution of these clouds and discuss how they depend on the explosion parameters.
To this end, we utilize radial profiles of density, kinetic and thermal energy density to compute the size, mass and virial parameter of the clouds.We note that our simulations are strictly without self-gravity.The timescale for self-gravity to have a qualitative effect on the evolution of the SNR is much longer than the simulation time.Nonetheless, it might have an effect on the evolution of the dense cloud, and in order to estimate the importance of self-gravity for the cloud and predict whether or not it might become self-gravitating and eventually form stars it is useful to look at the virial parameter.
We define the size of the cloud as the radius of the interface between the innermost radial bin, where the density falls below the ambient density.The cloud mass is defined as the integral of the density profile up to that radius and the virial parameter is defined as the ratio of the sum of kinetic and thermal energy and the modulus of the potential energy of the cloud.The kinetic and thermal energy are computed by integrating over the respective radial profiles and the potential energy is defined as Here we assume that the density profile within the cloud is flat and that all motion within the cloud is opposing the gravitational pull.Both of these assumptions approximately hold true (see e.g. Figure 4).
In Figure 8 we present the time evolution of the cloud's size, mass and virial parameter for the different models.
In panel a) the cloud radius is shown as a function of time.The cloud grows to a size of several to few tens of parsecs within 10 Myr.In lower density environments clouds grow larger, with little dependence on the explosion energy.
Panel b) shows the cloud's mass as a function of time.At the end of the simulation, the mass of the central cloud has reached a value of several thousand to about 2×10 4 solar masses, with only a small dependence on the density of the ambient medium.Clouds in high density environments become somewhat more massive.There is only little dependence on the explosion energy.
Panel c) shows the cloud's virial parameter as a function of time.The initial virial parameter scales with the ambient density roughly as where n 0 = n H /cm −3 , and decreases steadily approaching and dropping below α vir ∼ 1 for the high density runs, suggesting that indeed the clouds would become self-gravitating.There is little variation due to the explosion energy.
The masses and radii of the clouds formed by SN implosion are comparable to the values used in the initial conditions for the clouds in the STARFORGE (Grudić et al. 2022(Grudić et al. , 2023;;Farias et al. 2023) simulation suite, which studies the star formation from the collapse of a single giant molecular cloud.This suggests, that the thus formed clouds might indeed trigger another generation of star formation.

DISCUSSION
In the previous section we have described a new mechanism by which a cloud can form inside a radiative SNR due to its implosion as the shell pressure approaches that of the ISM.In the following, we will discuss some of the limitations of our model, the role of some of our model ingredients and the implications of our findings in the context of galaxy evolution.

Limitations
In this work, we have shown the existence of SN implosions and subsequent cloud formation in the center of the SNR, using a suite of hydrodynamic simulations of SN explosions in a uniform and stationary medium.Of course, a uniform and stationary ISM is a great simplification of the complexity of a realistic ISM.
In a more realistic model for the ISM, like the kind of turbulent, stratified box used in state-of-the-art ISM simulations (see e.g.Walch et al. 2015;Kim & Ostriker 2017), the existence of backflows as described in our work is a priori not guaranteed.Continuous or sufficiently frequent energy injection can keep the SNR overpressurized and prevent an implosion (Kim et al. 2017).Indeed, stellar populations are expected to explode SNe in regular intervals before they run out of fuel (see e.g.Leitherer et al. 1999).However, in the cases of runaway stars or populations hosting a sufficiently small number of massive stars this limitation does not apply and even in the case of stellar populations that remain active for a long time, eventually the bubbles are going to evacuate and cool off enough to make an implosion feasible.
Besides the importance of continuous driving, the role of the ambient medium cannot be ignored.In a more realistic description of the ISM, the ambient medium is highly structured due to the combined effect of turbulence, shear and stratification.SNe exploding in such an environment will follow the geometry of the ISM (Makarenko et al. 2023), as the shock wave can only slowly penetrate into dense structures, but will quickly fill out the volume filling low density medium, leading to highly amorphous SNR shapes (Kim & Ostriker 2015;Lancaster et al. 2021).In such a configuration, the SNR will reach pressure equilibrium at different times in different directions, leading to a displacement and deformation of the clouds formed in this way.If the momentum carried by the backflowing gas from different directions is not equal and opposite, the cloud would further end up with a net momentum leading to a drift.Similar asymmetries can follow from the interaction of the SNR with neighboring shocks, e.g.due to neighboring superbubbles (Breitschwerdt et al. 2000) or small-scale turbulence, which can locally contribute to the pressure opposing the shock expansion.In Appendix D we show, that asymmetries in the ambient pressure can indeed trigger an implosion locally.
Besides the limitations due to the environment, more complete physics might also qualitatively modify our conclusions.
As mentioned in section 3.5 the virial parameter of the clouds drops below unity a few Myr after their formation.While it seems plausible that self-gravity is too weak to have an important effect on the implosion mechanism and the crossing of the inflowing gas, once the cloud has formed, it might collapse and fragment due to its self-gravity.Kim & Ostriker (2015) have demonstrated that magnetic fields play only a subdominant role in SNR evolution.On the contrary, Gentry et al. (2019) show that magnetic fields suppress the growth of instabilities at the bubble-shell interface, which reduces mixing and delays radiative cooling in the case of multiple consecutive SNe.However, it is important to note that Kim & Ostriker (2015) used a mesh code, while Gentry et al. ( 2019) used a Lagrangian method.While both methods achieve comparable resolution in the shell, the mass resolution in the bubble is several orders of magnitude higher in mesh codes, which are thus likely less affected by spurious cooling of bubble gas.Nonetheless, contributions to the pressure from a magnetic field and cosmic rays can potentially alter the timescales for implosion and cloud formation or even prevent these processes, if they can maintain a high enough bubble pressure.Indeed, in the case of the hydromagnetic RTI Breitschwerdt et al. (2000) have shown, that the magnetic field has a stabilizing effect, that implies a characteristic length scale at which the instability can act.In this case the backflowing gas fragments into blobs of size similar to the fastest growing wavelength.Sharma et al. (2014) conducted 1D spherically symmetric simulations of SNe.They note that most of the heat losses occur in the unresolved layer between the bubble and the shell.The physical width of this layer is much too small to resolve, even with their 1D method, but they show that nonetheless the cooling losses are converged, even for moderate resolution.
Similarly, as discussed by Fierlinger et al. (2016), while the length scale associated with thermal conduction is much too small to be resolvable with current techniques, its effect is negligible.However, recent results by El-Badry et al. (2019) in the context of continuously driven superbubbles suggest that heat conduction does in fact play an important role for the transport of energy and mass across the bubble-shell interface.However, it is worth nothing that El-Badry et al. ( 2019) artificially enhanced the conduction rate to model turbulent mixing due to 3D instabilities, which makes a direct comparison difficult.Lancaster et al. (2021) simulate the expansion of a continuously driven wind bubble in a turbulent medium confirming that turbulence indeed enhances the mixing across the bubble-shell interface leading to catastrophic cooling losses in ideal hydrodynamics.Yet, the importance of these effects remains unclear in a picture where magnetic fields suppress the growth of the instabilities responsible for the mixing.Further studies of individual SNe with resolved heat conduction, magnetic fields and turbulence are required in order to settle this ongoing debate.

Role of the cooling model
In section 3.2 we have presented a model for the launching of the backflow.There are two ingredients of our model that are sensitive to assumed cooling physics.First, we have assumed that the temperature of the shell right after shell formation remains stable at 10 4 K, as the time it takes to cool beyond this temperature is much longer than the dynamical timescale.Second, we assume that the ISM is in cooling equilibrium, when we equate the shell pressure with the ISM pressure.Both of these assumptions invoke the assumed cooling physics, which may affect the resulting timescales and cloud properties.Kim et al. (2023) compare their detailed radiative transfer model to a range of commonly used cooling functions (see their Figure 17).They find that in these functions the equilibrium pressure at a given density may differ by up to three orders of magnitude between the models.In the Ramses cooling model utilized in this work the equilibrium pressure is indeed relatively high.As a consequence the SNR is expected to implode 100 times earlier with our cooling model, than e.g. with the model by Ploeckinger & Schaye (2020), which has an exceptionally low equilibrium pressure.Further complications like thermal instability and non-equilibrium effects (Katsuragawa et al. 2022), might also qualitatively alter our conclusions.
For a more detailed discussion, we refer the interested reader to appendix C, where we compare the results of simulations with different cooling models.The comparison indicates, that indeed a lower ambient pressure will delay the implosion and that the details of the cooling physics, may have a slight effect on the details of cloud formation, even at a comparable ambient pressure.

SN Implosion in the Literature
Despite the fact that the evolution of SNRs in a uniform medium has been studied in great detail for more than 40 years, to our knowledge there has been no mention of SN implosion.Here we discuss the various reasons for why this process might have remained unnoticed for so long.
Many authors (e.g.Chevalier 1974;Straka 1974;Thornton et al. 1998;Kim & Ostriker 2015) only focus on the transition to the radiative phase and would therefore not advance their simulations far enough to reach the implosion stage.Cioffi et al. (1988) set the pressure of the ambient medium to an artificially low value in order to maintain a strong shock, which in turn delays the implosion.Fierlinger et al. (2016) use one-dimensional hydrodynamic simulations to study the energy input from SNe in a uniform and stationary media with a range of densities and an initial temperature of 1000 K.They advance their models until the shock velocity approaches the sound speed of the ISM and thus in principle should have been able to see an implosion.However, the internal structure of the SNRs was out of their scope and thus they did not report any backflow.Gent et al. (2020) utilize three-dimensional hydrodynamic simulations of radiative SNRs to validate the PENCIL code.Even though they advance their models until the shock becomes sonic, as they mostly focus on benchmarking their code with previous work, they do not report any backflow.Breitschwerdt et al. (2000) use linear perturbation theory to describe a type of hydromagnetic RTI that, in the limit of vanishing magnetic field, is very similar in nature to the SN implosion described here.They show, that the interaction of two SNRs can lead to an inward flow of clouds originating from the interaction region.Our results, which appear to correspond to the same kind of instability, seem to confirm the linear prediction of Breitschwerdt et al. (2000).
To our knowledge there are no instances of SN implosions reported in observations.This might however simply attest to the fact, that radiative SNRs, in particular those close to merging with the ISM, are very dim and thus are often difficult to observe (Green 2019;Koo et al. 2020;Zhou et al. 2023).Moreover, the fact that the morphology of imploding SNRs differs qualitatively from traditional SNRs might have lead to a misclassification of imploding SNRs as something other than a SNR.

Implications for Galaxy Evolution
We have shown that under quite general conditions, old SNRs will implode and form a compact, massive, and potentially self-gravitating cloud in its center.Such a cloud could collapse and fragment under its own selfgravity to form stars, suggesting a novel mode of positive feedback.
Furthermore, given that such clouds would be highly enriched with the SN ejecta, and therefore with shortlived radionuclides (SLRs) like 26 Al, this provides an attractive pathway to the formation of planetary systems, where the heating due to SLRs plays an important role (Urey 1955).Indeed, it has been concluded by Forbes et al. (2021) that the enrichment with SLRs would have to occur prior to core formation, a condition that at face value is readily fulfilled by our proposed mechanism.However, further studies are necessary to investigate to what extent mixing due to small-scale turbulence might further dilute the imploding gas.
Besides the importance for star and planet formation we make a clear prediction for the lifetime of hot cavities in the ISM, which are filled shortly after the implosion is launched.This number is an important parameter in models for the multiphase structure of the ISM (McKee & Ostriker 1977;Wolfire et al. 2003;Draine 2011), which are used to estimate a wide variety of ISM properties.

CONCLUDING REMARKS
We have performed 3D hydrodynamic simulations of SNRs in a uniform, stationary medium with nonnegligible thermal pressure in order to study their evolution after shell formation.Our simulations reveal that radiative SNRs implode after the shell reaches pressure equilibrium with the ISM.The implosion leads to the formation of a compact, massive cloud that might soon become self-gravitating and that is highly enriched with SN ejecta.As we discuss, this novel mechanism of cloud formation provides attractive initial conditions for star and planet formation and might have some important implications for the theory of the ISM.
While the idealized setup is useful for understanding the underlying physical mechanism, understanding the role of SN implosion and subsequent cloud formation in a more realistic setup deserves further investigation.
We conclude that the dispersal and merging of SNRs with the ISM offers a wealth of hidden complexities, which deserve further study as they can help understand the physics of the ISM.
We are grateful to Guang-Xing Li for fruitful discussions.We thank the anonymous referee for their insightful comments and suggestions that helped to improve the quality of this work.Computations were performed on the HPC system Cobra at the Max Planck Computing and Data Facility.This research was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy -EXC 2094 -390783311.
Software: Julia v1.6.5 (Bezanson et al. 2014) In our fiducial simulation suite presented in the main body of the paper we have used AMR to reduce the numerical cost of the simulations.In our prescription we only refine cells, which have been sufficiently enriched by SN ejecta.We thus create a central refinement region that should expand at roughly the same rate as the shock.However, as the ejecta are physically confined behind the contact discontinuity, the refinement region might lag behind the shock and thus one might expect the shock to be slightly less refined than the rarefied gas behind it.
In SNRs the contact discontinuity and the shock are essentially at the same location and thus we expect this lag to be negligible, especially considering the presence of numerical diffusion, which acts to smear out the contact discontinuity.
In order to test whether the use of AMR can qualitatively modify our results in Figure 9 we compare the results of the models N1 n2 L13 and N1 n2 L13 noAMR.The two models are almost identical, but while in N1 n2 L13 we use our AMR prescription in N1 n2 L13 noAMR we use a static mesh with a grid spacing equal to the finest spacing in N1 n2 L13.
Figure 9 shows that the two models are essentially identical.Negligibly small differences arise due to small differences in the time when snapshots are written, since they are only written when the entire grid is synchronized.
We thus conclude that the adopted AMR technique does not affect our results in any meaningful way.

B. CONVERGENCE
In order to test, whether our results are converged, we ran the N1 n2 model at three different resolutions, denoted by N1 n2 L12, N1 n2 L13 and N1 n2 L14.The results are shown in Figure 10.Shown are the mass, momentum, energy, pressure, volume and ejecta content of the different gas components as a function of time.
Before shell formation, there are only minor differences between the models, which seem to stem from slightly different snapshot times and the analysis procedure.There is a slight dip in the mass, momentum and energy of the hot phase, in the lowest resolution run, which likely might have arisen from the SNR selection criterion (Figure 1).
The shell formation timescale, as well as the properties of the SNR at shell formation appear converged in line with the convergence criteria by Kim & Ostriker (2015).
Slight differences appear after shell formation.The pressure of the cold phase, during the PDS increases with resolution.Since the temperature of the shell during this phase is roughly fixed to about 10 4 K, the pressure is determined by the density of the shell, which depends on the width of the shell.In all runs, the shell during the PDS is resolved with only few resolution elements and thus indeed is not converged, as already noted by Sharma et al. (2014) who estimate that the width of the shell should be ∼ 0.001 R sf ∼ 10 −3 − 10 −2 pc.Despite the different PDS shell pressure, the PDS ends at a similar time in all runs and the shell pressure approaches a very similar evolution during the MCS and following phases.
In the lower resolution runs more of the ejecta tend to be incorporated in the shell right after shell formation.However, these differences are negligible by t = 0.1 Myr.
During the MCS, prior to implosion, there is already a small non-growing fraction of backflowing cold gas in the highest resolution run.This component is in pressure equilibrium with the ISM suggesting that it arises from small blisters in the shell's outer edge, which form, when the shell fragments due to thin shell shell instabilities which are unresolved in the lower resolution runs.
The implosion is launched at roughly the same time after ∼ 300−500 Myr and forms a central cloud (V hot → 0) after 1 Myr in all runs, indicating that these timescales are converged.The mass, momentum, thermal energy and volume of the cold components diverge at late times.Despite this non-convergence all SNRs exhibit the same qualitative features, and thus it does not affect our conclusions in any meaningful way.

C. COOLING MODELS
In order to investigate the role of the cooling function, we have rerun the model N1 n1 L11 with two different cooling models taken from Ploeckinger & Schaye (2020) integrated using an exact integration scheme (Townsend 2009;Zhu et al. 2017).For details of the implementation we refer the reader to Behrendt et al. (in prep.).The cooling functions in the models N1 n1 L11 Dust (Dust) and N1 n1 L11 PS20 (PS20) correspond to their models UVB dust1 CR0 G0 shield0 and UVB dust1 CR1 G1 shield1, respectively.The models' cooling-equilibrium curves in the P − n H -plane are shown in Figure 11.
Dust differs from the Ramses cooling model in that the equilibrium curve has a pronounced kink at around n H ∼ 1 cm −3 , while PS20 has a steep drop in pressure at around n H ∼ 0.01 cm −3 above which the pressure is 2-3 orders of magnitude below the pressure in the Ramses model.
In Figure 12 we show a comparison of the results for the different models.As expected the SNR evolution before shell formation is not affected by the cooling.The length of the PDS phase and the momentum boost during this phase are slightly increased for the Dust and PS20 models.Despite these differences, the timescales for implosion and cloud formation hardly differ between the Ramses and Dust models.On the contrary, for the PS20 model, since the equilibrium ISM pressure is several orders of magnitude lower, implosion is delayed by ≳ 3 Myr and does not lead to cloud formation within the simulated time.
These results indicate that it is indeed the ISM pressure, which controls the implosion, since all other quantities that might have an effect do not differ very much between the runs.

D. LOCALLY TRIGGERED IMPLOSIONS
In order to investigate, whether or not the implosion can be triggered by local pressure enhancements, we have rerun the model N1 n-1 L11 for 1.5 Myr with a cold (T ∼ 800 K), dense (n H ∼ 10 cm −3 ) cloud of Radius R cl = 15 pc centered at a distance of d cl = 100 pc from the explosion center.The pressure in the cloud is about an order of magnitude higher than that of the ambient medium.
As shown in Figure 7, N1 n-1 L11 does not implode until t = 1.75 ± 0.25 Myr and thus there should be no global implosion within the simulated time.However, if the implosion is indeed a local effect, the increased pressure within the cloud, should trigger a local implosion shortly after impact.
In Figure 13 we show slices of the radial mass flux, 1.5 Myr years after the explosion for the runs with and without a cloud.In the run without a cloud there is no implosion, while in the run with the cloud, there is a significant backflowing component behind the shock coming from the direction of the cloud.This confirms that the implosion can indeed be triggered locally.
We note, that the cloud is slowly expanding due to the pressure gradient relative to the ambient medium, leading to a radially inflowing component downstream of the shock.This flow is unrelated to the implosion, which is necessarily upstream.

Figure 2 .
Figure 2. Schematic overview of our proposed cloud formation mechanism.We show volume renderings of the radial mass flux in the model N1 n2 L14, during different phases of SNR late stage evolution.The physical scale differs between the frames.Opacity alpha is scaled with the logarithm of the density and is set to zero for densities in the range nH ∈ (95, 105) cm −3 in order to remove the foreground.The octant facing the observer has been made transparent.

Figure 3 .
Figure 3. Slices through the XY-plane of the model N1 n2 L14.The left and right columns display density and radial mass flux at different stages during the radiative phase of SNR evolution, respectively.The panels from top to bottom correspond to t = 0.006, 0.1, 1 and 10 Myr.Each panel showing density has a different color scale due to the large changes in dynamic range.The color scale for the density slices is logarithmic and asymmetrically centered around the ambient density nH, ISM = 100 cm −3 .The color scale for the radial mass flux is the same in all panels.

Figure 5 .
Figure 5.Time evolution of various quantities for model N1 n2 L13: (a) mass, (b) momentum, (c) energy, (d) pressure, (e) volume, (f) ejecta.Differently colored lines correspond to different gas components as described in Figure 1.The vertical dashed line marks the theoretical estimate of the shell formation time (Equation2).In panel (c) the thermal (solid) and kinetic (dashed) energy is shown.The horizontal green dashed line in panel (d) marks the equilibrium pressure of the ISM.The bubble pressure is initially very high (≳ 10 9 K cm −3 ) and thus outside of the axis limits.

Figure 6 .
Figure6.Rescaled version of Figure5for various models in media with different density.Solid, dashed and dotted lines correspond to the hot bubble, outflowing, and inflowing shell, respectively.The time is normalized by t n sf (see text).Quantities with subscript 'sf' correspond to the value of the whole SNR at t = t n sf .In panel (c) the thermal and kinetic energy are plotted with an opacity alpha of 1 and 0.5, respectively.

Figure 7 .
Figure 7. Launching timescale (top panels, a & b) and cloud formation timescale (bottom panels, c & d) as a function of explosion parameters.Left panels (a & c) show the timescales as a function of ISM density and right panels (b & d) as a function of explosion energy.Solid lines correspond to the model described in section 3.2.Data points are colored by the respective other explosion parameter.We added a small displacement (≤ 5 %) to the x-values in order to reduce the overlap of the markers.

Figure 8 .
Figure 8.Time evolution of the central cloud's radius (a), mass(b) and virial parameter (c) for the different models.Solid, dashed and dotted lines correspond to an explosion energy of E51 = 1.0,E51 = 5.0 and E51 = 14.0, respectively.Lines are colored by the ISM density.

Figure 11 .
Figure 11.Equilibrium pressure as a function of density for different cooling models.The brown curve corresponds to the default Ramses cooling, and the red and blue curves correspond to different models taken from Ploeckinger & Schaye (2020).

Table 1 .
Overview of the simulation suite