# Shocking interactions of supernova remnants with atomic and molecular clouds
- the interplay between shocks, thermal instability and gravity in the large
cloud regime

M. M. Kupilas,1 J. M. Pittard,1 C. J. Wareing,1 S. A. E. G. Falle2
1 School of Physics and Astronomy, University of Leeds, LS2 9JT, UK
2 Department of Applied Mathematics, University of Leeds, LS2 9JT, UK
<EMAIL_ADDRESS>
(Accepted 2022 April 04. Received 2022 March 17; in original form 2022 January
12)

###### Abstract

Using the adaptive mesh refinement code mg, we perform 3D hydrodynamic
simulations of a supernova-cloud interaction in the "large cloud regime". The
cloud is initially atomic and evolving due to the thermal instability (TI) and
gravity. We study interactions in a "pre-TI" and "post-TI" stage when cold and
dense clumps are present, and compare these results to idealised shock-cloud
scenarios in the "small cloud regime", and a scenario without shocks. On
aggregate, the supernova disruption is significantly weaker than that from an
idealised shock due to the supernova impact being instantaneous, and not
continuous. In both supernova-cloud interactions, we observe two shocks impact
the cloud, followed by the development of a weak 10 km s-1 upstream flow on
the cloud interface, and a global ambient pressure drop. When the cloud is
still atomic, it expands due to this drop. Additionally, the TI is triggered
at the front of the cloud, causing the formation of a cap-like structure with
clumps embedded inside. The upstream flow converges in this region, resulting
in a lobe-like cloud morphology. When the cloud is molecular, the transmitted
shock disrupts the inter-clump material and causes the clumps’ outer envelopes
to expand slightly and form tail-like morphologies. These effects are less
pronounced than those in our shock-cloud scenarios, and more pronounced that
those in our un-shocked scenario. After $\sim$ 3.5 Myrs, the effects from the
supernova decay and the cloud returns to an almost indistinguishable state
from an un-shocked cloud, in spite of the global ambient pressure drop. In
neither supernova-cloud scenario do we see any local gravitational collapse.

###### keywords:

Hydrodynamics – ISM: clouds – shock waves – supernova remnants – methods:
numerical

††pubyear: 2019††pagerange: Shocking interactions of supernova remnants with
atomic and molecular clouds - the interplay between shocks, thermal
instability and gravity in the large cloud regime–LABEL:lastpage

##  1 Introduction

The interstellar medium (ISM) is composed of multiple phases, ranging from a
hot ionized plasma (T $\gtrsim$ 106 K, n $\lesssim$ 0.01 cm-3) to a cold
molecular gas (T $\sim$ 10 – 20 K, n $\gtrsim$ 100 cm-3) (McKee & Ostriker,
1977; Wolfire et al., 1995, 2003). Stars form in molecular clouds, and thus to
convert ISM material into a star involves a web of processes occurring on time
and length scales spanning many orders of magnitude, ranging from galactic
scales down to the scale of the individual star. For a self-consistent theory
of star formation, details of the physics in the different phases across the
time and length scales must be connected and understood. This is far from
being the case, and many outstanding problems remain (e.g. see Krumholz,
2014).

In the context of molecular cloud formation, the interplay between the physics
of the thermal instability (TI, Parker, 1953; Field, 1965), gravity and
magnetic fields was studied in detail by Wareing et al. (2016a, hereafter
WPFVL16) in order to help elucidate the relative importance of these processes
at different stages. They performed idealised 2D and 3D simulations of an
initially quiescent cloud seeded with $\pm$10 per cent density perturbations
around an equilibrium state of n${}_{\textup{cl}}$ = 1.1 cm-3, in the warm,
neutral atomic unstable phase. They found that the TI alone was able to create
clumps, and filaments in the magnetic field case, whose properties connected
well with observations. The models were designed to incrementally study the
relative impact of different physics, and continued to increase in complexity
including feedback (Wareing et al., 2016b, 2017, 2018), and larger cloud
masses with higher resolution (Wareing & Falle, 2019, hereafter WFP19) and
magnetic fields (Wareing et al., 2021).

In Kupilas et al. (2021, Paper I) we looked at an extension of the 3D
hydrodynamical scenarios in WPFVL16 by including a shock. The simulations were
set up as a plane-parallel shock-cloud system commonly seen in the literature
(e.g. Stone & Norman, 1992; Klein et al., 1994; Nakamura et al., 2006; Van Loo
et al., 2007; Pittard et al., 2009; Pittard, 2011; Kinoshita et al., 2021;
Kanjilal et al., 2021) and the relative importance of the TI, gravity and the
shock impact was explored. Two scenarios were studied consisting of a shock-
cloud interaction with a "pre-TI" warm atomic cloud, and a "post-TI" molecular
cloud consisting of cold and dense (T $\sim$ 50 – 160 K, n $\sim\,$100 – 1000
cm-3) clumps embedded in a warm and diffuse (T $\sim$ 1000 K, n $\approx$ 0.8
cm-3) interclump gas.

In both cases, the constant impact due to the post-shock flow compressed the
cloud significantly and accelerated the global gravitational collapse of the
cloud. Evidence of local gravitational collapse was seen also, occurring on a
$\sim$ 5 Myr timescale in both scenarios.

Planar shocks like that in Paper I and the literature are usually
approximations of more complex flow patterns such as those resulting from
supernova remnants, a so-called "small-cloud approximation" (e.g. see Klein et
al., 1994). It is prudent therefore to investigate deviations from such an
approximation, and explore more realistic scenarios in a "large cloud regime",
where the idealised shock is replaced by a supernova explosion. Indeed, there
exists a history of studies with a vast range of initial conditions that
consider the impact of supernova explosions on molecular clouds and the wider
ISM (e.g. Wada et al., 2000; Joung & Mac Low, 2006; Rogers & Pittard, 2013,
2014; Padoan et al., 2016; Walch & Naab, 2015; Li et al., 2015; Seifried et
al., 2018; Lu et al., 2020). To our knowledge however, a systematic comparison
between the physical behaviour of an idealised shock-cloud and supernova-cloud
set-up has not yet been performed. This is the goal of the current paper, and
it is structured as follows. In Section 2 we present details of the model,
introduce the specific scenarios, define the small, medium and large cloud
regimes, and outline the computational details. In Section 3 we present our
results, followed by a brief discussion of physical implications in Section 4,
and we conclude our work in Section 5.

##  2 Methods

###  2.1 Numerical method

This work presents 3D hydrodynamical simulations of the interaction of a
supernova remnant with a cloud that is evolving due to the thermal instability
and gravity, and the interactions are compared to their analogous planar
shock-cloud and un-shocked scenarios studied in Paper I. All calculations were
performed using the finite-volume adaptive mesh refinement (AMR) code mg
(Falle, 1991, 2005; Hubber et al., 2013). The code solves the Euler equations
of hydrodynamics with free boundary conditions, is second order in space and
time, employs a hybrid cooling function and treats gravity by solving the
Poisson equation using the full-approximation multigrid technique. For the
reader interested in the details of the numerical scheme, they can be found in
Falle (1991), WPFVL16 and Paper I. We note however that in Paper I the
simulations were initially ran with the Godunov (GOD) Riemann solver (Godunov,
1959) until the gradients in the domain were too high for the simulations to
continue. At that point the solver was switched to Kurganov-Tadmor (KT)
(Kurganov & Tadmor, 2000). In the current paper, the KT solver is chosen from
the beginning.

###  2.2 Model

In this work we adopt the same initial model as scenario 3 in WPFVL16 and that
in Paper I. Namely, 17 000 M⊙ of diffuse material is initialised in a sphere
with radius r = 50 pc, which we refer to as the cloud (see figure 1 in Paper
I). The cloud lies at the origin of a Cartesian domain with a numerical extent
of -3 $<$ xyz $<$ 3\. In our model a unit length corresponds to 50 pc,
resulting in a domain with 300 pc extent on each side. 8 levels of AMR are
used, with the grid fixed to maximum resolution within a sphere with radius r
= 75 pc centered at the origin to accurately capture the behaviour of the
thermal instability (TI). The grid refines and de-refines dynamically outside
this region, with level 5 chosen as the minimum base level. For the supernova
scenarios, the region outside r = 75pc is then forced to be resolved down to
the 5th level to manage the computational cost associated with the grid
creating levels and cells as the remnant expands into regions that are of no
interest to the interaction. Note that for a fully refined grid, the maximum
number of grid cells would be 10243, however due to the AMR, $\lesssim$ 10 per
cent of this is used. The maximum resolution on the finest grid is $\Delta$x =
0.29 pc. A discussion of the dependence of the thermal instability on the grid
resolution and our arguments as to why it is acceptably resolved in our
simulations are presented in the appendix of Paper I.

