© The Authors 2025

1. Introduction

Advances in observational techniques over the last decades have made it possible to study the three-dimensional (3D) geometry of structures in the nearby galactic interstellar medium (hereafter ISM, e.g., Arenou et al. 1992; Lallement et al. 2019; Edenhofer et al. 2024). The Local Bubble (hereafter, LB Cox & Reynolds 1987; Linsky & Redfield 2021) is of particular interest in this context. It is characterized by a diffuse, X-ray emitting cavity, with a diameter of several hundred parsec, which we happen to be observing right from the center (Zucker et al. 2022; Yeung et al. 2024). The LB is believed to be a superbubble (SB) that had been evacuated due to collective feedback from massive stars, such as ionizing radiation (Linsky & Redfield 2021), stellar winds (Heiles 1998), and supernovae (hereafter SNe Breitschwerdt & de Avillez 2006; Wallner et al. 2021).

The geometry of SBs and supernova remnants (SNRs) provides a valuable tool for understanding phenomena such as galactic outflows, chemical enrichment, and star-formation, based on both observation and theory. Moreover, while the LB is the only SB whose 3D geometry has been studied in great detail thus far, novel techniques and a wealth of data will enable the study of many more Galactic SBs (Leike et al. 2020; Edenhofer et al. 2024). Despite the lack of 3D information, extragalactic observations also provide hints to the geometry of SBs (Watkins et al. 2023; Jiménez et al. 2024). To interpret this wealth of data, predictions from numerical simulations and analytical models for the geometry of SNRs and SBs are required.

Over the past five decades, the evolution of spherical SBs expanding into a uniform ISM has been studied in great detail (e.g., Chevalier 1974; Cioffi et al. 1988; Truelove & McKee 1999). While these efforts have provided useful intuition for the different processes dominating the dynamics of expanding SBs and shaped the theoretical methods used to describe their evolution (Kim & Ostriker 2015; Romano et al. 2024a), they lack the complexity needed to explore the physical processes governing the departure from spherical symmetry.

The processes that might deform SNRs are manifold. It was recognized early on that blastwaves expanding into a vertically stratified atmosphere are stretched out along the density gradient (Kompaneets 1960; Laumbach & Probstein 1969). Overall, SNRs have been found to preferentially expand into low-density channels, following the density structure of the ambient ISM, shaped by gravity and turbulence (Kim & Ostriker 2015; Ohlin et al. 2019; Makarenko et al. 2023; Lau & Bonnell 2025). Moreover, galactic shear might stretch out a SB along the direction of rotation (Tenorio-Tagle & Palous 1987; Bisnovatyi-Kogan & Silich 1995).

Observations of starburst galaxies reveal that many galaxies host galactic outflows (Xu et al. 2022). This suggests that vertical stratification plays an important role in shaping the geometry of SBs, provided they are powered by a sufficiently strong source. Studies of SBs in nearby star-forming galaxies report ellipsoidal geometries, aligned with the galactic rotation (Watkins et al. 2023), indicating that galactic shear might be at play. However, from the same observations, it has become clear that density structures, such as low-density channels and high-density filaments, align themselves in the same way (Xie et al. 2024), making it difficult to disentangle the role of shear and density structure in shaping the geometry of SNRs.

While these studies have shown the effectiveness of these various physical processes in deforming SNRs in isolation, there have been only a small number of works dedicated to addressing how they affect the geometry in concert (e.g., Jiménez et al. 2024, however, this study neglects radiative cooling). Indeed, most studies investigating the effect of stellar feedback in turbulent, stratified, and occasionally shearing media, focus on the collective effect stellar feedback has on the average properties of the multi-phase ISM and galactic outflows (e.g., de Avillez & Breitschwerdt 2005; Walch et al. 2015; Fielding et al. 2018; Kim & Ostriker 2017). However, a clear picture of how the different processes affecting SNR geometry compete remains unavailable.

In this paper, we present the Supernovae In a Stratified, Shearing ISM (SISSI) simulation suite to address this gap and enable a more comprehensive study of the phenomenology of SNRs. The SISSI project, which aims to evolve well-resolved SNRs in a realistic, but controlled environment, will enhance our theoretical understanding of the complex interaction of SNRs with their environment and provide future observational studies with new tools for disentangling the complex physics of SNRs in the galactic ISM.

The remainder of this paper is organized as follows. In Sects. 2 and 3, we describe the numerical and analysis methods and give a description of the SISSI simulation suite. In Sects. 4 and 5, we give an overview of the time evolution of our simulated sample of SNRs as well as an analysis of the geometry. We discuss our results in Sect. 6. Finally, we summarize our findings and present our conclusions in Sect. 7. In the appendix, we present the properties of the ISM of our simulated galaxy and provide some additional background to some of models and data used in our analysis.

2. Numerical methods

We modeled the evolution of SNRs embedded in an isolated disk galaxy, using the adaptive mesh refinement (AMR) code RAMSES (Teyssier 2002). It solves the system of hydrodynamic equations on a finite-volume, Cartesian grid using a second-order unsplit Godunov method (MUSCL scheme). The code reconstructs variables at the cell interfaces from the cell-centered values utilizing the HLLC Riemann solver with MinMod total variation diminishing scheme (Toro et al. 1994). RAMSES employs a conjugate gradient method and cloud-in-cell interpolation of particle contributions to solve the Poisson equation.

We related the gas pressure and internal energy using an adiabatic index of γ = 5/3. We implemented radiative cooling and heating based on the HEIKOU integration scheme (Behrendt et al., in prep.), based on the exact integration scheme (Townsend 2009; Zhu et al. 2017) and utilizing the UVB_dust1_CR0_G0_shield0 cooling table from Ploeckinger & Schaye (2020) at solar metallicity. We modeled the star formation by allowing gas with densities n_H > 100 cm⁻³ and temperatures T < 150 K to form star particles with m_⋆ = 10³ M_⊙ at the rate given by a local Schmidt-law (see, e.g., Katz 1992; Springel & Hernquist 2003; Shimizu et al. 2019; Oku et al. 2022), with ϵ_ff = 1%.

The simulation was separated into two stages. In the first stage, we relaxed an isolated disk galaxy into a quasi-steady state, where gravitational collapse and cooling were balanced by stellar-feedback-driven turbulence and heating. In the second stage, we turned off the feedback and zoomed in on the ISM in various locations, where we injected energy and mass to model the evolution of SNRs in a self-consistently generated galactic ISM.

2.1. Setup: Isolated disk galaxy

The SISSI galaxy is part of the AVALON galaxy formation and evolution project (Behrendt et al., in prep.), which utilizes the GALAXY COMPOSER package (Behrendt et al., in prep.) to generate the initial conditions of an isolated Milky-Way-like galaxy with galaxy parameters taken from Bland-Hawthorn & Gerhard (2016). The simulation domain is a cubical box with side length L = 48 kpc and outflow boundaries, subdivided into a coarse grid of 256 cubic cells, corresponding to a maximum cell size of Δx_max = 187.5 pc. Cells are refined up to an effective resolution of 2¹² (l_{max, ISM} = 12) or Δx_{min, ISM} ≈ 11.7 pc if they are larger than N_Jeans = 8 local Jeans lengths or if they contain a mass exceeding 20 (star) particle masses, which ensures that star-forming cells are Jeans-unstable. We modeled the influence of the stellar disk, bulge and dark matter halo as a static, axisymmetric background-potential. The gas was initially set up as a combination of a warm, isothermal disk in vertical hydrostatic equilibrium and a hot, diffuse uniform background.

During the initial relaxation stage, we modeled stellar feedback by injecting a thermal energy of 2 × 10⁵² erg and a mass of 200 M_⊙ into a single cell hosting a star particle 8 Myr after its formation. We avoided overcooling by flagging cells affected by stellar feedback with a passive scalar that disables cooling for the first ∼500 kyr after the feedback event.

We evolved the isolated disk galaxy for ≲500 Myr until it settled into a quasi-steady state, where gravitational collapse and cooling are balanced by feedback-driven heating and turbulence. We show a projection of the surface density of the ISM after the initial relaxation in Fig. 1. This figure depicts only a small cut-out of the simulation domain focusing on the galactic ISM. While the large box size is required to ensure a realistic galactic eco-system (galactic outflows and large-scale fountain-flows) and reduce numerical effects due to the domain boundaries, for this study, the details of the circumgalactic medium can be ignored, as our focus lies mainly on the central ∼8 kpc.

Fig. 1. Face-on (top) and edge-on (bottom) projection of the simulated galaxy at t = 0. We mark the explosion sites of the SNRs with star markers. Different marker colors correspond to the different passive scalars associated with the SN ejecta. The ISM in the inner ∼10 kpc is highly structured with filamentary outflows that reach several kpc above the midplane, while the ISM in the outskirts is rather smooth without any prominent vertical features.

2.2. Zoom-in: Treatment of SNRs

We flagged 30 star particles at three galactocentric radii R ∈ {2,4.5,8} kpc and z ∼ 0, spaced equidistantly in polar direction as SNR particles (see markers in Figs. 1 and 2). An overview of the local ISM properties in the selected regions is given in Appendix A.

Fig. 2. Initial vertical height of the explosion sites, grouped by galactocentric radius (star markers). Black dots denote the local galactic midplane; error bars the vertical scale height defined in the Appendix A. Radial coordinates, corresponding to R = 2,  4.5 and 8 kpc, were shifted for visibility. Even though the explosion sites were chosen to be close to z = 0, due to the warping of the disk, some of the SNRs are located outside the midplane.

We refined all cells within r_zoom, l = N_zoomΔx_l of an SNR particle up to a maximum zoom-in resolution of l_max = 18, corresponding to Δx_min ≈ 0.18 pc, where N_zoom = 15. We further relaxed the system for ≲50 kyr to avoid numerical artifacts due to the sudden refinement. Unless specified otherwise, we measured the time from the time of the snapshot at the end of this final relaxation step (t = 0). Starting from t = 0, each SNR particle injects N_SN SNe per injection; SN injections can happen every Δt_SN. The models differ only by the choice of N_SN and Δt_SN.

