# Evolution and distribution of superbubbles in simulated Milky Way-like
galaxies

Chengzhe Li1, Hui Li1,2, Wei Cui1, Federico Marinacci3,4, Laura V. Sales5 Mark
Vogelsberger6 and Paul Torrey7
1Department of Astronomy, Tsinghua University, Beijing 100084, China
2Department of Astronomy, Columbia University, New York, NY 10027, USA
3Department of Physics & Astronomy “Augusto Righi", University of Bologna, via
Gobetti 93/2, 40129 Bologna, Italy
4INAF, Astrophysics and Space Science Observatory Bologna, Via P. Gobetti
93/3, I-40129 Bologna, Italy
5Department of Physics and Astronomy, University of California, Riverside, CA,
92521, USA
6Department of Physics & Kavli Institute for Astrophysics and Space Research,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA
7Department of Astronomy, University of Florida, Gainesville, FL 32611, USA
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Stellar feedback plays a crucial role in regulating baryon cycles of a
galactic ecosystem, and may manifest itself in the formation of superbubbles
in the interstellar medium. In this work, we used a set of high-resolution
simulations to systematically study the properties and evolution of
superbubbles in galactic environments. The simulations were based on the
SMUGGLE galaxy formation framework using the hydrodynamical moving-mesh code
Arepo, reaching a spatial resolution of $\sim 4\,\rm pc$ and mass resolution
of $\sim 10^{3}\,\rm M_{\sun}$. We identified superbubbles and tracked their
time evolution using the parent stellar associations within the bubbles. The
X-ray luminosity-size distribution of superbubbles in the fiducial run is
largely consistent with the observations of nearby galaxies. The size of
superbubbles shows a double-peaked distribution, with the peaks attributed to
early feedback (radiative and stellar wind feedback) and supernova feedback.
The early feedback tends to suppress the subsequent supernova feedback, and it
is strongly influenced by star formation efficiency, which regulates the
environmental density. Our results show that the volume filling factor of hot
gas ($T>10^{5.5}~{}\mathrm{K}$) is about $12\%$ averaged over a region of 4
kpc in height and 20 kpc in radius centered on the disk of the galaxy.
Overall, the properties of superbubbles are sensitive to the choice of subgrid
galaxy formation models and can, therefore, be used to constrain these models.

###### keywords:

galaxies: evolution – galaxies: ISM – ISM: bubbles – methods: numerical

††pubyear: 2023††pagerange: Evolution and distribution of superbubbles in
simulated Milky Way-like galaxies–B

##  1 Introduction

Stellar feedback plays a critical role in the evolution of galaxies. The
absence of stellar feedback would result in excessive cooling of the
interstellar gas, and the formed stellar mass would be much larger than the
observed values (i.e. the overcooling problem, see recent reviews, e.g., Naab
& Ostriker, 2017; Vogelsberger et al., 2020). Stellar feedback, mainly from
massive stars, injects energy and momentum into the surrounding gas and
regulates the phases of the interstellar medium (ISM). Among possible stellar
feedback channels, supernovae and radiation are the most widely studied (e.g.
Stinson et al., 2013; Sales et al., 2014; Gupta et al., 2016). Stellar
feedback can transport metal-rich gas to the circumgalactic medium (CGM) or
even the intergalactic medium (IGM) through supernova driven galactic winds,
and is therefore an important part of baryon cycling in a galaxy ecosystem
(Tumlinson et al., 2017).

The most direct product of stellar feedback is the emergence of a multiphase
ISM structure, with the hot phase mainly corresponding to superbubbles and
galactic winds. Superbubbles are usually considered as large structures (about
$100~{}\mathrm{pc}$ in scale) filled with hot gas that are blown by successive
supernovae or stellar winds. The Solar System is located inside a superbubble
(i.e. the Local Bubble, Cox & Reynolds, 1987; McCammon et al., 1971), which is
mainly filled with highly ionized gas of temperature $T\sim
10^{6}~{}\mathrm{K}$ and electron density $n_{\rm{e}}\sim
0.005~{}\mathrm{cm^{-3}}$ (Snowden et al., 1998). Outside the Local Bubble
lies cooler gas. The amount of neutral or partially ionized gas is expected to
increase substantially on the boundary, which leads to a fairly accurate
measurement of the Local Bubble extension, ranging from $100\,\rm{pc}$ to
$200\,\rm{pc}$ for different sight lines (e.g. Lallement et al., 2003; Welsh
et al., 2010). Apart from the Local Bubble, other superbubbles were also
discovered with X-ray observations in the Milky Way and nearby galaxies, such
as the Large Magellanic Cloud (LMC). Chu & Mac Low (1990) searched for diffuse
X-ray emission around OB associations in the LMC and identified a number of
superbubbles with high significance ($>4\sigma$). Fainter superbubbles were
subsequently found (Chu et al., 1995; Dunne et al., 2001). Studies of other
X-ray superbubbles were presented in Wang & Helfand (1991b); Bomans et al.
(1994); Oey (1996); Burrows & Zhiyu (1996); Chu et al. (2003); Dunne et al.
(2003); Townsley et al. (2006); Pidopryhora et al. (2007); Tüllmann et al.
(2008, 2009); Kavanagh et al. (2012) and Krause et al. (2018), as well as of
superbubble candidates in Long et al. (2010); Sasaki et al. (2012); Sun et al.
(2012) and Kavanagh et al. (2020). The superbubbles that manage to break the
cold gas disk of a galaxy would manifest themselves observationally as H i
holes, when viewed face-on (e.g. Brinks & Bajaja, 1986; Boomsma et al., 2008;
Pokhrel et al., 2020). The size distribution of the observed H i holes largely
follows a power-law distribution with an index of about $-3$ (Oey & Clarke,
1997; Bagetakos et al., 2011). Recently, several works on PHANGS galaxies have
investigated some other superbubble-related structures in other wavelengths
such as molecular shells and H ii regions (e.g. Watkins et al., 2023; Egorov
et al., 2023). Also, the Local Volume Mapper (LVM Kollmeier et al., 2019) is
able to give a very detailed spectroscopic inspection of these structures by
the integral-field spectroscopic survey.

Based on observations, Weaver et al. (1977) constructed an analytical model
for individual self-similar interstellar bubbles, powered by continuous
injection of energy from stellar winds. This model was extended by Mac Low &
McCray (1988), who showed that successive supernovae can also produce
superbubbles and presented a very comprehensive scenario of the dynamical
structure and evolution of superbubbles. In complement to observations and
analytical models, numerical simulations have been carried out on stellar
feedback processes (e.g. de Avillez & Breitschwerdt, 2004; Gent et al., 2013;
Keller et al., 2014; Li et al., 2017; Kim et al., 2017; Marinacci et al.,
2019; Sartorio et al., 2021; Ali et al., 2022; Bieri et al., 2023), and
isolated superbubbles (e.g. Krause & Diehl, 2014; Kim & Ostriker, 2015; Yadav
et al., 2017; Fielding et al., 2018; El-Badry et al., 2019). Simulations by
El-Badry et al. (2019) suggested that both cooling at the interface and
evaporation of shell gas into the interior of the bubble are important in the
evolution of superbubbles. Yadav et al. (2017) performed three-dimensional
simulations of isolated superbubbles and found that most of the input energy
is lost by radiative cooling, and the superbubbles retain less than $10\%$ of
the input energy. Based on 1-d simulations, Nath et al. (2020) found that the
size of superbubbles has a power-law distribution with an index of about
$-2.7$.

Built upon the models for individual superbubbles, McKee & Ostriker (1977)
proposed a self-consistent ISM model (the so-called "three-phase model"),
consisting of cold neutral medium ($T\sim 10^{2}~{}\mathrm{K}$), warm medium
($T\sim 10^{4}~{}\mathrm{K}$) and hot medium ($T\sim 10^{6}~{}\mathrm{K}$).
The three-phase model is reasonably successful in explaining some of the
important observations related to multi-phase ISM (see McKee, 1990, for a
review). Cox & Smith (1974) suggested that a supernova event rate of about 1
per 50 years is sufficient for supernova remnants to overlap and form a
network of "hot tunnels", with the hot gas filling about $10\%$ of the
interstellar volume. On the other hand, McKee & Ostriker (1977) argued that
the overall ISM structure would largely follow the local ISM, implying that
the hot ISM occupies a large fraction of the Galactic disk and acts as a
background phase. In their three-phase ISM model, the volume filling factor of
hot gas can be larger than $70\%$. In observations, it is difficult to
directly measure the amount of hot gas in the Milky Way outside of the Local
Bubble. Attempts were made to estimate the filling factor of hot ISM in nearby
face-on galaxies (Cui et al., 1996), but also turned out to be challenging,
due to limited spatial and spectral resolutions. Recently, simulations show
that the filling factor of hot gas ranges from $17\%\,-\,44\%$ in the ISM of
galaxies (de Avillez & Breitschwerdt, 2004), depending on starburst
activities. However, previously with galactic-scale numerical simulations, no
statistical analysis of superbubbles and their evolution was conducted.

In this work, we study superbubbles systematically using state-of-the-art
galaxy simulations. We are able to resolve individual superbubbles, track
their evolution and perform a systematic study of superbubbles in high-
resolution hydrodynamical simulations based on the SMUGGLE model (Marinacci et
al., 2019). SMUGGLE is a comprehensive, physically-motivated star formation
and feedback model which can be widely used in detailed studies of galaxy
evolution. Li et al. (2020) presented a comprehensive study of giant molecular
clouds (GMCs) in Milky-Way mass galaxies and found the influence of various
model configurations (i.e. different feedback channels and different choices
of star formation efficiency). This model is also widely used in the studies
of young massive clusters (Li et al., 2022), black hole fueling (Sivasankaran
et al., 2022), Lyman-$\alpha$ properties (Smith et al., 2022) and feedback in
dwarf galaxies (Burger et al., 2022a, b; Jahn et al., 2023).

The outline of the paper is as follows. In Section 2, we describe the ISM
model adopted, the algorithms developed to identify superbubbles in 3-d maps
and H i holes in column density maps, and the methods used to find embedded
stellar associations and track the evolution of superbubbles. In Section 3,
the results about properties and evolution of superbubbles are presented. We
discuss the implications of the results and open questions in Section 4, and
briefly summarize the work in Section 5.