The cloud is seeded with $\pm$ 10 per cent random density variations around
n${}_{\textup{cl}}$ = 1.1 cm-3 in the thermally unstable phase. The pressure
for this is then P/k = 4800 $\pm$ 300 K cm-3. We note that we assume a mean
particle mass of 2.0 $\times$ 10-24 g. Just like in Paper I, all material in
the domain is initially quiescent and the cloud is embedded in a medium with a
density contrast $\chi$ = 50 ($\chi$ =
n${}_{\textup{cl}}$/n${}_{\textup{amb}}$) setting the ambient medium
n${}_{\textup{amb}}$ = 0.022 cm-3. An advected scalar $\alpha$ is set to 1 in
the cloud and 0 everywhere else, and the energy source term is switched on
only in regions where $\alpha$ $>$ 0.9. The value of the scalar and all the
other fluid states transition monotonically between the cloud and ambient
values over an interface region with thickness of $\sim$ 5 cells. At the
moment of supernova/shock injection, the surroundings are remapped to have
n${}_{\textup{amb}}$ = 0.0022 and T${}_{\textup{amb}}$ = 2.19 $\times$ 106 K
($\chi$ = 500) so that the energy source term can be switched on in the whole
domain whilst the surroundings remain adiabatic due to the large cooling times
of the order t${}_{\textup{cool}}$ $\gtrsim$ 100 Myr at such temperatures.

####  2.2.1 Scenarios

Five scenarios are considered in this paper: NoShock, S1, S2, SN1, SN2. The
first three, NoShock, S1 and S2 are effectively the cases presented in Paper
I. For self-consistency, the data from Paper I was not used and instead the
models were re-simulated using the KT solver for accurate comparisons to the
new SN1 and SN2 scenarios. Just like in Paper I, the NoShock scenario follows
the evolution of the initial cloud for a free-fall timescale of t $\approx$ 35
Myr without any additional disturbances. S1 is then the 12Shock analog, where
a planar shock is introduced at 11.78 Myr into the evolution of NoShock. SN1
then corresponds directly to S1, with the idealised planar shock being
replaced by a supernova explosion introduced at the exact same time. S2 is
then the 24Shock analog, however this time a shock is introduced at 26.5 Myr
into the evolution of NoShock instead of 24 Myr. The reason for this later
injection time is due to the KT solver being more diffusive than the GOD
solver, causing the thermal instability to take longer to trigger in the re-
simulated NoShock scenario. This causes the phase transition to occur $\sim$ 2
Myr later than in Paper I, and as such we choose a later time for
shock/supernova injection in the S2/SN2 scenarios.

Figure 1: Shown are snapshots of the z = 0 plane of the density logarithm of
(a): The cloud at injection time t${}_{\textup{inj}}$ = 11.78 Myr for SN1/S1
scenarios and (b): Cloud at injection time t${}_{\textup{inj}}$ = 26.5 Myr for
SN2/S2 scenarios. The location of the explosion can be seen in both snapshots.
Table 1: Summary of all the scenarios studied in this work. The columns show
the scenario name, the source of impact onto the cloud, the distance from the
blast origin to the nearest cloud edge R${}_{\textup{c}}$, the Mach number (M)
of the disturbance shock immediately upon impact, and the injection time
(t${}_{\textup{inj}}$) of the disturbance. Scenarios NoShock, S1 and S2 are
effectively NoShock, 12Shock and 24Shock from Paper I which have been re-
simulated for this study for self-consistent comparisons. Scenarios SN1 and
SN2 introduce supernova explosions at the same time as the shock in S1/S2.

Scenario |  Disturbance |  R${}_{\textup{c}}$ [pc] |  M (impact) |  t${}_{\textup{inj}}$ [Myr]
---|---|---|---|---
NoShock |  None |  N/A |  N/A |  N/A
S1 |  Shock |  $\infty$ |  1.5 |  11.78
SN1 |  Supernova |  50 |  $\approx$ 7* |  11.78
S2 |  Shock |  $\infty$ |  1.5 |  26.5
SN2 |  Supernova |  70 |  $\approx$ 7* |  26.5

  * *

Due to a 10 M⊙ and 1051 erg explosion located at the polar coordinate (r,
$\theta$) = (100 pc, 160∘) in the z = 0 plane. $\theta$ is measured from the
positive x axis.

For our idealised shock-cloud models we follow an assumption known as the
small cloud approximation (Klein et al., 1994, hereafter KMC94), a condition
that enables a physical set-up where an incoming shock onto a cloud has no
curvature in its structure, and the shock-driving pressure is time
independent. Implications of such an approximation can be understood in terms
of the time-dependence of the shock-driving pressure, discussed in KMC94, and
formally analysed in McKee et al. (1987, hereafter MHST87). In their theory of
weakly time-dependent interstellar shocks, MHST87 introduce the pressure
variation timescale t${}_{\textup{p}}$, and find that it approximates well as

$t_{\textup{p}}\simeq 0.1\frac{R_{\textup{c}}}{v_{\textup{b}}},$ (1)

where R${}_{\textup{c}}$ is the distance of the nearest part of the surface of
the cloud from the blast epicentre, and v${}_{\textup{b}}$ =
dR${}_{\textup{b}}$/dt is the velocity of the blast wave.

Two other timescales of interest are the shock-crossing timescale

$t_{\textup{sc}}=\frac{2r_{\textup{cl}}}{v_{\textup{b}}},$ (2)

i.e. the time taken for the external blast wave shock to sweep over the cloud,
and the cloud-crushing timescale

$t_{\textup{cc}}=\frac{\chi^{1/2}r_{\textup{cl}}}{v_{\textup{b}}},$ (3)

i.e. the time taken for the transmitted cloud shock with velocity
v${}_{\textup{b,cl}}$ = v${}_{\textup{b}}$/$\chi^{1/2}$ to traverse the cloud
radius r${}_{\textup{cl}}$ and crush the cloud.

Comparison of these timescales allows the definition of three regimes (see
e.g. KMC94 and references therein). The small cloud regime requires the cloud
to be sufficiently small (or blast to be sufficiently large) for
t${}_{\textup{p}}$ to satisfy $t_{\textup{cc}}\ll t_{\textup{p}}$, thus by
extension for r${}_{\textup{cl}}$ to satisfy

$r_{\mathrm{cl}}\ll 0.1\frac{R_{\mathrm{c}}}{\chi^{1/2}}.$ (4)

The medium cloud regime requires the cloud to be a size such that the blast
wave does not change significantly as it sweeps over the cloud, but does
change as the cloud is crushed. This requires t${}_{\textup{p}}$ to satisfy
$t_{\textup{cc}}\gtrsim t_{\textup{p}}\gtrsim t_{\textup{sc}}$, and by
extension for r${}_{\textup{cl}}$ to satisfy

$0.1\frac{R_{\textup{c}}}{\chi^{1/2}}\lesssim
r_{\textup{cl}}\,\lesssim\,0.05R_{\textup{c}}.$ (5)

Finally, the large cloud regime requires that the cloud is large enough (or
blast small enough) that the blast wave ages significantly as it sweeps over
the cloud, resulting in vastly weaker compression at the rear of the cloud
than the front, such that t${}_{\textup{p}}$ satisfies
$t_{\textup{sc}}\,>\,t_{\textup{p}}$ and r${}_{\textup{cl}}$ satisfies

$r_{\textup{cl}}>0.05R_{\textup{c}}.$ (6)

For the S1/S2 scenarios, an M = 1.5 shock is artificially imposed on the grid
and boundary cells lying within the post-shock region are fixed to the post-
shock values. The shock front is located at the polar coordinate of (r,
$\theta$) = (50 pc, 160∘) in the z = 0 plane, placing the shock immediately on
the cloud’s edge. The angle is measured counterclockwise from the positive
x-axis, and is chosen to minimise the effects of the Quirk instability, or the
so-called carbuncle phenomenon, that can cause perturbations along the shock
symmetry axis, with strongest perturbation occurring when aligned with the
grid (Quirk, 1997; Elling, 2009).

For our choice of parameters, S1/S2 scenarios can be considered to have
R${}_{\textup{c}}\,=\,\infty$, placing them in the small cloud regime as
initially assummed.