Per SN, each SNR particle distributes E_SN = 10⁵¹ erg of thermal energy and M_ej = 5 M_⊙ of ejecta mass evenly within a sphere of radius R_inj = 5Δx_min ≈ 0.92 pc centered on the SNR particle’s position. In addition each SNR particle injects one of two passive scalars Z_{ej, i}, corresponding to red and blue markers in Fig. 1, used to label the mass fraction of SN ejecta and distinguish between the ejecta of neighboring SNRs. We note that in Romano et al. (2024a), we demonstrated that the later evolution from the ST phase onward is correctly reproduced. This is in agreement with the results of Walch & Naab (2015), who found that different injection methods reproduce the ST phase equally well, provided short enough time steps are used. Other studies focusing on the numerical implementation of SN feedback further corroborate this conclusion (Dalla Vecchia & Schaye 2012; Kim & Ostriker 2015; Hopkins et al. 2018).

We refined polluted cells with Z_{ej, i} > 10⁻¹⁵ to at least l_{min, zoom} and even further, up to at most l_{max, zoom} if

$$\begin{array}{c}\hfill \mathrm{\Delta }x>0.1R_{\text{sf}}^{\text{KO15}}\left(n_{\text{H}}\right)=2.3{\left(\frac{n_{\text{H}}}{{\text{cm}}^{-3}}\right)}^{-0.42}\text{pc},\end{array}$$(1)

which roughly resembles the convergence criterion proposed by Kim & Ostriker (2015). We show an idealized refinement map in Fig. 3. We initially set the zoom-in resolutions to l_{min, zoom} = 14 and l_{max, zoom} = 18, corresponding to Δx_{max, zoom} ≈ 2.9 pc and Δx_{min, zoom} ≈ 0.18 pc, respectively, and reduced the resolution as the SNRs grow to keep the numerical cost at a manageable level (as shown in Fig. 4). With this refinement prescription, we were thus able to resolve the momentum generation during the Sedov-Taylor phase and the transition to the radiative phase for a single SN up to densities of n_{H, max} ∼ 430(Δx/Δx_{min, zoom})^−2.4 cm⁻³ (see, e.g., Kim & Ostriker 2015; Romano et al. 2024a).

Fig. 3. Refinement map produced with the refinement method outlined in Sect. 2.2 for the idealized situation of a diffuse bubble with a dense shell, designed to roughly resemble an SNR after shell formation. The solid-black, dashed-blue and dotted-green lines show the radial profiles of the refinement level, gas density, and ejecta fraction (scalar tracer field), respectively. The resolution is decreasing radially outward, levels off at lmin, zoom = 14 and increasing again to lmax, zoom = 18 inside the shell.

Fig. 4. Resolution in the zoom-in region as a function of time. Gray, red, and blue lines correspond to the three different runs, while different line styles correspond to different refinement parameters. The resolution was decreased between restarts of the simulation when the memory requirements became too large. The maximum resolution in the refinement regions around the central SNR particles was left untouched.

With our implementation of star-formation, the mass of star particles formed at higher resolution needs to be adjusted to ensure that stars are forming if (and only if) the cells are Jeans-unstable and fully refined. This condition can be satisfied by scaling m_⋆, l ∝ Δx_l.

2.3. Simulation suite: Overview

Our simulation suite consists of four different runs: a baseline simulation without SNe (N0) and three simulations with SNe labeled N1, N10, and N1×10, corresponding to (N_SN, Δt_SN) = (1, ∞), (10, ∞) and (1, 1 Myr), respectively. In particular, models N1 and N10 feature a single explosion event at t = 0, while it is only in N1×10 that there are subsequent explosion events every Δt_N1 × 10 = 1 Myr. In N0, no zoom-in was applied. To estimate the effect the refinement might have, we ran a fifth simulation labeled N0_zoom without SNe, but with N_zoom = 85.

3. Analysis

3.1. Classification of ISM components

To be able to meaningfully analyze the SNRs’ properties, we need to reliably differentiate between SNRs and the unperturbed ISM. Moreover, we classified different components of the SNRs, similarly to the approach of Romano et al. (2024a) for a single SNR in a uniform ISM.

We adopted the same method of using the passive scalars to flag cells belonging to an SNR. Neighboring SNRs inject different passive scalars, which enabled us to resolve ambiguities if the SNRs approach or even overlap. Naïvely, each SNR corresponds to the set of cells polluted with the respective scalar that are closest to its center (e.g., the corresponding SNR particle). However, in practice, since some SNRs get significantly larger than others, we find that this simple prescription would lead to a large number of cells being grouped incorrectly once the SNRs become too large. We avoided this problem by creating a weighted Voronoi-tesselation in face-on projection with cells centered at the position of the SNR particles and assigning weights, such that all polluted cells belonging to an SNR would end up within the corresponding cell. We assigned these weights by visual inspection.

As opposed to the case studied in Romano et al. (2024a), in this study, the ISM into which the SNRs are expanding is undergoing constant change. Thus, to study how the properties of the SNRs depend on the that of the ISM, we need to find an appropriate definition for the local ISM. Here, we define the local ISM as the contents of the smallest rectangular box, containing the entire SNR at all times. The unperturbed, local ISM, then corresponds to the contents of the local ISM without the SNRs, which necessarily shrinks as the SNRs grow. This leads to the slight bias, whereby once a region is swept up by the SNR, it ceases to contribute to the description of the unperturbed ISM. Nonetheless, the instantaneous state of the immediate surroundings of the SNRs describes the unperturbed ISM much more accurately than its state at a single point in time (i.e., at t = 0 or at the time of observation). We note that while we did not make use of this definition in the present study, future analyses using the SISSI simulations will include it (Romano et al., in prep.).

We further classified different components of the unperturbed ISM and the SNRs. For the SNRs, we followed the classification of Romano et al. (2024a). We distinguished between radially inflowing and outflowing shell and bubble components. The bubble corresponds to polluted, hot (T > 2 × 10⁴ K) or diffuse (n_H < 10⁻² cm⁻³) gas, while the shell corresponds to cold and dense, polluted gas. We decided whether the gas is in- or outflowing by measuring the radial velocity, measured from the center of mass of the SNR in the co-rotating, center-of-mass frame of each SNR.

For the unperturbed ISM, we distinguished between cold (T < 7 × 10³ K), warm (7 × 10³ K < T < 10⁵ K) and hot (10⁵ K < T) gas phases, which would be expected to coexist co-spatially in a turbulent medium with inhomogeneities driven by SN explosions and differential cooling (see e.g., McKee & Ostriker 1977; Cox 2005). The choice of 7 × 10³ K for the threshold between warm and cold gas, slightly less than the commonly used ∼10⁴ K, arises from the adopted cooling function, which produces persistent gas at this temperature as is shown in Appendix A. We also classified the stars within the ISM boxes based on whether they are old (i.e., formed before t = 0) or young. For the young stars, we further distinguished between stars that are formed from polluted or pristine gas.

3.2. Definition of polluted cells

We defined a cell as polluted if its passive scalar concentration exceeds some threshold value of Z_{ej, thr}. The choice of this threshold value is arbitrary and can systematically bias our results. If we choose a value of Z_{ej, thr} that is too low, we risk including gas that is only (slightly) polluted due to numerical noise, but that physically is not associated with the SNRs. On the other hand, if we choose a value that is too high, we risk missing parts of the SNRs.

In practice, it seems impossible to entirely prevent both effects from happening, so we aimed for a compromise and chose to give our results in terms of range of plausible values based on a slightly high and a slightly low threshold value. We first performed our analysis for a slightly low value Z_{thr, low} = 10⁻¹², comparable to the value used in Romano et al. (2024a). After defining the local ISM boxes, based on the SNRs defined by the choice of Z_{thr, low}, we defined $Z_{\text{thr, high}}^i\left(t\right)$ for each SNR and snapshot, by requiring that the total ejecta mass of cells with $Z_{\text{ej, i}}>Z_{\text{thr, high}}^i\left(t\right)$ just exceeds 99.99 percent of the total gas-phase ejecta-mass in the ISM box.

3.3. SNR geometry

We studied the dynamical evolution of the SNRs’ geometry by analyzing how their shape tensors evolve over time. We defined the shape tensor as

$$\begin{array}{c}\hfill S_{\mathit{ij}}=V_{\text{SNR}}^{-1}{\displaystyle {\int }_{\text{SNR}}}(||\mathbf{x}|{|}^2{\delta }_{\mathit{ij}}-x_i{x}_j){\text{d}}^3\mathbf{x},\end{array}$$(2)

which is the volume weighted inertia tensor, assuming a constant density of unity. By assuming an approximately ellipsoidal shape, we can define the three ellipsoidal radii, defined as

$$\begin{array}{c}\hfill r_i=\sqrt{2.5(\text{tr}\left(S\right)-2S_i)},\end{array}$$(3)

where S_i are the eigenvalues of S_ij and tr(S) the trace. We refer to the smallest, intermediate, and largest eigenvalue as the minor a, semi-major b and major c axis, respectively. We defined the effective size of an SNR as the geometric mean of the three eigenvalues

$$\begin{array}{c}\hfill r_{\text{eff}}={\left(abc\right)}^{1/3}.\end{array}$$(4)

To determine the alignment of the SNRs within the galaxy, we measured their pitch angle, α, and polar direction, cos(θ), for both the major and minor axes. The pitch angle is defined relative to the direction of galactic rotation, with α = 90° and α = −90° corresponding to the galactic center and anti-center, respectively. The magnitude of the polar direction is 0 (1) for directions parallel (perpendicular) to the galactic plane.

4. Time evolution of SNRs

4.1. Showcase: Supernovae in a relatively uniform medium