##  2 Methods

Table 1: Summary of the model configurations in the simulations we use. The "SN" ("Radiative and wind") column tells whether the supernova (radiative and stellar wind) feedback is considered. $\epsilon_{\rm ff}$ is the star formation efficiency defined in Equation (1). Name | $\epsilon_{\rm ff}$ | SN | Radiative and wind
---|---|---|---
SFE1 | 0.01 | Yes | Yes
SFE10 | 0.1 | Yes | Yes
SFE100 | 1 | Yes | Yes
SN | 0.01 | Yes | No
Rad | 0.01 | No | Yes

We analyze a set of hydrodynamic simulations of isolated Milky Way-like
galaxies presented in Li et al. (2020). The simulations were performed by the
moving-mesh hydrodynamic code Arepo (Springel, 2010) with a physically-
motivated star formation and feedback sub-grid model, SMUGGLE (Marinacci et
al., 2019). Star particles are formed in cold and dense molecular gas with a
rate $\dot{M}_{\star}$ defined as

$\displaystyle\dot{M}_{\star}$ $\displaystyle=\epsilon_{\rm ff}\;\frac{M_{\rm
gas}}{t_{\rm ff}},$ (1)

where $\epsilon_{\rm ff}$ is the local star formation efficiency, $M_{\rm
gas}$ is the mass of the gas cell, and $t_{\rm
ff}=\sqrt{\frac{3\pi}{32G\rho_{\rm gas}}}$ is the free-fall time of the gas
cell with density $\rho_{\rm gas}$. The simulations reach a minimum softening
length of $3.6\,\rm pc$ and mass resolution of about $1.4\times 10^{3}\,\rm
M_{\sun}$, which are thus able to spatially resolve the structure of
superbubbles around stellar associations. To investigate the time evolution of
superbubbles, we output the simulation snapshots every 1 Myr. We analyze
simulations with different sub-grid models and parameters, as summarized in
Table 1, to investigate how the properties of superbubble changes with galaxy
formation physics. The SFE1 run is the fiducial run, which considers all
stellar feedback channels with a star formation efficiency $\epsilon_{\rm
ff}=0.01$. The SFE10 and SFE100 runs change $\epsilon_{\rm ff}$ to 0.1 and 1,
respectively. The Rad run considers radiative and stellar wind feedback, while
the SN run considers only supernova feedback. We refer the readers to Li et
al. (2020) for more details on the initial conditions and simulation setup. We
identify and analyze the properties of superbubbles by applying the method
described in the next subsection to the simulation snapshots between 0.25 Gyr
and 1 Gyr.

###  2.1 3D superbubble identification

(a) SFE1

(b) SFE10

(c) SFE100

(d) SN

(e) Rad

Figure 1: Panel (a): Identification of superbubbles (left) and HI holes
(right) in the SFE1 run. The left panel shows the distribution of X-ray
emission (0.1 - 2.4 keV) at 0.4 Gyr for the SFE1 run in a face-on projection,
with an extension of $40\times 40\;\rm kpc$. White circles correspond to the
identified three-dimensional superbubbles, which are defined as hot, X-ray
luminous structures. The positions and radii of the circles represent the
projected location and the half-mass radii of the corresponding superbubbles.
The right panel shows the face-on HI surface density map at 0.4 Gyr for the
SFE1 run with an extension of $40\times 40\;\rm kpc$. Circles show the
identified two-dimensional HI holes. Panel (b)-(e): similar to panel (a) but
for SFE10, SFE100, SN and Rad run.

Superbubbles are structures filled with hot gas blown out by stellar winds
(Weaver et al., 1977) or supernovae (Mac Low & McCray, 1988). In the
simulations, we consider superbubbles as structures of hot gas, and exclude
structures that are outside the cold galactic disk and flow into the CGM.
Superbubbles are born from stellar feedback in the cold ($T<10^{3.5}$ K) disk.
If the feedback process is strong enough, it may break the cold disk and blow
out. Sometimes the blown-out gas may cut off the connection to the cold disk
if feedback has already stopped or the cold disk is disturbed by some
surrounding activities. We study hot gas within 2 kpc above and below the mid-
plane, a region mainly composed of the cold disk and surrounding warm
($10^{3.5}<T<10^{5.5}$ K) ionized gas. Beyond this region, the environment
consists mostly of very diffuse hot gas, which can interfere with our
superbubble finding algorithm.

We first select gas cells within $2~{}\mathrm{kpc}$ above and below the mid-
plane. We then find neighbours of each gas cell in order to prepare for the
identification of structures. In the framework of Arepo, the adaptive, dynamic
Voronoi mesh is defined by the Voronoi tessellation of mesh-generating points.
Therefore we only need to perform its topological dual Delaunay triangulation
to the three-dimensional simulation output, and find the connected mesh-
generating points to locate the neighbours of each Voronoi cell. After
neighbours are found, we take gas cells with temperature
$T>10^{5.5}~{}\mathrm{K}$ and density $\rho>10^{-27}~{}\mathrm{g/cm^{3}}$.
This kind of gas is expected to efficiently emit X-rays, providing a
convenient tracer for the identification of superbubbles. Neighbouring gas
cells satisfying the above selecting criteria will be assigned to a structure,
which is considered as a superbubble candidate. Finally we need to exclude
some unreliable structures to ensure our identification of superbubbles. We do
not consider structures containing less than 5 gas cells or with radii smaller
than $10~{}\mathrm{pc}$, and also exclude structures without stellar
associations inside as they probably belong to some fragmented superbubbles or
they are just parts of the hot, low-density background.

For each identified superbubble, we define its center as the center of mass
and its effective radius $r_{\rm{bubble}}$ as the half-mass radius of all the
cells within. We have also tried using an ellipsoid to represent each
superbubble. We estimate the semi-principal axes of the ellipsoid ($a$, $b$
and $c$ ) of each superbubble using the shape tensor (e.g. Zemp et al., 2011),
and calculate its effective radius as $R_{\rm ellipsoid}=(abc)^{1/3}$. It
turns out that the results from these two methods agree well with each other,
so we choose to use the half-mass radius to represent the size of superbubbles
in this work for simplicity. On the left part of panel (a)-(e) in Fig. 1, we
show the location and radius of each identified superbubble over the the X-ray
emissivity map in $0.1-2.4~{}\mathrm{keV}$ band from hot
($T>10^{5.5}~{}\mathrm{K}$) and optically thin ($n<0.3~{}\mathrm{cm^{-3}}$)
gas with the APEC model (Smith et al., 2001), which assumes collisional
ionization equilibrium. The size of superbubbles can range from 10 pc to 1
kpc, with a typical size of about several hundred pc.

###  2.2 Finding stellar associations located in superbubbles

In observations, the existence of OB associations offers important guidance to
the identification of superbubbles (e.g. Chu & Mac Low, 1990). Moreover, the
properties of stellar associations, such as their mass, provide valuable
constraints on the formation and evolution of the superbubbles. In our
simulations, a stellar association is denoted as a star particle which
represents a coeval simple stellar population (Marinacci et al., 2019)
following the initial mass function presented by Chabrier (2001). As gas is
denoted by Voronoi cells, we match each star particle to the gas cell where it
lies in, so we need to find the nearest mesh-generating point of that star
particle.

We use the k-dimensional tree (KD-Tree) in SciPy (Virtanen et al., 2020) to
find nearest points. We consider only stellar associations younger than $100$
Myr, as older stellar associations are not likely to contribute feedback any
more. Some gas structures found in Section 2.1 may not contain stellar
associations if they are parts of the hot background or belong to some
fragmented superbubbles. In either case, we do not regard them as
superbubbles. The mass of stellar associations have some important effects on
the properties of superbubbles, and the related results are presented in
Section 3.3.

###  2.3 Tracing the evolution of superbubbles

Figure 2: An example shows how we trace the superbubbles between snapshots.
Superbubbles $A$ and $B$ are identified in snapshot $t$, while $C$ and $D$ are
identified in snapshot $t+1$. The total shared stellar mass between $A$ and
$C$, $B$ and $C$, $B$ and $D$ are defined as $m_{\rm AC},m_{\rm BC}$ and
$m_{\rm BD}$. In this example we take $m_{\rm AC}>m_{\rm BC}>m_{\rm BD}$. $B$
is fragmented into two parts. One becomes $D$, and the other merges with $A$
into $C$.

Tracing the movement of gas is not trivial in a Voronoi mesh, where gas is
always flowing across neighbouring Voronoi cells. Additionally, refinement and
de-refinement will further complicate the analysis by removing or adding cells
to the existing mesh.

Here we take the star particles in each superbubble as tracers. We assume that
stellar associations will not have significant movements relative to the
corresponding superbubbles, at least they will not totally leave the original
superbubbles. This assumption seems quite reasonable and is also a basic
assumption in observations to identify superbubbles, where most of the
observed structures have their corresponding OB associations detected (e.g.
Chu & Mac Low, 1990).

In Section 2.1, we have already identified superbubbles in each snapshot of
the simulations, and what we need to do is to match the superbubbles with
their progenitors and descendants in former and latter snapshots. The match
follows a methodology similar to Hopkins et al. (2012), but is not exactly the
same, considering the features of our simulations and objectives. In general,
several nearby superbubbles may merge into a bigger one, and when a big
superbubble is dying, it may split into several smaller structures. For a
superbubble in snapshot $t$, we define its progenitors (descendants) as
superbubbles that share at least one star particle in snapshot $t-1$ ($t+1$).
The main progenitor (descendant) of a superbubble is defined as the progenitor
(descendant) containing the highest amount of shared stellar mass. In snapshot
$t$, a superbubble is born if there is no progenitor in snapshot $t-1$, and
dies if there is no descendant in snapshot $t+1$. This approach allows us to
establish connections among superbubbles across different snapshots, forming
an evolutionary network with diverse evolutionary branches.