For the SN1/SN2 scenarios, a supernova is injected at the polar coordinate of
(r, $\theta$) = (100 pc, 160∘) in the z = 0 plane on the grid located such
that the first impact on the cloud is felt at the same location as the first
shock impact in S1/S2. The injection volume is artificially refined to the
highest AMR level and subsequently chosen to fill a region with a radius of 5
cells, i.e. 1.5 pc, into which 10 M⊙ of mass and 1051 erg of thermal energy
are injected over a 500 yr period, roughly the time it would take for a
remnant to expand to the injection volume. Note that immediately upon impact,
the forward shock of the supernova remnant has M $\approx$ 7\. The clouds at
11.78 Myr and 26.5 Myr are shown in figure 1 where the region of supernova
injection can also be seen.

For SN1, R${}_{\textup{c}}$ $\approx$ 50 pc which gives 0.05R${}_{\textup{c}}$
= 2.5 pc. As in this scenario r${}_{\textup{cl}}$ = 50 pc, it places this
model firmly in the large cloud regime. For SN2, R${}_{\textup{c}}$ $\approx$
70 pc as the evolution results in the cloud shrinking slightly due to the
formation of clumps and gravitational contraction, reducing
r${}_{\textup{cl}}$ to a size of 30 pc. Nevertheless, 0.05R${}_{\textup{c}}$
evaluates to 3.5 pc, also placing this model firmly in the large cloud regime.
Thus our new models showcase interactions in two contrasting, and very
different regimes, with the large cloud regime showing a much more realistic
interaction. A summary of the model parameters is presented in Table 1.

Figure 2: Shown are pressure logarithm plots in the z = 0 plane of the SN1
scenario. The cloud edge cannot be seen initially in panel (a) and is thus
marked by the white dashed line. The time in each panel corresponds to time
since supernova injection t \- t${}_{\textup{inj}}$, with t${}_{\textup{inj}}$
= 11.78 Myr. The final snapshot shows the moment immediately prior to the
primary shock moving off the grid, and at the time the upstream flow around
the cloud is being established.

##  3 Results

In this section we present our results. First, we present the SN1 scenario and
discuss the evolution outside the cloud followed by a discussion of the
internal evolution. Secondly we present the SN2 scenario by again looking at
the external and internal evolution separately.

###  3.1 Interaction with atomic cloud - SN1

Throughout the discussion, comparisons are made to the planar shock scenario
S1, whose detailed evolution is presented in Paper I. In addition, where
appropriate, the NoShock scenario is referenced. The 5.16 Myr timescale, where
we choose to stop our analysis, corresponds to the time of the first S1
snapshot that showed evidence of local gravitational collapse.

####  3.1.1 Dynamics outside the cloud

To illustrate the dynamics outside the cloud we first show in figure 2 the
pressure profile of the domain in the z = 0 plane. The first 0.74 Myr are
shown, as the most interesting features in the ambient medium are seen during
this timescale. Figure 2(a) shows a snapshot immediately post supernova
injection. A spherical region with a 50 pc radius centered at the origin
contains 17 000 M⊙ of material, whose density is a factor $\chi$ = 500 greater
than the surroundings with n${}_{\textup{amb}}$ = 0.0022 cm-3. A region with a
radius r = 1.5 pc contains 10 M⊙ and 1051 erg of thermal energy, and as such
is strongly over-pressured with respect to its surroundings, with internal
pressures exceeding 1010 K cm-3, 7 orders of magnitude larger than the ambient
pressure of approximately 4800 K cm-3. A powerful explosion thus follows,
accelerating stationary material to speeds exceeding 104 km s-1.

Before describing the interaction, we highlight a number of noticeable
numerical effects. Firstly, the cloud in Fig. 2(a) has a radius r = 50 pc,
however only a smaller spherical region with a radius r $\approx$ 25 pc can be
seen. This corresponds to where the internal pressure has dropped slightly due
to the early action of the thermal instability. This behaviour is expected and
acceptable, however we would expect to see the whole cloud behave this way.
Instead what we see is a numerical effect on the boundary resulting from how
the source term was initially set up (i.e. switched off for regions where the
scalar $\alpha$ $<$ 0.9, with $\alpha$ = 1 inside the cloud and $\alpha$ = 0
outside). Material below $\alpha$ $<$ 0.9 does not cool, and as such a
gradient develops starting at the interface (where $\alpha$ $=$ 0.9) from the
value of the ambient pressure (P/k $\approx$ 4800 K cm-3) to the lowest value
of the cloud centre (P/k $\approx$ 4750 K cm-3), and accounts for why the full
cloud is not visible in panel (a). Where the edge actually is has been marked
with a dashed white line. Whilst this effect is noticeable, it accounts for a
maximum of 1 per cent fractional difference in the fluid variables, and does
not affect the dynamics. Additionally, an outwardly propagating expansion wave
can also be seen at the far side of the cloud, most clearly in panels (b) and
(c). This results from the mapping of the ambient medium from $\chi$ = 50 to
$\chi$ = 500, and amounts to a fractional difference in the fluid state of
$\sim$ 5 per cent.

Again, this has little effect on the dynamics especially since these
differences become overwhelmed by the supernova impact shortly thereafter.

The impact with the cloud occurs $\sim$ 40 kyrs after explosion. At this
point, $\sim$ 33 M⊙ of material has been swept up, approximately 3 $\times$
the mass of the ejecta, and as such the profile immediately prior to impact is
Sedov-Taylor-like. The shock has a speed $v_{\textup{s}}$ = 1087 km s-1 upon
impact, and given a pre-shock sound speed c = 160 km s-1, this amounts to an
impact shock with Mach number M = $v_{\textup{s}}/c$ $\approx$ 7, roughly a
factor $\sim$ 4.5 $\times$ stronger than the M = 1.5 planar shock in S1. A
reverse shock from the blast moves towards the epicentre with
v${}_{\textup{rs}}$ = 450 km s-1, and a contact discontinuity separates the
forward and reverse shocks. Note that a purely analytical Sedov-Taylor profile
does not have these two features - it is assummed that all waves have settled
into a steady state such that a self-similar analytical solution is
appropriate.

In our simulations the profile behind the forward shock matches well with the
analytical Sedov-Taylor profile.

The profile interacts with the cloud as the remnant passes over it,
transmitting a shock into the cloud and reflecting a shock back into the
ambient medium. The reflected shock propagates towards the blast epicentre,
and leaves behind a low pressure void due to a trailing expansion. Meanwhile,
the reverse shock from the initial blast reflects at the epicentre, and as it
propagates back outwards, interacts with the cloud-reflected shock, generating
an additional set of disturbances. This shock subsequently approaches the
cloud and causes a secondary impact. We refer to the cloud edge where this
impact occurs as the front of the cloud, the first shock as the primary shock
and the second shock as the secondary shock.

Fig. 2(b) shows a snapshot at 92 kyrs. The primary shock has passed almost
half of the cloud, and a clear distortion in its interior structure is seen in
the region that interacts with the reflected shock. The primary shock
diffracts around the cloud which further distorts its structure. A low
pressure region ahead of the cloud can be seen, which forms due to the
expansion of gas behind the reflected primary shock. The gas at the front of
the cloud continues to expand during the secondary shock approach, and as such
there is a strong and sharp pressure contrast between this region and its
surroundings - as evident in Fig. 2(c) which shows the moment preceding the
secondary impact. By Fig. 2(d) the secondary impact has occurred, resulting in
turbulent eddies forming at the front of the cloud and a set of reflected
waves propagating towards the left x boundary. The level of which this may
seed turbulence is left to a future work. Both shocks can now be seen sweeping
over the cloud too, and the primary shock has nearly converged at the rear.

By Fig. 2(e), the primary shock has converged and the secondary shock is
trailing close behind. As the primary shock converges, an interaction region
is generated with portions of the resulting waves propagating both downstream
and upstream, with the upstream waves interacting with the secondary shock as
it sweeps over the cloud.

By Fig. 2(f) both shocks have passed over the cloud, and are now at the domain
edge. The secondary shock has now also converged at the rear, resulting in a
secondary set of waves with both downstream and upstream components. Around
the edges of the cloud the upstream components generate a flow which then
itself converges at the front, exaggerating the compression from the initial
impact.