The case of stellar feedback in an ambient medium with solar metallicity and an ambient density of n_H ∼ 1 cm⁻³ has been widely studied (e.g., Kim & Ostriker 2015; Fierlinger et al. 2016; Oku et al. 2022; Romano et al. 2024a). In this section, we showcase the results of SNR #22, which happens to explode in a relatively uniform medium with an ambient density close to 1 cm⁻³ and compare its time evolution to that found in previous studies.

In Fig. 5, we show slices of the density field through the center of SNR #22 parallel to the xy-plane at various characteristic times for the different models. In each panel, the outline of the SNR is shown by orange lines, depicting contours of constant Z_ej, corresponding to Z_{thr, low} (dashed line) and Z_{thr, high} (solid line). As can be seen in the bottom row, corresponding to the N0 model, the density field is indeed rather uniform, but some collapse into a filamentary structure over several Myr is visible.

Fig. 5. Density slices through the central plane of the SNR #22 at various points in time for each model. Arrows are depicting the velocity field in the co-rotating center-of-mass frame of the local ISM. The various timescales correspond to different points in time for the different models. Red and blue arrows in the top-right corner of each panel indicate the directions of the galactic rotation and the galactic center, respectively. The dashed orange contour corresponds to the surface where Zej = Zthr, low, while the solid contour corresponds to Zej = Zthr, high. Since the various timescales are undefined for the model no_expl, we are using the same times as model N10. The SNe explode into a fairly homogeneous ISM, with a slowly collapsing, slight overdensity right where the SNe explode. At similar evolutionary stages the SNR is about twice as large in N10 compared to N1, with very similar geometry; Spheroidal with a slight elongation in the direction of rotation. On the other hand, the geometry in the model N1×10 qualitatively differs from the other models, with an elongated cavity normal to the rotational direction, due to the elliptical orbit (vR ∼ 20 km/s) of the explosion site. Only in the model N1, after 10 Myr a dense cloud, aligned with the SNR is forming in the center as predicted by Romano et al. (2024a).

At shell formation (first column), the SNRs are spherical with a slightly underdense central region and a thin, overdense shell. The time of shell formation and the SNRs’ sizes are in agreement with previous works (e.g., Kim & Ostriker 2015).

After shell formation, SNRs enter the so-called pressure-driven snowplow (PDS) phase, which ends once the pressure in the cavity drops below that of the shell (second column). At this time, the SNRs are spherical, with an increasingly underdense central region and a thin, overdense shell. The time at which the PDS phase ends and the SNRs’ sizes are in agreement with previous estimates (Romano et al. 2024a).

Romano et al. (2024a) showed that SNRs implode as they merge with their ambient medium. In their simulations, a SNR in an ambient medium with n_H ∼ 1 cm⁻³, such as the one considered here, began to implode after ∼1 Myr. In the third column, we show the SNRs right after the onset of implosion. By this time, the SNRs are slightly elongated, parallel to the collapsing filament, which is at a slight angle to the direction of galactic rotation. In all cases, the implosion occurs significantly later than our expectation based on previous work. In the model N1×10, the implosion seems to be coincident with the explosion happening at t ∼ 9 Myr and is no longer visible by t = 10 Myr. We ruled out this “implosion” as a false positive and caution that with our definition of the implosion timescale we cannot distinguish between brief moments of radially “inflowing” ejecta in N1×10 due to a displacement of the explosion sites and sustained inflows of cold gas from the shell.

After 10 Myr (fourth column), the SNRs in models N1 and N10 have been stretched out considerably in the direction of the collapsing filament. The implosion in N1 has reached the center and condensed into a growing, filamentary implosion cloud, as predicted by Romano et al. (2024a,b), which is stabilized by the rapid radiative dissipation of the energy carried by the implosion shocks colliding in a central region. In N10 the center of the SNR is still underdense indicating that the implosion has not yet reached the center. Meanwhile, in N1×10 the SNR is stretched out predominantly in radial direction following the wake of the explosion center, which happens to be drifting radially outward. The interior of the superbubble remains strongly underdense.

The volume traced by the dashed line corresponding to Z_{thr, low} tends to be slightly larger than the SNRs, particularly along the directions aligned with the Cartesian grid at early times, and at late times the direction of galactic rotation.

In Fig. 6 we show various global properties of the SNRs as a function of time and compare them to those of a single SN exploding into a uniform medium with an ambient density of n_H = 1 cm⁻³ taken from Romano et al. (2024a). We note that they used a different cooling function, leading to a slightly lower equilibrium pressure (see panel d).

Fig. 6. Various properties (mass, momentum, kinetic and thermal energy, pressure, and volume) of the SNR #22 as a function of time for the various models using the classification introduced in Sect. 3.1. Shaded regions correspond to the margin of uncertainty introduced by the choice of Zej, thr, while the lines correspond to the geometric average of the values obtained with high and low values of Zej, thr. For comparison, we show the time evolution of an isolated SNR at a similar ambient density (nH = 1 cm−3) taken from Romano et al. (2024a). As expected, the models N1 and N10 exhibit similar behavior and the model N1 also agrees quantitatively well with the isolated SNR. In the model N1×10, the SNR initially follows the model N1 and then after the onset of the consecutive SNe diverges reaching a comparable mass, momentum and size as the model N10 after 10 Myr. However, the fraction of thermal energy in the bubble is higher in N1×10 compared to N10, indicating more efficient hot phase generation. The differences due to the choice of Zej, thr are largest before shell formation and are the most pronounced in the mass and momentum of the bubble, indicating that ejecta are initially lagging behind the shock, but catch up once a cold shell forms.

All quantities are defined in essentially the same way as in Romano et al. (2024a), namely, extensive quantities are computed by summing up the contributions from all cells belonging to the respective gas phases and the pressure is calculated as a volume weighted average. However, due to the differential movement of the ambient medium, differences might arise in the momentum and kinetic energy due to the choice of inertial system, which we chose as the center-of-mass system after subtracting the galactic rotation, azimuthally averaged in linearly spaced, radial bins with a spacing of ΔR ≈ 47 pc.

In the model N1, all quantities except the pressure and the thermal energy agree with the isolated SNR for the first ∼2 Myr within uncertainties. At late times the radial momentum and kinetic energy drop more rapidly. The kinetic energy, eventually recovers and levels off at ∼1% of the injected explosion energy. However, these small differences might well be explained by the choice of the inertial system.

Model N10 appears to be a rescaled version of N1, in line with the idea that SNRs undergo a series of self-similar evolutionary stages. After 1 Myr, model N1×10 starts to diverge from N1. The amount of swept up mass, the total radial momentum, and kinetic as well as thermal energy of the SB grow to be quite similar to those of N10, at t = 10 Myr; albeit with large temporal variations in the distribution between the bubble and the shell. This indicates that these quantities are mostly sensitive to the total amount of injected energy, regardless of the exact interval between injections, in stark contrast to the geometry and mass distribution within the SNR, as shown in Fig. 5. Nonetheless, in model N1×10 the SNR is retaining more radial momentum, and sustaining a higher pressure, indicating that the subsequent radiative losses are determined by the multi-phase structure of the ISM at the injection site.

4.2. The full SISSI sample

We have shown that SNR #22 adheres well to the expectations from isolated SNRs in uniform ambient media. However, this might just have been a special case that cannot be applied to the whole sample. Thus, in this section we evaluate to what extent our full sample of SNRs follows the expectations from previous work.

In Fig. 7 we show various characteristic timescales as a function of ambient density, defined here as the ratio of the swept-up mass and the volume covered by each SNR at each point in time, for the different models and compare them with analytical results from previous work, described in Appendix B, shown as red and gray lines in the different panels.

Fig. 7. Various timescales as a function of ambient density for our simulated sample of SNRs. Uncertainties arise due to the finite spacing of the snapshots and due to the choice of Zej, thr. The timescales of shell formation, the end of the PDS phase and SNR implosion agree well with the predictions from models based on simulations of isolated SNRs (Kim & Ostriker 2015; Romano et al. 2024a). The timescale measuring the onset of star formation within the SNRs is within a factor of 5 of 10 percent of the free-fall timescale of the star-forming SNRs.

We find that the shell-formation timescale of the simulated SNRs matches the theoretical estimate Eq. (B.2), in line with the expectation that SNRs are hardly affected by the galactic environment during the early adiabatic expansion phase. The same holds true for the timescale for the end of the PDS (Eq. (B.5)), with some exceptions at very low densities.

Differences to the purely analytic picture become more apparent when comparing the timescale of implosion Eq. (B.7), where we assume σ₁ = 0.8 corresponding to an ambient pressure of $P_{\text{ISM}}\sim {10}^4{k}_{\text{B}}\text{K}{\text{cm}}^{-3}$, matching the pressure of the isobaric phase of the ISM (Fig. A.1). Here, we define t_launch slightly differently from Romano et al. (2024a), who defined t_launch as the time of the first snapshot when at least 0.1 M_⊙ are in the form of backflowing shell gas. In SISSI, this condition would be met at almost all times, due to the uncertainties in the selection of the SNR gas and the turbulent motion of the background medium. We thus restrict the criterion to the ejecta, and define t_launch as the earliest time when the backflowing part of the shell contains at least 2% of the ejecta, tagged by the respective scalar tracer. We find that while the bulk of SNRs is not too far from the analytical model, there is considerable scatter and a number of extreme outliers. Moreover, as we note in the discussion of Fig. 5, the interpretation of SNR implosion in the context of model N1×10 is somewhat unclear, as the interior pressure of the SBs tends to remain high.

We also show the time after which stars start to form from material polluted by SNe. There is no star-formation from polluted gas for the first 1 Myr, but within about a factor of 5 of 0.1 t_ff stars begin to form within the SNRs, where

$$\begin{array}{c}\hfill t_{\text{ff}}=\sqrt{\frac{3\pi }{32G\rho }}\sim 44.9n_0^{-0.5}\text{Myr},\end{array}$$(5)

