# Where do stars explode in the ISM? – The distribution of dense gas around
massive stars and supernova remnants in M33

Sumit K. Sarbadhicary Department of Physics, The Ohio State University,
Columbus, Ohio 43210, USA  Center for Cosmology & Astro-Particle Physics, The
Ohio State University, Columbus, Ohio 43210, USA  Department of Astronomy, The
Ohio State University, Columbus, Ohio 43210, USA  Jordan Wagner  Department of
Physics, The Ohio State University, Columbus, Ohio 43210, USA  Center for
Cosmology & Astro-Particle Physics, The Ohio State University, Columbus, Ohio
43210, USA  Department of Astronomy, The Ohio State University, Columbus, Ohio
43210, USA  Eric W. Koch Center for Astrophysics $\mid$ Harvard & Smithsonian,
60 Garden St., 02138 Cambridge, MA, USA  Ness Mayker Chen Center for Cosmology
& Astro-Particle Physics, The Ohio State University, Columbus, Ohio 43210, USA
Department of Astronomy, The Ohio State University, Columbus, Ohio 43210, USA
Adam K. Leroy Center for Cosmology & Astro-Particle Physics, The Ohio State
University, Columbus, Ohio 43210, USA  Department of Astronomy, The Ohio State
University, Columbus, Ohio 43210, USA  Natalia Lahén Max-Planck-Institute für
Astrophysik, Karl-Schwarzschild-Straße 1, D-85740 Garching, Germany  Erik
Rosolowsky Department of Physics, University of Alberta, Edmonton, AB T6G 2E1,
Canada  Kathryn F. Neugent NASA Hubble Fellow  Center for Astrophysics $\mid$
Harvard & Smithsonian, 60 Garden St., 02138 Cambridge, MA, USA  Chang-Goo Kim
Department of Astrophysical Sciences, Princeton University, Princeton, NJ
08544, USA  Laura Chomiuk Center for Data Intensive and Time Domain Astronomy,
Department of Physics and Astronomy, Michigan State University, East Lansing,
MI 48824, USA  Julianne J. Dalcanton Department of Astronomy, University of
Washington, Box 351580, Seattle, WA 98195-1580  Center for Computational
Astrophysics, Flatiron Institute, 162 Fifth Ave, New York, NY 10010, USA
Laura A. Lopez Center for Cosmology & Astro-Particle Physics, The Ohio State
University, Columbus, Ohio 43210, USA  Department of Astronomy, The Ohio State
University, Columbus, Ohio 43210, USA  Nickolas M. Pingel Department of
Astronomy, University of Wisconsin-Madison, 475 North Charter St., Madison,
WI, 53706-15821, USA  Remy Indebetouw University of Virginia Astronomy, Dept.
530 McCormick Road, Charlottesville, VA 22904-4325  National Radio Astronomy
Observatory, 520 Edgemont Road, Charlottesville, VA 22903-2475  Thomas G.
Williams Sub-department of Astrophysics, Department of Physics, University of
Oxford, Keble Road, Oxford OX1 3RH, UK  Elizabeth Tarantino Space Telescope
Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA  Department
of Astronomy, University of Maryland, College Park, MD 20742, USA  Jennifer
Donovan Meyer National Radio Astronomy Observatory, 520 Edgemont Road,
Charlottesville, VA 22903-2475  Evan D. Skillman Minnesota Institute for
Astrophysics, University of Minnesota, Minneapolis, MN 55455, USA  Adam
Smercina Department of Astronomy, University of Washington, Box 351580,
Seattle, WA 98195-1580  Amanda A. Kepley National Radio Astronomy Observatory,
520 Edgemont Road, Charlottesville, VA 22903-2475  Eric J. Murphy National
Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903-2475
Jay Strader Center for Data Intensive and Time Domain Astronomy, Department of
Physics and Astronomy, Michigan State University, East Lansing, MI 48824, USA
Tony Wong Department of Astronomy, University of Illinois, Urbana, IL 61801,
USA  Snežana Stanimirović Department of Astronomy, University of Wisconsin-
Madison, 475 North Charter St., Madison, WI, 53706-15821, USA  Vicente
Villanueva Department of Astronomy, University of Maryland, College Park, MD
20742, USA  Fabian Walter Max Planck Institute for Astronomy, Königstuhl 17,
D-69117 Heidelberg, Germany  National Radio Astronomy Observatory, Pete V.
Domenici Array Science Center, P.O. Box O, Socorro, NM 87801, USA  Juergen Ott
New Mexico Institute of Mining and Technology, 801 Leroy Place, Socorro, NM
87801, USA  National Radio Astronomy Observatory, Pete V. Domenici Array
Science Center, P.O. Box O, Socorro, NM 87801, USA  Jeremy Darling Center for
Astrophysics and Space Astronomy, Department of Astrophysical and Planetary
Sciences, University of Colorado, 389 UCB, Boulder, CO 80309-0389, USA  Julia
Roman-Duval Space Telescope Science Institute, 3700 San Martin Drive,
Baltimore, MD 21218, USA  Claire E. Murray Space Telescope Science Institute,
3700 San Martin Drive, Baltimore, MD 21218, USA  Department of Physics &
Astronomy, Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD
21218 Sumit K. Sarbadhicary<EMAIL_ADDRESS>
###### Abstract

Star formation in galaxies is regulated by turbulence, outflows, gas heating
and cloud dispersal – processes which depend sensitively on the properties of
the interstellar medium (ISM) into which supernovae (SNe) explode.
Unfortunately, direct measurements of ISM environments around SNe remain
scarce, as SNe are rare and often distant. Here we demonstrate a new approach:
mapping the ISM around the massive stars that are soon to explode. This
provides a much larger census of explosion sites than possible with only SNe,
and allows comparison with sensitive, high-resolution maps of the atomic and
molecular gas from the Jansky VLA and ALMA. In the well-resolved Local Group
spiral M33, we specifically observe the environments of red supergiants (RSGs,
progenitors of Type II SNe), Wolf-Rayet stars (WRs, tracing stars $>$30M⊙, and
possibly future stripped-envelope SNe), and supernova remnants (SNRs,
locations where SNe have exploded). We find that massive stars evolve not only
in dense, molecular-dominated gas (with younger stars in denser gas), but also
a substantial fraction ($\sim$45% of WRs; higher for RSGs) evolve in lower-
density, atomic-gas-dominated, inter-cloud media. We show that these
measurements are consistent with expectations from different stellar-age
tracer maps, and can be useful for validating SN feedback models in numerical
simulations of galaxies. Along with the discovery of a 20-pc diameter
molecular gas cavity around a WR, these findings re-emphasize the importance
of pre-SN/correlated-SN feedback evacuating the dense gas around massive stars
before explosion, and the need for high-resolution (down to pc-scale) surveys
of the multi-phase ISM in nearby galaxies.

Interstellar medium (847) — Stellar feedback (1602) — Massive stars (732) —
Supernovae (1668) — Radio astronomy (1338) — Millimeter astronomy (1061)

††facilities: VLA, ALMA, Mayall, UKIRT††software: CASA (McMullin et al.,
2007), astropy (Astropy Collaboration et al., 2022), matplotlib (Hunter,
2007), scipy (Virtanen et al., 2020), numpy (Harris et al., 2020), photutils
(Bradley et al., 2023).

##  1 Introduction

Feedback from supernovae (SNe) plays a central role in galaxy evolution. SN
explosions drive outflows (e.g. Strickland & Heckman, 2009; Heckman &
Thompson, 2017; Veilleux et al., 2020), turbulence and the multi-phase
distribution of the interstellar medium (ISM, e.g. McKee & Ostriker, 1977;
Joung & Mac Low, 2006; Ostriker & Shetty, 2011a; Hill et al., 2012; Kim et
al., 2013; Martizzi et al., 2016; Kannan et al., 2020; Rathjen et al., 2021;
Kim et al., 2022), accelerate cosmic rays (e.g. Ackermann et al., 2013;
Caprioli & Spitkovsky, 2014; Bykov et al., 2018), and disperse metals
synthesized in stars (see e.g. Walch et al., 2015; Andrews et al., 2017;
Telford et al., 2019; Li & Tonnesen, 2020, for recent simulations and
measurements). The balance of these mechanisms determines the star-forming
properties of the ISM (as shown e.g. in the recent simulations of Ostriker &
Kim, 2022, one of many examples).

While the importance of SNe in galaxy evolution is broadly recognized, many
aspects of the physics remain poorly understood. For example, the properties
of the ISM, which spans several orders of magnitude in density and
temperature, affect the spatial scales and timescales on which SN blastwaves
cool and subsequently share the thermal energy and momentum with the ambient
gas (Li et al., 2015; Kim & Ostriker, 2015a; Martizzi et al., 2015; Karpov et
al., 2020; Koo et al., 2020; Oku et al., 2022). As such, the accumulation of
thermal energy and momentum from SNe can significantly vary with the
environments in which they are seeded. Simulations have confirmed the
importance of the SN ambient density, finding large differences in the global
properties of the ISM such as star-formation, outflow rates and phase
distribution, depending on the correlation between SN locations and ISM
density peaks (Gatto et al., 2015; Walch et al., 2015; Girichidis et al.,
2016; Martizzi et al., 2016, e.g).

Unfortunately, the locations where stars explode depend on a number of poorly
understood physical properties of massive stars, such as their evolutionary
timescales (e.g., Zapartas et al., 2017; Eldridge & Stanway, 2022), potential
for explosion from different progenitor masses $\gtrsim$18 M⊙(Sukhbold et al.,
2016), single/binary-evolution-driven mass-loss (Langer, 2012; Smith, 2014),
3D kinematics from interactions in many-body environments (Oh & Kroupa, 2016;
Oey et al., 2018) or kicks from SN explosions in binary systems (Eldridge et
al., 2011; Renzo et al., 2019; Dallas et al., 2022), and the co-evolution
timescales with their parent molecular clouds (e.g., Kruijssen et al., 2019;
Chevance et al., 2020). While modern high-resolution simulations are
increasingly making great strides in capturing many of these intricate
processes and their effects on the ISM (e.g. Kim & Ostriker, 2017; Hu et al.,
2017; Gutcke et al., 2021), a comprehensive set of observations of the diverse
ambient environments of SNe is needed to serve as a useful guide for
simulations, but presently lacking.

In recent years, detailed mapping of the multi-phase ISM down to scales
($\sim$100 pc) of individual clouds and their co-eval stellar associations
have become possible with modern synthesis arrays like the Atacama Large
Millimeter Array (ALMA) and the Karl G. Jansky Very Large Array (VLA) –
effectively opening the window to _direct_ observations of the cloud-scale ISM
around SNe. Using the PHANGS survey of nearby galaxies within 20 Mpc (Leroy et
al., 2021a), Mayker Chen et al. (2022) (hereafter MC22) studied the molecular
gas environments of 63 historical SNe in 31 galaxies at physical resolutions
spanning 60-150 pc, and found that about 15$\%$ of stripped envelope (SE) SNe
(Type IIb, Ib, Ic) and 50$\%$ of Type II SNe explode at least 150 pc away from
the sites of giant molecular clouds (GMCs). This was a _direct_ confirmation
that a significant fraction (if not all) core-collapse SNe explode in low-
density environments outside of GMCs.

The potential to expand this analysis, while no longer limited by the lack of
cloud-scale resolution that was inaccessible to previous generation
extragalactic CO-surveys (e.g. Galbany et al., 2017), faces a more fundamental
limitation. SNe in individual galaxies are rare and require observing large
cosmological volumes to obtain a sizable sample (Li et al., 2011), precluding
the high spatial resolution needed to characterize their environments. While
the mm-wavelengths are short enough that ALMA can still resolve molecular
clouds out to 20 Mpc with longer baselines, the ability to observe the atomic
(H I) ISM – which is a significant component by both mass and volume in star-
forming galaxies (Kulkarni & Heiles, 1987; Dickey & Lockman, 1990; Walter et
al., 2008) – at arcsecond resolution out to the same distances is unfeasible,
at least until the advent of the next-generation VLA (Murphy et al., 2018) and
Square Kilometer Array (SKA). Thus while studies of the molecular gas around
SNe are promising in the nearby universe, the _total_ ISM density at the scale
of individual clouds around SNe restricts our observing volume to the nearest
galaxies, where again, SN counts become vanishingly small.

Here we propose a novel strategy for uncovering the ISM distribution near SNe
– by focusing instead on evolved massive stars that will likely explode within
a few Myr in Local Group galaxies. The approach has two major advantages.
Firstly, even though SNe are rare, we have access to high-completeness
catalogs of thousands of evolved massive stars like red supergiants (RSGs) and
Wolf-Rayet stars (WRs, Neugent & Massey, 2011; Massey et al., 2021b) that have
been meticulously constructed over decades and are strong candidates for
future core-collapse SNe (see Section 2). We also have access to supernova
remnants (SNRs, White et al., 2019) which trace locations of actual explosions
that have occurred in the past. This enables a much larger catalog of
‘explosion sites’ than achievable with a SN catalog while retaining
substantial spatial detail of the environments. Secondly, the proximity of
Local Group galaxies111While not the focus of this paper, we note that another
advantage in the Local Group is availability of detailed _star-formation
histories_ from observations of resolved stellar populations (Massey et al.,
2006; Cioni et al., 2011; Dalcanton et al., 2012; Nidever et al., 2017;
Lazzarini et al., 2022). We use them briefly in our analyses to obtain useful
insight into the distribution of our massive stars in certain sections. offers
not only cloud-scale observations of molecular gas (Engargiola et al., 2003;
Rosolowsky et al., 2007; Miura et al., 2012; Wong et al., 2011; Druard et al.,
2014; Caldú-Primo & Schruba, 2016), but _also_ the atomic gas phase (e.g.
Braun et al., 2009; Braun, 2012; Koch et al., 2018; Pingel et al., 2022), with
detailed information about kinematics, opacity and the cold/warm phases that
is not as reliably recovered in more distant targets.

We carry out our work with published and new datasets of M33, and organize our
paper as follows. In Section 2, we elaborate on our overall strategy of
observing the ISM around evolved massive stars. Section 3 describes the
published catalogs of RSGs, WRs and SNRs that we use, presents our new
observations of H I and CO with the VLA and ALMA. Section 4 describes our
results, including the relative spatial distribution of the stars and gas, and
the density distributions; we statistically verify these densities with multi-
wavelength star-formation maps. Section 5 discusses the implications of the
results, specifically the correlation between stellar age and ambient ISM
density, the significant fraction of massive stars in the low-density atomic
ISM, the need for higher ($\sim$pc) resolution molecular data to identify
cavities and substructures in the ISM around stars, the novelty of our data
for the purpose of validating feedback models, and comparison of our findings
with those in SN environment and SNR interaction studies.

##  2 Locating where supernovae explode with evolved massive stars

Figure 1: The spatial distribution of red supergiants (RSGs), Wolf-Rayet stars
(WRs) and supernova remnants (SNRs) with respect to the cold (H2 \+ H I) ISM
in the M33 ACA survey region (shown in the inset plot). Orange circles
represent RSGs with log $(L/L_{\odot})\geq 4.7$ (roughly corresponding to
$\gtrsim$12M⊙). Bluish circles show WRs and green stars show SNRs. H I is
shown in grey-scale with the density range shown in the bottom colorbar.
Contours of CO (2-1) emission (converted to $H_{2}$ surface densities) are
also shown at levels of 1.9, 7.6, 15.3 and 30.6 M⊙ pc-2. The prominent star-
forming regions in the ACA footprint – NGC 604, IC 142 and NGC 595 – are
specifically labelled (showing high WR concentrations). The black box shows a
region containing a WR and SNR with higher-resolution ALMA data in Figure 6.
The map shows the 50 pc-scale cold ISM environments where core-collapse SNe
will occur in the future (or have occured in the case of SNRs). Statistical
analysis of these environments and their implications are discussed in Section
4 and 5.