After Fig. 2(f), the ambient medium evolves due to existing pressure gradients
and the remaining flows on the grid. The ambient pressure in panel (f) can
already be seen to be lower than most of the cloud’s edge, and this continues
to drop as the model evolves. The portions of the cloud that can respond to
these changes, i.e. those that have been shocked, gradually expand in an
attempt to balance the ambient pressures and the post-shock pressures in the
cloud. This is halted however at the front of the cloud due to the ram
pressure impact of the flow that converges there. Thus the outer envelope of
the cloud expands everywhere but the front, resulting in a lobe-like
morphology with the front being pinched inwards.

As early as the first 1 Myr of evolution, considerable differences between the
SN1 and S1 scenario can be seen. Two main ones can be identified however which
exaggerate as the simulations evolve, and account for all of the differences
seen between our shock-cloud and supernova-cloud scenarios. Firstly, for an
idealised, infinitely extending planar shock, the flow behind it is
continuously replenished, resulting in disruption of the cloud throughout its
entire evolution. Secondly, after shock passage, the ambient pressure is
relatively time-independent, i.e after the medium is shocked and its pressures
exceed that of the cloud, it stays that way, and provides continuous
compression. These two features are present in all adiabatic, planar shock-
cloud systems, and thus a more realistic model like ours can usefully be
contrasted against them.

Figure 3: Shown are slices in the z = 0 plane showing the velocity (panels a
and c) and pressure (panels b and d) of SN1 (panels a and b) and S1 (panels c
and d) at t = 5.16 Myrs - the final snapshot presented in this work and the
timescale at which the S1 scenario witnessed local gravitational collapse.

The impact of this in our models can clearly be seen in the last snapshot at
5.16 Myrs in figure 3 which shows the pressure and velocity field in the z = 0
plane. In Fig. 3b the ambient thermal pressure can be seen to have decreased
down to $\sim$ 4000 K cm-3, approximately 80 % of its original value of 4800 K
cm-3. The final S1 pressure in comparison is $\sim$ 104 K cm-3, a factor of
2.5 larger than the final SN1 pressure, and has not dropped since the initial
shock passage. The velocity difference is even larger, with impact amplitudes
on the order of 10 km s-1 in SN1 (Fig. 3a) and 100 km s-1 in S1 (Fig. 3c),
resulting in ram pressure impact that is $\sim$ 100 times greater
($\rho|\textbf{u}|^{2}|_{\textit{\tiny{S1}}}/\rho|\textbf{u}|^{2}|_{\textit{\tiny{SN1}}}$,
$|\textbf{u}|^{2}$ = $u_{x}^{2}+u_{y}^{2}+u_{z}^{2}$) in S1 than in SN1. Thus
in aggregate, the supernova has a profoundly weaker compressive effect on the
cloud, and is mainly responsible for creating external conditions that couple
strongly to the cloud behaviour and result in cloud expansion.

We note that it is not entirely clear to what extent boundary effects impact
the ambient evolution. Many waves are expected due to back and forth
reflections within the remnant (e.g. Cioffi et al., 1988). As seen in Fig. 2,
the shocks move off the grid soon after 0.75 Myrs and so these reflections
cannot occur. We would not however expect the fluid state to be significantly
altered by their presence, and 1D spherically symmetric tests of supernova
explosions with equivalent parameters as our models suggest this to be true.
In fact, in our 1D tests, once the reverse shock reflects at the blast
epicentre, it is unable to catch up to the forward shock to reflect backwards
once more. An additional concern is the 26 km s-1 flows near the left x
boundary. Boundary effects are likely responsible for these flows as they do
not align with the spherical symmetry of the initial blast. However as they
are not near the cloud and flowing off the grid, we do not expect them to
affect the cloud. Finally, the internal pressure of any supernova remnant
drops as it expands. As we see such a pressure drop occur in our models, we
are confident that our simulations are representative of what would happen if
the whole remnant was captured until the latest times considered.

####  3.1.2 Dynamics inside the cloud

Figure 4: Panels (a – h) show density slices in the z = 0 plane for the SN1
scenario. The logarithm of the number density is shown, and the time on each
panel corresponds to the time elapsed since supernova injection t \-
t${}_{\textup{inj}}$, where t${}_{\textup{inj}}$ = 11.78 Myr. Panels (i), (j)
and (k) respectively show the maximum density, phase fractions and energy
fractions in the cloud, showing data for SN1 (red lines), S1 (blue lines) and
NoShock (black lines) for direct comparison at equivalent times. Panel (j)
shows mass fractions in cold (T $<$ 160 K, solid line), unstable (160 $<$ T
$<$ 5000 K, dashed line) and warm (5000 $<$ T $<$ 10 000 K, dot-dashed line)
regimes. Panel (k) shows the fractions of the thermal energy
E${}_{\textup{thermal}}$ (solid line) and kinetic energy
E${}_{\textup{kinetic}}$ (dashed line) out of E${}_{\textup{tot}}$ =
E${}_{\textup{thermal}}$ \+ E${}_{\textup{kinetic}}$ in the cloud.
Gravitational energy is ignored.

To describe the evolution inside the cloud, we first show snapshots of the
logarithm of the density through the z = 0 plane in figure 4. Additionally, in
Fig. 4(i – k) we respectively compare the maximum density, fractions of
material existing in warm (10 000 $>$ T $>$ 5000 K), unstable (5000 $>$ T $>$
160 K) and cold (T $<$ 160 K) thermal regimes, and the energy partition
between thermal and kinetic energies. Comparisons are drawn between the
NoShock (black line), SN1 (red line) and the S1 scenarios. Here, only cloud
material is traced, i.e. where the passively advected scalar $\alpha$ $>$ 0.9.
Note that in our discussions, material lying on the equilibrium curve is
referred to as belonging to a phase (e.g. thermally unstable phase). When
talking about material within a temperature bracket (i.e. all material,
equilibrium + non-equilibrium), we refer to this as belonging to that
particular regime (e.g. thermally unstable regime).

The impact of the primary shock is felt strongest on the cloud front. As the
cloud and ambient medium are initially in pressure balance, the transmitted
shock has roughly equivalent strength in the cloud as it does in the
surroundings, resulting in an M $\approx$ 7 shock (M = vs/c where
v${}_{\textup{s}}$ $\approx$ 50 km s-1 and c $\approx$ 7 km s-1) travelling
into the cloud and thermalising a thin layer of material - as seen in Fig.
4(a). Just like for a planar shock, the cloud causes the supernova shock to
diffract and reduce its strength, thus reducing the strength of the
transmitted shock. Its strength is weakened further as the remnant expansion
itself slows down, which is not an effect that occurs with a planar shock due
to a continuously driven flow. Thus when the primary shock reconnects at the
rear, much of its strength has been lost, with the convergence having little
effect on amplifying the shock impact. This is in stark contrast to the S1
scenario, where the strongest compression was seen at the back of the cloud
due to the converging flow.

By Fig. 4(b), the shock has travelled over the whole cloud and the transmitted
shock is moving inwards from all sides. Where the impact was the strongest,
the thermalised state experiences an increase in cooling, causing the material
to be compressed by the surroundings as its temperature and pressure drops.
This effect is only evident at the front of the cloud, resulting in the
formation of the single cap-like structure that after 2 – 3 Myrs appears
filament-like in the density slices shown. The shape of this structure is
distorted by the converging flow, and the pinching of the front causes the
structure to accelerate faster at the centre than the sides. This is an effect
more vividly seen after $\sim$ 3 Myrs (Fig. 4e) and is exaggerated by the
expansion of the surrounding envelope. The remaining evolution is gradual,
where the cloud continues to expand and the cold structure continues to be
distorted by the converging flow. The shock inside the cloud continues to slow
down, and by Fig. 4(h) it has still not converged at the centre, in contrast
to S1 which saw the cloud shock converge on this timescale.