is the free-fall timescale. Importantly, in many cases this star formation does not appear to be triggered within the SNRs; rather, it is the continued star-formation in pre-existing star-forming regions that are swept up and enriched by the SNRs. A more detailed analysis of the potential triggering of star-formation in SISSI is out of the scope of this work and will be the focus of future publications (Romano et al., in prep.).

In Fig. 8, we show the effective size (as defined in Sect. 3.3) as a function of ambient density for the different models at various characteristic points in time. Red, gray and blue lines depict the expected sizes, building on the theoretical models described in Appendix B. In the panel corresponding to the last snapshot at 10 Myr, the orange and blue square corresponds to the effective size of the LB derived from the data products of Edenhofer et al. (2024) as described in a companion paper (Romano et al., in prep.).

Fig. 8. SNR size at various characteristic points in time as a function of ambient density for our simulated sample of SNRs. Uncertainties arise due to the finite spacing of the snapshots and due to the choice of Zej, thr. An orange and blue square depicts the effective radius of the LB derived from the 3D dust maps of Edenhofer et al. (2024) in the panel corresponding to t = 10 Myr. Error bars are smaller than the marker and thus not shown. The radii at shell formation, the end of the PDS phase and at SNR implosion agree well with the predictions from models based on simulations of isolated SNRs (Kim & Ostriker 2015; Romano et al. 2024a) for sufficiently large ambient densities. At low densities, the sizes tend to exceed the model predictions. After 10 Myr the SNRs are about twice the expected size.

We find that overall SNR sizes are in line with theoretical expectations during the stages of SNR evolution before merging with the ISM, namely, before t = t_launch, but start growing larger than expected at later times. SNRs in low-density environments n_H ≲ 0.1 cm⁻³ start to diverge from the theoretical expectation by t_imp ≳ 1 Myr.

After 10 Myr, all SNRs are about twice the expected size, indicating the need for better models of old SNRs in a shearing, stratified ISM. Interestingly, the LB is on the smaller end of the sizes for simulated SNRs in similar density media, even though it is expected to be older, namely, t_LB ∼ 14 Myr (Zucker et al. 2022; Breitschwerdt & de Avillez 2006). We further discuss this point and its implications in a companion paper (Romano et al., in prep.).

In Fig. 9 we show the momentum input per SN at the end of the PDS stage and compare it to the theoretically expected value, assuming a momentum enhancement after shell-formation of ∼20%, slightly lower than the ∼50% reported by Kim & Ostriker (2015). The momentum input in the denser regions n_H ≳ 0.1 cm⁻³ roughly follows the theoretical expectation, with little scatter. In contrast to lower-density regions, where the momentum per SN drops off with large scatter. This behavior is likely due to the large size (≳100 pc) of these SNRs, leading to more frequent energy dissipation due to interactions with high-density structures.

Fig. 9. Outward radial momentum per SN at the end of the PDS phase as a function of ambient density. Uncertainties arise due to the finite spacing of the snapshots and due to the choice of Zej, thr. At sufficiently high densities, nH ≳ 0.1 cm−3, the momentum input offers a good match to the prediction of simulations of isolated SNRs (Kim & Ostriker 2015; Romano et al. 2024a). On the other hand, at lower densities, the momentum per SN is often lower than expected, with large error bars.

5. Geometry of simulated SNRs

In the previous section, we show that while young SNRs are well described by the theory based on models in a uniform, stationary medium, the models start to fail on longer timescales, namely, ≳1 Myr. To obtain some clues as to what might be causing these differences, we study their geometry, which reveals a preferential alignment that might point us towards the governing physical processes.

5.1. The shape phase-space

In Fig. 10, we show the trajectories of the SNR #22 for the different explosion models in the shape phase-space, defined by the minor-to-major and semi-major-to-major axis ratio. By definition, at t = 0 the SNR starts as a perfect sphere (a/c = b/c = 1) and by deformation through various processes might evolve to become increasingly prolate (b/c < 2/3) or oblate (a/c → 0 and b/c > 2/3). In purple, we also show the trajectory of a shearing sphere as described in Appendix C. The sphere of radius r₀ = 100 pc is initially located at a galactocentric radius of R₀ = 8 kpc and is rotating at a constant rotation velocity, matching that measured in the simulation.

Fig. 10. Evolutionary tracks of the SNR #22 in the shape phase-space for the various explosion models. Uncertainties due to the choice of Zej, thr are shown as shaded regions. In different parts of the phase space the SNRs are either spherical (S), oblate-spheroidal (OS), prolate (P), or oblate (O). The SNR starts out as a perfect sphere and becomes increasingly prolate over time. The ratio of the two minor axes remains close to one and never falls below 2/3. In the model N10 the SNR remains spherical for longer compared to N1. Similarly the consecutive SN explosions in the model N1×10 restore spherical symmetry.

In the models N1 and N10, the trajectory in the shape phase-space is smooth, with almost constant minor-to-semi-major axis-ratio a/b > 2/3 and ever decreasing a/c, namely, the SNRs are becoming increasingly prolate. In model N1, a/b is slightly larger than in N10, namely, the SNR is slightly more prolate. The simulated SNRs are significantly more deformed than the shearing sphere, which after 10 Myr is still quite spherical (a/c ∼ 0.75, b/c ∼ 0.85).

The trajectory in the model N1×10 has a kink, corresponding to the onset of further explosions, which deform the SNR in chaotic ways, ultimately leading to a more spherical shape. The final shape is similar to that of the shearing sphere.

In Fig. 11, we show the locations of our sample of SNRs in the shape phase-space at various characteristic points in time. Markers are colored by the ambient density. In the panel corresponding to the last snapshot at 10 Myr, we compare the shape of the LB derived from the data products of Edenhofer et al. (2024), shown as an orange and blue square, to our simulated sample.

Fig. 11. Distribution of SNR shapes at various characteristic points in time for the different explosion models. Uncertainties due to the choice of Zej, thr are represented by error bars. Different regions are labeled as in Fig. 10. An orange and blue square depicts the shape of the LB derived from the 3D dust maps of Edenhofer et al. (2024) in the panel corresponding to t = 10 Myr. Error bars are smaller than the marker and thus not shown. At shell formation the SNRs tend to be spherical, with SNRs in lower-density environments being somewhat less spherical. At the end of the PDS stage and the onset of the implosion, most SNRs are still spherical or oblate spheroids with a/b ≳ 0.5 − 0.67. However, some of the lower-density SNRs are already quite asymmetric falling into the prolate and oblate category. The SNRs in the N10 model tend to be somewhat more spherical. At 10 Myr, the majority of SNRs are asymmetric. In dense environments, SNRs tend to be quite asymmetric, with low a/c ∼ 0.2. The model N1 tends to have lower a/b ∼ 0.5 compared to the other explosion models, which tend to have a/b ≳ 2/3.

The three panels corresponding to shell-formation, the end of the PDS phase and the onset of the implosion reveal that most SNRs remain close to spherical throughout the main stages of SNR evolution, with a/b ≳ 2/3 and b/c ≳ 2/3. In contrast, SNRs in very low-density ambient media n_H ≲ 10⁻² cm⁻³ already begin to deviate from spherical symmetry before shell formation, likely due to their older age and larger size at the same evolutionary stage, indicating that they are likely tracing a more anisotropic environment than their high-density counterparts.

After 10 Myr, the trend is reversed. The SNRs with the highest ambient densities are deformed the most, exhibiting highly anisotropic shapes a/c ∼ 0.2 with a wide range of geometries 1/3 ≲ b/c ≲ 1. Only 4 SNRs remain spherical, with most SNRs being slightly prolate and some oblates. The LB has a usual shape for an SNR with its ambient density, being slightly prolate with a/b ∼ 0.8 and a/c ∼ 0.5.

5.2. Alignment of SNRs within the galaxy

In the previous subsection, we show that the simulated SNRs evolve towards increasingly anisotropic geometries, suggesting that they might expand more in certain directions than others. To check whether there are any preferential directions, in Fig. 12 we show the time-span weighted distribution of the pitch angles and the magnitude of the polar directions as defined in Sect. 3.3 of the minor axis (oblates, orange) and major axis (prolates, purple). We also show the alignment of the minor- and major-axes of the LB derived from the data products of Edenhofer et al. (2024), shown as an orange and a purple square, respectively. The 1σ-confidence intervals are shown as shaded regions in the one-dimensional histograms.

Fig. 12. Time-span-weighted histograms showing the distribution of SNR directions for asymmetric SNRs. Distributions of oblate (prolate) SNRs are colored orange (violet). For oblate (prolate) SNRs, we show the direction of the minor (major) axis. We show the polar direction normal to the disk plane and the pitch angle relative to the direction of galactic rotation. We do not differentiate between directions above or below the disk. Positive pitch angles point between the galactic center and negative angles point towards the galactic outskirts. Orange and purple squares depict the directions of the minor and major axes of the LB, respectively, derived from the 3D dust maps of Edenhofer et al. (2024). The 1σ-ranges for the LB are also indicated as shaded areas in the one-dimensional histograms. Oblate SNRs tend to point vertically out of the disk and towards the galactic outskirts with a typical pitch angle of αminor ≲ −50 deg. On the other hand, prolate SNRs tend to lie within the disk plane |cos(θ)| ≲ 0.5 pointing slightly towards the galactic center αmajor ∼ 10 − 60 deg.

We find that most of the time, the minor axis of the oblate SNRs is directed perpendicular to the disk plane, with a broad distribution of negative pitch angles, centered around α_oblate ∼ −60°. In contrast, in the case of prolate SNRs, the polar direction of the major axis is broadly distributed, with most of the weight lying below |cos(θ)| ≲ 0.5, corresponding to the directions within the galactic plane. The distribution of pitch angles has three peaks around α_major ∼ 15° , 25° and 50°, in line with the expectations for structures deformed by shear (Appendix C, see also the alignment of underdense substructure in Fig. 1).

While the LB is slightly prolate, its minor axis points in a direction in agreement with that of oblate SNRs. On the other hand, its major axis is pointing slightly towards the galactic outskirts and is slightly more perpendicular than the bulk of our sample of SNRs.