To trace the locations where stars will explode in M33, we use catalogs of
WRs, RSGs and SNRs. WRs and RSGs are evolved stages of OB stars, with expected
post-main-sequence lifetimes $\lesssim 1-4$ Myr ($\sim$3-10% of the main-
sequence lifetime of OB stars, e.g., Meynet & Maeder, 2005; Georgy et al.,
2013; Meynet et al., 2015), and SNRs are markers of recent SNe, with lifetimes
of only 104-105 yrs after SN (Sarbadhicary et al., 2017, e.g.). All three of
these tracers are expected to be spatially and temporarily ‘near’ the location
of their associated SNe. For typical RSG and WR lifetimes, we expect stars to
migrate $\sim$$(\mathrm{10pc})(v_{*}/10\mathrm{km/s})(t_{*}/1\mathrm{Myr})$
from their current location, assuming a characteristic velocity dispersion of
10 km/s. This migration of even less of a concern for SNRs given their even
shorter lifetimes. We now discuss these 3 types of source in a bit more
detail.

The reason for using RSGs is relatively straightforward – they have been
directly identified as the progenitors of Type II SNe, which form the majority
category of core-collapse SNe in flux-limited surveys (Smartt, 2009, 2015).

The use of WRs as SN locators deserves some discussion. WRs were classically
considered progenitors of stripped-envelope (SE) SNe, which form the second
dominant category of core-collapse (CC) SNe after Type II (Smartt, 2009,
e.g.), and are known to generate powerful winds to blow off their outer
H-envelopes (Conti et al., 1988; Woosley et al., 2002; Crowther, 2007). This
picture however has been called into question owing to the faintness of SE-SN
progenitors compared to LMC WRs (e.g. Eldridge et al., 2013), low ejecta
masses based on narrow light curves (e.g. Drout et al., 2011; Taddia et al.,
2015; Lyman et al., 2016; Taddia et al., 2018), SN rates in tension with
single stellar models (Smith et al., 2011; Graur et al., 2017), and decreasing
explodability of higher mass stars (Pejcha & Thompson, 2015; Sukhbold et al.,
2016, e.g.), leading to suggestions that SE-SN progenitors are mostly lower-
mass He-stars stripped by binary interactions. On the other hand, at least
$\sim$10% of SE-SNe (Taddia et al., 2019; Karamehmetoglu et al., 2022),
particularly those associated with long-gamma-ray bursts (e.g. see Crowther,
2007; Smartt, 2015; Roy, 2021, and references therein) and strong
circumstellar interaction (e.g. Van Dyk et al., 2018; Xiang et al., 2019;
Horesh et al., 2020; Gal-Yam et al., 2022; Kool et al., 2021; Pellegrino et
al., 2022) are still consistent with WR-like progenitors.

Given the ongoing research about SE-SN progenitors, we justify our use of WR
tracers as follows. Firstly, to first-order, our WRs at least trace the
locations of the most massive stars in M33. For reference, at LMC-like
metallicities, which is comparable to the outer portions of M33 (Magrini et
al., 2007), WRs are expected to be $>$40 Msun ($>$20 M⊙) based on single
(binary) evolution models (Eldridge et al., 2017a). In the inner regions of
M33, where the metallicity is higher, WRs form at slightly lower masses, but
are still among the most massive stars in the surrounding environment.
Secondly, by definition, WRs are tracers of where significant pre-SN feedback
is taking place, since WRs are closely associated with OB stars, and their
winds carry substantial mechanical energy to influence their environments.
Lastly, WRs and SE-SNe appear to evolve in environments with very similar
star-formation rates (e.g., Kangas et al., 2017; Maund, 2018), so we expect
the explosion sites to be similar at our resolution.

Complementing WRs and RSGs are M33 SNRs, which denote the locations where
stars have already exploded in the past $10^{4}-10^{5}$ years (and thus
already cleared their ambient gas out to some radius, unlike the massive stars
which will explode in the future). In particular, the sample of SNRs in M33 is
among the most complete population known in the Local Group, with substantial
observational programs spanning optical to radio wavelengths that have been
dedicated to the galaxy (Long et al., 2010; Lee & Lee, 2014; Long et al.,
2018; White et al., 2019).