Many significant differences between the S1 and SN1 scenarios can be
identified and can all be explained by this fact that the continuously
replenished flow in S1 provides constant ram pressure impact and constant
thermal compression, and the remnant does not. Because of this, in S1 the
period between 1 – 4 Myrs saw the formation of a cold thin shell spanning most
of the cloud’s edge which was continuously impacted by the incident flow. This
flow distorted the shell, triggering Rayleigh-Taylor and Vishniac
instabilities at the front, with the sides being distorted by the Kelvin-
Helmholtz instability. Where the flow converged at the rear, the shell became
winged-like, with a protruding needle forming in the centre that was
subsequently accelerated against the motion of the cloud. The surrounding flow
maintained high thermal pressures, and in combination with the strong ram
pressure, the cloud was continuously compressed to create conditions where
gravity could take over, and collapse individual structures. Whereas here,
instead of a cold thin shell forming around the edges, the only cold region is
found in the small layer at the front. This forms due to an instantaneous
impact with no subsequent continuous flow. As such, the ram pressure is
smaller by a factor of $\sim$ 100 and thermal pressure by a factor of $\sim$
2.5 (see Section 3.1.1). None of the hydrodynamical instabilities present in
S1 are triggered, and the continuous drop in pressure means that the cloud
expands rather than contracts. The only other source of compression - the
secondary shock, has no effect on the material in the cold layer and only
somewhat affects the warm material at the interface. We do not expect to see
star formation take place on the same timescales as S1. The stronger and more
instantaneous impact of the supernova can be seen in the maximum density (Fig.
4i), subsequently followed by a rapid decline in the next 200 kyrs. The phase
and energy fractions (Fig. 4j,k) also show this behaviour, with the kinetic
energy specifically showing a rapid decay akin to that observed in the studies
of Seifried et al. (2018). S1 overshoots SN1 in the kinetic energy after
$\sim$ 300 kyrs, maximum density after $\sim$ 1 Myr, and cold mass fraction
after $\sim$ 2 Myr. The sustained compression and acceleration of material in
the S1 scenario subsequently results in the steady increase of these
properties, with the final maximum density growing to values of $\sim$ 104
cm-3, 2 orders of magnitude higher than SN1. At n = 104 cm-3, gravity is
dominating, and since 70 per cent of the S1 material is in the cold regime
(Fig. 4j), the probability of star formation is high. Indeed our results in
Paper I showed evidence of local gravitational collapse. At a maximum n = 100
cm-3 however and with less than 5 per cent of material in the cold regime,
clearly in SN1 the conditions for self-gravity to dominate and cause local
gravitational collapse are not created. We note however that the cold mass
fraction does appear to increase in Fig. 4(j) (although the maximum density
decreases) suggesting that the structure is not in the final stage of
formation.

####  3.1.3 Formation of cold material

Figure 5: Panels (a – l) show the mass distribution of cloud material
represented in pressure-density space for the SN1 scenario. Isotherms
differentiating the hot (T $>$ 10 000 K), warm (10 000 $>$ T $>$ 5000 K),
unstable (5000 $>$ T $>$ 160 K) and cold (T $<$ 160 K) thermal regimes are
shown on each panel. Presented is a period of 2.95 Myrs from
t${}_{\textup{inj}}$ = 11.78 Myrs, after which material settles into a
relatively steady state where the cold material evolves primarily due to
gravity, and the warm envelope expands into the surroundings.

Figure 5 shows the mass distribution of cloud material in pressure-density
space, along with the equilibrium curve and isotherms for the hot (T $>$ 10
000 K), warm (10 000 $>$ T $>$ 5000 K), unstable (5000 $>$ T $>$ 160 K) and
cold (T $<$ 160 K) thermal regimes. Note that the data shown is only for
material where the advected scalar $\alpha$ $>$ 0.9, i.e. cloud material.

Panel (a) shows the moment before any impact. Note that the remapped material
exists where $\alpha$ $<$ 0.9, and as such it is not seen in the panel. Panel
(b) clearly shows the compression of cloud gas due to the transmitted shock,
and the heating of unstable material to the warm regime. The low-pressure,
low-density gas that results from the reflection of the primary shock has now
mixed with the cloud and is traced by the scalar $\alpha$, and thus can be
seen to exist in the hot regime. In panel (c), a clear rise in pressure and
density can be seen with this material, which is an effect that can be
credited to the impact of the secondary shock. In the warm regime, a drop in
density and a broader pressure distribution can also be seen as material
begins to cool towards thermal equilibrium, a pattern that continues in panels
(d) and (e). What is most striking is the behaviour of the gas over the next
1.5 Myr. In the S1 scenario, post-shock gas is seen to cool from the warm
regime directly to the cold phase, completely passing the unstable phase. Here
we instead see material cool from the warm regime back to the unstable phase,
and then migrate towards the cold phase from the unstable phase - in a fashion
that resembles action of the thermal instability. This is exactly what was
seen in the NoShock scenario presented in Paper I (in particular their figs.
4(a – e)). In NoShock, the thermal instability formed cold and dense clumps
embedded in warm intercloud gas. We would therefore expect to see such
structures form here, which indeed appears to be the case. What looks like a
filament in the density slices shown in figure 4(c – h) is in fact a cross-
section through a broader cap-like structure, with what looks like clumps
embedded in a diffuse warm envelope, features better seen in projection.

Figure 6: Shown are logarithm plots of column density projections along two
lines of sight for the SN1 scenario. The dashed line in panel (a) corresponds
to the line of sight for panel (b), and vice versa in panel (b). The length of
each line corresponds to the projection depth. Note that the axes in the image
do not correspond to the original axes of the simulation, and instead form the
axes of the projected image centered on the centre of mass of the cold
structure (T $<$ 160 K).

In Fig. 6 we show column density projections at t = 5.16 Myr. Two lines-of-
sight (LOS) are considered, with the image in panel (a) being a projection
perpendicular to the z = 0 plane (LOS shown in panel b), and (b) being along
the axis of symmetry as determined by the blast (LOS shown in panel a). Here,
the length of the lines represent the projection depth, such that only
material within the volume as defined by the depth is projected. The axes of
the resultant projections are subsequently shown in the centre of mass frame
of the cold material.

In a projection perpendicular to the x – y plane (Fig. 6a) the singular
filamentary-like identity seen in Fig. 4 is mostly lost, and the projection
reflects a morphology consistent with a conic, cap-like structure that has
been shaped by the convergent flow. In the projection shown in Fig. 6(b) we
see that the filamentary identity is completely lost and within the surface
there appears to be a small complex of clumps separated by distances typical
of structures created by the thermal instability ($\sim$ 3-5 pc, see WPFVL16
and Paper I). The clumps themselves have average densities n $\approx$ 80
cm-3, temperatures T $\approx$ 30 – 50 K and sizes r $\approx$ 1 – 2 pc. We
thus see that the supernova shock has triggered accelerated thermal evolution
in the localised area of the ‘front’ of the cloud, closest to the supernova,
and on realistic timescales formed the TI-driven structure expected from such
a phase transition (e.g. see Falle et al., 2020). As mentioned, the cold mass
fraction is gradually increasing, suggesting that these clumps are not in
their final stages of formation, however no star formation is observed.

###  3.2 Interaction with molecular cloud - SN2

Just like SN1, we present the SN2 scenario by focussing first on the dynamics
outside the cloud and secondly on the dynamics inside the cloud. Throughout
our discussion we compare and contrast this behaviour to the idealised planar
shock scenario S2, the un-shocked scenario NoShock, and the previously
presented scenario SN1.

####  3.2.1 Dynamics outside the cloud

We illustrate the external behaviour by again showing snapshots of the
pressure logarithm in the z = 0 plane (figure 7) and a snapshot at 5.16 Myrs
showing pressure and velocity profiles in the same plane for both SN2 and S2
(figure 8). The full evolution is presented on the same timescale of 5.16 Myrs
which just like for S1 is the timescale of the first snapshot in S2 that
showed evidence of local gravitational collapse.

Figure 7: Shown are pressure logarithm plots in the z = 0 plane of the SN2
scenario. The time in each panel corresponds to the time since supernova
injection t-t${}_{\textup{inj}}$, with t${}_{\textup{inj}}$ = 26.5 Myr. Just
like the SN1 scenario, the final snapshot shows the moment immediately prior
to the primary shock moving off the grid and the development of the upstream
flow.

Fig. 7(a) shows the moment once all of the 10 M⊙ and 1051 erg are injected
into the domain. A supernova explosion follows, and the remnant propagates out
into the surroundings sweeping up mass and thermalising material as it does
so. Just like in SN1, a Sedov-Taylor-like profile develops with a forward
shock, reverse shock, and a contact discontinuity separating the two. The
reverse shock then reflects at the blast epicentre, sending a secondary shock
propagating back towards the cloud.

The subsequent evolution is very similar to SN1, however subtle differences
are worth pointing out. Firstly, the cloud no longer has a smooth edge,
resulting in multiple wave reflections that distort both the primary and
secondary shocks. As the primary shock sweeps over the cloud and interacts
with the inhomogeneities, sets of oblique shocks form which distort the
profile further. These are seen to persist for up to 180 kyrs (Fig. 7c), and
are destroyed after the passage of the secondary shock. Secondly, the cloud is
smaller than in SN1, resulting in shorter timescales for the passage of the
shocks. As can be seen in Fig. 7(d) - the primary shock has converged, and the
panel captures the moment prior to secondary shock convergence. Contrasting
this to SN1 (Fig. 2d), the primary shock has yet to converge at this time. In
spite of these differences, in the next $\sim$ 4 Myrs the ambient medium
behaves almost identically to the SN1 scenario: a low velocity, laminar
upstream flow develops that converges at the front of the cloud, and the
thermal pressure drops globally.