5.3. Deformation timescale

In the previous subsections we have found that our sample of simulated SNRs evolving into the shearing, stratified ISM of the SISSI galaxy grow increasingly anisotropic, assuming a geometry that aligns with the sheared structure of the galaxy. The time it takes for an initially spherically symmetric structure, such as an SNR, to become deformed hints at the processes governing deformation overall. To this end, we define the deformation timescale as the time at which the minor-to-major ratio drops below a/c = 2/3.

The shearing-sphere model (Appendix C) indicates that shear can deform a spherical structure within a few percent of the orbital timescale. To test, whether shear alone is enough to explain the deformation of the SNRs we show the deformation timescale as a function of galactocentric radius in Fig. 13. Markers are colored based on the shape classification of the SNRs.

Fig. 13. Deformation timescale as a function of galactocentric radius. Uncertainties arise due to the finite spacing of the snapshots and due to the choice of Zej, thr. For visibility, markers are slightly shifted around their respective radii Rgal = 2, 4.5 and 8 kpc, with the same shift used for the same SNR, but different explosion model. Typical deformation timescales are on the order of a few percent of the orbital timescale at each radius, slightly shorter than what is expected from deformation by galactic shear alone.

The majority of the SNRs are deformed within ≲1% of the orbital timescale, with several SNRs being deformed much before even a thousandth of an orbit. More spherical SNRs, namely, SNRs that are classified neither as prolate or oblate, tend to have longer deformation timescales, more plausibly explicable by shear alone. We note that there are more oblate SNRs observed at larger galactocentric radii. Overall, the deformation of the SNRs is too rapid to be explained by shear alone.

Another likely relevant source of deformation are preexisting density anisotropies in the ambient ISM, which imprint onto the geometry of the SNRs as they expand into them (Makarenko et al. 2023). We quantify the degree of spatial variation in the density field, by measuring its relative variation at t = 0 across spatial scales, by averaging over nested ISM patches of side-lengths ℓ = 0.2, 0.5, 1, 1.5 and 2 kpc, corresponding to the scatter in Fig. A.2.

In Fig. 14, we show the deformation timescale as a function of the above-defined density dispersion. We find a steep decline in the deformation timescale with increasing density dispersion ∝(δρ/ρ)⁻³, in qualitative agreement with the expectation that indeed anisotropies in the density distribution might be dictating the geometry of SNRs in a turbulent ISM. Since the dense structures in the ISM themselves are subject to differential rotation, they can be stretched out considerably by galactic shear over timescales that are much longer than the age of the SNR and thus imprint a relatively larger degree of anisotropy than the expansion of an SNR subject to shear alone.

Fig. 14. Deformation timescale as a function of density dispersion. Uncertainties arise due to the finite spacing of the snapshots and due to the choice of Zej, thr. The deformation timescale is roughly ∝(δρ/ρ)−3 with significant scatter, qualitatively in line with the expectation that SNRs are deformed earlier in more anisotropic media.

6. Discussion

In the previous sections, we describe the evolution of the geometry of SNRs expanding into an ISM structured by the complex interplay of gravity, galactic rotation, and turbulence. In the following, we discuss some of the limitations of our simulations and how our results compare to observations of SNRs, along with the implications of our findings to the study of galaxy evolution and the structure of the ISM.

6.1. Limitations

The SISSI simulation suite aims to simulate the evolution of SNRs in a realistic galactic environment. Of course, simulating a realistic galactic ISM is challenging and the numerical prescriptions and sub-grid physics involved can greatly influence the phase structure and the morphology of the galaxy as a whole. A discussion of the quality of the simulated ISM of the SISSI galaxy, which is part of the AVALON simulation suite, focusing on the modeling of various aspects of galaxy evolution related to the structure of the ISM is out of the scope of this work, and will be presented elsewhere (Behrendt et al., in prep.).

Highly resolved simulations of SNRs, alongside an entire galaxy are computationally challenging and the computational resources to resolve large patches of the ISM with our maximum zoom-in resolution, are out of reach for currently available computing hardware. Therefore, we had to resort to the refinement strategy outlined in Sect. 2.2, which might itself introduce numerical artifacts. In the case of a uniform medium and without gravity, Romano et al. (2024a) have shown that our method of a co-evolving refinement region does not greatly affect the evolution of SNRs. However, we do find some differences in the star-forming activity and the partition of energy between the models N0 and N0_zoom, suggesting that the properties of the background ISM might differ substantially, based on the resolution. We accounted for this fact by relaxing the initial conditions for 50 kyr; however, it remains unclear whether there is an optimal relaxation duration, given the differences between the high- and low-resolution ISM do not seem to reach an asymptotic state. Instead, this could be attributed to chaos, due to unresolved gravitational collapse coupled to stochastic star-formation.

More detailed modelling of physical processes, such as cosmic rays, magnetic fields, thermal conductivity, and non-equilibrium radiation chemistry, as well as more detailed stellar models that include sources of early stellar feedback, can influence the dynamics as well as the geometry of SNRs (e.g., Gentry et al. 2019; Makarenko et al. 2023; Guo et al. 2025; Diesing & Gupta 2025). In the present work, we opted for a lightweight physics model to lower the computational cost and allow for a higher resolution. However, future efforts involving more detailed physics models will certainly be worthwhile.

6.2. Observations of SNR geometry

As SNRs evolve, so do the wavelengths of light in which they can be observed. Young SNRs are usually observed in the optical, infrared and X-ray (e.g., Fesen et al. 2023; Kobashi et al. 2024; De Looze et al. 2024) and their geometry is dictated by the explosion mechanism as well as their immediate surroundings. These SNRs are usually fairly close to spherical symmetry, justifying our spherically symmetric injection of energy and ejecta mass.

Once SNRs enter the ST phase, they are extremely hot and bright in X-rays with diffuse X-ray emission coming from their center (Khabibullin et al. 2023; Reynolds & Borkowski 2024). In this evolutionary stage, most SNRs are very close to spherically symmetric, though interactions with nearby clouds can lead to asymmetric features (Chi et al. 2024) in agreement with our sample of simulated SNRs, where the majority of SNRs remain close to spherically symmetric, with the exception of those exploding in low-density conditions, which are likely affected by interactions with clouds and low-density channels.

Galactic SNRs are usually only observed until shortly after they enter the radiative stage, as they quickly become too faint to be observed. Observed radiative SNRs, tend to be quite spherically symmetric (Paylı et al. 2024), with just a few exceptions attributed to interactions with nearby density structures (Arias et al. 2024). This finding is in agreement with our simulations, which indicate that most SNRs remain spherically symmetric between shell-formation and the end of the PDS phase.

Older SNRs are usually too faint to be observed directly. However, once they become large enough, they might be indirectly observed by looking for large cavities in the dust distribution. Such large cavities are routinely observed both in our Galaxy (Pelgrims et al. 2020; Zucker et al. 2022; Li et al. 2022; Verma et al. 2023) as well as in nearby galaxies (Watkins et al. 2023; Sánchez-Cruces & Rosado 2023; Li et al. 2024). At first glance, observed SBs exhibit a wide variety of shapes and orientations; however, due to the scarcity of detailed analyses of their shapes and orientation with respect to galactic structure, it is difficult to say to what extent the observed sample agrees or disagrees with our simulated sample.

Fortunately, thanks to the recent 3D dust map made available by Edenhofer et al. (2024), we are in a position to study the geometry of the LB, a SB believed to be excavated by the SN explosions of ∼𝒪(10) massive stars within the last ∼𝒪(10) Myr (Breitschwerdt & de Avillez 2006; Wallner et al. 2021; Zucker et al. 2022). We find that, while the orientation of the LB is slightly unusual for an SB its age, its shape fits right into the range of shapes that we report for our sample of simulated SNRs after 10 Myr. Moreover, we find that its effective size of R_eff ∼ 212.3 ± 1.0 pc is on the lower end of sizes obtained after 10 Myr. In a companion paper focusing on this aspect (Romano et al., in prep.), we further discuss these implications, particularly with respect to current age estimates of the LB.

6.3. Implications and future directions

Our numerical simulations show that the geometry of evolved SNRs is changed due to the complex interplay of the expanding shell and a variety of environmental factors. Since different processes affect the geometry of the SNR on different timescales, it might be possible to disentangle their contributions and deepen our understanding of the underlying physical processes shaping the internal structure of galaxies.

By leveraging novel analysis techniques (Edenhofer et al. 2024), and increasingly detailed observations (Gaia Collaboration 2023) it will soon be possible to study the geometry of an ever growing sample of galactic SNRs. Already, the data products of Edenhofer et al. (2024) can be used to study the neighboring known SBs, such as the Per-Tau SB (Bialy et al. 2021) and GSH 238+00+09 (Heiles 1998), which could provide additional hints to the assembly of structures in the solar neighborhood.

While we have focused on SNR geometry in this work, there are many more aspects of SNR phenomenology that can be addressed by the SISSI simulations. In future studies we aim to investigate the role of SNRs in driving interstellar turbulence, their coupling to and potential driving of galactic outflows as well as triggered star-formation.

7. Concluding remarks

We introduce the SISSI simulation suite, featuring 3D hydrodynamic zoom-in simulations of SNRs embedded in the realistic, self-consistently generated ISM of an isolated, Milky-Way-like galaxy. Our aim in this work is to deepen our understanding of various aspects of SNR physics. We focus on the geometry of the SNRs and show how it can be used as a useful observational diagnostic for understanding the various environmental effects, affecting the evolution of the system. Here, we summarize our most important findings:

1. The dynamics of young SNRs (≲1 Myr) are well described by standard analytical models. However, these models become less accurate for SNRs exploding in low-density (≲0.1 cm⁻³) environments, likely due to the large size of the SNRs, which increases the likelihood of interactions with both high- and low-density structures, such as clouds and channels. At later times, environmental effects start to dominate the dynamics of all simulated SNRs, leading to SNRs that are ≳2 times larger than predicted by models of SNRs in uniform, stationary media (see Fig. 8).