Fig. 2 shows an example, in snapshot $t$ we identify superbubbles $A$ and $B$,
while in snapshot $t+1$ we identify superbubbles $C$ and $D$. Superbubbles are
linked with an arrow if they have some star particles in common. The total
stellar mass in common between $A$ and $C$, $B$ and $C$, $B$ and $D$ are
defined as $m_{\rm AC},m_{\rm BC}$ and $m_{\rm BD}$. In this example we take
$m_{\rm AC}>m_{\rm BC}>m_{\rm BD}$. Superbubble $B$ is fragmented into two
parts. One becomes superbubble $D$, and the other merges with superbubble $A$
into superbubble $C$. Therefore, superbubble $C$ has two progenitors $A$ and
$B$, while $A$ is its main progenitor. Superbubble $B$ has two descendants $C$
and $D$, while $C$ is its main descendant. As superbubble $A$ is the main
progenitor of superbubble $C$ and also superbubble $C$ is the main descendant
of superbubble $A$, the evolutionary branch from $A$ to $C$ is defined as a
main branch. We note that we do not consider superbubble $D$ as a newly born
superbubble, because it is just a piece of the previously fragmented
superbubble. Also, our main results are based on comparing the influence of
different simulation configurations on superbubble properties. Using the
identical analysis method on different simulation runs provides us a simple
yet robust way of identifying superbubbles. Basic results on the size
evolution of the main-branch superbubbles are shown in Section 3.4. A brief
comparison to the analytical model will be discussed in Section 4.5.

###  2.4 2D H i hole identification

Superbubbles triggered by feedback from their natal stellar associations can
sometimes break the disk of cold gas above and below the mid-plane and
puncture holes that are visible in H i observations. In observations, we
should note first that several studies rely on identification of H i holes "by
eye" (e.g. Bagetakos et al., 2011; Pokhrel et al., 2020), while there are some
attempts trying to identify H i holes automatically. For instance, Thilker et
al. (1998) attempted identification using an expanding model, which is based
on the superbubble model given by Weaver et al. (1977). H i holes were
identified with spatial and velocity information in the observed H i data. The
algorithm was refined in Mashchenko et al. (1999) to consider simulation
results of three-dimensional H i shells. This kind of algorithms are only able
to study structures with regular shapes, and do not work well when structures
become non-spherical or incomplete. A model independent hole-searching
algorithm was developed by Ehlerová & Palouš (2005). It was first used to find
high-velocity clouds in de Heij et al. (2002). The basic idea of this
algorithm is to find local minima (or maxima) as the seeds of structures, then
expand outward to take in brighter (or fainter) neighbour pixels. Based on
similar ideas, here we adopt a tool, astrodendro, to identify H i holes.

astrodendro 111http://www.dendrograms.org/ (Rosolowsky et al., 2008) is a
dendrogram-based algorithm, which is widely used in finding extended sources,
such as giant molecular clouds (GMCs) and star-forming regions, for both
observational (e.g. Cosens et al., 2022) and simulations (e.g. Li et al.,
2020) datasets. astrodendro turns the input map into a dendrogram, namely a
tree that represents the hierarchy of the structures on the input map. The
algorithm stops at a value $\sigma_{\rm base,min}$, which is mainly used to
account for the background noise, and values below this threshold will not be
considered. Another parameter $\sigma_{\rm delta,min}$ indicates a minimum
significance for the structures in order to exclude "spikes" induced by
fluctuations. The last parameter, $N_{\rm pix,min}$, represents a minimum
pixel number of a structure to be identified. This is to avoid some spurious
structures induced by the limitation of spatial resolution.

To identify H i holes in simulation snapshots, we first generate the surface
density map of the neutral hydrogen gas projected along the z-axis, namely a
face-on projection of the gaseous disk. We create the H i surface density map
for the central $40\times 40$ $\rm kpc^{2}$ region for simulations snapshots
between 0.25 Gyr and 1 Gyr for all five runs. The map is resolved with
$10000\times 10000$ pixels, with the size of each pixel corresponding to 4 pc,
which is approximately the limit of the adaptive softening of the simulations.
astrodendro is developed to find peaks, but in our work we seek to find holes
(valleys) in HI surface density ($\sigma_{\rm\text{H\,{i}}}$) map. Therefore,
we take the inverse of our HI surface density map so that HI holes become
peaks in the inverse map ($\frac{1}{\sigma_{\rm\text{H\,{i}}}}$). In practice,
the value of HI surface density varies over six orders of magnitudes, so we
feed the values of surface density map in logarithms ($-{\rm
log}\,\sigma_{\rm\text{H\,{i}}}$) into astrodendro. We set the threshold
$\sigma_{\rm base,min}=-19.5$ in the $-{\rm log}\,\sigma_{\rm\text{H\,{i}}}$
map, which means a surface density upper limit of $10^{19.5}\,\rm cm^{-2}$ in
the $\sigma_{\rm\text{H\,{i}}}$ map. $\sigma_{\rm delta,min}$ is the minimum
height difference between peaks, and it is also the minimum height of one peak
to be recognized. Parameter $\sigma_{\rm delta,min}=0.5\;\rm dex$ in the
$-{\rm log}\;\sigma_{\rm\text{H\,{i}}}$ map means that we find peaks and their
minimum heights are 0.5 dex. Therefore, in the $\sigma_{\rm\text{H\,{i}}}$ map
we find valleys (local minima), and their minimum heights are $10^{0.5}(\sim
3.162)$ times the local minima. Finally, we take $N_{\rm pix,min}=10$ in this
work to avoid potential spurious structures.

The right part of panel (a)-(e) in Fig. 1 show the identified H i holes. In
general, H i holes are located between the filamentary structures. Comparing
to the left part of panel (a)-(e) of Fig. 1, some H i holes have co-spatial
superbubble counterparts, but some do not. The relationship between
superbubbles and H i holes is complicated, and will be discussed in detail in
Section 4.1.

##  3 Results

Table 2: Catalogue of X-ray identified superbubbles, organized by host galaxy.
For each superbubble, we summarize its name, radius $r$, X-ray luminosity
$L_{\rm X}$ and the corresponding OB association. References for each observed
superbubble are also given.

NO. | Name |  $r$ (pc) | log $L_{\rm X}$ (erg / s) ${\rm a}$ | OB Association | Reference | Host Galaxy
---|---|---|---|---|---|---
1 | Local Bubble | 165.0 | 36.08 |  UCL / LCC ${\rm b}$ |  Snowden et al. (1990); | MW
Zucker et al. (2022)
2 | Scorpius–Centaurus superbubble | 120.0 | 35.30 | USco ${\rm b}$ | Krause et al. (2018)
3 | Cygnus superbubble | 225.0 | 36.70 (36.64) | Cygnus OB2 | Cash et al. (1980)
4 | Orion-Eridanus superbubble | 113.0 | 37.61 | Orion OB1 | Burrows & Zhiyu (1996)
5 | N44 ${\rm c}$ (DEM 152 $\rm d$ ) | 27.4 | 36.15 | LH 47 $\rm e$ , LH 48 |  | LMC
6 | N51D (DEM 192) | 53.0 | 35.52 | LH 51, LH 54 |  Chu & Mac Low (1990);
7 | N57A (DEM 229) | 50.0 | 34.87 | LH 76 | Wang & Helfand (1991a)
8 | N70 (DEM 301) | 55.0 | 35.26 | LH 114 |
9 | N9 (DEM 31) | 54.5 | 35.00 (35.00) | LH 6 | Chu et al. (1995)
10 | N30A (DEM 105) | 63.2 | 35.28 (35.28) | LH 37
11 | N30C (DEM 106) | 34.2 | 34.95 (34.95) | LH 38
12 | DEM137 | 85.0 | 35.15 (35.15) | LH 43
13 | N11 shell 1 (DEM 34) | 61.2 | 37.17 (37.90) | LH 9 | Dunne et al. (2001)
14 | N103B (DEM 84) | 60.0 | 34.85 (35.00) | NGC 1850
15 | N105 (DEM 86) | 36.7 | 35.40 (35.62) | LH 31
16 | N144 (DEM 199) | 47.4 | 34.74 (34.95) | LH 58
17 | N154 (DEM 246) | 73.5 | 35.52 (35.75) | LH 81, LH 87
18 | N158 (DEM 269) | 54.8 | 35.85 (36.00) | LH 101, LH 104
19 | N160 (DEM 284) | 82.2 | 35.80 (35.95) | LH 103
20 | N206 (DEM 221) | 45.0 | 35.09 (35.31) | LH 66, LH 69
21 | 30 Doradus C | 47.4 | 36.00 (36.19) | LH 90
22 | N185 (DEM 25) | 43.0 | 35.26 | - | Oey (1996)
23 | N186 CE (DEM 50) | 50.0 | 35.62 | -
24 | LMC-2 $\rm f$ | 500.0 | 37.30 (37.33) | - | Wang & Helfand (1991b)
25 | NGC 604 | 122.5 | 36.16 (36.28) | - | Tüllmann et al. (2008) | M33
26 | IC 131 | 100.0 | 36.09 (36.43) | - | Tüllmann et al. (2009)

  * *

We note that some superbubbles are identified in other wavelengths and do not
have clear identifications in X-ray, such as the Ophiuchus superbubble
(Pidopryhora et al., 2007) and the Perseus-Taurus superbubble (Bialy et al.,
2021).

  * ${\rm a}$

The energy range used for the calculation of luminosities may vary between
different observations. Corrected values for 0.1 - 2.4 keV band are shown in
brackets for observations with available spectral fitting models.

  * ${\rm b}$

There are three subgroups of the Scorpius–Centaurus association: Upper
Centaurus Lupus (UCL), Lower Centaurus Crux (LCC) and Upper Scorpius (USco).

  * ${\rm c}$

Following the name used in Henize (1956).

  * ${\rm d}$

Named after Davies et al. (1976).

  * ${\rm e}$

Named after Lucke & Hodge (1970).

  * ${\rm f}$

Following the name used in Meaburn (1980).

###  3.1 X-ray luminosity – size distribution of superbubbles

Figure 3: X-ray luminosity – size distribution of superbubbles from
simulations and observations. Each colored dot shows an identified superbubble
in the SFE1 run, with its color representing the embedded stellar association
mass. Dashed lines show the number density contours of 68%, 95% and 99.7%. The
Local Bubble and the eROSITA Bubble are denoted as red star and square. The
X-ray identified superbubbles in the Milky Way, LMC and M33 are labeled with
(I), while some candidates in M33, M31 and M101 with (C). A few examples of
central diffuse emission (CDE) in other galaxies are shown for comparison.