Figure 8: Shown are plots in the z = 0 plane showing the velocity (panels a
and c) and pressure (panels b and d) of SN2 (panels a and b) and S2 (panels c
and d) at t \- t${}_{\textup{inj}}$ = 5.16 Myrs - the final snapshot presented
in this work and the timescale at which the S2 scenario witnessed local
gravitational collapse.

The final state of this behaviour is seen in figure 8, showing the same order
of magnitude differences between the velocity and thermal pressures, giving a
thermal pressure difference of a factor $\sim$ 2.5 lower in SN2 than S2, and
ram pressure difference a factor $\sim$ 100 lower (on the cloud edge). The
most surprising behaviour in this scenario is the state of the SN2 cloud
itself, showing that although the cloud becomes immersed in under-pressured
surroundings just like in SN1, the cloud has not expanded. The cloud in fact
continues to contract (albeit not as fast as S2), which is clearly a
gravitational effect. In this scenario therefore, the ambient behaviour is
only weakly coupled to the cloud.

####  3.2.2 Dynamics inside the cloud

Figure 9: Panels (a – h) show density slices in the z = 0 plane for the SN2
scenario. The logarithm of the number density is shown, and the time on each
panel corresponds to the time elapsed since supernova injection t \-
t${}_{\textup{inj}}$, where t${}_{\textup{inj}}$ = 26.5 Myr. Panels (i), (j)
and (k) respectively show the maximum density, phase fractions and energy
fractions in the cloud, showing data for SN2 (red lines), S2 (blue lines) and
NoShock (black lines) for direct comparison at equivalent times. Panel (j)
shows mass fractions in cold (T $<$ 160 K, solid line), unstable (160 $<$ T
$<$ 5000 K, dashed line) and warm (5000 $<$ T $<$ 10 000 K, dot-dashed line)
regimes. Panel (k) shows the fractions of the thermal energy
E${}_{\textup{thermal}}$ (solid line) and kinetic energy
E${}_{\textup{kinetic}}$ (dashed line) out of E${}_{\textup{tot}}$ =
E${}_{\textup{thermal}}$ \+ E${}_{\textup{kinetic}}$ in the cloud.
Gravitational energy is ignored.

We now discuss the dynamics inside the cloud and refer to the density
logarithm snapshots and statistics in figure 9. Again we show the maximum
density (i), mass fractions in the cold, unstable and warm thermal regimes
(j), and the fractions of thermal and kinetic energies (k).

This scenario begins at 26.5 Myr. The evolution of the material inside the
cloud has so far been dominated by the effects of the thermal instability
which has generated a complex of cold and dense clumps embedded in a warm
diffuse gas. The evolution of the thermal instability in the NoShock case was
covered in section 3.1 of Paper I. To summarise, the initial grid-scale
inhomogeneities initialised around the thermally unstable equilibrium smooth
out and seed the instability within the first few thousand years of evolution.
The cloud remains quiescent for $\sim$ 15 Myr, corresponding to a period that
sees a growth of the inhomogeneities, enhancing the pressure and density
differences within the cloud. A critical point is reached within the next
$\sim$ 5 Myrs where the cooling rates increase such that cold material cools
further, causing it to be further compressed by the over-pressured
surroundings. This continues to cool until temperatures reach a stable cold
equilibrium, resulting in densities increasing from n $\sim$ 1 cm-3 to n $>$
100 cm-3 in a period of $\sim$ 3 Myr. Gravity subsequently increases these
densities to $\sim$ 1000 cm-3 prior to the global collapse of the cloud at
$\sim$ 35 Myr.

Following phase transition and prior to the global collapse, the state of the
cloud is as a complex of cold and dense (T $\sim$ 50 – 160 K, n $\sim$ 100 –
1000 cm-3) clumps embedded in a warm and diffuse (T $\approx$ 5000 K, n
$\approx$ 0.8 cm-3) gas. Low-density, diffuse (T $\approx$ 700 K, n $\approx$
5 cm-3) structures resulting from thermal flows connect the clumps, and the
clumps appear somewhat elongated towards one another, serving as evidence that
these clump-connecting structures are possible regions of filament formation
(e.g. seen in the high resolution hydrodynamical simulations of WFP19). As
mentioned, the first detailed study of this behaviour was conducted by WPFVL16
and for extra details and comparisons between non-gravitational simulations
and ones including magnetic fields - we refer the reader there. We now
describe the impact of the supernova on material inside the cloud.

Upon impact, the primary shock has an equivalent strength to that in SN1 with
an approach velocity of $\sim$ 1000 km s-1. The interclump material is
impacted first and accelerated towards the cloud interior. By Fig. 9(b) the
shock has interacted with the clumps, and many individual shock-clump
interactions ensue. The clumps are accelerated less than their surroundings,
thus interclump material is seen to be pushed inwards between the clumps -
most apparent in panels (c – d). Note that by panel (d), the upstream flow
around the cloud contributes to this effect.

As the shock interacts with individual clumps, they appear to slightly
elongate towards the cloud centre and develop tail-like morphologies, an
effect due to their outer envelopes being stripped in the direction of the
post-shock flow and potentially due to impact with the material swept up by
the shock itself (Pittard, 2011). The material connecting the clumps is
disturbed also and contributes to this tail-like appearance, again orienting
radially towards the centre. Note that the irregular interface and presence of
clumps distorts the shock front, and as such its position is not obvious in
the panels. However by panel (d) it is evident the shock has interacted with
outer clumps around all of the cloud’s edge, as they all appear to have these
elongated morphologies from their stripped envelopes.

Tail-like features were seen both in NoShock and S2, but to very different
extents. In the NoShock scenario, material was accelerated radially inwards by
gravity only, with clumps accelerated slower than the interclump material. The
resulting velocity gradients were quite weak, although large enough to
slightly ablate the clump envelopes, resulting in similar, albeit much less
pronounced tail-like morphologies compared to S2/SN2. This is seen in both
Paper I and WPFVL16 (who also observe that equivalent models in 2D exaggerate
this effect due to motion being restricted to a plane). In the S2 scenario,
this effect is much more pronounced, as after shock passage the clumps become
entrained in a continuous flow. Although the initial impact is weaker than in
SN2, the post-shock flow is in total much more disruptive, with clump material
being continuously ablated. The constant flow results in continuous
acceleration of the interclump material, resulting in some clumps breaking
away from the cloud. These clumps become fully entrained in the external flow,
which further exaggerates the formation of elongated tails, and destroys some
clumps completely. In contrast, no clumps are seen to break away from the
cloud in the SN2 scenario.

In addition to seeing clumps elongate, they also appear to increase in size
and coalesce into larger cold regions. This appears to be due to their outer
envelopes expanding. As the shock enters the cloud, its post shock pressure
drops due to the combination of increased cooling in the cloud, gas attempting
to match the decreasing external pressure, and gas expansion effects
characteristic of Sedov-Taylor-like evolution. Due to these effects, the
clumps become embedded in material with lower pressure than the initial
interclump material, and their outer envelopes expand, with the remainder of
the material bound together by gravity. The cooling is strongest at the front,
and as such we see the majority of clump expansions occur there (see Fig. 9d –
f). Thus this increase in size is not to be mistaken by a gain in clump mass -
also apparent by the lack of deviation from the cold mass fraction NoShock
(Fig. 9j).

Note that what we could also be seeing is thermal-instability induced material
forming at the front in a manner similar to SN1. However, it is difficult to
find explicit evidence that would differentiate this from the individual clump
expansion, as the effects of the thermal instability would be seen at the
front (similar to SN1), which is where the clumps expand the most. In
addition, any migration across the phase diagram akin to figure 5 is obscured
as the existence of the two-phase medium occupies exactly that portion of the
diagram. The only possible evidence is a short-lived increase of cold material
seen in Fig. 9(j) which occurs on the order of the cooling timescale of the
post-shock gas. However this is negligible ($\lesssim$ 1% increase, $\sim$ 100
M⊙), and returns almost immediately to the NoShock fraction. Thus in this
scenario we conclude that once the clumps have formed, the supernova impact
does not trigger any further thermal instability.