2. Overall, SNRs tend to be deformed greatly on timescales shorter than a few percent of the orbital timescale. In environments with larger density fluctuations, SNRs get deformed earlier, as they tend to follow the geometry of dense structures, which are deformed by shear and aligned with the galactic rotation.

3. The deformation of SNRs has preferred directions. The minor axis of oblate SNRs tends to be aligned with the galactic poles, with a slight tilt towards the galactic outskirts, suggesting that the vertical expansion of these SNRs is stalled by the graviational pull of the galactic disk. The polar angle of the major directions is broadly distributed, slightly favoring directions in the galactic plane, with pitch angles peaked between ∼15° and ∼50° slightly pointing towards the galactic center. This is in agreement with the expectation from alignment due to galactic shear.

4. The LB has a typical geometry for a SB of its age and size; however, it appears slightly small compared to the size of the SNRs in the SISSI sample at 10 Myr. This suggests that previous estimates of the age based on idealized, 1D models of SB expansion might need to be revised. The LB might be just old – and large – enough to be affected by its galactic environment. This makes it a unique laboratory for studying the expansion of SBs at the interface between local ISM physics and galactic dynamics.

We conclude that SNR geometry offers a novel observational tool for understanding the complex physics of galaxies and their impact on galactic substructure, leveraging the full potential of recent high-quality observations.

Acknowledgments

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 systems Cobra and Viper 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. LR thanks the developers of the following software and packages that were used in this work: JULIA v1.10.0 (Bezanson et al. 2017), MATPLOTLIB v3.5.1 (Hunter 2007), MERA v1.4.4 (Behrendt 2023), RAMSES v19.10 (Teyssier 2002), and HEALPIX v2.3.0 (Tomasi & Li 2021).

References

Appendix A: The ISM of the SISSI galaxy

The SNRs of the SISSI simulations expand into a complex environment which differs between different regions of the ISM. To help interpret the role of the ISM, we start by describing the properties of this environment.

Figure A.1 shows the T − n_H phase diagram of the gas in the SISSI galaxy at t = 0. In the ISM most of the gas is concentrated in two distinct gas phases: Warm neutral gas at T ∼ 7 × 10³ K for densities in the range ∼10⁻² − 10² cm⁻³, and colder gas at a constant pressure of $P\sim {10}^4{k}_{\text{B}}\text{K}{\text{cm}}^{-3}$ for denser gas above n_H ∼ 1 cm⁻³. For gas above n_H ∼ 1 cm⁻³ the warm phase is unstable and cools after ∼1 Myr. Cold (T ≲ 10² K) and dense (n_H ≳ 10² cm⁻³) gas is star-forming and, thus, it is steadily being consumed in this process. The galactic ISM is surrounded by a hot, diffuse circumgalactic medium (CGM), corresponding to a roughly adiabatic phase with T_CGM ∼ 10⁷ (n_H/(10⁻³ cm⁻³))^γ − 1 with a maximum density of ∼10⁻³ cm⁻³ at a temperature of ∼10⁷ K. It is in pressure balance with the cold ISM.

Fig. A.1. Temperature-density phase-diagram of the global ISM. The galactic ISM features a stable phase at P/kB ∼ 104 K cm−3 for nH ∼ 1 - 100 cm−3, along with an unstable phase at T ≲ 104 K, which becomes stable below nH ∼ 1 cm−3.

Locally, the ISM properties might not adhere to this simple picture. In Fig. A.2 we show the initial properties of the local ISM at the various explosion sites. All quantities are derived from density, and vertical velocity profiles averaged over ISM patches of side lengths ℓ = 0.2, 0.5, 1, 1.5,  and 2 kpc parallel to the galactic plane centered around the explosion site. The error bars indicate how much a quantity varies with scale.

Fig. A.2. Corner plot showing the distribution of ISM properties (Galactocentric radius R, midplane density nH, mid, velocity dispersion σ, vertical scale height zscale and vertical velocity gradient ∂zvz) at the SNR locations averaged over various length scales. We average over quadratic apertures with side lengths L = 0.2,  0.5,  1.0,  1.5 and 2.0 kpc. As expected the mean of nH, mid decreases with R, with large scatter at R > 2 kpc. Compared to the ICs, the mean density profile has steepened. The velocity dispersion is roughly constant σ ∼ 10 km/s throughout the disk with considerable scatter. The disk scale height follows the scaling behavior predicted by vertical hydrostatic equilibrium of the gas $z_{\text{scale}}\propto \sigma /n_{\text{H}}^{0.5}$, but with a slightly lower normalization, likely due to the presence of stars. The vertical velocity gradient on average is zero with increasing scatter towards larger R. The maximum values reveal a preferred scale of $\left|{\partial }_z{v}_z\right|\sim \sigma /z_{\text{scale}}\sim t_{\text{ff}}^{-1}$ indicating that these expansion and contraction motions correspond to gravitational breathing modes of the disk. The velocity gradients peak around nH, mid ∼ 1 cm−3, namely, the breathing of the disk coincides with the presence of multi-phase gas (see Fig. A.1).

We define the vertical scale height as half the distance between the ∼12 % and the ∼88 % mass percentiles of the density profile, roughly matching the definition of a sech²(z/z_scale)-profile, as would be appropriate for a single-component, isothermal disk. While the simulated galaxy does neither consist of only a single component, nor is it isothermal, our definition still yields a reasonable definition for an effective gas scale-height. Correspondingly, we define the midplane density as

$$\begin{array}{c}\hfill n_{\text{H, mp}}=\frac{1}{2\text{tanh}\left(1\right)z_{\text{scale}}}{\displaystyle {\int }_{z_0-z_{\text{scale}}}^{z_0+z_{\text{scale}}}}{n}_{\text{H}}\left(z\right)\text{d}z,\end{array}$$(A.1)

where z₀ is the midplane, right between the ∼12 % and the ∼88 % mass percentiles of the density profile. The velocity dispersion is defined as the average of the three components of the velocity dispersion vector, namely, σ² = (∑_{i ∈ {x,y,z}}σ_i²)/3, within all vertical bins within the midplane, namely, within z₀ ± z_scale.

We find a diverse range of midplane densities spanning over two orders of magnitude. The densities roughly follow the radial trend of the initial conditions, albeit with considerable scatter and a steepening towards R_gal = 2 kpc. Thus, we expect SNRs at larger galactic radii to grow bigger, with more variation between regions.

The velocity dispersion is roughly constant throughout the sample of regions with a typical value of σ ∼ 10 km s⁻¹, though with larger spatial variations in some regions. Therefore, SNRs would be expected to merge with the ISM at around the same time in all regions.

The measured scale heights indicate that, while the overall trend follows the expected scaling from dynamical equilibrium considerations in a single-component disk, it is more compact due to the dominant gravitational potential of the stellar disk. SNRs will start to be affected by vertical stratification once their size grows similar to this scale height, indicating that these effects might become important earlier for SNRs in higher-density regions.

In some regions we find that the mean vertical velocity is increasing (decreasing) linearly as a function of height with a midplane vertical velocity gradient on the order of σ/z_scale. These motions appear to be strongest around n_{H, mp} ∼ cm⁻³, indicating that gas at this density can be thermally unstable and is driven towards lower (higher) densities (∂_zv_z > 0 (< 0)). We interpret these motions as disk breathing-modes around the dynamical equilibrium, with an expected period of about a free-fall timescale, much longer than the dynamical timescale of an expanding SNR. Thus we expect the velocity gradients to be frozen-in during the lifetime of an SNR. SNRs expanding into a positive velocity gradient will grow faster as they sweep-up co-expanding material, while SNRs expanding into a collapsing region will be slowed down.

Appendix B: Analytic theory of SNR evolution in a uniform medium

In this section, we briefly review the analytic theory for the dynamics of radiative SNRs and SBs (see e.g., Kim & Ostriker 2015; Oku et al. 2022; Romano et al. 2024a). We consider the case of spherical expansion driven by point-explosions with explosion energy E_SN = 10⁵¹ E₅₁ erg into a uniform medium with hydrogen number density n_H = n₀ cm⁻³, solar metallicity and pressure P = μ n_H σ², where μ = 1.4 is the mean atomic weight and $\sigma =10{\sigma }_1\text{km}{\text{s}}^{-1}$ is the sound speed, which in a supersonically turbulent medium such as the ISM, might be replaced with the turbulent velocity dispersion.

Since here we are mostly interested in the dynamics of old radiative SNRs we skip the dynamics of the initial ejecta dominated expansion and start directly with that of adiabatic expansion; the so-called Sedov-Taylor (ST) phase (Sedov 1959). The internal structure of the Sedov-Taylor blastwave is described by a similarity solution with similarity parameter ξ = r/(E_SNt³/ρ)^1/5, with ξ₀ ≈ 1.15167 at the position of shock radius. During the ST phase, the radially outward momentum increases as a function of time and is given by (Kim & Ostriker 2015)

$$\begin{array}{c}\hfill p_{\text{ST}}=2.21\times {10}^4{E}_{51}^{4/5}{n}_0^{1/5}{t}_3^{3/5}{\text{M}}_{\odot }\text{km}{\text{s}}^{-1},\end{array}$$(B.1)

where t = t₃ kyr = t₆ Myr.

The ST phase ends, once radiative cooling becomes dominant and a thin shell forms right behind the shock front, after (Kim & Ostriker 2015)

$$\begin{array}{c}\hfill t_{\text{sf}}\sim 44E_{51}^{0.22}{n}_0^{-0.55}\text{kyr},\end{array}$$(B.2)

at which point the SNR has a size of

$$\begin{array}{c}\hfill R_{\text{sf}}=22.6E_{51}^{0.29}{n}_0^{-0.42}\text{pc},\end{array}$$(B.3)

and a momentum of