While OB star catalogs are also available in the Local Group galaxies, we use
their evolved counterparts to identify SN-related environments instead for a
few reasons. One is that OB stars are in the H-burning main-sequence phase
with typical durations about 10 times longer than the post-main-sequence
phases where RSGs and WRs manifest (Ekström et al., 2012). As a result, for
the same velocities, OB stars would travel 10 times farther before exploding,
so their explosion sites have a larger ‘error’ circle from the present
location, compared to RSGs and WRs (for reference, a star moving at a velocity
$v_{*}$ will move
$(\mathrm{10pc})(v_{*}/10\mathrm{km/s})(t_{*}/1\mathrm{Myr})$. Additionally
(and arguably the bigger issue) is that existing OB catalogs rely on detailed
spectroscopic follow-up for classification, and existing catalogs in these
galaxies are believed to be largely incomplete given the number counts of WRs
and RSGs (Massey et al., 2016), so currently the only feasible way to probe
the locations where massive stars explode is with their evolved counterparts.

Finally, we note that our measured ambient densities on 50 pc scales should be
representative of the densities at the time of explosion in most regions since
the turbulent sound-crossing timescale, for a typical $\sigma=10$ km/s, is
about 5 Myrs, longer than the expected time left for the evolved massive stars
to explode ($<$0.5 Myrs for stars $>$30 M⊙, and $\sim$0.5-4 Myrs for 8-30
M⊙RSGs, Ekström et al., 2012). SN blastwaves however can alter the gas
distribution on a much shorter timescale due to their faster speeds. We return
to this point in the discussion in Section 5.1.

##  3 DATA

###  3.1 Stellar Catalogs

####  3.1.1 Wolf-Rayet Stars

We use the M33 WR catalog of Neugent & Massey (2011). The study cataloged 206
WRs within 0.32 deg2 ($\sim$ 8.5 kpc radius) of M33 using observations from
the 4 m Kitt Peak Mayall telescope and Mosaic CCD. Candidate WRs were
identified with crowded-field photometry and image subtraction with narrowband
filter images centered on C III $\lambda$4650 and He II $\lambda 4686$ lines
(typically seen in WRs), and then spectroscopically confirmed using the 6.5 m
MMT Hectospec instrument. Objects were classified as WN (strong lines of He,
N) and WC (strong lines of He, C), with further subclassification (e.g.
WN3-WN8) based on relative line strengths of He, C, N and O ionization states
as per convention (Conti et al., 1988).

The WR sample is expected to have high completeness according to Neugent &
Massey (2011), who estimated a 95% completeness based on comparisons with
previous studies in the Local Group (Massey & Johnson, 1998; Neugent et al.,
2018). Neugent & Massey (2023) further confirmed the high completeness of the
M33 WRs; the faintest WR star in the catalog is still about $\sim$1 mag
brighter than their nominal sensitivity limit. We also note from calculations
in Section 4.1 that even in regions like NGC 604 that are expected to have
high extinction, the observed number of WRs are consistent with expectations
from the underlying star-formation history in the region (Lazzarini et al.,
2022). The MMT spectroscopy also ensured high sample purity, removing cases
where previously-identified WRs in NGC 604 and NGC 595 (Drissen et al., 2008)
were revealed to be O-type stars with WR-like emission lines from the
spectroscopy of Neugent & Massey (2011). We do not subdivide our WR sample
into WC/WN (as provided in the catalog) since these subcategories have similar
masses, and we are primarily interested in the densities where stars in the
general mass range of WRs explode.

Of the 206 WRs, we restrict our catalog to 119 objects within the footprint of
our CO data (described in Section 3.2).

####  3.1.2 Red Supergiant Stars

We use the RSG catalog constructed by Massey et al. (2021b) from a UKIRT
$J-K_{s}$ band survey of the inner 3 deg2 of M33. Foreground stars were
separated from M33 stars using _Gaia_ proper motion and parallax measurements,
and RSGs were separated from asymptotic giant branch (AGB) stars and
background galaxies using cuts in $J-K_{s}$ colorspace, constructed to
optimize RSG selection based on previous tests in the Magellanic Clouds (Yang
et al., 2019; Neugent et al., 2020). The final catalog consists of 7088 RSGs,
of which we use the 812 RSGs inside our CO footprint. Physical quantities like
luminosity and temperature were also reported by Massey et al. (2021b) after
correcting for extinction and reddening, and appear to be in qualitative
agreement with predictions of stellar evolution models (Ekström et al., 2012).

####  3.1.3 Supernova Remnants

We use the SNR sample in M33 compiled in White et al. (2019), which is the
latest and most complete sample of SNRs in the galaxy. The SNRs were
discovered using optical-line emission data, specifically those with [S
II]/H$\alpha\gtrsim 4$ (to separate from photoionized nebulae) using data from
4m Mayall Telescope and 0.9m Burrell Schmidt telescope at KPNO ratios (Long et
al., 2010; Lee & Lee, 2014). These candidates were further spectroscopically
confirmed as SNRs with 6.5 MMT (Long et al., 2018), Chandra X-ray observations
(Long et al., 2010) and sensitive JVLA observations at 1.4 and 5 GHz taken in
White et al. (2019), covering roughly the inner 0.85 deg2 (0.35 deg2 at 5 GHz)
of M33. Of the 215 SNRs in the sample, we use the 71 inside the footprint of
the CO map.

Table 1: Surface densities of detectable H2 and H I around different categories of stars, and the corresponding fraction of stars with H2/HI non-detections, as measured in Section 4.2. The first row are values for all pixels in M33, and subsequent rows denote pixel values of stars in the different categories considered in the paper. Columns from left to right: Category of star, sample size ($N$), the median gas surface density ($\langle\Sigma\rangle_{\mathrm{med}}$), the 16th-84th percentile range of surface densities ($\langle\Sigma\rangle_{\mathrm{16-84}}$), and non-detection fractions ($f^{\mathrm{non}}$) of H2, H I and the total (H2 \+ H I) gas column respectively. All median and percentile calculations are calculated for values higher than the respective completeness limits.  Category | $N$ | $\langle\Sigma_{\mathrm{H_{2}}}\rangle_{\mathrm{med}}$ | $\langle\Sigma_{\mathrm{H_{2}}}\rangle_{\mathrm{16-84}}$ | $f^{\mathrm{non}}_{\mathrm{H_{2}}}$ | $\langle\Sigma_{\mathrm{HI}}\rangle_{\mathrm{med}}$ | $\langle\Sigma_{\mathrm{HI}}\rangle_{\mathrm{16-84}}$ | $f^{\mathrm{non}}_{\mathrm{HI}}$ | $\langle\Sigma_{\mathrm{tot}}\rangle_{\mathrm{med}}$ | $\langle\Sigma_{\mathrm{tot}}\rangle_{\mathrm{16-84}}$ | $f^{\mathrm{non}}_{\mathrm{tot}}$
---|---|---|---|---|---|---|---|---|---|---
|  | (M⊙ pc-2) | (M⊙ pc-2) |  | (M⊙ pc-2) | (M⊙ pc-2) |  | (M⊙ pc-2) | (M⊙ pc-2) |
M33 | – | 6.8 | 2.5–20.5 | 0.75 | 7.5 | 3.6–12.4 | 0.11 | 8.9 | 4.5–18.0 | 0.17
WRs | 119 | 15.5 | 3.8–33.8 | 0.45 | 9.7 | 5.0–14.0 | 0.08 | 16.7 | 6.5–39.8 | 0.13
RSGs (all) | 812 | 6.7 | 2.4–18.9 | 0.72 | 6.8 | 3.5–11.4 | 0.11 | 8.5 | 4.5–17.7 | 0.18
RSGs ($>10^{5}$ L⊙) | 98 | 11.6 | 4.8–30.5 | 0.56 | 8.4 | 4.8–13.5 | 0.09 | 11.2 | 5.4–26.2 | 0.12
RSGs ($10^{4.7}-10^{5}$ L⊙) | 178 | 6.4 | 2.4–17.5 | 0.68 | 6.6 | 3.2–11.3 | 0.1 | 9.2 | 4.3–18.9 | 0.18
RSGs ($<10^{4.7}$ L⊙) | 536 | 4.7 | 2.1–18.2 | 0.77 | 6.5 | 3.4–11.1 | 0.12 | 7.9 | 4.3–15.3 | 0.18
SNRs | 71 | 7.4 | 4.5–33.8 | 0.68 | 8.8 | 5.9–13.8 | 0.04 | 10.8 | 6.4–18.3 | 0.04

###  3.2 ISM maps

####  3.2.1 21 cm Atomic (H I) Map

We use two H I maps of M33, both obtained from the VLA and using short-spacing
observations from the GBT (Lockman et al., 2012). The coarser 20″ (80 pc)
resolution map is the 13-pt mosaic from Koch et al. (2018, project ID:
14B-088) using C configuration observations, which has a spectral resolution
of $\sim$0.2 km s-1 and rms noise of 2.8 K per channel. The data are feathered
with the GBT H I data from Lockman et al. (2012) with a single dish scaling
factor of unity (Koch et al., 2018).

The second map uses new VLA B configuration observations (project ID: 17B-162)
combined with the C configuration from Koch et al. (2018). We briefly describe
these data here, and refer readers to further details in E. Koch et al. (in
preparation). We use the inner 7 pointings from the C-configuration
observations and combine these with new B-configuration observations of the
same pointings. We image the data in 1.2 km s-1 channels with natural
weighting which yields a beam size of 8″(32 pc). The rms noise per channel is
6.2 K per 1.2 km s-1 channel. We use the same feathering approach as described
above with the C-configuration mosaic to incorporate the GBT H I short-spacing
information.

In this work, we use the integrated intensity (moment 0) maps calculated from
both spectral-line data cubes following the signal masking approach from Koch
et al. (2018). Briefly, we use a minimum $2\sigma$ threshold on the data, and
keep only regions that contain a $>4\sigma$ peak and span 10 consecutive
spectral channels, as appropriate for the typical H I line widths sampled with
$1$ km s-1 channel. Following the findings from Koch et al. (2021), we assume
the H I is generally optically thin on these spatial scales ($30\mbox{--}80$
pc) and convert to the atomic gas surface density in M⊙ pc-2 using

$\frac{\Sigma_{\mathrm{HI}}}{\mathrm{M}_{\odot}\ \mathrm{pc}^{-2}}=0.0196\
\mathrm{cos}(i)\left(\dfrac{I_{\mathrm{HI}}}{\mathrm{K\ km\ s^{-1}}}\right)\\\
$ (1)

where $i=55^{\circ}$ (Koch et al., 2018) is the inclination angle of M33. The
noise minimally varies across our field of interest, and so we adopt a median
3$\sigma$ completeness limit of 0.94 M⊙ pc-2. The cubes also provided a map of
the line-width of HI (km/s), which enters our calculations of the scale height
and total gas volume density, as discussed in Appendix A.

####  3.2.2 ACA CO(2-1)-traced H2 map

We trace the molecular gas in M33 using ALMA Band 6 12CO(2–1) observations
obtained by the 7-m array in the Atacama Compact Array (ACA; project IDs:
2017.1.00901.S; 2019.1.01182.S). The complete observations, including CO
isotopologues, are explained in detail in E. Koch et al. (in preparation).
Briefly, the 12CO(2–1) observations are taken in 15 individually-observed
mosaics that we image separately using a deconvolution approach similar to the
PHANGS-ALMA imaging pipeline (Leroy et al., 2021b). We then use linear
mosaicking to combine each individually observed field into one complete
mosaic covering the entire galaxy (see Figure 1). Because of the northern
location of M33 with respect to the ALMA site, the synthesized ACA beam is
highly elongated in some observations (e.g., up to a factor of 3 between the
major and minor axes). The map we use here is convolved to the coarsest common
round beam, giving us a final resolution of 12″($\sim$49 pc). We note that
this matches the resolution of the IRAM 30-m CO(2–1) map (Druard et al.,
2014), which maps a far larger region of M33. We verified that we recover
consistent fluxes between the maps, demonstrating that no significant diffuse
CO component is filtered out of the ACA data (Koch et al. in prep.). However,
the ACA data have a factor of 4 times finer spectral resolution (0.7 km s-1)
and a factor $\sim 2$ higher sensitivity (25 mK per 0.7 km s-1 channel).

We compute the CO(2–1) integrated intensity after calculating a signal mask
with an approach similar to the “strict” method from Leroy et al. (2021b),
which we note does mildly differ from the approach used for the H I but has
been well-tested for CO emission that is clumpier with narrower line widths
than the H I. As such, the signal masking approaches used for H I and CO are
appropriate given the expected properties (diffuse, extended vs. clumpy) of
each tracer. We convert the CO intensity to a mass surface density using,

$\frac{\Sigma_{\mathrm{H_{2}}}}{\mathrm{M}_{\odot}\
\mathrm{pc}^{-2}}=\alpha_{\mathrm{CO}}R_{21}^{-1}\mathrm{cos}(i)\left(\dfrac{I_{\mathrm{CO}}}{\mathrm{K\
km\ s^{-1}}}\right)\\\ $ (2)

where we assume $\alpha_{\mathrm{CO}}=4.3\
\mathrm{M}_{\odot}\mathrm{pc}^{-2}(\mathrm{K\ km\ s^{-1}})^{-1}Z^{-1.6}$ is
the empirical conversion factor between CO intensity and H2 surface density,
$Z$ is the metallicity normalized to the solar value,and $R_{21}=0.65$ is the
CO (2-1)/(1-0) line-ratio (Sun et al., 2018) . Since $\alpha_{\mathrm{CO}}$
can vary with metallicity, we also checked our results with a variable
$\alpha_{\mathrm{CO}}$ prescription, namely one where $Z$ varies with
galactocentric radius based on the measured metallicity gradient in M33
(Bresolin, 2011). The difference in $\Sigma_{\mathrm{H_{2}}}$ between our
variable and constant $\alpha_{\mathrm{CO}}$ is $\sim$30-40%, much smaller
than the orders of magnitude variation in $I_{\mathrm{CO}}$ and
$I_{\mathrm{HI}}$. We therefore continue using the constant
$\alpha_{\mathrm{CO}}$ prescription for simplicity. We assume $Z=0.6$ Z⊙,
appropriate for the inner part of M33, based on nebular abundances in the
inner 5 kpc of M33 (Rogers et al., 2022). From our noise cubes, we obtain a
median 3$\sigma$ noise (assuming detection over 3 channels) of 0.16 K km s-1,
and therefore a median completeness limit of our H2 maps of 1.37 M⊙ pc-2.

Figure 2: Cumulative distribution functions showing H2, H I, and total (H2 \+
H I) surface densities. Distributions are shown for the locations of WRs (blue
lines), SNRs (green lines) and RSGs (red lines), while grey histograms denote
the distribution of densities of the full M33 region, representing what would
be expected for a random distribution. Vertical dashed lines denote the
completeness limit of each dataset. Colored circles show the median of the
distribution above the completeness limit. The distributions show a higher
probability of more massive stars (e.g. the WRs) to evolve in denser gas.
Figure 3: Cumulative distribution functions showing H2, H I and total gas
surface densities (as in Figure 2). Here only the RSGs are plotted, dividing
them into 3 luminosity bins, roughly tracing mass ranges of 8-12 M⊙(peach
lines), 12-15 M⊙(orange lines) and $>$15 M⊙(red lines). Grey histograms denote
the distributions for the full M33 region. Colored circles show the median of
the distribution above the completeness limit.

####  3.2.3 High-resolution ALMA CO(2–1)-traced H2 map

We also make use of a single ALMA field observed around a WR and SNR in close
proximity (Figure 1, box region) to highlight the finer substructure in cloud
environments of SNe that can be revealed at much higher resolution.
Specifically, we look at the region around SNR L10-045 from ALMA Band 6 data
published in Sano et al. (2021, project ID: 2018.1.00378.S). The area observed
is a single ALMA 12-m pointing. These data resolve 0.46″($\sim 2$ pc) scales
in 12CO(2-1) and so resolve the filamentary structure within this single GMC.
The data we show here are re-processed, using the ALMA delivered QA2 data from
the archive with imaging using the PHANGS-ALMA pipeline (Leroy et al., 2021b).
We image the data using Brigg’s weighting with robust=0.5 in $1.2$ km s-1
channels. The corresponding rms is 38 mK per $1.2$ km s-1 channel. We signal
mask and generate moment maps using the PHANGS-ALMA pipeline (Leroy et al.,
2021b); specifically we use the integrated intensity map with strict signal
masking.

###  3.3 Galaxy simulation

In Section 5.4 we will compare our measured ambient densities (cm-3) around
WRs and RSGs with the ambient densities of SNe in a hydrodynamical simulation
of a dwarf starburst galaxy system in Lahén et al. (2020). We provide here a
brief summary of the simulation, and refer the reader to the paper for
details. The simulation follows the merger of two identical gas-rich dwarf
galaxies of initial stellar masses of $2\times 10^{7}$ M⊙, at a baryonic mass
resolution of $\sim 4$ M⊙ and sub-pc spatial resolution. We chose this
simulation for our example because the high baryonic resolution ($\sim$4 M⊙)
leads to well-resolved ISM environments (Steinwandel et al., 2020), the
starburst system results in a wide range of ISM densities and star formation
rates, and because each SN event can be directly traced back to a single
progenitor star. The simulation was performed using SPHGal (Hu et al. 2017,
and references therein), an improved version of the smoothed particle
hydrodynamics code GADGET-3 (Springel, 2005), which uses a simplified chemical
network to model the non-equilibrium cooling and heating processes within the
ISM including H2 formation and carbon chemistry. Stellar feedback processes
driving the multiphase ISM include the interstellar far-ultraviolet radiation
field, photoionizing radiation, core-collapse supernovae and winds of
asymptotic giant branch stars, all modelled according to the stochastically
sampled IMF (for details see Hu et al. 2017; Lahén et al. 2020). Massive stars
($>4$ M⊙) are realized as single particles in the simulation, allowing us to
measure the densities around each SN progenitor right before the explosion.

We will compare both the “local” and “average” ambient densities of SNe in the
simulation with the observations. Here “local” refers to the density
immediately surrounding the SN particles, which we obtain from the weighted
average density of the 100 closest gas particles. However, the local densities
of stars in observations are not directly measurable. To better compare the
simulation with observations, we obtain the column-averaged densities by
projecting the gas distribution in each simulation snapshot onto a grid along
the z-axis, at a resolution of 50 pc per pixel similar to our observation, and
then dividing the surface densities by a scale height of 100 pc. The gas
density maps are recorded in 1 Myr intervals, and the SN densities are
measured in the pixels that correspond to the locations of the evolved massive
stars ($8$–$50$ M⊙) in the snapshot prior to the SN explosion, i.e. up to 1
Myr earlier. We record densities for all SN progenitors produced during the
pre-starburst phase spanning 70-140 Myrs in the simulation when the two
galaxies are relatively quiescent (see Fig. 5 in Lahén et al. 2020). We expect
this phase to be a better match to M33, which is currently not undergoing a
starburst. We divide these SN progenitors into those with initial stellar mass
of 8-30 M⊙ and $>$30 M⊙ to directly compare with the RSG and WR stars.

##  4 ANALYSIS

With the datasets described in Section 3, we quantify differences in the ISM
properties of the different populations by comparing the spatial distributions
of RSGs, WRs and SNRs with respect to the H I and H2 distribution within the
ACA area (Figure 1). Surface densities of H2 and H I at the pixel locations of
stars are summarized in Table 1 and shown in Figures 2 and 3.

###  4.1 Spatial distribution of stars and gas

The relative distribution of the cold (H2 \+ H I) ISM and stars at $\sim$50 pc
spatial resolution is shown in detail in Figure 1, and the corresponding
density ranges are summarized in Table 1. Locations of the stars and SNRs are
shown as colored circles/stars, and the H I and H2 distributions are shown as
greyscale and contours (respectively) in units of surface density.

The H I distribution fills most of the survey area in Figure 1 compared to H2,
and shows the characteristic frothy structure typically seen in nearby
galaxies (e.g Stanimirovic et al., 1999; Walter et al., 2008; Pingel et al.,
2022). In comparison, the H2 contours are clumpier, coinciding with dense
molecular clouds in active star-forming regions. The detailed structure of the
H I and H2 of M33 will be covered in a forthcoming publication (Koch et al, in
prep). We note that M33 is about an order-of-magnitude lower in stellar mass
and average molecular surface density compared to the population of nearby,
star-forming galaxies in PHANGS analyzed by MC22 (Leroy et al., 2021a). Our
work is therefore complementary to MC22 in that it probes SN environment
densities in a lower-mass dwarf spiral galaxy.

While the stars generally follow the gas distribution, we note the stark
contrast in the relative spatial distribution of the gas and WRs, compared to
SNRs and RSGs. Dense concentrations of WRs in Figure 1 are particularly
prominent in major star-forming regions such as NGC 604, NGC 595, and IC 142,
as well as the numerous H II regions interspersed in the southern arm. This WR
distribution is a result of the abundant young ($<$10 Myr) stellar populations
in these regions. For example, the total stellar mass formed in the last
$\sim$10 Myrs in NGC 604 is $\left(1.06_{-0.42}^{+0.67}\right)\times 10^{5}$
M⊙ based on the PHATTER survey (Lazzarini et al., 2022, note that maps measure
star formation histories in $\sim$100-pc-sized cells). The production rate of
WRs is roughly 1.3 WR per 104 M⊙(Eldridge et al., 2017b; Dorn-Wallenstein &
Levesque, 2018a), which gives an expected number of WRs in NGC 604 to be
$14^{+8}_{-6}$, consistent with the 11 observed in the region. In contrast,
the 10 Myr stellar mass of IC 142 is about 7 times lower, so only about 2 WRs
are expected in the region, consistent with Figure 1. NGC 595 on the other
hand has the same number of WRs as NGC 604 but with a factor of 2 smaller
stellar mass, likely due to a more complex dependence on the underlying
metallicity and star-formation history (Drissen et al., 1993). This tendency
of WRs to concentrate in actively star-forming regions is likely due to their
higher zero-age main-sequence masses ($>$30 M⊙) and leads to a stronger
correlation of the WR population with dense molecular gas compared to RSGs
($\gtrsim$8-30 M⊙) and SNRs, which we discuss in more detail in Section 4.2.

###  4.2 ISM densities at the location of stars

We now compare the surface density distributions222The density measurements
will be made publicly available after acceptance and publication of the
manuscript. Please feel free to contact us if you would like to use the
dataset beforehand. of gas tracers at the locations of different categories of
stars in Figure 2. We also list their median, 16th-84th percentile range
(corresponding to a 1$\sigma$ interval), and fraction of non-detections (i.e.,
pixel values below the gas surface density completeness limit) in Table 1. The
median and percentiles are estimated only using values above the completeness
limits.

To visualize the distributions, we plot the cumulative distribution function
(CDF) of H2 and H I surface densities for WRs, RSGs and SNRs in Figures 2 and
3. The CDF captures the same information as a standard histogram or kernel
density estimation without the need to adopt bins or a smoothing kernel. It
allows easy comparison between our samples, which have very different overall
sample size. The CDFs flatten to 1 at the highest pixel value of the maps,
i.e., 100% of the pixels are below that maximum value. On the other end, the
CDFs flatten to a limiting fraction value as they reach the completeness
limits of the maps (vertical dashed lines). One should read this limiting
value as the fraction of that object sample without a detection in the
associated gas type. This limiting value is quite high for the H2 distribution
of M33 ($\sim 75\%$), reflecting the fact that the molecular gas is
concentrated into individual dense clouds, and as a result CO emission covers
only a small fraction of the total area in M33 (Figure 1). On the other hand,
atomic (H I) gas has a much larger filling fraction, and nearly 90% of the
pixels have some H I detection.

Notably, we find WRs coincident with H2 column densities on average 3 times
higher than that of all pixels in M33 and those containing RSGs, and almost 2
times higher than pixels with SNRs (Table 1). Only 4% of the RSG pixels and
8.5% SNR pixels have densities $\geq 20$ M⊙ pc-2 (typical for molecular clouds
in M33, Rosolowsky et al., 2007) while 21% of the WR pixels are above that
value. Also, only 45% of the WRs are located in H2 column densities below our
detection limit (i.e. $<1.4$ M⊙ pc-2), compared to 68% and 72% for RSGs and
SNRs respectively (we will frequently refer to this as the ‘non-detection’
fraction).

These results suggest that WRs are much more closely associated with dense gas
than the bulk RSG and SNR population, which is consistent with the idea that
WRs arise from more massive progenitors that have shorter lifetimes, and thus
more likely to be found near their birth clouds. A significant fraction of WRs
however ($\sim$45%), appear outside any detectable molecular gas, which has
important implications for feedback that will be discussed later in Section 5.

Unlike H2, most of our stars are evolving in pixels with some detectable
atomic gas. The majority of the pixels in our ACA region have H I detections
down to densities $\sim$$0.94$ M⊙ pc-2. Of the 45% of WRs (53) without a
significant H2 detection, 83% (44 out of 53) are in pixels with detectable
atomic gas. The statistics are also similar for RSGs and SNRs: about 85% RSGs
with H2 non-detections (499 out of 587) and 94% SNRs with H2 non-detections
(45 out of 48) are in pixels with detectable atomic gas. About $\sim$10% of
the stars in all categories are H I non-detections, but as shown with more
sensitive low-resolution H I data from Koch et al. (2018) in Appendix D, there
is still some atomic gas at the locations of these stars, just with surface
densities below our sensitivity limit of $\Sigma_{HI}\lesssim 1$ M⊙ pc-2.

The resulting _total_ surface density of neutral gas
($\Sigma_{HI}+\Sigma_{H_{2}}$), which would be of interest in simulations, is
shown in the bottom left panel of Figure 2, clearly showing a mix of the H I
and H2 distribution properties. Almost 80-90% of stars in the three categories
are evolving in column densities $\gtrsim$2 M⊙ pc-2, but WRs are still found
in environments almost twice as dense as RSGs and SNRs that have lower-mass
progenitors than WRs.

The correlation between progenitor mass and ISM densities is more explicity
shown in Figure 3. Here again we show the H I and H2 CDFs similar to Figure 2
but only for RSGs, divided into three luminosity (and therefore roughly
progenitor mass) bins: $\mathrm{log}(L/L_{\odot})=4.2-4.7$,
$\mathrm{log}(L/L_{\odot})=4.7{-}5$, and $\mathrm{log}(L/L_{\odot})>5$. Based
on comparison with stellar evolution models shown in Figure 8 of Massey et al.
(2021b), these ranges roughly correspond to zero-age main-sequence masses of
$8-12$M⊙, 12–15M⊙, and $>$15M⊙. From here on, we will refer to these RSGs by
their mass ranges. The lower limit of $\mathrm{log}(L/L_{\odot})=4.2$ also
helps to avoid contamination from bright AGB stars. We find that, similar to
the results in Figure 2, a greater fraction of higher mass RSGs are evolving
in denser H2 and H I gas, with the differences in H I gas distributions being
less pronounced than H2. Notably, only 56% of RSGs $>15$M⊙ have H2 non-
detections, compared to the 77% of $8-12$M⊙ RSGs having H2 non-detections.
This is similar to how WRs have a smaller fraction of H2 non-detections than
the full RSG sample.

####  4.2.1 Statistical significance of the distributions

Figures 2 and 3 suggest systematic differences between the gas distributions
found at the locations of different evolved stars. But are these differences
statistically significant?

We first check this with the two-sample Kolmogorov-Smirnov (KS) Statistic,
which tests the similarity between two unequal-sized distributions by
measuring the maximum distance between the sample CDFs. While the KS test has
weak sensitivity to differences at the tails of distributions, it is quite
sensitive to the difference in medians. We compare the H I and H2 surface
density CDFs of WRs, the three subcategories of RSGs, and SNRs with that of
the bulk H2 and H I distribution of M33 using the KS test, and assess whether
the null hypothesis (that the locations of massive stars are drawn from the
bulk H2 and H I distribution) can be rejected using $p$-values from the tests
(we use a significance level of 0.05 for this purpose). 333We implement this
using the scipy.stats package, which returns the KS statistic and two-sided
$p$ value, which are listed in Table 2.

Table 2: Summary of KS statistic tests for different star categories, specifically comparing their H2 and H I distributions with the bulk distribution of M33 as described in Section 4.2.1. Columns from left to right: Name of stellar category, KS Statistic for H2 surface density distribution, $p$-value for H2 KS test, KS Statistic for H I surface density distribution, and $p$-value for H I KS test. Boldfaced numbers emphasize $p$-values less than or equal to 0.05.  Category | KS (H2) | $p$ | KS (H I) | $p$
---|---|---|---|---
WRs | 0.33 | $\mathbf{9.2\times 10^{-12}}$ | 0.20 | $\mathbf{1.6\times 10^{-4}}$
RSGs (all) | 0.03 | 0.37 | 0.06 | $\mathbf{3.4\times 10^{-3}}$
RSGs |  |  |  |
– ($>15$ M⊙) | 0.22 | $\mathbf{1.5\times 10^{-4}}$ | 0.11 | 0.17
– ($12-15$ M⊙) | 0.08 | 0.24 | 0.08 | 0.15
– ($<12$ M⊙) | 0.04 | 0.22 | 0.09 | $\mathbf{4.3\times 10^{-4}}$
SNRs | 0.12 | 0.24 | 0.25 | $\mathbf{3.0\times 10^{-4}}$

The results are summarized in Table 2. For H2, we find that the WRs have the
largest difference between its CDF and the bulk M33 H2, and the null
hypothesis can be rejected with high confidence, at a significance level
$p\ll$0.05. A similar result is also found for the RSGs $>$15 M⊙, but not for
RSGs $<$15 M⊙ and SNRs, which appear at least marginally consistent with
having been drawn from the same parent population as the overall distribution
within our field of view. For H I, the WRs and SNRs have the largest CDF
difference, and the null hypothesis is rejected at the 0.05 level. The H I CDF
of high-mass RSGs, similar to H2, have the largest difference of the three RSG
categories with the bulk distribution, but it is not statistically significant
in this case. The H I of low-mass RSGs on the other hand appear to have a
statistically significant difference, but smaller than the high-mass ones. A
close inspection of Figure 3 shows that the low-mass RSGs are exploding in
slightly lower densities than both the bulk H I and H2 of M33. We also re-
checked this with the Anderson-Darling (AD) test, which is more sensitive to
the tails of the distributions. We find that the above cases where $p<0.05$
with the KS test also have $p<0.05$ with the AD test. The AD test returns two
cases implying statistical significance at $p<0.05$ where the KS test does
not: the H2 CDFs of the 12-15 M⊙RSGs, and SNRs.

These results provide an initial quantitative support (which we investigate
further in the next Section 4.3) that higher mass progenitors (as traced by
WRs and RSGs $>$15 M⊙) are indeed evolving in denser H2 pixels than the lower
mass stars (RSGs $<$15 M⊙), and their densities are not just random sampling
from the bulk distribution of M33. The same holds for the H I densities of
WRs, though the difference between the H I densities of RSGs $>15$M⊙and M33 is
not statistically significant. Visually, these results relate back to the
spatial distributions in Figure 1, where we see WRs appear more frequently
within H2 clouds marked by the contours, and sometimes well within giant
molecular clouds like in NGC 604 and NGC 595, compared to the RSGs.

###  4.3 Comparison with the general stellar population

Figure 4: Distribution of mock stars (purple) randomly drawn from different
stellar population tracers – uniform random (top left), 3.4 $\mu$m near-IR
tracing bulk stellar mass (top right), the 24 $\mu$m-corrected far-ultraviolet
(FUV) tracing the $<100$ Myr star-formation (bottom left) and 24
$\mu$m-corrected H$\alpha$ tracing the $<10$ Myr star-formation (bottom
right). Each map shows a population of $N=500$ randomly-drawn stars. The maps
are themselves shown in grey-scale. Orange contours denote H2 clouds with
densities $>$4.3 M⊙ pc-2 from the ACA survey. The figure visualizes where
stars from different progenitor populations occur, with their corresponding
densities shown in Figure 5. Figure 5: Cumulative density histograms of the
observed and modeled populations of RSGs, WRs, and SNRs from Section 4.3. Dark
grey histogram in each panel shows the bulk H2 distribution in M33 in the ACA
area. Solid colored line denotes the H2 surface density of the star category.
Colored shaded region shows the 5th-95th percentile region of the mock
populations drawn according to Section 4.3.

Are the above distributions therefore a consequence of the progenitor
properties of WRs, RSGs and SNRs? We test this by creating mock populations of
these stars from different star-formation tracer maps and comparing their
density distributions with what is observed. The primary goal is not to
determine the age distribution of RSGs, WRs and SNRs (which has been done
elsewhere), but rather to determine the degree to which their observed ISM
distribution is _astrophysical_ in nature, as opposed to a random alignment.
This will also assess the stochasticity of the CDFs due to the finite sample
sizes.

To create the mock populations, we use maps of tracers of star-formation on
different timecales. The probability ($p_{i}$) of a star occurring in pixel
$i$ is weighted by the pixel luminosity in the tracer map (listed below) with
respect to the total map luminosity, i.e., $p_{i}=L_{i}/\sum L_{i}$444Note
that by definition we are drawing stars in a light-weighted manner, not mass-
weighted, meaning the drawn populations are biased towards younger ages
instead of sampling the underlying star-formation history, which requires more
detailed modeling. Thus, in general, we are tracing progressively lower
mass/older stellar populations as we go from H$\alpha$ to FUV to near-IR
light, but the resulting stars have a light-weighted mass-range, which may not
accurately reflect the specific mass range of our RSGs, WRs and SNRs.. For
each category, we generate populations from four different tracers as follows:

  1. 1.

H$\alpha$+24$\mu$m – The H$\alpha$ emission line traces stellar populations
younger than 10 Myr, with a mean age of 3 Myr (Kennicutt & Evans, 2012). We
use _MIPS_ 24 $\mu$m maps (Verley et al., 2007) to correct the continuum-
subtracted H$\alpha$ map of M33 from Hoopes & Walterbos (2000) for
extinction555The native maps in emission measure (EM, cm-6 pc) units were
converted to flux density ($S_{\mathrm{H}\alpha}$, ergs s-1 cm-2 arcsec-2)
units using the relation $EM=4.858\times
10^{17}T_{4}^{0.9}S_{\mathrm{H}\alpha}$ (https://www.astronomy.ohio-
state.edu/pogge.1/Ast871/Notes/Rayleighs.pdf), where we assume temperature
$T_{4}=T/10^{4}\mathrm{K}=1$.. All maps are convolved and regridded to the ACA
12″resolution. We follow the extinction-correction prescriptions in Belfiore
et al. (2023) to carry this out (see Appendix B for details).

  2. 2.

FUV+24$\mu$m – Stars younger than 100 Myr account for about 90% of UV
emission, with a mean light-weighted age of about 10 Myr (Kennicutt & Evans,
2012). We use _GALEX_ FUV maps (Gil de Paz et al., 2007); both convolved and
regridded to the ACA 12″resolution, and extinction-corrected using MIPS 24
$\mu$m using the prescriptions in Belfiore et al. (2023), as discussed in the
Appendix B.

  3. 3.

Near-IR – A population tracing the bulk stellar mass of the disk, which we can
treat as a more realistic ‘random’ population, i.e., stars follow the actual
light profile of the galaxy. Analyses of resolved stellar populations in
nearby galaxies have shown that near-IR bands are primarily dominated by red-
giant and asymptotic-giant branch (AGB) stars in stellar populations with ages
up to a few Gyrs ($\gtrsim$2M⊙; Dalcanton et al., 2012), although some
contribution from evolved massive stars such as RSGs, red core-He-burning and
high-mass AGB stars are also expected (Melbourne et al., 2012; Melbourne &
Boyer, 2013). We use the bulk stellar mass map computed from the WISE W1 band
(3.4 $\mu$m) of M33 to trace locations of stars for this case (Leroy et al.,
2019). We refer the reader to Appendix A for details of the conversion from W1
band map to stellar mass. We use the common astrometric and beam-matched 7.5″
maps assembled in the z0MGs survey (Leroy et al., 2019) and convolve to the
12″ resolution and pixel scale of the ACA map.

  4. 4.

Random – A completely random population, assuming each pixel has the same
probability of hosting a star (i.e. $p_{i}=1/\sum L_{i}$). This is the
simplest control case – the positions of stars in our map cannot be 100%
random, since they will follow the galactic substructure, but this gives a
basic ‘control’ sample for comparison to the other populations.

The spatial distribution of the mock stars is shown in Figure 4. Stars drawn
from “FUV” and “H$\alpha$” maps (iii and iv) show more concentration along the
star-forming and H II regions, with a small fraction appearing in the inter-
cloud and inter-arm regions. As expected, the concentration of stars in H II
regions is much higher when drawn from the H$\alpha$ maps than the FUV maps.
In contrast, stars drawn from “Random” and “near-IR” are nearly uniformly
distributed throughout the ACA region, with the near-IR population slightly
more weighted towards the M33 center. However, even some of these random/near-
IR stars fall within H2 contours, which demonstrates how a fraction (seen by
the limiting fractions in Figure 5) of our samples could _appear_ to be
associated with molecular clouds by chance alignment at our 50 pc scales even
if they are not physically associated.

For each tracer (i–iv) and star category (WRs, RSGs and SNRs), we generate 500
mock populations each of size $N$ (where $N$ is the observed sample size of
the star category given in Table 1). We then compare the median, 5th and 95th
percentile range of H2 density of these randomly drawn stellar populations
from each tracer to the observed RSG, WR and SNR distributions. Note that we
only compare with H2, since the differences between the H I or H I +H2 CDFs of
the three categories are less prominent than H2 alone. Figure 5 shows the
comparison of the H2 CDFs between our randomly drawn mock populations and the
observed stars. The shaded (2$\sigma$) indicates the variance due to repeat
sampling, with larger variance seen for smaller samples. The CDFs of the mock
populations behave similar to the observed ones, flattening to 1 at the
maximum pixel values, and to a limiting value at pixel values approaching the
completeness limit of the CO map.

We find in Figure 5 that the H2 CDF of the WRs is most consistent with the
mock CDFs drawn from the H$\alpha$ maps, and in tension with the other
tracers. This is consistent with the results in Section 4.2.1 and Table 2,
which showed that H2 CDF of the WRs is statistically different from the bulk
M33 H2 distribution. The mock CDFs also reproduce a fraction of stars with no
detectable H2, similar to the WRs. This indicates that H2 distributions of WRs
and the general $<$10 Myr stellar population in M33 are statistically similar,
with a substantial fraction not coincident with molecular gas. The mock CDFs
in Figure 5 do indicate that stars drawn from H$\alpha$+24$\mu$m are
associated with slightly denser H2 than the WRs. In fact, as shown in Appendix
B and Figure 12, CDFs of stars drawn from H$\alpha$ without 24$\mu$m
correction are more consistent with the WR CDF, although both these cases are
consistent with WRs within the shaded error regions. We return to this point
in Section 5.

In contrast to WRs, the H2 CDFs of the RSGs are more consistent with the FUV
and near-IR-drawn populations. Specifically, the $>$15 M⊙ RSGs are best-
reproduced by the FUV light, while 12–15 M⊙ RSGs are reproduced partly by the
near-IR light and partly by purely random sampling. This is consistent with
the fact that FUV emission is weighted by younger stars than those
contributing to 3.4$\mu$m emission, and further supports that lower mass stars
are less likely to be associated with H2. Interestingly, the 12–15 M⊙ bracket
is more consistent with random/near-IR than the younger (FUV/H$\alpha$)
tracers. We expect 12–15 M⊙ stars to have ages of $\sim$10-20 Myrs (Ekström et
al., 2012), while the dominant contribution at $\lambda>1$ $\mu$m is from old
stellar populations with light-weighted ages approaching a few Gyrs (assuming
star-formation histories typical of local universe galaxies, see Conroy,
2013). This likely implies that by the time $<$15 M⊙ stars are evolving well
beyond the main sequence, they are almost unassociated with H2 gas, and any
spatial correlation with H2 is statistically indistinguishable from that of
the bulk older stellar population. The tension between the $>$15 M⊙ RSGs and
H$\alpha$ in contrast with the WRs is likely due to WRs having higher
progenitor masses than RSGs. Stellar evolution models predict that at higher
masses (and metallicities), the effects of binarity, rotation and mass-loss
limits the production of cool RSGs with respect to WRs (Eldridge et al.,
2017b; Massey et al., 2021a).

The analysis in this section also provides some context for understanding the
H I and H2 CDFs of SNRs. We see from Figures 2–5 that the H I and H2 CDFs of
SNRs are more similar to the lower mass RSGs (9-15 M⊙) than the WRs or high-
mass RSGs. Notably, the SNRs are inconsistent with being drawn purely from the
H$\alpha$-emitting population. These results are overall consistent with their
progenitor masses measured from the surrounding star-formation
histories666Note that ‘progenitor masses’ derived this way are a statistical
translation of the age distribution of stars around the SNRs into a single
stellar-mass. They are not a direct measurement of the progenitor mass, which
is difficult to obtain except for young ejecta-dominated SNRs, or ones with
light echoes. The masses can have uncertainties exceeding a factor of 2
depending on the stellar age distribution. However, this still remains the
_only_ way to obtain reasonable estimates of progenitor masses of individual
extragalactic SNRs.(Jennings et al., 2014; Díaz-Rodríguez et al., 2018;
Koplitz et al., 2023). The progenitor mass distribution peaks at $\sim$10 M⊙,
with a small tail towards 30–40 M⊙. According to Koplitz et al. (2023), about
40% of the SNRs having masses $<$8 M⊙, which may indicate SNe Ia or delayed
core-collapse SN origins (Zapartas et al., 2017). We do find in Appendix C and
Figure 13 that a 40% Type Ia contribution (assuming they are drawn from the
near-IR map) produces a slightly better match to the observed H2 CDF of SNRs
in Figure 5. In contrast, all the WRs come from high-mass ($>$30 M⊙)
progenitors, so they are more likely than the general SNR population to be
near dense H2 gas. A secondary effect could be that some of the SNRs have
cleared out H2 within their current shock radius ,though this is difficult to
confirm at our 50 pc resolution, which is larger than the mean diameters of
M33 SNRs (Lee & Lee, 2014; White et al., 2019). As mentioned in Section 5.3,
we have ongoing $\sim$pc-scale ALMA observations of M33, which we will use in
future papers to better confirm the fate of the gas within the SNR volume.

Overall, the analysis in this section drives home two important points – that
the density distributions of our stars in Figure 2 and 3 are indeed
statistically unique and not just due to chance alignment with the ISM
(consistent with Section 4.2.1), and that their observed density distributions
are _astrophysical_ in nature, consistent with the densities in which their
general progenitor populations would also be evolving, with some stochasticity
due to finite sample sizes.

##  5 Where do stars explode in the ISM? – The general picture

In this section, we weave the above results into a general picture of where
stars are exploding and the physical mechanisms driving these trends, the
implications for subgrid models of feedback in hydrodynamical simulations, and
comparison with existing literature on SNe and SNR environments.

###  5.1 Correlation between ISM density and progenitor age

Figure 6: CDFs of the average cold (H I + H2) gas volume densities for WRs,
SNRs and RSGs, derived according to Appendix A. Grey CDF denotes the bulk cold
density distribution in M33. See Section 5.1 for details. Figure 7: Average
cold (HI + H2) gas volume density, sampled at 48 pc resolution, vs progenitor
mass range of red supergiants and Wolf-Rayet stars. The grey shading indicates
the 16th-84th percentile range. This gives a rough idea of the progressively
increasing density expected around massive stars with younger ages.

One of the major implications of our results is that there is a clear
progenitor age/mass dependence on the ISM densities where stars explode. WRs,
which are the most massive (and shortest-lived) of the three categories, have
the highest sample fraction (55%) evolving in dense, molecular gas-rich
environments. More luminous (and thus more massive and shorter-lived) RSGs
also have a higher fraction (44%) of stars evolving in dense, molecular-gas
dominated environments than less luminous RSGs. The progenitor mass
distribution of SNRs peaks at 10 M⊙(Section 4.3) which is much lower than WRs,
and there may be some Type Ia SNRs with lower mass progenitors. We therefore
see a larger sample fraction of SNRs with non-detections than WRs (Figures 2,
3). Similar trends are also observed for these stars in relation to age-
sensitive stellar population tracers as discussed in Section 4.3.

The simplest physical interpretation of this trend is that more massive stars
have shorter lifetimes, and are thus more likely to explode close to their
birth clouds, or before the clouds have been dissipated. We caution however
that this fraction of stars correlated with dense gas at our 50 pc resolution
should be treated as an upper limit to the actual number of stars that will
interact with dense gas. While the H2 cloud densities on these spatial scales
may not be significanty modified by turbulence before our WRs and RSGs explode
(Section 2), they may get modified in regions of high gas density and dense
stellar clustering (e.g. NGC 604) by powerful pre-SN feedback and frequent SN
episodes. Thus some of the massive progenitors in such dense regions may
explode in evacuated regions (which can be better examined with higher-
resolution CO data as discussed in Section 5.3), while others might still end
up interacting with dense gas (indeed many SNRs have been observed to be
interacting with dense molecular gas, Section 5.5).

We can also express the surface density distributions in units of volume
densities to make them more accessible to simulations. In contrast to surface
densities however, which are directly measured from the line fluxes, volume
densities require some assumption about the line of sight distribution of gas,
and thus are somewhat model-dependent. Any line of sight through the turbulent
ISM disk of a star-forming galaxy will coincide with bubbles and
inhomogeneities that are not captured in projection, but simulations have
shown such disks evolve to an equilibrium state where the gravitational
pressure of dark matter, stars, and gas is balanced by thermal and turbulent
pressure primarily from stellar feedback (Ostriker et al., 2010; Ostriker &
Shetty, 2011b; Kim et al., 2011). Assuming such an equilibrium in M33, we
derive H I and H2 scale heights and volume densities in our ACA region in
Appendix A in each pixel given the locally measured stellar mass, gas surface
densities, and velocity dispersions. These should be treated as densities
averaged over 50-pc beams.

Figure 6 shows the volume densities of three star categories. All stars have
average ambient densities between $0.1-50$ cm-3. The general trends from
Figures 2 and 3 are still noticeable. WRs once again evolve in the highest
range of densities compared to the SNRs and RSGs. In Figure 7, we show the
average densities for WRs and RSGs, specifically the median and the 16th-84th
percentile range vs their approximate range of single stellar progenitor
masses. We also subdivide the RSG sample into the three mass bins of $8-12$
M⊙, $12-15$ M⊙and $>15$M⊙ as described in Section 4.2. While each category
evolves in an order-of-magnitude range of densities, one can clearly see a
systematic increase in average densities around the stars with higher
progenitor mass.

###  5.2 Significant fraction of future SNe outside molecular clouds

####  5.2.1 General statistics

Another key implication from Figures 2, 3 and Section 4.2 is that a
significant fraction of future core-collapse SNe will occur in low-density
inter-cloud regions, dominated by atomic gas. Only about 30-40% of the RSGs
are in pixels that contain H2, and the rest in H I-dominated pixels. A similar
fraction of low-density environments are also inferred for SNRs, which likely
originate from both massive and low-mass progenitors. The results are
consistent with MC22, who also showed that almost 50% Type II SNe (primarily
arising from red and yellow supergiant stars) were not coincident with
detectable molecular clouds.

We note however that the shocks from future SN explosions of these stars in
the intercloud region can still travel tens of pc and interact with nearby H2
gas clouds. This association can be seen somewhat by eye in Figure 1, where
many SNRs are not coincident with, but are near a molecular cloud. In terms of
feedback, the ability of SNe to disperse dense gas reduces with distance from
the molecular cloud (Hennebelle & Iffrig, 2014; Iffrig & Hennebelle, 2015),
although the edge interaction can still accelerate cosmic rays (Sano & Fukui,
2021) and/or compress gas potentially to star-forming densities (Cosentino et
al., 2022), but with a lower angular filling fraction than a fully embedded
spherical shock.

We therefore assess the H2 non-detection fraction around stars as a function
of a search radius (called ‘aperture radius’, following the nomenclature of
MC22) in Figure 8. Specifically, we plot the sample fraction containing no
detectable H2 pixels within circles of increasing radii from the central star.
To distinguish from the H2 distribution of stars uncorrelated with molecular
clouds, we also show in Figure 8 the aperture non-detection fraction for the
“Random” population introduced in Section 4.3.

Figure 8 highlights the importance of spatial resolution for characterizing
the physical association of SNe with molecular gas (a point we will return to
in Section 5.3). The sample fractions of non-detection are the highest for
each category at the ACA map pixel scale ($\sim$4 pc), and then slowly
decreases as we increase the aperture radius. At $\gtrsim$300 pc, every star
has some detectable H2 pixel. Figure 8 also reaffirms that the distribution of
non-detection fraction vs aperture radius is different from the case of stars
randomly aligned with respect to molecular clouds. Similarly, MC22 showed in
their Figure 2 and Table 2 that their H2 detection fraction around core-
collapse SNe increases from 60% at 150 pc resolution to 94% at 1 kpc
resolution. Coincidentally, this result bears similarity to the findings that
CO and star formation are strongly spatially correlated on large scales, but
becomes de-correlated on the scales of individual giant molecular clouds (e.g.
Schruba et al., 2010; Kruijssen et al., 2019).

We then compare the sample fractions of non-detection vs radii in Figure 8
with the typical cooling radius expected from a 1051 erg SN explosion in
atomic gas (n = 1 cm-3). The SNR cooling radius is an important resolution
criterion in ISM simulations to identify whether SN feedback is implemented
primarily in the form of thermal energy or momentum (Hopkins et al., 2013;
Martizzi et al., 2015; Kim & Ostriker, 2015b). By the time the SNR has
expanded to its cooling radius, which typically marks the end of the adiabatic
Sedov-Taylor stage, it has deposited $\gtrsim$50% of its final momentum into
the surrounding gas, and has begun quickly radiating away thermal energy.

Kim & Ostriker (2015b) expressed the cooling radius ($R_{c}$; or the shell
formation radius) fitted to numerical simulations of SNRs from 1051 erg SN
explosions as

$R_{c}=\begin{cases}(22\ \mathrm{pc})(Z^{2}+0.04^{2})^{-0.061}n_{1}^{-0.43}\\\
\qquad\qquad\qquad\qquad\qquad\quad\text{(homogenous ISM)}\\\\[8.5359pt] (30\
\mathrm{pc})(Z^{2}+0.04^{2})^{-0.061}n_{1}^{-0.46}\\\
\qquad\qquad\qquad\qquad\qquad\quad(\text{inhomogenous ISM})\\\ \end{cases}$
(3)

where $Z$ is the ISM metallicity in units of $Z_{\odot}$ and
$n_{1}=n/1\mathrm{cm}^{-3}$ is the average gas number density (which is
characteristic of the median densities inferred from Figure 6). The
metallicity dependence is taken from Kim et al. (2023). The cooling radius is
larger in an inhomogeneous ISM as the SNR forward shock escapes more quickly
through the larger volume-filling low-density channels. We assume a typical
atomic gas density of 1 cm-3 at Z = 0.6 Z${\odot}$, which gives $R_{c}=23.4$
pc (homogenous) and $R_{c}=31.9$ pc (inhomogeneous), shown as grey shaded
regions in Figure 8.

We find that 10-20% WRs, 21-34% SNRs, and 30-50% RSGs have no detectable H2
within the typical cooling radii shown in Figure 8. In other words, the
blastwaves from the explosions of these stars in the inter-cloud media,
assuming typical ISM densities of $\sim$1 cm-3, are unlikely to affect nearby
molecular gas before the onset of significant thermal energy losses in the
shock. They may of course still interact if the inter-cloud densities are much
smaller (e.g. $R_{c}\approx 63-170$ pc for the homogeneous case with
$n_{1}=0.1-0.01$ at $Z=0.6$ Z⊙).

Figure 8: Fraction of stars in the samples of WRs, RSGs and SNRs without any
detectable pixels of H2 (traced by CO 2-1) within an aperture of given search
radius. For comparison with stars uncorrelated with molecular gas, we also
show the non-detection fraction of a mock population drawn randomly in the ACA
area (the “Random” case from Figures 4 and 5 ) with thin solid and dashed
lines, with solid showing the median, and dashed showing the 5th and 95
percentile bounds. The shaded region denotes the cooling radii for a supernova
shock evolving in density = 1 cm-3 and metallicity = 0.6 $Z_{\odot}$ for a
uniform ISM (dark shaded region) and inhomogeneous ISM (light shaded region;
Eq. 3). The horizontal scales show the maximum projected distance that a star
would travel within its lifetime (vertical ticks) for velocities = 1, 10 and
100 km/s. Note that the vertical positioning of the scales are arbitrary.

####  5.2.2 Long lifetimes, runaways, cloud dispersal?

While the H2 non-detection fraction of RSGs and SNR progenitors could be
explained by them simply drifting away from their birth clouds over their long
lifetimes (upto $\sim$40 Myrs), it is intriguing why a substantial fraction of
WRs have no detectable H2 despite their much shorter single stellar ages of a
few Myrs (Eldridge et al., 2017b; Dorn-Wallenstein & Levesque, 2018b). Since
we do not have detailed 3D kinematic information for the WRs (see Neugent &
Massey, 2014, for a discussion of the complexities of radial velocity
measurements of WRs), we make use of the separation information in Figure 8
and some published observations to assess possible scenarios.

One possibility is that these WRs had escaped their birth clouds at high
velocities acquired from either dynamical interactions or supernovae (believed
to be the cause of runaway OB stars; Fujii & Portegies Zwart 2011). We assume
the separation between the WR and natal cloud is $d=v_{*}t_{\mathrm{wr}}$,
where $v_{*}$ is the relative velocity between the WR and the cloud, and
$t_{\mathrm{wr}}$ is the age of the WR progenitor. In a single stellar
population at LMC-like metallicity, WRs have an average age of
$t_{\mathrm{wr}}=4$ Myrs, while in a 100% binary population,
$t_{\mathrm{wr}}=10$ Myrs, with the majority being within 16 Myrs (Eldridge et
al., 2017b; Dorn-Wallenstein & Levesque, 2018b). Therefore, the 20% of WRs in
Figure 8 that have no detectable H2 within atleast 32 pc (the assumed cooling
radius in the homogenous case in Section 5.2.1) would need $v_{*}>8$ km/s if
they were single stars, and $v_{*}>3.2$ km/s if they were in an interacting
binary system777We use ‘$>$’ because the actual velocity needed to produce the
projected separation may be larger if the velocity vector is not in the plane
of sky; in addition; the WR may not have originated in the nearest H2-detected
pixel, but farther out.. The 5% WRs in Figure 8 that are even more remote
(with no detectable H2 within $>$100 pc) would require $v_{*}>25$ km/s
(single) and $v_{*}>10$ km/s (binary). These velocities are not unlike those
predicted by models of walkaway and runaway stars produced by dynamical
ejections and supernova (e.g., Oh & Kroupa, 2016; Renzo et al., 2019; Dorigo
Jones et al., 2020). Such runaways have been considered in simulations as a
way to boost outflows in low-mass galaxies (e.g. Andersson et al., 2020;
Steinwandel et al., 2020). A non-trivial number of runaway OB stars are known
to exist in the Galaxy and Magellanic Clouds, though the estimated fraction of
the total population varies in the literature (e.g. Blaauw, 1961; de Wit et
al., 2005; Gvaramadze et al., 2012; Lamb et al., 2016). Therefore it is not
implausible that the H2 non-detection WRs are a result of them being displaced
from their birth locations at some walkaway/runaway-star-like velocity.

It is also possible that the parent clouds of WRs were quickly destroyed by
pre-SN and prior-SN feedback. Molecular clouds have been estimated to disperse
on 1-5 Myr timescales after emergence of massive stars, based both on
observations of giant molecular clouds and H$\alpha$-emitting stellar
populations of nearby galaxies (e.g. Chevance et al., 2020) as well as
detailed high-resolution simulations of turbulent molecular clouds (e.g.
Grudić et al., 2022). Nearly 79% of the WRs with H2 non-detection (i.e. 42 out
of 53 WRs) are coincident with known OB associations according to Neugent &
Massey (2011), who cross-matched their WR sample with the OB associations of
Humphreys & Sandage (1980). The H2 detected WR population also has a similar
($\sim$80%) fraction of WRs coincident with OB associations (50 out of 63
WRs). So it is possible that the WR population with H2 non-detections are
still within their parent OB association, and since the typical stellar
velocity dispersions in OB associations is $\sim$4 km/s (Mel’nik & Dambis,
2017; Wright et al., 2022), these WRs may not have acquired the velocities
estimated in the previous paragraph to produce the observed cloud-star
separations in Figure 8, which leaves the cloud destruction scenario by pre-
SN/prior-SN feedback a more plausible explanation. The remaining 21% of WRs
with neither detectable H2 nor identified within any known OB associations may
be stronger candidates for walkaways/runaways.

Figure 9: Zoom-in region of the Wolf-Rayet/SNR pair from Figure 1,
demonstrating the benefit of CO observations with increasing spatial
resolution. (Left) Small patch of the ACA CO(2-1) + 21 cm map at 12″resolution
from Figure 1, more centered on the WR/SNR pair. Symbols denoting RSGs
(orange), SNRs (green) and WRs (purple) are the same as Figure 1. The black
box is the $\sim$(122 pc)2 region detailed in the smaller panels. (Top middle)
Same low-resolution (12″$\approx$48 pc) CO (2-1) map as the left panel, but
zoomed in. The WR and SNR are denoted. (Top right) High-resolution
(0.45″$\approx$1.8 pc) ALMA CO map of the same region (Section 3.2.3), showing
more detailed substructure of the molecular ISM. The WR is clearly located in
a vacant cavity in the ISM. (Bottom middle) Continuum-subtracted [S
II]/H$\alpha$ from narrowband images taken by Massey et al 2007. The
colorscale is split at a value of 0.4, commonly used to distinguish
photoionized (bluer) from shock-heated (redder) regions. We clearly see the
SNR region as a shock-heated nebula, with the WR embedded in a more
ionization-dominated region. (Bottom Right) RGB image showing broadband $U$
(blue), V (green) and narrowband H$\alpha$ (red) from Massey et al (2006) and
Massey et al (2007). The Wolf-Rayet star is clearly visible as a hot star, and
the SNR nebula is visible below it. The higher-resolution CO map reveals the
detailed substructure around the WR, particularly its evolution inside a
cavity presumably driven by powerful pre-SN feedback.

###  5.3 Cavities and substructures on scales $\ll$ 50 pc

Although we are finding younger massive progenitors closer to dense gas, a
caveat is the spatial resolution of 50 pc. While this is higher than
extragalactic SN studies (e.g., MC22), the 50 pc beam still only captures a
relatively large-scale mean gaseous environment, and misses substantial
structure and inhomogeneity on smaller scales.

As a case study, we show in Figure 9 a region of M33 (the black box from
Figure 1) containing a giant molecular cloud, a WR (J013335.73+30.3629.1 in
Neugent & Massey, 2011) and an SNR (L10-045 in White et al., 2019), where we
also have pc-scale CO (2-1) data from the ALMA archive (Section 3.2.3). The
low-resolution ACA image shows the WR evolving in detectable H2 and close to
the peak H2 emission of the cloud, but the higher resolution image shows the
WR inside a $\sim$10 pc radius cavity. The peak emission comes from a
neighboring cloud complex likely produced by ongoing cloud collisions (Sano &
Fukui, 2021). We also see a similar cavity in the H$\alpha$ image (bottom
right panel), where the WR star is surrounded by an H$\alpha$ shell of similar
extent as the H2. The H$\alpha$ shell also shows [S II]/H$\alpha$ $<0.4$,
indicating it is photoionization-dominated as opposed to shock-heated
(Matonick & Fesen, 1997). This proves that it is a real region of gas deficit,
most likely encased in a photoionized shell sweeping up surrounding CO-
emitting H2. Similar structures have been observed around WRs in the Galaxy
(e.g., Baug et al., 2019).

A $\sim$10 pc-radius cavity could have been produced by a previous SN
explosion or pre-SN feedback. In the first of two possibilities, a 1051 erg SN
explosion will sweep up ambient gas until the shock velocity is reduced to the
turbulent velocity dispersion ($\sigma$) of the ISM. Based on Eq 20 in
Martizzi et al. (2015), the maximum radius of the swept-up momentum-driven SNR
shell is,

$R_{\mathrm{SNR}}=(12.5\ \mathrm{pc})\ n_{2}^{-0.39}\sigma_{5}^{-0.33}$ (4)

for an ambient gas with density $n_{2}=n_{0}/10^{2}$ cm-3, and turbulent
velocity dispersion $\sigma_{5}=\sigma/5$ km s-1, where we normalized to 5
km/s which is the average value inside the inner 5 kpc of M33 (Koch et al.,
2019). Eq (4) shows that a SN explosion with typical energy output can easily
produce a 10 pc-radius cavity inside molecular gas with typically observed
densities and turbulence.

If a SN has not yet exploded in the region however, pre-SN feedback from the
WR and its preceding O-star phase could also have also carved out the present-
day cavity. The WR was not identified as part of any known OB association by
Neugent & Massey (2011), and no catalogued star-cluster is in the cavity
(Sarajedini & Mancone, 2007; Fan & de Grijs, 2014; Johnson et al., 2022), so
it is likely the WR star would be the primary power source888We also note that
no other spectroscopically-confirmed O/B star from Massey et al. (2016) is in
the vicinity, though completeness of the O-star catalog is still an issue..
Thermal pressure from ionized gas has been shown to dominate the expansion of
H II regions and dispersal of molecular gas (e.g. Lopez et al., 2014; Olivier
et al., 2021; Chevance et al., 2020), and this may be evident in the H$\alpha$
shell lined with the H2 cavity in Figure 9. The expansion of a D-type
ionization front around OB stars in a uniform ambient medium with density
$n_{0}$ can be expressed as in Hosokawa & Inutsuka (2006) Eq (36) as

$R_{\mathrm{II}}(t)=R_{\mathrm{st}}\left(1+\frac{7}{2\sqrt{3}}\frac{c_{s}t}{R_{\mathrm{st}}}\right)^{4/7}$
(5)

where $R_{\mathrm{st}}$ is the Stromgern sphere given by

$R_{\mathrm{st}}=\left(\frac{3Q_{0}}{4\pi\alpha_{B}n_{0}^{2}}\right)^{1/3}$
(6)

where $\alpha_{B}=3.023\times 10^{-13}$ cm3 s-1 is the case-B recombination
coefficient, $Q_{0}$ is the ionizing luminosity (s-1), $c_{s}\approx(11.6\
\mathrm{km/s})(T/10^{4}\mathrm{K})^{0.5}$is the sound speed inside the H II
region. Based on Eq (5), stars $>$30 M⊙with typical $Q_{0}=10^{49}-10^{50}$
s-1 (Sternberg et al., 2003) can easily produce ionized bubbles $>$10 pc
inside gas with densities $n_{0}=10^{2}$ cm-3 within $t\approx 0.4-1$ Myrs.
These numbers are also consistent with results from more detailed numerical
simulations (Kim et al., 2018, 2021, 2023). We can therefore conclude that
even if a SN had not exploded in this region, the expansion of an H II region
driven by ionizing radiation from the WR star or its preceding O-star phase
could have also easily created this cavity.

It is also possible that the cavity was produced by strong stellar winds,
given the typical mass-loss rates in O-stars being in the range of
$10^{-7}-10^{-5}$ M⊙yr-1 at velocities $>$103 km s-1, and around
$10^{-5}-10^{-4}$ M⊙yr-1 for WRs (Smith, 2014). According to Castor et al.
(1975) and Weaver et al. (1977), the mechanical luminosity from this wind can
drive an idealized adiabatic wind bubble of radius

$R_{b,\mathrm{ad}}=(8.45\ \mathrm{pc})\
\dot{M}_{w,-6}^{1/5}v_{w,3}^{2/5}n_{2}^{-1/5}t_{6}^{3/5}$ (7)

where $\dot{M}_{w,-6}=\dot{M}_{w}/10^{-6}$ M⊙yr-1 is the mass-loss rate,
$v_{w,3}=v_{w}/10^{3}$ km s-1 is the wind speed, and $t_{6}=t/10^{6}$ yr is
the age of the bubble. Again, the blast-wave extent in Eq (7) is similar to
the observed cavity size for fiducial values.

On the other hand, cooling at the interface of the shocked wind and swept up
shell can limit the expansion, as was investigated in detail by Lancaster et
al. (2021a, b), who derived the corresponding momentum-driven bubble radius

$R_{b,\mathrm{mom}}=(3.55\ \mathrm{pc})\
\dot{M}_{w,-6}^{1/4}v_{w,3}^{1/4}n_{2}^{-1/4}t_{6}^{1/2}\\\ $ (8)

In this case, one would require larger wind luminosity (e.g. increasing
$\dot{M}_{w,-6}>10$ and $v_{w,3}>4$) than the adiabatic case to produce our
observed cavity, assuming same cloud density and driving age.

The relevant takeaway, regardless of the specific feedback history or
mechanism, is that the WR is in a cavity that was only revealed at higher
resolution. Had the WR exploded in the future directly in the dense gas in the
absence of a cavity, the thermal energy of the blast wave would have been
efficiently radiated away on 6–12 pc scales (Eq 3), limiting the spatial
extent of its impact. In the case of explosion inside a cavity, where we
assume the gas density is now quite low ($\lesssim$10-2 cm-3), the cooling
radius is about 150 pc, which means the blast wave will easily retain its
total energy until it impacts the surrounding cloud material.

Such small-scale substructures could also be pervading other peaks of gas
density seen in Figure 1 where massive WRs and RSGs are coincident, but
smeared out by our 12″resolution. This further highlights the need for mapping
the cold, dense ISM in these galaxies at pc-scale resolution in order to
unveil the true distribution of high and low-density channels of ambient gas
around massive stars through which their future SN blast-waves will eventually
evolve.

###  5.4 A new observational anchor for supernova feedback models

Figure 10: Example comparison of the H I + H2 gas volume density ($n_{tot}$)
of RSGs (_left_ panels) and WRs (_right_) at 50 pc resolution with the ambient
SN densities from a high-resolution hydrodynamical simulation of a dwarf
starburst galaxy (Lahén et al., 2020). Simulated densities were obtained for
all SNe during the quiescent, pre-starburst phase between 70-140 Myrs (See Fig
5 in Lahén et al., 2020), and divided into 8-30 M⊙and $>30$M⊙mass ranges to
better compare with RSGs and WRs respectively. The observations, i.e. WRs and
RSGs, are shown as colored histograms. The ‘local’ theoretical densities (i.e.
measured in the immediate region around the SN) are shown in the top panel as
thin dashed histograms, while the ‘average’ theoretical densities (measured
from 2D-projected gas maps) are shown in thick histograms. Bottom panel zooms
into the observed and the average density histograms, which are more directly
comparable. The plot shows the feasibility of comparing ambient densities of
simulations and observations as an independent way to assess the ISM and
feedback physics (see Section 5.4 for a detailed discussion).

Stellar feedback models are widely recognized as a major source of systematic
uncertainty in the properties of galaxies predicted by numerical simulations
(Naab & Ostriker, 2017; Vogelsberger et al., 2020). The largest simulations of
galaxies with mass-resolutions of $\sim$105-107 M⊙(e.g., Vogelsberger et al.,
2014; Schaye et al., 2015; Khandai et al., 2015; Feng et al., 2016; Davé et
al., 2016) use subgrid models of SN feedback that involve significant
abstraction, such as imparting hydrodynamically-decoupled ‘kicks’ to
neighboring gas particles of SNe (e.g., Mihos & Hernquist, 1994; Springel &
Hernquist, 2003; Okamoto et al., 2010; Vogelsberger et al., 2013) and
artificial suppression of cooling in SN-heated gas (e.g., Thacker & Couchman,
2000; Sommer-Larsen et al., 2003; Stinson et al., 2006). Zoom-in simulations
of individual galaxies reach higher baryonic resolution ($\lesssim$103 M⊙, or
$\lesssim$10 pc) at the cost of simulated volume, and can more explicitly
account for the thermal energy and momentum budget of SNRs, as well as massive
stellar winds, UV radiation and H II regions (Hopkins et al., 2011, 2012,
2013, 2018). Simulations of stratified ISM disks (e.g., Hill et al., 2012;
Walch et al., 2015; Kim & Ostriker, 2015a; Martizzi et al., 2016; Kim &
Ostriker, 2017; Kim et al., 2017, 2020; Rathjen et al., 2021) and low-mass
dwarf galaxies (e.g., Hu et al., 2017; Smith et al., 2018; Hu, 2019; Emerick
et al., 2019; Steinwandel et al., 2020; Lahén et al., 2020; Smith et al.,
2021; Gutcke et al., 2021; Andersson et al., 2022) achieve the highest
resolutions ($\sim$ few M⊙, $<$pc), which can further spatially-resolve
processes like the formation of individual stars, the formation and
interaction of blast waves, turbulent multi-phase ISM and outflows. These
simulations have particularly highlighted the importance of mass-resolution
for resolving SN thermal feedback, and pre-SN feedback and SN clustering in
pre-processing the SN environments. Overall, simulations across the resolution
ladder reinforce the pivotal role of stellar feedback, and the need for
capturing all the relevant physics in simulations.

Here we propose that high angular resolution observations of individual
massive stars and their ISM environments can be a novel constraint on stellar
feedback models for future simulations, particularly as a guide for the
seeding strategy of SNe.

Such comparisons can be done as shown in Figure 10, where as a preliminary
demonstration we compare the observed densities of our WRs and RSGs from
Figure 6 (colored histograms) to the ambient SN densities from a high-
resolution ($\sim$4 M⊙) simulation of a dwarf starburst galaxy from Lahén et
al. (2020) mentioned in Section 3.3. Since our observed densities are derived
from the projected ISM densities divided by a scale height (Section 5.1 and
Appendix A), they are best compared with similarly obtained average SN
densities from the simulation, as described in Section 3.3 (solid black
histograms). For the sake of discussion, we show both the local (thin-dashed)
and averaged (thick-solid) densities around simulated SNe in Figure 10. We
notice firstly that the local densities are much lower than the projected
densities. This is a common prediction of high-resolution galaxy simulations
where majority of stars explode in low densities carved by photoionization and
SN clustering, and due to stellar motion away from the birth clouds (e.g. Hu
et al., 2017; Lahén et al., 2020; Rathjen et al., 2021; Andersson et al.,
2022). Unfortunately, the immediate local densities are not measurable in
observations where we always see the ISM and stars in projection. In addition,
our ISM tracers of H I and H2 are mainly sensitive to gas densities $>0.1$
cm-3 (since our main goal in this paper was to measure the quantity of dense
gas around SNe). In future we can include recombination line and diffuse X-ray
maps of galaxies to estimate the proportions of the ionized and hot ISM gas
along the line of sight of stars.

Nevertheless, the histogram of the average densities are still the most direct
observations possible of where stars are exploding in the ISM and, as shown in
Figure 10, provide some interesting lines of comparison with simulations. The
range of average densities of the simulated SNe and the observed stars are
somewhat similar, falling in the range of 0.1-100 cm-3. This rough agreement
with observations is encouraging given that this is one of the highest-
resolution simulations of an entire galaxy, such that the feedback-driven ISM
and star-formation is effectively generated from first-principles with little
to no reliance on subgrid models.

However, there are also subtle differences between the observed and simulated
densities. The peak of the RSG density is a factor of 3 smaller than the 8-30
M⊙stars, and the bimodal WR density (with one peak at $\sim$ 1 cm-3 associated
with primarily atomic densities, and other at $\sim$ 10 cm-3 associated with
primarily molecular gas) contrasts with the unimodal $>$30 M⊙distribution in
the simulation. Interestingly, the bimodality is seen in the local ISM
densities in Figure 10, but more prominently for the 8-30 M⊙stars instead of
the $>$30 M⊙stars.

Exactly why these densities don’t agree is unclear presently but can be
understood in more detailed future investigations. On one hand, ambient
densities of SNe can depend significantly on assumptions in the pre-SN
feedback and clustering models in simulations (e.g. Gentry et al., 2017; Naab
& Ostriker, 2017; Smith et al., 2018; Rathjen et al., 2021), runaways (e.g.
Andersson et al., 2020) and choice of numerical models (e.g. Hu et al., 2022).
On the other hand, not all aspects of our specific simulation-observation
comparison in Figure 10 are optimal. The simulated dwarf galaxy system is
about a factor of 100 lower in stellar mass than M33 (Corbelli, 2003), and is
undergoing some interaction (though we exclude the starburst phase during the
second encounter), so the overall ISM phase distribution and porosity could be
quite different from M33. Inclusion of the warm and hot ISM tracers in the
observations as mentioned earlier can provide a fairer comparison of the low-
density tail with simulations. And finally, as mentioned in Fig 9, higher-
resolution observations will better clarify whether some of the WRs near high
density peaks are located in low-density cavities even in projection.

We conclude this section encouraging further comparisons of this kind to
validate the ISM and feedback physics in simulations. It will likely be most
beneficial to compare simulations and observations of similar galaxies, so
from the observations side – we will expand our analysis to a wider range of
galaxies, from the rest of M33 and other dwarf galaxies in our Local Group, to
the nearby PHANGS targets that fully cover the star-forming main-sequence
(Leroy et al., 2021a).

###  5.5 Comparison with SN and SNR studies

While our paper was focused specifically on the ISM properties and
implications for feedback, we note the parallels of our work with past studies
of star-formation rates around SNe, which show that SE-SNe (Type Ib/c, Type
IIb) are more likely to be near bright H II regions than Type II or Type Ia
SNe (e.g. Kelly et al., 2008; Anderson & James, 2008; Anderson et al., 2012;
Galbany et al., 2016a, b, 2018; Kuncarayakti et al., 2018; Cronin et al.,
2021). The trend has been attributed to an increasing progenitor mass range
going from Type II to SE progenitors, hence their association with younger
star-forming populations. Similar results were also found with molecular gas,
which correlates with recent star-formation rate, in kpc-scale maps of
galaxies (Galbany et al., 2017) as well as the GMC scale studies of MC22. Our
work also shows that RSGs, which are expected to be progenitors of Type II
SNe, are less closely associated with star-forming H2 gas than WRs. Similar to
our Figure 5, Figures 10 and 11 in MC22 show that Type Ib/c’s can be
reproduced by distributions drawn from the H2 density peaks, while Type II SNe
closely follow the stellar distribution traced by near-IR maps. A large
fraction (53%) of the MC22 sample also has CO non-detections similar to our
findings in Figures 2 and 3. Zoom-in images of the SNe in Figures 8 in MC22
display the wide variety of environments apparent even at 150 pc, with SNe
evolving in bright dense gas peaks, some in regions of lower density, clouds
and general gas deficit, and some in outright voids in the ISM.

We also note the parallels of our work with Kangas et al. (2017), who compared
the locations of massive stars (including RSGs and WRs) on H$\alpha$ images of
LMC and M33, and compared the correlations with that observed for SNe vs
H$\alpha$ observed in nearby galaxies. They find that WRs and SE-SNe have
similar and stronger spatial correlation with H$\alpha$ emission (and thus
recent star-formation) compared to RSGs and Type II SNe.

Overall, it is encouraging that the results with massive star progenitors and
SNe converge, validating the use of these massive stars as a local analog of
SN sites, where we have substantial detail compared to distant galaxies,
allowing us to deconstruct their environmental details. We believe higher
resolution studies of both our Local Group stars, and PHANGS SNe promise to
provide even more subtle details involving high mass/short-lifetime
progenitors.

A final comparison we make is with the statistics of known SNRs interacting
with molecular gas. Many well-studied interacting SNRs exist in our Galaxy
such as IC443, W44, W28 and W51C, visible in broadened spectral line and
continuum emission across the EM spectrum (e.g. see Slane et al., 2015, and
references therein). These provide the most direct evidence that SNe interact
with molecular clouds. The exact fraction of interacting SNRs in our Galaxy
however is highly uncertain due to the unknown classifications of most SNRs,
and unknown completeness of the highly heterogeneous collection of Galactic
SNRs. To our knowledge, the most systematic study of interacting SNRs was done
by Kilpatrick et al. (2016), who looked for broadened CO emission towards a
subset of SNRs. They found 17/46 SNRs ($\sim$37%), or 7/15 ($\sim$47%) if
restricted to SNRs with known CC classification, showing potential evidence of
interaction. This is roughly consistent with our 39% SNRs in a detected H2
pixel, though it is much lower if we consider SNRs within H2 pixels within a
cooling radius (66-80%). While this may suggest the interacting SNRs in our
Galaxy may be an undercount, we are only going by positional coincidence. It
is possible the SNRs may be in the foreground or background of the cloud.
These SNRs would be interesting targets for follow-up with more sensitive
high-resolution ALMA observations, where signatures of interaction could be
inferred from moment data and their correlation with optical line/X-ray
luminosities. We leave this for a future work.

##  6 Conclusion

With ambient density being a decisive parameter for SN thermal and momentum
feedback, we have embarked on a much-needed detailed observational campaign to
map the gas distribution around SNe by employing a novel and unique strategy –
targeting evolved massive stars that are soon to explode near their present
location, and their cold ISM environments. This not only yields an order-of-
magnitude larger sample of ‘explosion sites’ of stars than possible with a SN
survey in the nearby universe, but also exploits the proximity of Local Group
galaxies (in our case M33), which has not only cloud-scale molecular gas maps,
but also atomic gas maps, giving estimates of the H I, H2 and total cold ISM
densities in the 50 pc environments of massive stars.

Specifically, we use published catalogs of evolved massive stars such as red
supergiants (RSGs) and Wolf-Rayet (WRs) stars that have been connected to SN
explosions. These stars will have lifetimes of roughly 0.1–1 Myr, and thus
explode within a few pcs of their current location. While RSGs plausibly trace
locations of future Type II SNe, WRs to first order trace where the most
massive ($>$30 M⊙) and youngest stellar populations are located and may
explode in the future as stripped-envelope SNe and gamma-ray bursts. We also
include supernova remnants (SNRs), verified by high-quality multi-wavelength
data, as a tracer of where stars have already exploded.

To trace the cold ISM, we use new ALMA ACA maps of CO(2-1)-traced H2 and 21 cm
atomic H I in the inner projected $\sim 3\times 5$ kpc2 of M33 (E. Koch et
al., in prep), complete down to $\sim$1–2M⊙ pc-2, and at 12″ resolution
($\sim$50 pc). We then assess the relative cumulative distribution function of
their densities, and assess the physical and statistical nature of these
distributions using multi-wavelength star-formation maps.

Our investigation revealed the following about where stars explode in the ISM:

  1. 1.

We find an increasing tendency of higher mass progenitors to explode in denser
gas, as evidenced by 55% of WRs exploding in pixels with detectable H2 with a
median density of 15.5 M⊙ pc-2, compared to 27-32% for RSGs and SNRs at median
densities of $\approx 7$ M⊙ pc-2 (Figures 2,3, 7, Table 2). This trend is also
observed within RSGs: more luminous RSGs with log$(L/L_{\odot})>5$,
corresponding to masses $>$15 M⊙, have a higher fraction of stars (44%)
exploding in denser gas (with a median surface density of 11.6 M⊙ pc-2) than
the bulk RSG population.

  2. 2.

A significant fraction of massive stars ($\sim$45% for WRs, 56-77% for the
RSGs, depending on mass) are located in pixels with no detectable H2. This is
similar to the result of Mayker Chen et al. (2022) who found that 50% of Type
II SNe explode outside GMCs at 150-pc resolution. Of our H2 non-detected
population, 83-94% of WRs, RSGs and SNRs are in pixels with detectable H I,
implying they will likely explode in the lower density, atomic-gas dominated
intercloud media. Analysis of the star-molecular cloud separation distances in
Section 5.2.2 indicate that pre-SN and prior-SN clearing of H2 clouds is the
likely cause (particularly for the WRs), though stellar motion away from
clouds (particularly for lower mass progenitors) cannot be ruled out.

  3. 3.

For the first time, we have measured _both_ molecular and atomic gas densities
around SN progenitors at 50 pc resolution (Figures 2+3). The H I or H I + H2
densities of WRs, RSGs and SNRs are not as different as H2 alone, due to the
uniform distribution of atomic gas. The general trends though are similar –
WRs explode at the highest median H I surface densities of 9.7 M⊙ pc-2
followed by SNRs (8.8 M⊙ pc-2). and then RSGs (6.8 M⊙ pc-2).

  4. 4.

The density distribution around WRs and RSGs $>$15 M⊙ are statistically
distinct from the bulk gas distribution of M33, and these densities can
reproduced by stellar populations drawn from H$\alpha$ and far-ultraviolet
maps respectively (Figure 4 and 5). Ambient densities of the lower mass RSGs
and SNRs are marginally different from the bulk densities, and appear to
follow the near-infrared light in M33 tracing bulk stellar mass.

  5. 5.

We discovered a 20 pc-scale cavity around a WR from an archival $\sim$pc-
resolution ALMA field that was not apparent at 50 pc resolution (Figure 9).
The cavity is also apparent as an H$\alpha$-shell around the WR. A cavity
could have been produced by a previous SN explosion, ionizing radiation or
stellar winds from the WR/O-star. This revelation only at $\sim$pc-scale
further highlights the importance of pre-SN feedback/clustered SNe in
preprocessing SN environments, and the value of high angular-resolution
observations for revealing such ISM substructure.

  6. 6.

Our measured densities can be directly compared to the average ambient
densities of SNe in hydrodynamical simulations of galaxies (an example
comparison is done with the high-resolution dwarf galaxy system by Lahén et
al. (2020) in Figure 10). The comparison reveals that the simulated and
observed densities span a similar range, but with subtle differences in the
shapes of the density distributions. We will continue expanding such
observations to galaxies spanning a range of stellar masses and star-formation
rates in the local universe, and encourage the use of this data to validate
the feedback models in simulations for _any_ range of baryonic resolutions.

## Acknowledgements

We acknowledge useful discussions with Drs. Thorsten Naab, Eve Ostriker and
Todd Thompson on stellar feedback models. LC grateful for support from the
National Science Foundation in the form of grants AST-2107070 and AST-2205628.
V. V. acknowledges support from the scholarship ANID-FULBRIGHT BIO 2016 -
56160020, funding from NRAO Student Observing Support (SOS) - SOSPADA-015, and
partial support from NSF-AST1615960. The National Radio Astronomy Observatory
is a facility of the National Science Foundation operated under cooperative
agreement by Associated Universities, Inc. This paper makes use of the
following ALMA data: 2017.1.00901.S; 2019.1.01182.S. ALMA is a partnership of
ESO (representing its member states), NSF (USA) and NINS (Japan), together
with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in
cooperation with the Republic of Chile. The Joint ALMA Observatory is operated
by ESO, AUI/NRAO and NAOJ Support for K.F.N was provided by NASA through the
NASA Hubble Fellowship grant (HST-HF2-51516) awarded by the Space Telescope
Science Institute, which is operated by the Association of Universities for
Research in Astronomy, Inc., for NASA, under contract NAS5-26555. NL
acknowledges the computational resources granted by CSC – IT Center for
Science Ltd. in Finland and The Max Planck Computing and Data Facility (MPCDF)
in Garching, Germany. SS and NMP acknowledge support provided by the
University of Wisconsin - Madison Office of the Vice Chancellor for Research
and Graduate Education with funding from the Wisconsin Alumni Research
Foundation, and the NSF Award AST-2205630. This research made use of
Photutils, an Astropy package for detection and photometry of astronomical
sources (Bradley et al. 2023)

##  Appendix A Calculating the spatially-varying volume density

Figure 11: Variation in the H I, H2, and stellar mass surface densities and H
I velocity dispersion (_top_), scale height (_middle_) and total H I + H2
volume density (_bottom_) versus deprojected radius of M33. Dashed lines
represent the approximate extents of the ACA region in Figure 1 in radius
units. The dotted line in the top panel represents an exponential disk model,
assumed in the text. Shaded regions represent the 16th-84th percentile range,
while solid lines represent the median. See Appendix A for details.

Given measurements of $\Sigma_{HI}$ and $\Sigma_{H_{2}}$, the average total (H
I + H2) density ($n_{\mathrm{tot}}$, in units of atoms per cm-3) can be
obtained as

$n_{\mathrm{tot}}=\frac{1}{2\mu
m_{p}}\left(\frac{\Sigma_{HI}}{h_{HI}}+2\frac{\Sigma_{H_{2}}}{h_{H_{2}}}\right),$
(A1)

where the factor of 2 with $\Sigma_{H_{2}}$ accounts for 2 atoms of H per
molecule, $\mu=1.4$, and $m_{p}=1.6\times 10^{-24}$ g. The underlying
assumption is that both gas phases are uniformly distributed along each
column, characterized by scale heights $h_{H_{I}}$ and $h_{H_{2}}$. While
obviously the real feedback-driven ISM is inhomogeneous and turbulent, this
simple parameterization of the line-of-sight gas gives an easy-to-interpret
average density per pixel as a function of the observed projected surface
density in order to enable first-order comparisons with simulations.

The scale height of H I ($h_{HI}$) can be estimated assuming the disk of M33
is in a state of dynamical equilibrium, i.e., the gravitational weight of the
inner disk due to stellar and gaseous components is balanced by thermal and
turbulent pressure of the diffuse gas driven by star-formation and stellar
feedback (Ostriker et al., 2010; Ostriker & Shetty, 2011b; Shetty & Ostriker,
2012). Numerical simulations have verified that within an orbital time, star-
forming disks reach such quasi-steady equilibrium properties due to self-
regulation of star formation by feedback (e.g., Kim et al., 2011, 2013;
Ostriker & Kim, 2022), and the resulting star-formation rates and related
observables are broadly consistent with observations (e.g., Gallagher et al.,
2018; Sun et al., 2020, 2023). This formalism therefore provides a unique way
to obtain the _spatially-varying_ scale height and line-of-sight volume
density of the atomic ISM for a nearly face-on galaxy like M33, as a function
of locally measured quantities such as the gas and stellar mass surface
density, kinematics, and star-formation rates.

For our purpose, we closely follow the arguments in Utomo et al. (2019) to
obtain the scale height, with a few minor modifications. We assume that the
vertical pressure in the H I gas in M33 is related to its velocity dispersion
and density by

$P_{HI}=\rho_{HI}\ \sigma_{HI}^{2}.$ (A2)

In a disk under dynamical equilibrium, the midplane pressure is given by,

$P_{\mathrm{tot}}=f_{d}\frac{\pi
G\Sigma^{2}}{4}\left[(2-f_{d})+\Biggl{\\{}(2-f_{d})^{2}+\frac{32\sigma^{2}\rho_{\mathrm{sd}}}{\pi^{2}G\Sigma^{2}}\Biggr{\\}}^{1/2}\right]$
(A3)

in the form given in Eq (35) of Kim et al. (2011). Here
$\Sigma=\Sigma_{HI}+\Sigma_{H_{2}}$ is the total gas surface density, and
$f_{d}$ is the diffuse gas fraction, i.e., $f_{d}=\Sigma_{HI}/\Sigma$. The gas
velocity dispersion $\sigma$ is assumed to be equal to $\sigma_{HI}$. The
total stellar + dark-matter density is $\rho_{sd}$, but since we are only
interested in the inner few kpc of M33, we can approximate
$\rho_{sd}\approx\rho_{*}$, i.e., only the stellar density (Corbelli, 2003).

The stellar volume density, $\rho_{*}$, can be computed from a near-IR map of
M33. Similar to H I and H2, we assume $\rho_{*}\approx\Sigma_{*}/2h_{*}$,
where $h_{*}$ is the stellar scale height obtained by assuming a flattening
ratio between the stellar scale lengths and heights, $l_{*}/h_{*}=7.3$ (Kregel
et al., 2002), where we assume the scale length $l_{*}\approx 2$ kpc based on
photometric observations in M33 (Barker et al., 2007), and also in Figure 11.
This gives $h_{*}=0.27$ kpc. We compute $\Sigma_{*}$ from the _WISE_ W1 (3.4
$\mu$m) map following the conversion given in Leroy et al. (2019)

$\frac{\Sigma_{*}}{\mathrm{M}_{\odot}\mathrm{pc}^{-2}}\approx
330\left(\frac{M/L}{0.5}\right)\left(\frac{I_{3.4\mu m}}{\mathrm{MJy}\
\mathrm{sr}^{-1}}\right),$ (A4)

where $M/L$ is the mass-to-light ratio in the 3.4 $\mu$m range. The $M/L$ is
expected to correlate with the underlying age of the stellar population, and
thus vary across the galaxy (Stanway & Eldridge, 2018). We account for this
based on the empirical relation with the specific star-formation rate given in
Eq (24) in Leroy et al. (2019) as

$M/L=\begin{cases}0.5&\text{if}\ Q<a\\\ 0.5+b(Q-a)&\text{if}\ a<Q<c\\\
0.2&\text{if}\ Q>c\end{cases}$ (A5)

Here $Q=\mathrm{log_{10}}(L_{\mathrm{FUV}}/L_{W1})$ is the ratio between the
FUV and W1 luminosities, roughy equivalent to the specific star-formation
rate, $a=-2.5$, $b=-0.167$ and $c=-0.7$. Both FUV and W1 maps were convolved
to the resolution and pixel scale of the ACA map before the conversions were
applied.

With the above quantities in place and assuming
$\rho_{HI}=\Sigma_{HI}/2h_{HI}$, we obtain the scale height $h_{HI}$ by
equating Eq (A2) and (A3), giving

$h_{HI}=\frac{1}{2}\frac{\Sigma_{HI}\sigma_{HI}^{2}}{P_{\mathrm{tot}}}$ (A6)

providing a spatially-varying H I scale height as a function of the local
measured atomic, molecular and stellar mass surface densities, and the atomic
line-width maps. We have verified from simulations of Ostriker & Kim (2022)
that these scale height equation provide a reasonable estimate (within 10-20%)
of the average density at 50 pc resolution. Further accuracy can be obtained
in future work by including the mass of hot gas from X-ray observations.

The molecular scale height, $h_{H_{2}}$, however is less trivially obtained as
the gas is primarily locked in gravitationally-bound clouds that are not
necessarily in pressure equilibrium with the disk, and equations of
hydrostatic balance (similar to Eq (A3)) have been shown to over-estimate the
actual scale height of the clumpy molecular disk (Jeffreson et al., 2022). For
simplicity, we assume $h_{H_{2}}=50$ pc, which is consistent with that
expected for the Milky Way (Jeffreson et al., 2022). This may be slightly
smaller than the expected value for M33, which has a much smaller stellar mass
(and thus smaller gravitational potential) than the Milky Way. We found that
$n_{\mathrm{tot}}$ obtained by assuming $h_{H_{2}}=50$ pc is on average only
about 40% larger than the limiting case of $h_{H_{2}}=h_{HI}$. This is much
smaller than the orders of magnitude variation in densities expected across
the disk in the ACA region, so we believe $h_{H_{2}}=50$ pc is a reasonable
value for M33.

Figure 12: The H2 surface density CDFs similar to Figure 5 (dashed line
represents median, shaded region represents 95% interval), except showing the
different 24 $\mu$m extinction-corrections for FUV (_left_) and H$\alpha$
(_right_) as discussed in Section B. “Only FUV/H$\alpha$” refers to the maps
without extinction correction, “Const” C${}^{FUV,H\alpha}_{24}$ refers to a
constant value of the extinction-correction coefficient, and “Variable”
C${}^{FUV,H\alpha}_{24}$ refers to the coefficient scaled by the WISE W1-band
maps according to Belfiore et al (2023) to account for dust-heating by old
stellar populations. For simplicity, since the focus is on differences between
the prescriptions, we only show the high-mass RSGs on the FUV panel and WRs on
the H$\alpha$ panel, since these were most consistent with the respective
wavelengths in Figure 5.

Figures 11 shows the variation in the relevant quantities in the above
equations as a function of deprojected radius of M33 (note that the
calculations in the paper are pixel-by-pixel, not radial; we just show these
quantities vs radius to give a succinct idea of their variation). Deprojection
is done assuming the center of M33 is at RA(2000) = 01h:33m:50.89657s,
Dec(2000)= +30∘:39′:36.630403′′, an inclination angle $55^{\circ}$, and
position angle of $200^{\circ}$ (Koch et al., 2018). The profiles are plotted
in radial bins of $\Delta r=100$ pc. Within each bin, for the stellar and gas
mass densities, we estimate a median, 16th, and 84th percentile density,
multiply by the area covered by the _detected_ pixels to get the total mass,
and divide this mass by the total area inside the $\Delta r$ to get the radial
surface density. For the volume density and scale height, we show their median
and 16th-84th percentile range of their pixel values within each $\Delta r$.

In particular, we highlight the profiles of the scale height and total volume
density in Figure 11. The H I scale height as shown in the middle panel
increases with distance from the center because $\Sigma_{*}$, which provides
the dominant contribution to the gravitational potential of the disk within 5
kpc, drops with radius as $\sim$$\mathrm{e}^{-r/l_{*}}$, while the diffuse gas
properties that provide outward pressure (i.e., $\Sigma_{HI}$, $\sigma_{HI}$)
remain relatively flat. The scatter in the scale height is primarily due to
region-to-region variation in $\Sigma_{HI}$, $\Sigma_{H_{2}}$, $\Sigma_{*}$,
and $\sigma_{HI}$, which we are able to capture with this above approach (as
opposed to assuming a constant value throughout the disk). The final panel in
Figure 11 shows the total resulting volume density (assuming h${}_{H_{2}}$=50
pc) vs deprojected radius from Eq (A1). The scatter in density folds in the
scatter from the above quantities, but overall the density decreases with
distance from the center. The density within 1 kpc of the center has
particularly large scatter and is higher than at outer radii due to more
significant contribution from molecular gas (whose surface density also drops
off exponentially with radius as seen in the top panel). Both the scale height
and the density values are roughly consistent with values in the Milky Way
(Dickey & Lockman, 1990) and other edge-on galaxies (e.g., Patra, 2019;
Randriamampandry et al., 2021; Zheng et al., 2022), as well as simulations of
the multiphase ISM in realistic galaxy environments (e.g., Hill et al., 2012;
Hopkins et al., 2012; Kim & Ostriker, 2017).

##  Appendix B Extinction-corrected FUV and H$\alpha$ maps

We extinction-corrected the FUV and H$\alpha$ maps using MIPS 24 $\mu$m, as
prescribed in Belfiore et al. (2023) (hereafter B23). We first convolve and
regrid all maps to 12″resolution and 1″ pixel scale as the ACA maps, and
convert them to units of luminosity. We then convert them to units of hybrid
star-formation rate as given in Eq (4) as

$SFR_{\mathrm{FUV+24}}=C_{\mathrm{FUV}}L_{\mathrm{FUV}}+C^{\mathrm{FUV}}_{24}L_{24}$
(B1)

where SFRFUV+24 is the 24$\mu$m corrected FUV-based star-formation rate in
units of M⊙yr-1, LFUV is the FUV luminosity in erg s-1 and L24 is the 24$\mu$m
MIPS luminosity in erg s-1, $C_{\mathrm{FUV}}=10^{-43.42}$ (erg s-1)(M⊙
yr-1)-1 given in Table 2 of B23 and derived in Leroy et al. (2019), and
$C^{\mathrm{FUV}}_{24}$ is the 24$\mu$m correction factor for FUV. Similarly
for H$\alpha$ we follow Eq (5) in B23 as

$SFR_{\mathrm{H\alpha+24}}=C_{\mathrm{H\alpha}}L_{\mathrm{H\alpha}}+C^{\mathrm{H\alpha}}_{24}L_{24}$
(B2)

where the terms carry the same meaning as above, except for H$\alpha$ maps. We
use $C_{\mathrm{H\alpha}}=10^{-41.26}$ (erg s-1)(M⊙ yr-1)-1 as given in Table
2 of B23, derived in Calzetti et al. (2007).

The correction factors C${}^{FUV}_{24}$ and C${}^{\mathrm{H\alpha}}_{24}$ were
adjusted by B23 to more accurately account for dust heating by lower-mass
stellar populations. We use the calibration form of Eq 7 in B23

$\mathrm{log}\ C^{i}_{24}=\begin{cases}\mathrm{log}\ C_{\mathrm{max}}\\\
\quad+\ a_{1}(\mathrm{log}\ Q-\mathrm{log}\ Q_{\mathrm{max}})&\text{if}\
Q<Q_{\mathrm{max}}\\\ \mathrm{log}\ C_{\mathrm{max}}&\text{if}\
Q>Q_{\mathrm{max}}\\\ \end{cases}$ (B3)

where the subscript $i$ refers to H$\alpha$ or FUV, and $Q$ refers to the
quantities log$(L_{FUV}/L_{W1})$ and log$(L_{H\alpha}/L_{W1})$, as already
defined in Eq A5. Based on Table 3 in B23, for H$\alpha$, we adopt
$a_{1}=0.45$, log $C_{\mathrm{max}}=-42.88$, and log $Q_{\mathrm{max}}=-1.52$.
Similarly, for FUV, we adopt $a_{1}=0.23$, log $C_{\mathrm{max}}=-42.73$, and
log $Q_{\mathrm{max}}=0.6$. We note that B23 estimated these scaling relations
for log $C^{i}_{24}$ in the WISE W4 (22 $\mu$m) band, which was different from
measurements in MIPS 24 $\mu$m by 0.08 dex. As a result, we add 0.08 to log
$C^{i}_{24}$ in Eq B3 since we are working with MIPS 24$\mu$m maps.

Figure 13: The H2 surface density CDFs for mock population of stars drawn from
near-IR and FUV maps meant to replicate SNRs with some combination of core-
collapse and Type Ia SNe. This combination is represented by the ‘Ia
fraction’, the fraction of the population drawn from near-IR (rest drawn from
FUV). See Section C for details.

Figure 12 shows that CDFs exhibit small variations, mostly within error-bars
of each other, for the different prescriptions discussed above. Including the
24$\mu$m correction, the median fraction of the H$\alpha$-drawn star sample
occurring in detectable H2 is higher by 8% and the fraction of FUV-drawn stars
occurring in H2 is higher by 7%, compared to the maps without 24 $\mu$m
correction. The higher fraction is due to the infra-red capturing more
obscured regions in M33, which are expected to be associated with denser gas
and younger, more obscured stellar populations. Interestingly, the WR
population in Figure 12 is more consistent with the H$\alpha$-population
without extinction correction. This indicates that the spatial distribution of
WRs is similar to the unobscured, ionizing OB stars, and reinforces that we
are not missing a significant fraction ($\lesssim$7%) of WRs in highly
obscured gas clouds. For completeness, we also show the CDFs for FUV and
H$\alpha$ maps with a constant correction factor, as has been traditionally
used in the literature. For this purpose, we used the constant values given in
Table 2 of B23, specifically C${}^{H\alpha}_{24}=10^{-42.78}$ (Calzetti et
al., 2007) and C${}^{FUV}_{24}=10^{-42.83}$ (Hao et al., 2011; Leroy et al.,
2019). The difference in the corresponding CDFs lies in between the no-
extinction and variable coefficient cases as seen in Figure 12.

As explained in Section 4.3, these differences are not consequential to the
main science result, that the densities of the observed RSG and WR population
are not randomly drawn, and are roughly consistent with their progenitor
population.

##  Appendix C The contribution of Type Ia SNe in the SNR densities

We assess how the observed CDF of H2 densities of SNRs in the last row of
Figure 5 compares with model CDFs of H2 densities drawn from a combination of
SF-tracer maps that represents some proportion of SNe Ia (with low-mass
progenitors) and core-collapse SNe (with high mass progenitors) in the sample.
We know from Koplitz et al. (2023) that the Type Ia fraction is potentially
40% of the SNRs in the ACA region based on star-formation histories. We repeat
the experiment in Section 4.3 and Figure 5 for SNRs, except now assuming that
some fraction of each randomly-drawn sample of SNRs comes from the near-IR map
(representing old progenitors of Type Ia SNe) and the rest from the FUV map
(representing the massive progenitors of core-collapse SNe).

Figure 14: Same as Figure 2, but using the lower resolution C-configuration H
I data.

The results are shown in Figure 13. While all three models are consistent with
each other within the shaded error regions, we do find the median model CDF of
H2 densities, where we draw 40% of the sample from the near-IR to represent
Type Ia SNe, is slightly more consistent with the observed CDF (purple) than
the two extreme scenarios, i.e. all SNRs are Type Ia (i.e. drawn from near-IR)
or all are core-collapse (drawn from FUV). This points to some Type Ia
contribution to the observed H2 density CDF of SNRs. Revisiting this
experiment with the full M33 SNR population (which we can utilize once the
full CO map of M33 is available) will help reduce the shaded error region, and
put stronger constraints on the Ia contribution to the observed H2 densities.

##  Appendix D Lower resolution HI data

Figure 14 shows the same result as the middle panel of Figure 2, but with the
lower-resolution C-configuration ($\sim$21″$\sim$86 pc) H I data from Koch et
al. (2018). The increased surface brightness sensitivity of the lower
resolution data leads to a lower fraction of non-detections in the 3
categories of stars, specifically 2%, 3% and 0% for WRs, RSGs and SNRs
respectively.

## References

  * Ackermann et al. (2013) Ackermann, M., Ajello, M., Allafort, A., et al. 2013, Science, 339, 807, doi: 10.1126/science.1231160
  * Anderson et al. (2012) Anderson, J. P., Habergham, S. M., James, P. A., & Hamuy, M. 2012, MNRAS, 424, 1372, doi: 10.1111/j.1365-2966.2012.21324.x
  * Anderson & James (2008) Anderson, J. P., & James, P. A. 2008, MNRAS, 390, 1527, doi: 10.1111/j.1365-2966.2008.13843.x
  * Andersson et al. (2020) Andersson, E. P., Agertz, O., & Renaud, F. 2020, MNRAS, 494, 3328, doi: 10.1093/mnras/staa889
  * Andersson et al. (2022) Andersson, E. P., Agertz, O., Renaud, F., & Teyssier, R. 2022, arXiv e-prints, arXiv:2209.06218.  https://arxiv.org/abs/2209.06218
  * Andrews et al. (2017) Andrews, B. H., Weinberg, D. H., Schönrich, R., & Johnson, J. A. 2017, ApJ, 835, 224, doi: 10.3847/1538-4357/835/2/224
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Barker et al. (2007) Barker, M. K., Sarajedini, A., Geisler, D., Harding, P., & Schommer, R. 2007, AJ, 133, 1125, doi: 10.1086/511185
  * Baug et al. (2019) Baug, T., de Grijs, R., Dewangan, L. K., et al. 2019, ApJ, 885, 68, doi: 10.3847/1538-4357/ab46be
  * Belfiore et al. (2023) Belfiore, F., Leroy, A. K., Sun, J., et al. 2023, A&A, 670, A67, doi: 10.1051/0004-6361/202244863
  * Blaauw (1961) Blaauw, A. 1961, Bull. Astron. Inst. Netherlands, 15, 265
  * Bradley et al. (2023) Bradley, L., Sipőcz, B., Robitaille, T., et al. 2023, astropy/photutils: 1.9.0, 1.9.0, Zenodo, Zenodo, doi: 10.5281/zenodo.8248020
  * Braun (2012) Braun, R. 2012, ApJ, 749, 87, doi: 10.1088/0004-637X/749/1/87
  * Braun et al. (2009) Braun, R., Thilker, D. A., Walterbos, R. A. M., & Corbelli, E. 2009, ApJ, 695, 937, doi: 10.1088/0004-637X/695/2/937
  * Bresolin (2011) Bresolin, F. 2011, ApJ, 730, 129, doi: 10.1088/0004-637X/730/2/129
  * Bykov et al. (2018) Bykov, A. M., Ellison, D. C., Marcowith, A., & Osipov, S. M. 2018, Space Sci. Rev., 214, 41, doi: 10.1007/s11214-018-0479-4
  * Caldú-Primo & Schruba (2016) Caldú-Primo, A., & Schruba, A. 2016, AJ, 151, 34, doi: 10.3847/0004-6256/151/2/34
  * Calzetti et al. (2007) Calzetti, D., Kennicutt, R. C., Engelbracht, C. W., et al. 2007, ApJ, 666, 870, doi: 10.1086/520082
  * Caprioli & Spitkovsky (2014) Caprioli, D., & Spitkovsky, A. 2014, ApJ, 783, 91, doi: 10.1088/0004-637X/783/2/91
  * Castor et al. (1975) Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157, doi: 10.1086/153315
  * Chevance et al. (2020) Chevance, M., Kruijssen, J. M. D., Hygate, A. P. S., et al. 2020, MNRAS, 493, 2872, doi: 10.1093/mnras/stz3525
  * Cioni et al. (2011) Cioni, M. R. L., Clementini, G., Girardi, L., et al. 2011, A&A, 527, A116, doi: 10.1051/0004-6361/201016137
  * Conroy (2013) Conroy, C. 2013, ARA&A, 51, 393, doi: 10.1146/annurev-astro-082812-141017
  * Conti et al. (1988) Conti, P. S., Underhill, A. B., Jordan, S., et al. 1988, O stars and Wolf-Rayet stars, Vol. 497
  * Corbelli (2003) Corbelli, E. 2003, MNRAS, 342, 199, doi: 10.1046/j.1365-8711.2003.06531.x
  * Cosentino et al. (2022) Cosentino, G., Jiménez-Serra, I., Tan, J. C., et al. 2022, MNRAS, 511, 953, doi: 10.1093/mnras/stac070
  * Cronin et al. (2021) Cronin, S. A., Utomo, D., Leroy, A. K., et al. 2021, ApJ, 923, 86, doi: 10.3847/1538-4357/ac28a2
  * Crowther (2007) Crowther, P. A. 2007, ARA&A, 45, 177, doi: 10.1146/annurev.astro.45.051806.110615
  * Dalcanton et al. (2012) Dalcanton, J. J., Williams, B. F., Lang, D., et al. 2012, ApJS, 200, 18, doi: 10.1088/0067-0049/200/2/18
  * Dallas et al. (2022) Dallas, M. M., Oey, M. S., & Castro, N. 2022, ApJ, 936, 112, doi: 10.3847/1538-4357/ac8988
  * Davé et al. (2016) Davé, R., Thompson, R., & Hopkins, P. F. 2016, MNRAS, 462, 3265, doi: 10.1093/mnras/stw1862
  * de Wit et al. (2005) de Wit, W. J., Testi, L., Palla, F., & Zinnecker, H. 2005, A&A, 437, 247, doi: 10.1051/0004-6361:20042489
  * Díaz-Rodríguez et al. (2018) Díaz-Rodríguez, M., Murphy, J. W., Rubin, D. A., et al. 2018, ApJ, 861, 92, doi: 10.3847/1538-4357/aac6e1