At later stages of evolution, many of the clumps appear to collide and merge.
This is interesting, as we would expect this effect to raise the maximum
density in the cloud. Inspecting individual clumps shows that densities do
increase in merging regions, however these are small and do not exceed the
maximum. There does appear to be a slight increase in the maximum density for
a brief period (see Fig. 9i) however again this returns to NoShock values, and
eventually decreases further. In S2, we can see that at 5.16 Myrs densities
are attained that are two orders of magnitude larger ($\sim$ 105 cm-3) than in
SN2 ($\sim$ 103 cm-3). Once again, this can be attributed to the fact that the
S2 cloud was continuously compressed, resulting in local gravitational
collapse, and an acceleration of the global collapse of the cloud.

It is interesting therefore to observe that a supernova blast is unable to
create such conditions. In fact, looking at the density snapshots, and
examining the extracted data from the model - after the first 3.5 Myrs the
cloud is almost indistinguishable from the NoShock scenario. Once the initial
impact of the supernova is seen, the cloud settles into a steadily, globally
collapsing state with gravity dominating the evolution. Nevertheless, even
with more than 50 per cent of the cloud material being cold and dense
(increasing up to $\sim$ 90 per cent), a single supernova event is unable to
aid gravity in collapsing individual structures prior to the global collapse
of the cloud. Thus we conclude that no extra star formation is to be expected
from this interaction. Further work is warranted to see if such behaviour is
seen when a supernova interacts with TI formed clouds with different masses,
and ones simulated at higher resolution that results in the hydrodynamic
formation of filaments (as seen in WFP19).

##  4 Discussion

One would expect that all else being equal, the closer a blast source is to an
object, the stronger the impact. Indeed, studying the impact of supernovae on
molecular cloud velocity dispersions, Seifried et al. (2018) found that
doubling R${}_{\textup{c}}$ decreased the maximum impact by roughly 40%. It is
interesting therefore to observe that shocks significantly disrupt clouds in
the small cloud regime, while they do not have much impact on clouds in the
large cloud regime, given that the small cloud regime is often employed as an
approximation to supernova remnants. It is prudent to ask how reflective such
simulations are of realistic astrophysical systems.

To illustrate such considerations with our clouds, a density contrast $\chi$
$\approx$ 500 evaluates equation (4) to r${}_{\textup{cl}}$ $\ll$
0.0047R${}_{\textup{c}}$, requiring a cloud to have a radius that is
approximately 0.5 % the size of the blast radius R${}_{\textup{c}}$ to qualify
for the small cloud regime. For a cloud with r${}_{\textup{cl}}$ = 50 pc this
requires a minimum R${}_{\textup{c}}$ $\approx$ 11 kpc, which is unrealistic
for a single supernova. In fact, using equation (39.31) in Draine (2010) to
calculate the "fade-away radius" - the theoretical length-scale at which the
shock Mach number falls below unity and a supernova remnant mixes with the
ISM, we get R${}_{\textup{fade}}$ $\approx$ 200 pc for our 10 M⊙, 1051 erg
supernova exploding into a medium with n${}_{\textup{amb}}$ = 0.0022 cm-3.
Using R${}_{\textup{fade}}$ to approximate the maximum distance
R${}_{\textup{c}}$ that would lead to an interaction, according to equation
(4), only clouds with r${}_{\textup{cl}}$ $\ll$ 1 pc would fall into the small
cloud regime. Clearly, our S1/S2 models are not representative of a supernova
like that seen in SN1/SN2 exploding far away.

Whilst a single supernova may not drive the kind of wind seen in S1/S2, a
superbubble due to a cluster of supernovae might, opening up the parameter
space to exploring regimes between a single supernova, to whole clusters of
supernovae (e.g. Padoan et al., 2016; Fielding et al., 2018). For medium and
large cloud regimes, it is clear that a constraint of the parameter space
requires a supernova rate of 1 per t${}_{\textup{p}}$ (equation 1) as a
necessary condition in a clustered environment. Given the size of our clouds,
it is likely that our S1/S2 scenarios are thus more reflective of the kinds of
clouds embedded in galactic winds and superbubbles.

To speculate further however, given typical wind speeds (500 – 1500 km s-1)
and molecular phase speeds (50 – 300 km s-1) (see e.g. Strickland & Heckman,
2009; Rupke, 2018) it’s unclear at what region of the parameter space these
models occupy. For instance, at 5.16 Myrs when evidence of local gravitational
collapse was seen and the simulations were stopped, the cloud’s centre of mass
velocities were only $\sim$ 6 km s-1. Therefore it does not appear likely that
our clouds could be accelerated to speeds close to 50 km s-1 before the cloud
is destroyed by the external wind or subsequent feedback.

With this in mind, it is nevertheless possible that in a scenario with
feedback, a powerful wind due to a star with mass $>$ 100 M⊙ breaks out (e.g.
Rogers & Pittard, 2013; Wareing et al., 2017) and mass loads the external
flow. In such a scenario, warm cloud material could be expelled out of the
cloud and subsequently cool and condense, perhaps due to the thermal
instability, forming fast moving clumps and resulting in a system with a hot
ionized component and a cold molecular component (see e.g. Zhang et al., 2017;
Gronke & Oh, 2018). However at this stage this behaviour is very speculative,
and further exploration of the parameter space is needed to draw any more
conclusions. This is beyond the scope of this paper, and is left for future
work.

##  5 Conclusion

In this work, we study the interaction of a supernova remnant with a cloud
that is evolving due to the thermal instability (TI) and gravity, and draw
comparisons to idealised shock-cloud scenarios and a scenario without shocks.
A total of five 3D hydrodynamical simulations were performed: an un-shocked
scenario NoShock, two shock-cloud scenarios in the small cloud regime, S1/S2,
respectively analogous to 12Shock/24Shock in Paper I, and two supernova-cloud
scenarios in the large cloud regime, SN1/SN2, looking at impacts at equivalent
times to S1/S2. Just like Paper I, the disturbance is introduced at a "pre-TI"
(S1/SN1) and "post-TI" (S2/SN2) stage, and the models are evolved for a
timescale corresponding to the first snapshots when the shock-cloud scenarios
S1/S2 showed evidence of local gravitational collapse.

For our chosen parameters, the Mach number of the supernova shock is roughly a
factor $\sim$ 4.5 $\times$ stronger at impact than that of the planar shock.
However on aggregate the supernova is significantly less disruptive than the
sustained impact from the idealised post-shock flow. In S1/S2, the post-shock
flow is constantly replenished and approaches the cloud at $\sim$ 100 km s-1.
This provides a ram pressure that consistently throughout the simulation is a
factor $\sim$ 100 $\times$ larger than from the $\sim$ 10 km s-1 flows that
develop around the cloud in the SN1/SN2 scenarios. Additionally, the
constantly replenished flow means that the ambient thermal pressure is
maintained at higher values than those in the cloud, resulting in continuous
compression throughout the entire simulation. In contrast, the profile within
the remnant is Sedov-Taylor-like. As such the impact of the blast on the cloud
decays rapidly, the cloud is not met with a continuous flow, and the ambient
thermal pressure drops below that of the cloud in less than 1 Myr. These facts
accounts for all of the differences seen between the supernova/shock models,
and we observe the following particular outcomes:

  1. 1.

SN1/SN2: The cloud is impacted by two shocks - a primary shock that is the
forward shock from the initial blast, and a secondary shock resulting from the
reverse shock that reflects at the blast epicentre and propagates back
outwards and towards the cloud. The impact of the secondary shock on each
cloud is weak, and its influence is mainly felt at the cloud – ambient
interface.

  2. 2.

SN1/SN2: The primary and secondary shocks sweep over the cloud and generate
complex wave patterns that propagate both downstream and upstream. The
upstream components generate a flow on the cloud interface that converges at
the front.

  3. 3.

SN1/SN2: The primary shock loses much of its strength as it sweeps over the
cloud. When it converges at the rear of the cloud, the convergence has
negligible effect in amplifying the impact. In fact, as the model evolves, the
strongest ram pressures occur at the front of the cloud where the upstream
interface flow converges. This is in stark contrast to S1/S2, where the
strongest compression was at the rear.

  4. 4.

SN1: The lower ambient pressure compared to that of the cloud (ram + thermal)
results in a constant and gradual cloud expansion. Where the upstream flow
converges, the expansion is prevented - giving the cloud a lobe-like
morphology.

  5. 5.