$$\begin{array}{c}\hfill p_{\text{sf}}=2.17\times {10}^5{E}_{51}^{0.93}{n}_0^{-0.13}{\text{M}}_{\odot }\text{km}{\text{s}}^{-1}.\end{array}$$(B.4)

Right after shell formation, the interior of the SNR is still hot and at a higher pressure than the shell, which has a temperature of about T_shell ∼ 10⁴ K and is highly compressed relative to the ambient medium χ ∼ 10. During this so-called pressure-driven snowplow (PDS) stage, the SNR expands R ∝ t^2/7, leading to a slight enhancement of the radial momentum ∝t^1/7. The PDS ends when the pressure in the bubble becomes comparable to the pressure in the shell (Romano et al. 2024a) after

$$\begin{array}{c}\hfill t_{\text{PDS}}\sim 0.15E_{51}^{0.27}{n}_0^{-0.44}\text{Myr},\end{array}$$(B.5)

corresponding to a size of

$$\begin{array}{c}\hfill R_{\text{PDS}}\sim 32.1E_{51}^{0.3}{n}_0^{-0.39}\text{pc}.\end{array}$$(B.6)

Depending on the details of radiative cooling and incorporation of mass from the bubble into the shell, the radially outward momentum of the SNR is boosted by up to about ∼50 % during the PDS.

Once the pressure in the interior of the bubble has dropped, the SNR expands solely due to its inertia R ∝ t^1/4. During this so-called momentum-conserving snowplow (MCS) phase, the pressure of the shell is proportional to the shock velocity P_shell ∝ t^−3/2. During this stage, the back of the shell is unstable and a reflected shockwave or implosion is driven into the interior of the SNR once the pressure of the shell becomes comparable to that of the ambient medium (Romano et al. 2024a) after

$$\begin{array}{c}\hfill t_{\text{launch}}\sim 0.5E_{51}^{0.27}{n}_0^{-0.44}{\sigma }_1^{-4/3}\text{Myr},\end{array}$$(B.7)

at which point the SNR has a size of

$$\begin{array}{c}\hfill R_{\text{launch}}\sim 43.6E_{51}^{0.3}{n}_0^{-0.39}{\sigma }_1^{-1/3}\text{pc}.\end{array}$$(B.8)

These expressions differ from those derived by Romano et al. (2024a), due to the differences in our model for the ambient pressure.

There are various different models for the case of an SB, driven by subsequent SN explosions (e.g., El-Badry et al. 2019; Oku et al. 2022). The basic assumption in these models is that, if the age of the SB is greater than the average time between SN explosions t ≫ Δt_SN = Δt₆ Myr, the expansion can be approximately described by that of a wind with a constant mechanical luminosity $L=E_{\text{SN}}/\mathrm{\Delta }{t}_{\text{SN}}={10}^{45}{L}_{45}\text{erg}{\text{yr}}^{-1}$. The dynamics of a radiative SB depend on the efficiency of energy dissipiation, for instance, due to radiative cooling, facilitated by thermal conduction and turbulent mixing of the hot interior and the cold shell. When energy injection dominates over dissipation, the expansion can be described by that of an energy-driven wind (Weaver et al. 1977; El-Badry et al. 2019). In contrast, if cooling losses dominate, the expansion is effectively momentum-driven (Lancaster et al. 2021; Oku et al. 2022; Lancaster et al. 2024).

In the model N1x10, Δt_SN = 1 Myr ≫ t_cool, corresponding to the momentum-driven regime. The size of a momentum-driven SB is given by Eq. 26 of Oku et al. (2022)

$$\begin{array}{c}\hfill R=40t_6^{1/2}{L}_{45}^{0.23}{n}_0^{-0.28}\text{pc}.\end{array}$$(B.9)

Appendix C: Shearing-sphere model

To model the deformation by shear, here we derive the simplest possible model: starting from t = 0, a sphere of radius r₀, centered in the galactic midplane at a galactocentric radius R₀, is subjected to differential rotation with a constant rotation speed V_rot, corresponding to an angular frequency of Ω(R) = V_rot/R. Correspondingly, the orbital timescale is t_orb(R) = 2π Ω⁻¹(R).

We parameterize the surface of the sphere using the polar angles θ and φ at t = 0

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_t(\theta ,\varphi ;r_0)=\left(\begin{array}{c}R\text{cos}({\varphi }_0+\mathrm{\Omega }t)\\ R\text{sin}({\varphi }_0+\mathrm{\Omega }t)\\ r_0\text{cos}\left(\theta \right)\end{array}\right),\end{array}$$(C.1)

where

$$\begin{array}{c}\hfill R^2=R_0^2+2R_0{r}_0\text{cos}\left(\phi \right)\text{sin}\left(\theta \right)+r_0^2{\text{sin}}^2\left(\theta \right)\end{array}$$(C.2)

is the galactocentric radius of the point on the surface of the sphere and the initial azimuthal angle ϕ₀ is defined by

$$\begin{array}{cc}\hfill R\text{cos}\left({\varphi }_0\right)& =R_0+r_0\text{sin}\left(\theta \right)\text{cos}\left(\phi \right)\hfill \end{array}$$(C.3)

$$\begin{array}{cc}\hfill R\text{sin}\left({\varphi }_0\right)& =r_0\text{sin}\left(\theta \right)\text{sin}\left(\phi \right).\hfill \end{array}$$(C.4)

Due to the differential rotation, parts of the sphere that are at a larger galactocentric radius lag behind and the parts that are at a smaller R advance ahead, leading to deformation. It can be shown that in spite of the deformation, the volume remains constant.

We measure the geometry as defined, in Sect. 3.3, namely, by computing the shape tensor

$$\begin{array}{c}\hfill S_{\mathit{ij}}=V^{-1}{\displaystyle {\int }_V}(\left||\mathrm{\Delta }{\mathrm{\Phi }}_t\right|{|}^2{\delta }_{\mathit{ij}}-\mathrm{\Delta }{\mathrm{\Phi }}_{t,i}\mathrm{\Delta }{\mathrm{\Phi }}_{t,j})\text{d}V,\end{array}$$(C.5)

where ΔΦ_t = Φ_t(θ,ϕ;r) − Φ_c(t) is the coordinate vector of a point within the shearing ball, relative to the volume-weighted center

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\text{c}}\left(t\right)=V^{-1}{\displaystyle {\int }_V}{\mathrm{\Phi }}_t(\theta ,\varphi ;r)\text{d}V.\end{array}$$(C.6)

Since the polar points are always co-rotating with the center, the semi-major axis ${\widehat{e}}_{\text{semi-major}}={\widehat{e}}_z$ and b = r₀. The major axis evolves from α_{major, 0} ∼ 45° towards α_{major, ∞} ∼ 0°, and necessarily α_minor = α_major + 90°.

For the sake of a better intuition of the model, in Figs. C.1 and C.2 we show the shape phase-space trajectories as well as the time evolution of the pitch angle for various shearing spheres with different r₀/R₀. We show the time evolution over half an orbit. The trajectories in shape phase-space differ only marginally between different sized spheres for r₀ ≲ 0.1 R₀. However, at later times larger spheres tend to have slightly larger pitch angles and minor-to-major axis ratios.

Fig. C.1. Same as Fig. 10 for the shearing sphere model for spheres with different sizes, evolved for 0.5 torb(R0−r9). The phase-space trajectories of the spheres with different sizes are almost identical.

Fig. C.2. Time evolution of the pitch angle of the major axis of the shearing sphere for spheres with different sizes. The pitch angle starts off near 45° and decays over time. Larger spheres tend to have slightly larger pitch angles.

In Fig. C.3 we show the dependence of the deformation timescale on the size of the sphere for r₀ < 0.1 R₀. We find a weak linear dependence on the size, which is approximately fit by

Fig. C.3. Deformation timescale as a function of size. Spheres are deformed greatly after ∼6.5 percent of $t_{\text{orb}}\left(\overline{R}\right)$, where $\overline{R}=p{R}_0+(1-p)(R_0-r_0)$ is a characteristic radius lying between R0 − r0 and R0. Our linear fit yields p ∼ 0.65.

$$\begin{array}{c}\hfill t_{\text{deform}}\approx 0.065t_{\text{orb}}\left(\overline{R}\right),\end{array}$$(C.7)

where $\overline{R}=p{R}_0+(1-p)(R_0-r_0)$ is a characteristic radius lying between R₀ − r₀ and R₀. Our linear fit yields p ∼ 0.65.

We do not show any results for r₀ > 0.1 R₀ due to the large difference in orbital timescales, which leads to rapid deformation and even winding of the part of the sphere with R ≪ R₀, while the rest has hardly moved.

All Figures

Fig. 1. Face-on (top) and edge-on (bottom) projection of the simulated galaxy at t = 0. We mark the explosion sites of the SNRs with star markers. Different marker colors correspond to the different passive scalars associated with the SN ejecta. The ISM in the inner ∼10 kpc is highly structured with filamentary outflows that reach several kpc above the midplane, while the ISM in the outskirts is rather smooth without any prominent vertical features.

In the text

Fig. 2. Initial vertical height of the explosion sites, grouped by galactocentric radius (star markers). Black dots denote the local galactic midplane; error bars the vertical scale height defined in the Appendix A. Radial coordinates, corresponding to R = 2,  4.5 and 8 kpc, were shifted for visibility. Even though the explosion sites were chosen to be close to z = 0, due to the warping of the disk, some of the SNRs are located outside the midplane.

In the text

Fig. 3. Refinement map produced with the refinement method outlined in Sect. 2.2 for the idealized situation of a diffuse bubble with a dense shell, designed to roughly resemble an SNR after shell formation. The solid-black, dashed-blue and dotted-green lines show the radial profiles of the refinement level, gas density, and ejecta fraction (scalar tracer field), respectively. The resolution is decreasing radially outward, levels off at lmin, zoom = 14 and increasing again to lmax, zoom = 18 inside the shell.

In the text