X-ray observations provide direct information on superbubbles. Here, we
summarize the key properties of the observed superbubble candidates in X-ray,
and briefly compare them with those identified in our simulations. However,
due to limitations in spatial resolution and sensitivity in the soft X-ray
band, the observed sample size is relatively small, with over half of the
identified superbubbles located in the LMC. In Table 2, We have compiled the
name, radius, luminosity of each superbubble identified with confidence in
X-ray. Most of these bubbles have co-spatial OB associations detected. We
should note that in some of the observations the superbubbles are considered
elliptical in shape. In those cases, we compute the effective radius as
$\sqrt{ab}$, where a and b are the semi-major and semi-minor axes of the
ellipse. Generally, the observed superbubbles have radii ranging from about 10
pc to several hundred pc and X-ray luminosities from
$10^{34}~{}\mathrm{erg\,s^{-1}}$ to $10^{37}~{}\mathrm{erg\,s^{-1}}$. It is
important to note that these values were derived from data taken with
different observatories, so the energy ranges used for computing luminosity
may also be different. To address this issue, we have re-computed luminosity
for the energy range of 0.1 - 2.4 keV, for cases in which spectral fitting
results are available, and the results are shown in perentheses. Additionally,
different studies may adopt different algorithms to account for the
foreground. These modelling uncertainties are expected to influence the
luminosity values to some extent, but are unlikely to qualitatively affect the
trend. Some superbubble-like candidates also exist in the nearby galaxies. For
reference, we summarize the published results on those that were found in M31,
M33 and M101 in Table 3.

As described in Section 2.1, our identification of superbubbles is in
principle able to find some kpc-scale giant bubbles breaking out from the
galactic center, which have structures similar to the Fermi Bubbles (Su et
al., 2010) and the eROSITA Bubbles (Predehl et al., 2020). As mentioned
earlier, such diffuse emissions can arise from stellar feedback or AGN
feedback, or a combination of both. They can be classified as superbubbles if
the central super massive black hole is not active. For comparison, Table 4
show the published results (likely incomplete) on the diffuse emission from
the central region of a galaxy. These structures are usually large (exceeding
$400~{}\mathrm{pc}$ in radius) and bright (with luminosity reaching
$10^{38}~{}\mathrm{erg\,s^{-1}}$ or higher in the 0.1 - 2.4 keV band), with
large variation in size and luminosity, presumably due to their diverse
galactic environment and intrinsic properties.

Fig. 3 shows a luminosity-size diagram for the superbubbles identified in the
fiducial SFE1 run, along with identified X-ray superbubbles and candidates. We
also show our compilation of galactic central diffuse emissions (CDEs, yellow
asterisks). The Local Bubble and the eROSITA Bubbles are highlighted as red
star and red square, respectively. Note that we calculate the X-ray luminosity
of simulated superbubbles using the APEC model (Smith et al., 2001) provided
by AtomDB (Foster et al., 2012), which assumes optically-thin gas,
consequently we have excluded gas with density $n>0.3~{}\mathrm{cm^{-3}}$,
which lies in the optically-thick regime. For a better comparison with the
observations, we specifically employ the half-luminosity radius, as opposed to
the half-mass radii, to represent the size of superbubbles in Fig. 3. The
identified superbubbles in our simulations vary greatly in size, spanning from
10 pc to several kpc. The luminosity $L_{\rm X}$ varies from $10^{34}$ to
$10^{40}~{}\mathrm{erg\,s^{-1}}$, showing a positive correlation with
superbubble size and the mass of stellar association. Perhaps, larger mass of
stellar association generally imply stronger feedback, which leads to
superbubbles of larger size and higher luminosity. The X-ray luminosity-size
distribution of superbubbles in our simulations is largely consistent with
observations.

###  3.2 Size distribution of superbubbles and H i holes

Figure 4: Left panel: Size distribution of the identified 2D H i holes in
different runs. Right panel: Size distribution of the identified 3D
superbubbles in different runs.

As described in Section 2.1 and 2.4, we identify superbubbles from the 3D
distribution of gas cells, and H i holes from two-dimensional H i maps. The
size distribution of superbubbles is presented on the right panel of Fig. 4,
showing complexity and variation across different runs. Interestingly, in the
Rad run, where supernova feedback is not considered, bubbles of hot gas can
still form due to stellar winds and the radiation emitted by young massive
stars and star clusters. When comparing the SN run to the Rad run, we find
that superbubbles generated by supernova feedback tend to be much larger in
size compared to those resulting from radiative feedback and stellar winds.
Both the Rad run and the SN run exhibit a single-peak size distribution, with
the Rad run peaking at tens of parsecs which represents the size of wind-
blowing bubbles, and the SN run peaking at several hundred parsecs which
represents the size of SN-driven bubbles. The remaining three runs exhibit
double-peak distributions. Drawing on the findings from the Rad and SN runs,
we infer that one is dominated by radiative and stellar wind feedback (at a
smaller size), while the other peak by supernovae (at a larger size). By
comparing SFE1, SFE10, and SFE100, we find that higher star formation
efficiency results in a higher proportion of large superbubbles. We will see
in Section 3.3 and Fig. 6, that the larger the SFE, the larger the size of
superbubbles dominated by early feedback, so more bubbles originated from
early feedback will blend into the larger peak, making it appear like the
larger peak is more prominent. The detailed impact of star formation
efficiency will be discussed in Section 4.2.

In principle, lower SFE should in general form smaller stellar associations,
and should peak at a smaller size because of weaker feedback. However, as we
will see from Section 3.3 that supernova-dominated (SN-dominated) superbubbles
are not very sensitive to $\epsilon_{\rm ff}$. Moreover, in the simulations,
the mass resolution of a star particle is $10^{3}\;\rm M_{\odot}$. Star
particle with $10^{3}\;\rm M_{\odot}$ can generate about 6.5 SN in 10 Myr, and
lead to a bubble with minimum radius of about 155 pc at 10 Myr for typical
ambient density (see, e.g. the ambient density estimated in the SFE1 run as
shown in Section 4.2), so the place of larger peak might be limited by the
resolution of star particle mass. With higher mass resolution, possibly we can
see more SN-dominated superbubbles with smaller size.

The left panel of Fig. 4 shows the size distributions of H i holes in
different simulation runs. The size distributions in runs with supernovae
feedback are similar, ranging from $10~{}\mathrm{pc}$ to a kpc, with a peak at
a radius of approximately $100~{}\mathrm{pc}$. However, in the Rad run, most
of the identified H i holes are considerably larger. This difference can also
be seen from the H i surface density maps in Fig. 1. They are probably low
density voids between spiral arms/filaments and are not related to stellar
feedback. It is similar to the observed HI voids in grand-design galaxies
(like M51) (e.g. Schuster et al., 2007). Early feedback alone is not powerful
enough to frequently destroy the filaments and push gas to the inter-arm
regions. This is also consistent with the two-point correlation function of
giant molecular clouds analysed in Li et al. (2020).

The distribution of H i hole size is very different from that of superbubbles,
and shows only one peak in any run. The lack of the peak at smaller size might
reflect the fact that H i is not very sensitive to early feedback, because the
processes are not as powerful as supernovae in pushing away the surrounding
gas. Theoretically, the size distribution of H i holes could be expressed as
$N(R)\propto R^{1-2\beta}$ when assuming a power-law mechanical luminosity
function $\phi(L)\propto L^{-\beta}$ for OB associations and a constant star
formation rate (Oey & Clarke, 1997). Typically in observations, the size
distribution is fit with a power-law only for the larger sizes (i.e.,
$R_{\rm{hole}}>100~{}\mathrm{pc}$), taking into account that the spatial
resolution of observations may not be sufficient to resolve H i holes smaller
than $100~{}\mathrm{pc}$ (Oey & Clarke, 1997; Bagetakos et al., 2011). It is
important to note that the larger size part does not seem to follow a simple
power-law distribution in either simulations or observations. Furthermore, we
find that the probability of finding H i holes smaller than
$100~{}\mathrm{pc}$ decreases even when the spatial resolution is sufficient
to resolve them in our simulations. This tendency seems in line with the
observations (Oey & Clarke, 1997; Bagetakos et al., 2011). The peak position
is still limited by the resolution of star particle mass, as we mentioned
above, even though the spatial resolution of the simulations is much higher
than the observations.

###  3.3 Effects of stellar association mass in superbubbles

Figure 5: Distribution of superbubble density and stellar association mass for
the SN (orange) and Rad (green) run. The horizontal dashed line shows the
classification criterion $\rho_{\rm bubble}=10^{-25.223}{\rm g\;cm^{-3}}$.
Figure 6: The effects of stellar association mass on the superbubble size (top
panel) and density (middle panel). Two types of superbubbles (supernova-
dominated and early feedback-dominated) are shown separately. Different runs
are indicated with different colors. Solid lines show the median value of size
and density at a particular mass of stellar association. Shaded regions
corresponding to the $40^{\rm{}th}-60^{\rm{}th}$ percentiles. The bottom panel
shows the distribution of stellar association mass in different runs.

A higher mass of stellar associations typically results in stronger stellar
feedback processes that contribute to the formation of superbubbles. This
factor is crucial in studying the properties of superbubbles in both
simulations and observations. Moreover, the influence of stellar association
mass reveals distinct variations across different runs, and can offer some
valuable insights into the effects of star formation efficiency.

Firstly, the distribution of stellar association mass is shown as the bottom
panel of Fig. 6. The SN run has more massive stellar associations than other
runs, but once early feedback is considered, the mass of stellar association
is lowered. Therefore, early feedback plays a very important role on
decreasing the mass of stellar associations, by suppressing the formation of
star particles in the early stage, prior to supernova explosions. We will
discuss more about this in Section 4.3.

In the SFE1, SFE10 and SFE100 run, there exist two types of superbubbles,
which are dominated by supernova and early feedback. Fig. 5 shows the
distribution of superbubble density and stellar association mass for SN and
Rad runs. We can see that density of superbubble $\rho_{\rm bubble}$ is a
pretty good property to classify these two types of superbubbles. We use
logistic regression and find $\rho_{\rm bubble}=10^{-25.223}{\rm g\;cm^{-3}}$
is a good criterion with 96% accuracy. We applied this criterion onto all
other runs and show in Fig. 6 the results on the effects of stellar
association mass for these two types of superbubbles. We note that the
classification of these two types of superbubbles is not perfect. In Fig. 6,
the larger the $\epsilon_{\rm ff}$, the larger the size of superbubbles
dominated by early feedback, so more bubbles originated from early feedback
will blend into the supernova-dominated type.