SN1: The lack of continuous flow means that the Rayleigh-Taylor, Vishniac and
Kelvin-Helmholtz instabilities are not triggered.

  6. 6.

SN1: The TI is triggered behind the primary shock. This forms a cap-like
structure with clumps embedded inside with properties expected from such a
phase transition. The maximum densities in these structures do not exceed 100
cm-3, however the cold mass fraction is increasing, suggesting the structure
is not in its final stages of formation. Nevertheless, we do not expect star
formation to occur on the same timescales as S1.

  7. 7.

SN2: The passage of the shock through the interclump medium results in the
formation of tail-like structures due to clump envelopes being ablated. This
does not happen to the same extent as in S2, e.g. no clumps are seen to break
away from the parent cloud and get entrenched in an external flow. This effect
is however more pronounced than that due to flows resulting only from
gravitational acceleration of material in the NoShock scenario.

  8. 8.

SN2: After shock passage, the clumps appear to increase in size. This is not
due to an increase in mass, and instead due to expansion. This happens because
the inter-clump pressure drops below that of the clumps after shock passage.

  9. 9.

SN2: Clumps are seen to coalesce and merge. This however does not increase the
maximum density in the cloud or cause any local gravitational collapse.

  10. 10.

SN2: Although the exterior pressure drops below that of the cloud within the
first 1 Myr, the cloud does not expand as it is gravitationally dominated.
After $\sim$ 3.5 Myr, the cloud is almost indistinguishable from the NoShock
scenario. Thus, apart from the trivial global collapse of the cloud, we do not
expect to see any additional star formation.

## ACKNOWLEDGEMENTS

We thank the anonymous referee for her/his helpful comments that led to the
further development of the discussion section, and the physical implications
of the models.

MMK acknowledges support from the Science and Technology Facilities Council
(STFC) for PhD funding (STFC DTP grant ST/R504889/1). JMP and SAEGF
acknowledge STFC for project funding (STFC Research Grant ST/P00041X/1). CJW
thanks Leeds Institute for Fluid Dynamics and the School of Chemical and
Process Engineering for the visiting research fellowship. The calculations
herein were performed on ARC3, part of the High Performance Computing
facilities at the University of Leeds, UK. Additional simulations were
performed using the DiRAC Data Intensive service at Leicester, operated by the
University of Leicester IT Services, which forms part of the STFC DiRAC HPC
Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding
via STFC capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC
Operations grant ST/R001014/1. DiRAC is part of the National e-Infrastructure.
MMK would like to thank John Holden for fruitful discussions.

## Data Availability

The dataset associated with this article will be made available in the
Research Data Leeds Repository at https://doi.org/10.5518/1130

## References

  * Cioffi et al. (1988) Cioffi D. F., McKee C. F., Bertschinger E., 1988, ApJ, 334, 252
  * Draine (2010) Draine B. T., 2010, Physics of the interstellar and intergalactic medium. Vol. 19, Princeton University Press
  * Elling (2009) Elling V., 2009, Acta Mathematica Scientia, 29, 1647
  * Falle (1991) Falle S. A. E. G., 1991, MNRAS, 250, 581
  * Falle (2005) Falle S. A. E. G., 2005, in , Adaptive Mesh Refinement-Theory and Applications. Springer, pp 235–253
  * Falle et al. (2020) Falle S., Wareing C., Pittard J., 2020, MNRAS, 492, 4484
  * Field (1965) Field G. B., 1965, ApJ, 142, 531
  * Fielding et al. (2018) Fielding D., Quataert E., Martizzi D., 2018, Monthly Notices of the Royal Astronomical Society, 481, 3325
  * Godunov (1959) Godunov S. K., 1959, Sbornik: Mathematics, 89, 271
  * Gronke & Oh (2018) Gronke M., Oh S. P., 2018, Monthly Notices of the Royal Astronomical Society: Letters, 480, L111
  * Hubber et al. (2013) Hubber D., Falle S. A. E. G., Goodwin S., 2013, MNRAS, 432, 711
  * Joung & Mac Low (2006) Joung M. R., Mac Low M.-M., 2006, ApJ, 653, 1266
  * Kanjilal et al. (2021) Kanjilal V., Dutta A., Sharma P., 2021, MNRAS, 501, 1143
  * Kinoshita et al. (2021) Kinoshita S. W., Nakamura F., Wu B., 2021, ApJ, 921, 150
  * Klein et al. (1994) Klein R. I., McKee C. F., Colella P., 1994, ApJ, 420, 213
  * Krumholz (2014) Krumholz M. R., 2014, Physics Reports, 539, 49
  * Kupilas et al. (2021) Kupilas M. M., Wareing C. J., Pittard J. M., Falle S. A. E. G., 2021, MNRAS, 501, 3137, Paper I
  * Kurganov & Tadmor (2000) Kurganov A., Tadmor E., 2000, J. Comp. Phys., 160, 241
  * Li et al. (2015) Li M., Ostriker J. P., Cen R., Bryan G. L., Naab T., 2015, ApJ, 814, 4
  * Lu et al. (2020) Lu Z.-J., Pelkonen V.-M., Padoan P., Pan L., Haugbølle T., Nordlund Å., 2020, ApJ, 904, 58
  * McKee & Ostriker (1977) McKee C. F., Ostriker J. P., 1977, ApJ, 218, 148
  * McKee et al. (1987) McKee C. F., Hollenbach D. J., Seab G. C., Tielens A., 1987, ApJ, 318, 674
  * Nakamura et al. (2006) Nakamura F., McKee C. F., Klein R. I., Fisher R. T., 2006, ApJSS, 164, 477
  * Padoan et al. (2016) Padoan P., Pan L., Haugbølle T., Nordlund Å., 2016, ApJ, 822, 11
  * Parker (1953) Parker E. N., 1953, ApJ, 117, 431
  * Pittard (2011) Pittard J. M., 2011, MNRAS, 411, L41
  * Pittard et al. (2009) Pittard J. M., Falle S. A. E. G., Hartquist T. W., Dyson J. E., 2009, MNRAS, 394, 1351
  * Quirk (1997) Quirk J. J., 1997, in , Upwind and High-Resolution Schemes. Springer, pp 550–569
  * Rogers & Pittard (2013) Rogers H., Pittard J., 2013, MNRAS, 431, 1337
  * Rogers & Pittard (2014) Rogers H., Pittard J., 2014, MNRAS, 441, 964
  * Rupke (2018) Rupke D. S., 2018, Galaxies, 6, 138
  * Seifried et al. (2018) Seifried D., Walch S., Haid S., Girichidis P., Naab T., 2018, ApJ, 855, 81
  * Stone & Norman (1992) Stone J. M., Norman M. L., 1992, ApJ, 390, L17
  * Strickland & Heckman (2009) Strickland D. K., Heckman T. M., 2009, The Astrophysical Journal, 697, 2030
  * Van Loo et al. (2007) Van Loo S., Falle S. A. E. G., Hartquist T., Moore T., 2007, A&A, 471, 213
  * Wada et al. (2000) Wada K., Spaans M., Kim S., 2000, ApJ, 540, 797
  * Walch & Naab (2015) Walch S., Naab T., 2015, MNRAS, 451, 2757
  * Wareing & Falle (2019) Wareing C. J., Falle S. A. E. G.and Pittard J. M., 2019, MNRAS, 485, 4686, WFP19
  * Wareing et al. (2016a) Wareing C. J., Pittard J. M., Falle S. A. E. G., Van Loo S., 2016a, MNRAS, 459, 1803, WPFVL16
  * Wareing et al. (2016b) Wareing C., Pittard J., Falle S. A. E. G., 2016b, MNRAS, 465, 2757
  * Wareing et al. (2017) Wareing C. J., Pittard J. M., Falle S. A., 2017, MNRAS, 470, 2283
  * Wareing et al. (2018) Wareing C., Pittard J., Wright N., Falle S. A. E. G., 2018, MNRAS, 475, 3598
  * Wareing et al. (2021) Wareing C., Pittard J., Falle S. A. E. G., 2021, MNRAS, 500, 2831
  * Wolfire et al. (1995) Wolfire M. G., Hollenbach D., McKee C. F., Tielens A., Bakes E., 1995, ApJ, 443, 152
  * Wolfire et al. (2003) Wolfire M. G., McKee C. F., Hollenbach D., Tielens A., 2003, ApJ, 587, 278
  * Zhang et al. (2017) Zhang D., Thompson T. A., Quataert E., Murray N., 2017, Monthly Notices of the Royal Astronomical Society, 468, 4801