Fig. 4. Resolution in the zoom-in region as a function of time. Gray, red, and blue lines correspond to the three different runs, while different line styles correspond to different refinement parameters. The resolution was decreased between restarts of the simulation when the memory requirements became too large. The maximum resolution in the refinement regions around the central SNR particles was left untouched.

In the text

Fig. 5. Density slices through the central plane of the SNR #22 at various points in time for each model. Arrows are depicting the velocity field in the co-rotating center-of-mass frame of the local ISM. The various timescales correspond to different points in time for the different models. Red and blue arrows in the top-right corner of each panel indicate the directions of the galactic rotation and the galactic center, respectively. The dashed orange contour corresponds to the surface where Zej = Zthr, low, while the solid contour corresponds to Zej = Zthr, high. Since the various timescales are undefined for the model no_expl, we are using the same times as model N10. The SNe explode into a fairly homogeneous ISM, with a slowly collapsing, slight overdensity right where the SNe explode. At similar evolutionary stages the SNR is about twice as large in N10 compared to N1, with very similar geometry; Spheroidal with a slight elongation in the direction of rotation. On the other hand, the geometry in the model N1×10 qualitatively differs from the other models, with an elongated cavity normal to the rotational direction, due to the elliptical orbit (vR ∼ 20 km/s) of the explosion site. Only in the model N1, after 10 Myr a dense cloud, aligned with the SNR is forming in the center as predicted by Romano et al. (2024a).

In the text

Fig. 6. Various properties (mass, momentum, kinetic and thermal energy, pressure, and volume) of the SNR #22 as a function of time for the various models using the classification introduced in Sect. 3.1. Shaded regions correspond to the margin of uncertainty introduced by the choice of Zej, thr, while the lines correspond to the geometric average of the values obtained with high and low values of Zej, thr. For comparison, we show the time evolution of an isolated SNR at a similar ambient density (nH = 1 cm−3) taken from Romano et al. (2024a). As expected, the models N1 and N10 exhibit similar behavior and the model N1 also agrees quantitatively well with the isolated SNR. In the model N1×10, the SNR initially follows the model N1 and then after the onset of the consecutive SNe diverges reaching a comparable mass, momentum and size as the model N10 after 10 Myr. However, the fraction of thermal energy in the bubble is higher in N1×10 compared to N10, indicating more efficient hot phase generation. The differences due to the choice of Zej, thr are largest before shell formation and are the most pronounced in the mass and momentum of the bubble, indicating that ejecta are initially lagging behind the shock, but catch up once a cold shell forms.

In the text

Fig. 7. Various timescales as a function of ambient density for our simulated sample of SNRs. Uncertainties arise due to the finite spacing of the snapshots and due to the choice of Zej, thr. The timescales of shell formation, the end of the PDS phase and SNR implosion agree well with the predictions from models based on simulations of isolated SNRs (Kim & Ostriker 2015; Romano et al. 2024a). The timescale measuring the onset of star formation within the SNRs is within a factor of 5 of 10 percent of the free-fall timescale of the star-forming SNRs.

In the text

Fig. 8. SNR size at various characteristic points in time as a function of ambient density for our simulated sample of SNRs. Uncertainties arise due to the finite spacing of the snapshots and due to the choice of Zej, thr. An orange and blue square depicts the effective radius of the LB derived from the 3D dust maps of Edenhofer et al. (2024) in the panel corresponding to t = 10 Myr. Error bars are smaller than the marker and thus not shown. The radii at shell formation, the end of the PDS phase and at SNR implosion agree well with the predictions from models based on simulations of isolated SNRs (Kim & Ostriker 2015; Romano et al. 2024a) for sufficiently large ambient densities. At low densities, the sizes tend to exceed the model predictions. After 10 Myr the SNRs are about twice the expected size.

In the text

Fig. 9. Outward radial momentum per SN at the end of the PDS phase as a function of ambient density. Uncertainties arise due to the finite spacing of the snapshots and due to the choice of Zej, thr. At sufficiently high densities, nH ≳ 0.1 cm−3, the momentum input offers a good match to the prediction of simulations of isolated SNRs (Kim & Ostriker 2015; Romano et al. 2024a). On the other hand, at lower densities, the momentum per SN is often lower than expected, with large error bars.

In the text

Fig. 10. Evolutionary tracks of the SNR #22 in the shape phase-space for the various explosion models. Uncertainties due to the choice of Zej, thr are shown as shaded regions. In different parts of the phase space the SNRs are either spherical (S), oblate-spheroidal (OS), prolate (P), or oblate (O). The SNR starts out as a perfect sphere and becomes increasingly prolate over time. The ratio of the two minor axes remains close to one and never falls below 2/3. In the model N10 the SNR remains spherical for longer compared to N1. Similarly the consecutive SN explosions in the model N1×10 restore spherical symmetry.

In the text

Fig. 11. Distribution of SNR shapes at various characteristic points in time for the different explosion models. Uncertainties due to the choice of Zej, thr are represented by error bars. Different regions are labeled as in Fig. 10. An orange and blue square depicts the shape of the LB derived from the 3D dust maps of Edenhofer et al. (2024) in the panel corresponding to t = 10 Myr. Error bars are smaller than the marker and thus not shown. At shell formation the SNRs tend to be spherical, with SNRs in lower-density environments being somewhat less spherical. At the end of the PDS stage and the onset of the implosion, most SNRs are still spherical or oblate spheroids with a/b ≳ 0.5 − 0.67. However, some of the lower-density SNRs are already quite asymmetric falling into the prolate and oblate category. The SNRs in the N10 model tend to be somewhat more spherical. At 10 Myr, the majority of SNRs are asymmetric. In dense environments, SNRs tend to be quite asymmetric, with low a/c ∼ 0.2. The model N1 tends to have lower a/b ∼ 0.5 compared to the other explosion models, which tend to have a/b ≳ 2/3.

In the text

Fig. 12. Time-span-weighted histograms showing the distribution of SNR directions for asymmetric SNRs. Distributions of oblate (prolate) SNRs are colored orange (violet). For oblate (prolate) SNRs, we show the direction of the minor (major) axis. We show the polar direction normal to the disk plane and the pitch angle relative to the direction of galactic rotation. We do not differentiate between directions above or below the disk. Positive pitch angles point between the galactic center and negative angles point towards the galactic outskirts. Orange and purple squares depict the directions of the minor and major axes of the LB, respectively, derived from the 3D dust maps of Edenhofer et al. (2024). The 1σ-ranges for the LB are also indicated as shaded areas in the one-dimensional histograms. Oblate SNRs tend to point vertically out of the disk and towards the galactic outskirts with a typical pitch angle of αminor ≲ −50 deg. On the other hand, prolate SNRs tend to lie within the disk plane |cos(θ)| ≲ 0.5 pointing slightly towards the galactic center αmajor ∼ 10 − 60 deg.

In the text

Fig. 13. Deformation timescale as a function of galactocentric radius. Uncertainties arise due to the finite spacing of the snapshots and due to the choice of Zej, thr. For visibility, markers are slightly shifted around their respective radii Rgal = 2, 4.5 and 8 kpc, with the same shift used for the same SNR, but different explosion model. Typical deformation timescales are on the order of a few percent of the orbital timescale at each radius, slightly shorter than what is expected from deformation by galactic shear alone.

In the text

	Fig. 14. Deformation timescale as a function of density dispersion. Uncertainties arise due to the finite spacing of the snapshots and due to the choice of Zej, thr. The deformation timescale is roughly ∝(δρ/ρ)−3 with significant scatter, qualitatively in line with the expectation that SNRs are deformed earlier in more anisotropic media.
In the text

	Fig. A.1. Temperature-density phase-diagram of the global ISM. The galactic ISM features a stable phase at P/kB ∼ 104 K cm−3 for nH ∼ 1 - 100 cm−3, along with an unstable phase at T ≲ 104 K, which becomes stable below nH ∼ 1 cm−3.
In the text

Fig. A.2. Corner plot showing the distribution of ISM properties (Galactocentric radius R, midplane density nH, mid, velocity dispersion σ, vertical scale height zscale and vertical velocity gradient ∂zvz) at the SNR locations averaged over various length scales. We average over quadratic apertures with side lengths L = 0.2,  0.5,  1.0,  1.5 and 2.0 kpc. As expected the mean of nH, mid decreases with R, with large scatter at R > 2 kpc. Compared to the ICs, the mean density profile has steepened. The velocity dispersion is roughly constant σ ∼ 10 km/s throughout the disk with considerable scatter. The disk scale height follows the scaling behavior predicted by vertical hydrostatic equilibrium of the gas $z_{\text{scale}}\propto \sigma /n_{\text{H}}^{0.5}$, but with a slightly lower normalization, likely due to the presence of stars. The vertical velocity gradient on average is zero with increasing scatter towards larger R. The maximum values reveal a preferred scale of $\left|{\partial }_z{v}_z\right|\sim \sigma /z_{\text{scale}}\sim t_{\text{ff}}^{-1}$ indicating that these expansion and contraction motions correspond to gravitational breathing modes of the disk. The velocity gradients peak around nH, mid ∼ 1 cm−3, namely, the breathing of the disk coincides with the presence of multi-phase gas (see Fig. A.1).

In the text

	Fig. C.1. Same as Fig. 10 for the shearing sphere model for spheres with different sizes, evolved for 0.5 torb(R0−r9). The phase-space trajectories of the spheres with different sizes are almost identical.
In the text

	Fig. C.2. Time evolution of the pitch angle of the major axis of the shearing sphere for spheres with different sizes. The pitch angle starts off near 45° and decays over time. Larger spheres tend to have slightly larger pitch angles.
In the text

	Fig. C.3. Deformation timescale as a function of size. Spheres are deformed greatly after ∼6.5 percent of $t_{\text{orb}}\left(\overline{R}\right)$, where $\overline{R}=p{R}_0+(1-p)(R_0-r_0)$ is a characteristic radius lying between R0 − r0 and R0. Our linear fit yields p ∼ 0.65.
In the text