The top panel of Fig. 6 illustrates the influence of stellar association mass,
$M_{\rm{star}}$, on the size $r_{\rm{bubble}}$ of the superbubbles. For a
particular stellar association mass, superbubbles dominated by supernova
feedback tend to be larger than those dominated by early feedback. A positive
correlation between $r_{\rm bubble}$ and $M_{\rm star}$ can be seen for early
feedback-dominated superbubbles, and also for supernova-dominated superbubbles
with stellar association mass larger than $10^{5}\;{\rm M_{\odot}}$. However,
the relation seems to become flat for supernova-dominated superbubbles when
stellar association mass is lower than $10^{5}\;{\rm M_{\odot}}$. This is
possibly limited by the mass resolution of gas ($10^{3}\;{\rm M_{\odot}}$) and
the density criterion when classifying superbubbles (density lower than
$10^{-25.223}{\rm g\;cm^{-3}}$). Together with the minimum cell number we set,
which is 5, when identifying superbubbles in Section 2.1, we can get a lower
limit of radius as 178 pc (log 178 $\sim$ 2.25) for supernova-dominated
superbubbles.

The middle panel of Fig. 6 provides information on the effects of stellar
association mass on the density of superbubbles (estimated as the averaged
density of its gas cells). Firstly, supernova-dominated superbubbles keep low
densities (approximately $10^{-26}\,{\rm g\;cm^{-3}}$) across a wide range of
stellar association masses, from $10^{3}\,\mathrm{M_{\sun}}$ to
$\mathrm{10^{7}\,M_{\sun}}$. This may correspond to a typical density for
overlapping supernovae. In contrast, early feedback-dominated superbubbles
show significantly higher densities (above $10^{-24}\,{\rm g\;cm^{-3}}$),
likely due to the relatively weaker dynamic effects of radiative and stellar
wind feedback compared to supernovae. As stellar association mass increases,
radiative and stellar wind feedback gain more strength to disperse the
surrounding interstellar medium (ISM), resulting in lower superbubble
densities. It is worth noting that the density $10^{-24}\,{\rm g\;cm^{-3}}$
largely corresponds to the typical density of the main interstellar medium
(ISM) component, specifically phase (a) as shown in Fig. 10. This indicates
that early feedback-dominated superbubbles with lower stellar association
masses primarily form thermally by the heating of photoionization, together
with radiation pressure and stellar winds dispersing the original surrounding
dense clouds to a density similar to the environment ISM density. Since early
feedback-dominated superbubbles lack the significant power to sweep away the
surrounding gas, they are primarily confined to the cold disk.

Fig. 6 shows that star formation efficiency $\epsilon_{\rm ff}$ does not have
noticeable effect on supernova-dominated superbubbles, while superbubbles
dominated by early feedback generally become larger and more diffuse when
$\epsilon_{\rm ff}$ is higher. This indicates that $\epsilon_{\rm ff}$ mainly
influences the processes of radiative feedback and stellar winds. In general,
higher $\epsilon_{\rm ff}$ leads to lower environmental density. Radiative
feedback is very sensitive to density $\rho_{\rm amb}$ and therefore this kind
of superbubbles will be very different in size. Section 4.2 will discuss this
in more detail.

###  3.4 Evolution of superbubbles

Figure 7: The evolution of superbubble size (solid lines) in different runs,
which are shown in different colors, with a temporal resolution of 1 Myr.
Shaded regions correspond to the $40^{\rm{}th}-60^{\rm{}th}$ percentiles.

Based on the method described in Section 2.3, we use the star particles (each
of which represents a stellar association) to link superbubbles in different
snapshots and construct an evolutionary network including merging and
splitting of superbubbles between snapshots. In order to study the size
evolution of superbubbles qualitatively, for simplicity, we consider only main
branches of the evolutionary network. As defined in Section 2.3, the main
progenitor (descendant) of a superbubble in snapshot $t$ is the superbubble
containing the highest amount of shared stellar mass in snapshot $t-1$
(snapshot $t+1$). The evolution of main-branch superbubble size in each run is
presented in Fig. 7, with a temporal resolution of 1 Myr.

We can see that superbubbles in the SN run are larger than those in the other
runs at the beginning. In the SN run, the lack of early feedback means the
molecular clouds have a lot more time to grow larger. That will lead to more
massive stellar associations and stronger feedback will happen on the birth of
superbubbles.

Comparison across the SFE1, SFE10 and SFE100 runs reveals that star formation
efficiency $\epsilon_{\rm ff}$ impact only the first few million years in the
growth of superbubbles. This is consistent with the result shown in Section
3.3, which is that $\epsilon_{\rm ff}$ mainly influences the processes in
early feedback. In the beginning of superbubbles’ life, they might be
dominated by early feedback. Once supernova feedback starts, they may transit
into supernova-dominated superbubbles, and therefore, the influence of
$\epsilon_{\rm ff}$ gradually wanes. The possible transition from early
feedback-dominated type to supernova-dominated type is also why it is
difficult to separately study the evolution of these two types of
superbubbles. Further discussion about the effects of star formation
efficiency is in Section 4.2.

###  3.5 Volume filling factor of the multi-phase gas

Figure 8: Spatial distribution of volume filling factor of different ISM
phases. The interstellar gas is categorized into cold ($T<10^{3.5}$ K), warm
($10^{3.5}~{}\mathrm{K}<T<10^{5.5}$ K) and hot ($T>10^{5.5}$ K) phases. Left:
Radial profiles of hot gas volume filling factor in different runs. Right:
Radial filling factor profiles of cold (dash), warm (dashdot) and hot (solid)
gas in the SFE1 run, with an extension to 20 kpc from the galactic center.

In the SFE1 run, considering the gas within a cylinder centered at the
galactic center, with a radius of $20~{}\mathrm{kpc}$ and a total height of
$4~{}\mathrm{kpc}$, the overall volume filling factor of hot gas
($T>10^{5.5}~{}\mathrm{K}$) is found to be $12\%$, which is similar to the
findings of de Avillez & Breitschwerdt (2004). In fact, the spatial
distribution and filling factor of the multi-phase gas are not uniform. We
divide the interstellar gas into cold ($T<10^{3.5}$ K), warm
($10^{3.5}~{}\mathrm{K}<T<10^{5.5}$ K) and hot ($T>10^{5.5}$ K) phases (a more
detailed phase diagram is shown in Appendix B). To investigate the radial
variation of the filling factor, we maintain a constant thickness of $4$ kpc
for the cylinder and divide it into multiple "rings" with different
galactocentric radii. The non-uniform distribution of hot gas primarily arises
from the non-uniform distribution of young stellar associations, as they are
typically associated with stellar feedback.

The left panel of Fig. 8 shows the volume filling factors of hot gas in
different "rings" and runs. The volume filling factor of hot gas peaks at
$80\%$ in the galactic center and gradually decreases with radial distance,
without much variation across different runs. The right panel of Fig. 8
presents the radial profiles of the filling factors for different gas phases
in the SFE1 run. Cold gas occupies only a small fraction of space, with hot
gas predominantly filling the volume near the galaxy’s center, while the warm
gas phase dominates in outer regions.

##  4 Discussion

###  4.1 Relationship between superbubbles and H i holes

Feedback from stellar associations result in a multi-phase structure of their
surrounding environment. The enormous ejected energy and momentum give birth
to hot, X-ray emitting superbubbles. The neutral gas in the vicinity is heated
and ionized, leading to the formation of warm-gas shells (which have been
observed, e.g., Dunne et al., 2001), beyond which lies the neutral
environment. Feedback with sufficient power can break through the disk of
neutral gas above and below the mid-plane and puncture holes that are visible
in the view of H i. Therefore, one might envision the co-spatial existence of
stellar associations, X-ray superbubbles and H i holes. In reality, however,
the conditions are much more complex.

The right part of panel (a)-(e) in Fig. 1 show the identified H i holes.
Comparing to the left part of panel (a)-(e) of Fig. 1, we can see that some H
i holes have co-spatial superbubble counterparts, but some do not. This is not
very surprising because they trace two different gas phases. Firstly, from
Section 3.2 we know that H i gas is not to be sensitive to early feedback,
therefore small superbubbles dominated by early feedback are not powerful
enough to generate H i holes. While stellar winds can in principle blow away
the surrounding gas and reduce the density of H i, it is less effective
compared to supernova explosions. Moreover, an H i hole may be destroyed by
disturbances in the disk. On the other hand, superbubbles die when continuous
stellar feedback ceases, and the hot gas subsequently cools down. After the
death of a superbubble, the cavity blown out by stellar wind and supernovae
may still remain there as a H i hole. However, as the cavity loses pressure
support from hot gas, it will easily shrink away or be disturbed by other
surrounding activities. Finally, as mentioned in Section 3.2, some identified
H i holes in our simulations are probably just low density voids between
spiral arms/filaments and are not related to stellar feedback.

###  4.2 Effects of star formation efficiency

Not all of the gas in giant molecular clouds will form stars in a free-fall
time (Krumholz & Tan, 2007). The star formation efficiency $\epsilon_{\rm
ff}$, as defined in Equation (1), determines the fraction of star forming gas
and governs the local star formation rate. Based on the filling factor
profiles in Section 3.5, we see that $\epsilon_{\rm ff}$ has minimal influence
on the overall distribution of hot gas. However, properties of superbubbles
show a strong dependence on $\epsilon_{\rm ff}$, as demonstrated in Section
3.2, 3.3, and 3.4. There must be some fundamental effects of star formation
efficiency.

Superbubbles are basically formed through the interaction between stellar
associations and the surrounding interstellar medium, so the properties of
superbubbles will be regulated by both the environment and the energy injected
by stellar feedback. For the SFE1, SFE10 and SFE100 runs, we estimate the
average time interval between successive supernova events $\Delta t_{\rm SN}$
and the average density of the ambient medium $\rho_{\rm amb}$ (estimated as
the average density within the range from $R_{\rm bubble}$ to $2R_{\rm
bubble}$). We find their $\Delta t_{\rm SN}$ are 0.12, 0.15 and 0.12 Myr, and
their $\rho_{\rm amb}$ are $1.51\times 10^{-1},\;2.35\times 10^{-2}\;{\rm
and}\;6.08\times 10^{-3}\;\rm{}M_{\sun}\;pc^{-3}$. This result indicates that
simulation runs with different $\epsilon_{\rm ff}$ have similar supernova
rates, but larger $\epsilon_{\rm ff}$ will lead to smaller density of the
environment. This probably reveals differences in both the star forming stage
and the early feedback stage. In the star forming stage, larger star formation
efficiency means more efficient consumption of the cold gas, so the
environment density will be lower. In the early feedback stage, radiative
feedback and stellar winds heat the surrounding gas and reduce its density.
Photoionization is very sensitive to the environment density $\rho_{\rm amb}$,
as the Stromgren radius is proportional to $\rho_{\rm amb}^{-2/3}$.
Consequently, the pre-processing of the surrounding gas becomes more
pronounced.

###  4.3 Interaction between early feedback and supernova feedback

While the SFE1 run considers both supernova feedback and early feedback, its
results are not just simply the sum of the SN and Rad runs, and may give us
some inspections about the interaction between early feedback and supernova
feedback.

We note that early feedback, to some extent, suppresses the impact of
supernova feedback. Given the same star formation efficiency in the SN run and
the SFE1 run, the SN run exhibits a greater abundance of huge superbubbles
(with radius larger than $10^{2.4}$ pc) in Fig. 4. This can be attributed to
the absence of early feedback, resulting in more massive stellar associations.
As we see in the bottom panel of Fig. 6, The SN run has more massive stellar
associations than other runs, and early feedback plays a very important role
on decreasing the mass of stellar associations. For a given star-forming gas,
early feedback starts at an early stage right after the formation of the first
star particle. Therefore, it rapidly reduces the formation of star particles
while gas clouds are still collapsing. Without early feedback, clouds keep
collapsing and forming stars until the first SN explosion several million
years later. The effect of early feedback on suppressing star cluster
formation have also been shown in (e.g. Smith et al., 2021; Andersson et al.,
2024). Therefore, it consequently suppresses the power of supernova feedback
and the formation of those huge superbubbles.

###  4.4 Difference with observed X-ray luminosity - size distribution

In section 3.1 we have presented the X-ray luminosity - size distribution of
superbubbles from our fiducial run as well as observations in Fig 3. While
they are largely consistent, we note that the observed distribution seems to
be a little biased towards smaller size. Firstly, the environment and history
of various galaxies could be very different. Moreover, observations are
actually limited by the sensitivity of instruments. For a superbubble at a
given luminosity, larger in size means more diffuse and thus it has lower
surface brightness and harder to be observed. As far as we learn from Fig 3,
there appears to remain some faint superbubbles unexplored by current
observations, so the sample size of observed superbubbles may hopefully
increase in the future with advancements in more sensitive instruments.

###  4.5 Comparison with analytical evolution model

Figure 9: The evolution of superbubble size (solid line) in the SN run, with a
temporal resolution of 1 Myr. Shaded region corresponds to the
$40^{\rm{}th}-60^{\rm{}th}$ percentiles. The dashed line represents the
analytical evolution profiles for the SN run according to Equation (2).

We compare the size evolution of superbubbles in the SN run with the
analytical model given by Weaver et al. (1977) and Mac Low & McCray (1988), in
which self-similar superbubbles are generated by successive supernovae in a
uniform environment. In the pressure-driven snowplow phase, the radius $R$ of
a superbubble can be expressed as

$\displaystyle R$ $\displaystyle=\left(\frac{125}{154\pi}\right)^{1/5}L_{\rm
SN}^{1/5}\;\rho_{\rm amb}^{-1/5}\;t^{3/5}$ (2) $\displaystyle=26.7\;{\rm
pc}\;\left(\frac{E_{\rm SN}}{10^{51}\rm erg}\right)^{1/5}\;\left(\frac{\Delta
t_{\rm SN}}{\rm Myr}\right)^{-1/5}\;\left(\frac{\rho_{\rm amb}}{\rm
M_{\sun}pc^{-3}}\right)^{-1/5}\;\left(\frac{t}{\rm Myr}\right)^{3/5}$

where $t$ is the evolutionary time of a superbubble. The model has three
parameters: the injected energy of each supernova event $E_{\rm SN}$ (it is
set to $10^{51}\rm erg$); the average time interval between successive
supernova events $\Delta t_{\rm SN}$; and the average density of the ambient
medium $\rho_{\rm amb}$. To compare with the analytical model, we estimate
$\Delta t_{\rm SN}$ and $\rho_{\rm amb}$ for each identified superbubbles. In
the evolutionary network in our simulations, we record the number of supernova
events in each time step and calculate the average time interval of
supernovae. The average supernova rate $N_{\rm SN}$ in the SN run is about
$25.3\;\rm{}Myr^{-1}$, and we derive $\Delta t_{\rm SN}$ as the inverse of
$N_{\rm SN}$, which gives $\Delta t_{\rm SN}=0.04\;\rm Myr$. $\rho_{\rm
amb}=1.91\times 10^{-3}\;\rm{}M_{\sun}\;pc^{-3}$ is estimated as the average
density within the range from $R_{\rm bubble}$ to $2R_{\rm bubble}$.

Substituting the parameters into Equation (2), we can obtain predictions for
the SN run using the analytical model, and the result is shown in Fig. 7.
Evolution of superbubbles in the SN run is actually different from the
analytical model. In the analytical model, it is assumed that supernovae
explode one by one, with a uniform time interval $\Delta t_{\rm SN}$. However,
we found that the number of supernovae is usually larger in the beginning and
gradually decreases in our simulations. Moreover, while the analytical model
assumes an unchanged uniform environment, in our simulations the expansion
superbubbles will experience pressure confinement from the shocked surrounding
materials, and this suppression is more significant as superbubbles grow
larger. We also note that the zero-age point we defined in our simulation is
also different from that of the analytical model. In the analytical model, it
is assumed that a superbubble is born immediately when stellar feedback
begins. However, we define the birth of a superbubble at the time when we
first found it in our simulations.

###  4.6 Caveats

Some caveats should be noted. First, the simulations were made for isolated
$L_{*}$ galaxies. In comparison with observations, however, we included all
superbubbles including those in dwarf galaxies like LMC and M33, due to the
limited amount of X-ray observations. Second, the zoom-in simulations were
made without considering the effects of inflows from the large scale
structures into galaxies or AGN feedback. Although it is likely to lead to an
inaccurate account of baryon cycles in galaxy ecosystems, it is not expected
to significantly affect the properties of the ISM, which is the focus of this
work. Contribution from ultraluminous X-ray sources (ULXs, Oskinova et al.,
2019; Pinto & Walton, 2023) is not considered as well. Third, we have found
that early feedback processes (radiative feedback and stellar winds) can
impact the size and evolution of superbubbles. However, full radiative
transfer is not implemented in the simulations; the effects of radiative
feedback (including photoheating and radiation pressure) are calculated from
an effective subgrid model (McKee & Krumholz, 2010). Fourth, the
classification of two types of superbubbles in Section 3.3 is not perfect,
especially for high $\epsilon_{\rm ff}$ runs. Fifth, in this work,
superbubbles are identified as X-ray emitting bubbles, while the actual
structures of superbubbles may be more complex. The cooler outskirts of
superbubbles would emit at longer wavelengths. Finally, we note that our
automated H i hole identifying method is not perfect and might include some
local minima which are not related with feedback. However, the purpose of our
study is to find the effects of different SFE and feedback processes on
superbubbles and HI holes. The choice of parameters for HI hole identification
is not expected to affect our main result that SFE does not have much impact
on the distribution of HI hole size.

##  5 Summary

We use hydrodynamical simulations of isolated Milky Way-sized galaxies based
on the SMUGGLE model in Arepo to study the properties and evolution of
superbubbles. We are able to resolve and identify 2D H i holes and 3D
superbubbles, and construct evolution tracks of superbubbles using stellar
associations as tracers. We conduct statistical analysis of superbubbles and
take direct comparisons with both distribution in observations and evolution
in analytical models. Here, we summarize our main results.

  * •

The results from the fiducial run show a positive correlation between the
X-ray luminosity and size of superbubbles. The correlation is in good
agreement with observations of nearby galaxies.

  * •

The size distribution of H i holes shows only one peak in all runs, and is not
sensitive to radiative and stellar wind feedback.

  * •

The superbubbles show a double-peaked size distribution. The size of
superbubbles dominated by supernova feedback is generally larger than that of
early feedback (radiative and stellar wind feedback) dominated superbubbles.
The size distribution of superbubbles is different from that of H i holes. In
addition, not all superbubbles coincide spatially with H i holes.

  * •

Star formation efficiency does not have significant influence on the overall
distribution of hot gas but has significant impact on the properties of
superbubbles. Larger star formation efficiency will lead to a higher
proportion of large superbubbles with radius more than 100 pc. Star formation
efficiency mainly influences the processes in early feedback, by regulating
the environmental density.

  * •

Early feedback plays an important role in suppressing supernova feedback and
the mass of stellar associations. Through the process of heating and ionizing
the surrounding environment of newly born stars, early stellar feedback
actively restrains the formation of additional stars. Consequently, fewer star
particles are formed compared to the scenario without early feedback,
resulting in a decreased number of supernovae and thus a reduced supernova
feedback.

  * •

The average volume filling factor of the hot gas is about 12 %, but its
distribution is not uniform across the galaxy. The filling factor reaches more
than 80 % at the galactic center and decreases towards the outer disk, where
the warm gas mainly lies. The filling factor of the cold gas is a very low
throughout the galaxy.

## Acknowledgements

We thank the anonymous referee for the suggestions that have significantly
improved the manuscript, Volker Springel for allowing us access to the Arepo
code, and Yu-Ning Zhang, Rui Huang and Jiejia Liu for very useful discussion
and comments. This work was supported in part by the National Natural Science
Foundation of China through Grant No. 11821303 & 12373006, by the Ministry of
Science and Technology of China through Grant 2018YFA0404502 and
2023YFB3002502, and by China Manned Space Program through its Space
Application System. We also acknowledge the Tsinghua Astrophysics High-
Performance Computing (TAHPC) platform for providing computational and data
storage resources that have contributed to this work. The identification of 2D
H i holes in Section 2.4 made use of astrodendro, a Python package to compute
dendrograms of Astronomical data 222http://www.dendrograms.org/. LVS
acknowledges financial support from grants NSF-CAREER-1945310 and NSF
AST2107993.

## Data Availability

The data that support the findings of this study are available from the
corresponding author, upon reasonable request.

## References

  * Ali et al. (2022) Ali A. A., Bending T. J. R., Dobbs C. L., 2022, MNRAS, 510, 5592
  * Andersson et al. (2024) Andersson E. P., Mac Low M.-M., Agertz O., Renaud F., Li H., 2024, A&A, 681, A28
  * Bagetakos et al. (2011) Bagetakos I., Brinks E., Walter F., de Blok W. J. G., Usero A., Leroy A. K., Rich J. W., Kennicutt R. C. J., 2011, AJ, 141, 23
  * Bialy et al. (2021) Bialy S., et al., 2021, ApJ, 919, L5
  * Bieri et al. (2023) Bieri R., Naab T., Geen S., Coles J. P., Pakmor R., Walch S., 2023, MNRAS, 523, 6336
  * Bomans et al. (1994) Bomans D. J., Dennerl K., Kurster M., 1994, A&A, 283, L21
  * Boomsma et al. (2008) Boomsma R., Oosterloo T. A., Fraternali F., van der Hulst J. M., Sancisi R., 2008, A&A, 490, 555
  * Bregman (1980) Bregman J. N., 1980, ApJ, 236, 577
  * Brinks & Bajaja (1986) Brinks E., Bajaja E., 1986, A&A, 169, 14
  * Burger et al. (2022a) Burger J. D., Zavala J., Sales L. V., Vogelsberger M., Marinacci F., Torrey P., 2022a, Phys. Rev. Lett., 129, 191103
  * Burger et al. (2022b) Burger J. D., Zavala J., Sales L. V., Vogelsberger M., Marinacci F., Torrey P., 2022b, MNRAS, 513, 3458
  * Burrows & Zhiyu (1996) Burrows D. N., Zhiyu G., 1996, in Zimmermann H. U., Trümper J., Yorke H., eds, Roentgenstrahlung from the Universe. pp 221–224
  * Cash et al. (1980) Cash W., Charles P., Bowyer S., Walter F., Garmire G., Riegler G., 1980, ApJ, 238, L71
  * Chabrier (2001) Chabrier G., 2001, ApJ, 554, 1274
  * Chu & Mac Low (1990) Chu Y.-H., Mac Low M.-M., 1990, ApJ, 365, 510
  * Chu et al. (1995) Chu Y.-H., Chang H.-W., Su Y.-L., Mac Low M.-M., 1995, ApJ, 450, 157
  * Chu et al. (2003) Chu Y.-H., Guerrero M. A., Gruendl R. A., García-Segura G., Wendker H. J., 2003, ApJ, 599, 1189
  * Cosens et al. (2022) Cosens M., et al., 2022, ApJ, 929, 74
  * Cox & Reynolds (1987) Cox D. P., Reynolds R. J., 1987, ARA&A, 25, 303
  * Cox & Smith (1974) Cox D. P., Smith B. W., 1974, ApJ, 189, L105
  * Cui et al. (1996) Cui W., Sanders W. T., McCammon D., Snowden S. L., Womble D. S., 1996, ApJ, 468, 102
  * Davies et al. (1976) Davies R. D., Elliott K. H., Meaburn J., 1976, Mem. RAS, 81, 89
  * Dunne et al. (2001) Dunne B. C., Points S. D., Chu Y.-H., 2001, ApJS, 136, 119
  * Dunne et al. (2003) Dunne B. C., Chu Y.-H., Chen C. H. R., Lowry J. D., Townsley L., Gruendl R. A., Guerrero M. A., Rosado M., 2003, ApJ, 590, 306
  * Egorov et al. (2023) Egorov O. V., et al., 2023, ApJ, 944, L16
  * Ehlerová & Palouš (2005) Ehlerová S., Palouš J., 2005, A&A, 437, 101
  * El-Badry et al. (2019) El-Badry K., Ostriker E. C., Kim C.-G., Quataert E., Weisz D. R., 2019, MNRAS, 490, 1961
  * Fielding et al. (2018) Fielding D., Quataert E., Martizzi D., 2018, MNRAS, 481, 3325
  * Foster et al. (2012) Foster A. R., Ji L., Smith R. K., Brickhouse N. S., 2012, ApJ, 756, 128
  * Gent et al. (2013) Gent F. A., Shukurov A., Fletcher A., Sarson G. R., Mantere M. J., 2013, MNRAS, 432, 1396
  * Gupta et al. (2016) Gupta S., Nath B. B., Sharma P., Shchekinov Y., 2016, MNRAS, 462, 4532
  * Hartwell et al. (2004) Hartwell J. M., Stevens I. R., Strickland D. K., Heckman T. M., Summers L. K., 2004, MNRAS, 348, 406
  * Heckman et al. (1995) Heckman T. M., Dahlem M., Lehnert M. D., Fabbiano G., Gilmore D., Waller W. H., 1995, ApJ, 448, 98
  * Henize (1956) Henize K. G., 1956, ApJS, 2, 315
  * Hopkins et al. (2012) Hopkins P. F., Quataert E., Murray N., 2012, MNRAS, 421, 3488
  * Jahn et al. (2023) Jahn E. D., et al., 2023, MNRAS, 520, 461
  * Kavanagh et al. (2012) Kavanagh P. J., Sasaki M., Points S. D., 2012, A&A, 547, A19
  * Kavanagh et al. (2020) Kavanagh P. J., et al., 2020, A&A, 637, A12
  * Keller et al. (2014) Keller B. W., Wadsley J., Benincasa S. M., Couchman H. M. P., 2014, MNRAS, 442, 3013
  * Kim & Ostriker (2015) Kim C.-G., Ostriker E. C., 2015, ApJ, 802, 99
  * Kim et al. (2017) Kim C.-G., Ostriker E. C., Raileanu R., 2017, ApJ, 834, 25
  * Kollmeier et al. (2019) Kollmeier J., et al., 2019, in Bulletin of the American Astronomical Society. p. 274
  * Krause & Diehl (2014) Krause M. G. H., Diehl R., 2014, ApJ, 794, L21
  * Krause et al. (2018) Krause M. G. H., et al., 2018, A&A, 619, A120
  * Krumholz & Tan (2007) Krumholz M. R., Tan J. C., 2007, ApJ, 654, 304
  * Lallement et al. (2003) Lallement R., Welsh B. Y., Vergely J. L., Crifo F., Sfeir D., 2003, A&A, 411, 447
  * Lehner & Howk (2010) Lehner N., Howk J. C., 2010, ApJ, 709, L138
  * Li et al. (2017) Li M., Bryan G. L., Ostriker J. P., 2017, ApJ, 841, 101
  * Li et al. (2020) Li H., Vogelsberger M., Marinacci F., Sales L. V., Torrey P., 2020, MNRAS, 499, 5862
  * Li et al. (2022) Li H., Vogelsberger M., Bryan G. L., Marinacci F., Sales L. V., Torrey P., 2022, MNRAS, 514, 265
  * Long et al. (2010) Long K. S., et al., 2010, ApJS, 187, 495
  * Lucke & Hodge (1970) Lucke P. B., Hodge P. W., 1970, AJ, 75, 171
  * Mac Low & McCray (1988) Mac Low M.-M., McCray R., 1988, ApJ, 324, 776
  * Marinacci et al. (2019) Marinacci F., Sales L. V., Vogelsberger M., Torrey P., Springel V., 2019, MNRAS, 489, 4233
  * Mashchenko et al. (1999) Mashchenko S. Y., Thilker D. A., Braun R., 1999, A&A, 343, 352
  * McCammon et al. (1971) McCammon D., Bunner A. N., Coleman P. L., Kraushaar W. L., 1971, ApJ, 168, L33
  * McKee (1990) McKee C. F., 1990, in Blitz L., ed., Astronomical Society of the Pacific Conference Series Vol. 12, The Evolution of the Interstellar Medium. pp 3–29
  * McKee & Krumholz (2010) McKee C. F., Krumholz M. R., 2010, ApJ, 709, 308
  * McKee & Ostriker (1977) McKee C. F., Ostriker J. P., 1977, ApJ, 218, 148
  * Meaburn (1980) Meaburn J., 1980, MNRAS, 192, 365
  * Naab & Ostriker (2017) Naab T., Ostriker J. P., 2017, ARA&A, 55, 59
  * Nath et al. (2020) Nath B. B., Das P., Oey M. S., 2020, MNRAS, 493, 1034
  * Oey (1996) Oey M. S., 1996, ApJ, 467, 666
  * Oey & Clarke (1997) Oey M. S., Clarke C. J., 1997, MNRAS, 289, 570
  * Oskinova et al. (2019) Oskinova L. M., Bik A., Mas-Hesse J. M., Hayes M., Adamo A., Östlin G., Fürst F., Otí-Floranes H., 2019, A&A, 627, A63
  * Pidopryhora et al. (2007) Pidopryhora Y., Lockman F. J., Shields J. C., 2007, ApJ, 656, 928
  * Pietsch et al. (1998) Pietsch W., Trinchieri G., Vogler A., 1998, A&A, 340, 351
  * Pinto & Walton (2023) Pinto C., Walton D. J., 2023, arXiv e-prints, p. arXiv:2302.00006
  * Pokhrel et al. (2020) Pokhrel N. R., Simpson C. E., Bagetakos I., 2020, AJ, 160, 66
  * Predehl et al. (2020) Predehl P., et al., 2020, Nature, 588, 227
  * Rosolowsky et al. (2008) Rosolowsky E. W., Pineda J. E., Kauffmann J., Goodman A. A., 2008, ApJ, 679, 1338
  * Sales et al. (2014) Sales L. V., Marinacci F., Springel V., Petkova M., 2014, MNRAS, 439, 2990
  * Sartorio et al. (2021) Sartorio N. S., Vandenbroucke B., Falceta-Goncalves D., Wood K., 2021, MNRAS, 500, 1833
  * Sasaki et al. (2012) Sasaki M., Pietsch W., Haberl F., Hatzidimitriou D., Stiele H., Williams B., Kong A., Kolb U., 2012, A&A, 544, A144
  * Schuster et al. (2007) Schuster K. F., Kramer C., Hitschfeld M., Garcia-Burillo S., Mookerjea B., 2007, A&A, 461, 143
  * Sivasankaran et al. (2022) Sivasankaran A., et al., 2022, MNRAS, 517, 4752
  * Smith et al. (2001) Smith R. K., Brickhouse N. S., Liedahl D. A., Raymond J. C., 2001, ApJ, 556, L91
  * Smith et al. (2021) Smith M. C., Bryan G. L., Somerville R. S., Hu C.-Y., Teyssier R., Burkhart B., Hernquist L., 2021, MNRAS, 506, 3882
  * Smith et al. (2022) Smith A., et al., 2022, MNRAS, 517, 1
  * Snowden et al. (1990) Snowden S. L., Cox D. P., McCammon D., Sanders W. T., 1990, ApJ, 354, 211
  * Snowden et al. (1998) Snowden S. L., Egger R., Finkbeiner D. P., Freyberg M. J., Plucinsky P. P., 1998, ApJ, 493, 715
  * Springel (2010) Springel V., 2010, MNRAS, 401, 791
  * Stinson et al. (2013) Stinson G. S., Brook C., Macciò A. V., Wadsley J., Quinn T. R., Couchman H. M. P., 2013, MNRAS, 428, 129
  * Su et al. (2010) Su M., Slatyer T. R., Finkbeiner D. P., 2010, ApJ, 724, 1044
  * Summers et al. (2003) Summers L. K., Stevens I. R., Strickland D. K., Heckman T. M., 2003, MNRAS, 342, 690
  * Summers et al. (2004) Summers L. K., Stevens I. R., Strickland D. K., Heckman T. M., 2004, MNRAS, 351, 1
  * Sun et al. (2012) Sun W., Chen Y., Feng L., Chu Y.-H., Chen C. H. R., Wang Q. D., Li J.-T., 2012, ApJ, 760, 61
  * Thilker et al. (1998) Thilker D. A., Braun R., Walterbos R. M., 1998, A&A, 332, 429
  * Townsley et al. (2006) Townsley L. K., Broos P. S., Feigelson E. D., Brandl B. R., Chu Y.-H., Garmire G. P., Pavlov G. G., 2006, AJ, 131, 2140
  * Tüllmann et al. (2008) Tüllmann R., et al., 2008, ApJ, 685, 919
  * Tüllmann et al. (2009) Tüllmann R., et al., 2009, ApJ, 707, 1361
  * Tumlinson et al. (2017) Tumlinson J., Peeples M. S., Werk J. K., 2017, ARA&A, 55, 389
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Vogelsberger et al. (2020) Vogelsberger M., Marinacci F., Torrey P., Puchwein E., 2020, Nature Reviews Physics, 2, 42
  * Wang & Helfand (1991a) Wang Q., Helfand D. J., 1991a, ApJ, 373, 497
  * Wang & Helfand (1991b) Wang Q., Helfand D. J., 1991b, ApJ, 379, 327
  * Wang et al. (2005) Wang Q. D., Whitaker K. E., Williams R., 2005, MNRAS, 362, 1065
  * Watkins et al. (2023) Watkins E. J., et al., 2023, A&A, 676, A67
  * Weaver et al. (1977) Weaver R., McCray R., Castor J., Shapiro P., Moore R., 1977, ApJ, 218, 377
  * Welsh et al. (2010) Welsh B. Y., Lallement R., Vergely J. L., Raimond S., 2010, A&A, 510, A54
  * Yadav et al. (2017) Yadav N., Mukherjee D., Sharma P., Nath B. B., 2017, MNRAS, 465, 1720
  * Zemp et al. (2011) Zemp M., Gnedin O. Y., Gnedin N. Y., Kravtsov A. V., 2011, ApJS, 197, 30
  * Zucker et al. (2022) Zucker C., et al., 2022, Nature, 601, 334
  * de Avillez & Breitschwerdt (2004) de Avillez M. A., Breitschwerdt D., 2004, A&A, 425, 899
  * de Heij et al. (2002) de Heij V., Braun R., Burton W. B., 2002, A&A, 391, 159

##  Appendix A Catalogue of superbubble candidates and central diffuse
emissions

Table 3: Superbubble candidates in M33, M31 and M101. Corrected luminosities for 0.1 - 2.4 keV band are shown in brackets. NO. | Name | $r$ | log $L_{\rm X}$ | Reference
---|---|---|---|---
|  | (pc) | (erg / s) |
M33
1 | L10-012 | 26.0 | 35.08 | Long et al. (2010)
2 | L10-026 | 36.5 | 34.69 | Long et al. (2010)
3 | L10-043 | 42.5 | 35.08 | Long et al. (2010)
4 | L10-089 | 46.0 | 35.04 | Long et al. (2010)
5 | L10-092 | 50.5 | 34.62 | Long et al. (2010)
6 | L10-117 | 33.0 | 35.28 | Long et al. (2010)
7 | L10-119 | 16.5 | 35.48 | Long et al. (2010)
8 | L10-122 | 55.5 | 35.11 | Long et al. (2010)
M31
9 | K20 region 2 | 324.0 | 36.43 (36.50) | Kavanagh et al. (2020)
10 | K20 region 3 | 300.0 | 35.76 (35.85) | Kavanagh et al. (2020)
11 | SPH11-263 | 47.0 | 35.60 | Sasaki et al. (2012)
12 | SPH11-1234 | 25.0 | 36.90 (36.91) | Sasaki et al. (2012)
M101
13 | NGC 5461 | 516.2 | 38.46 (38.69) | Sun et al. (2012)
14 | NGC 5462 | 663.3 | 38.23 (38.50) | Sun et al. (2012)
Table 4: Diffuse emission from central region of galaxies. Corrected luminosities for 0.1 - 2.4 keV band are shown in brackets. NO. | Name | $r$ | log $L_{\rm X}$ | Reference
---|---|---|---|---
|  | (pc) | (erg / s) |
1 | eROSITA | 7000.0 | 39.01 (39.45) | Predehl et al. (2020)
Bubbles
2 | NGC 3079 | 1030.8 | 39.99 (39.99) | Pietsch et al. (1998)
3 | IC 10 | 408.0 | 37.90 (38.10) | Wang et al. (2005)
4 | NGC 1569 | 948.7 | 38.70 (38.70) | Heckman et al. (1995)
5 | NGC 4449 | 979.8 | 38.96 (39.09) | Summers et al. (2003)
6 | NGC 5253 | 530.0 | 38.60 (38.83) | Summers et al. (2004)
7 | NGC 4214 | 524.4 | 38.48 (38.83) | Hartwell et al. (2004)

We list in Table 3 the information of superbubble candidates in some nearby
galaxies (M33, M31 and M101). They are not identified with confidence, but
their size and luminosity are still valuable due to the very small sample size
of confirmed superbubbles. They are also valuable targets for future
observations and maybe some of these candidates can be confirmed as
superbubbles in the future. In principle, our superbubble identification
method is able to find some kpc-scale structures which are just breaking out
from the galactic center. This kind of source is usually very extensive and
bright. Table 4 gives some examples of central diffuse emissions (CDEs) from
the Milky Way (i.e. the eROSITA Bubbles) and other galaxies.

##  Appendix B ISM-CGM phase structure

Figure 10: Phase diagram of gas for the SFE1 run at 0.4 Gyr in the form of a
two-dimensional histogram in density and temperature space, with the color of
each bin represents the mass fraction of gas at the corresponding density-
temperature phase. Different phases are presented as different regions
separated by dashed black lines.

As mentioned in Marinacci et al. (2019), the SMUGGLE model is capable to self-
consistently generate a multi-phase structure of the ISM. Fig. 10 shows the
phase diagram of temperature and density in SFE1 run at 0.4 Gyr. The diagram
is divided into multiple regions, largely following the approach of Marinacci
et al. (2019), with some necessary modifications made according to our study
of superbubbles. We consider the interstellar medium (ISM) as the gas within
the galactic disk, corresponding to phase (a) to (d). Most of the gas retains
a warm state as phase (a), with a temperature of about $10^{4}~{}\mathrm{K}$.
Some of the gas will cool down into cold phase, collapse into dense clumps and
form giant molecular clouds (GMCs) in phase (b), with some molecular cores
evolve dense enough (i.e., $n>100\,\rm{}cm^{-3}$) and become star forming gas
in phase (c). Phase (d1) and (d2) both represent the gas emerging from stellar
feedback, with (d1) resulting mainly from radiative and stellar wind feedback
while (d2) is dominated by supernova feedback. After stars are formed,
radiation from stars may ionize the surrounding gas as early feedback, heat
the gas to a temperature of about $1.7\times 10^{4}~{}\mathrm{K}$ and push the
gas away through radiation pressure and stellar winds. This gas primarily
occupies phase (d1). The gas in phase (d2) will be even hotter and more
diffuse due to supernova events. In this work, we mainly focus on the
superbubbles, which largely correspond to the region (d1) and (d2) with
temperatures higher than $10^{5.5}\;\rm K$ (indicated by the dashed blue
line). We note that the warm gas (with temperature lower than $10^{5.5}\;\rm
K$) in phase (d1) and (d2) largely corresponds to the observed H ii regions.

Gas in phase (d2) may eventually flow out of the galaxy and into the
circumgalactic medium (CGM), which is largely referred to as phase (e). The
outflowing gas will undergo an expansion process, gradually decreasing the
temperature and becoming more diffuse. Some of the gas in phase (e) is
possible to fall back into the ISM as phase (a), and this kind of falling gas
is probably the origin of the observed high velocity clouds (HVCs, e.g.
Bregman, 1980; Lehner & Howk, 2010). Intense stellar feedback in the galactic
center can occasionally generate a massive outflow and inject the gas to
distances comparable to, or even larger than, the virial radius. Some of the
outflow gas may become too distant from the galaxy and cannot fall back,
instead it will gradually mix with the gas in phase (f), which represents the
background gas environment. It is important to note that the simulations we
employ do not account for the cosmological environment and there is no gas
inflow from the IGM, so the CGM consists of only the outflow gas from the ISM
and it might be considered as a lower limit to the presence of hot and diffuse
gas.