# PHANGS-JWST First Results: A Statistical View on Bubble Evolution in NGC 628

Elizabeth J. Watkins, Ashley T. Barnes, Kiana Henny, Hwihyun Kim, Kathryn Kreckel, Sharon E. Meidt, Ralf S. Klessen, Simon C. O. Glover, Thomas G. Williams, Benjamin W. Keller, et al.

► To cite this version:

Elizabeth J. Watkins, Ashley T. Barnes, Kiana Henny, Hwihyun Kim, Kathryn Kreckel, et al. PHANGS-JWST First Results: A Statistical View on Bubble Evolution in NGC 628. The Astrophysical journal letters, 2023, 944, 10.3847/2041-8213/aca6e4. insu-04198461

HAL Id: insu-04198461

https://insu.hal.science/insu-04198461

Submitted on 7 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PHANGS–JWST First Results: A Statistical View on Bubble Evolution in NGC 628

Elizabeth J. Watkins¹ ID, Ashley T. Barnes² ID, Kiana Henny³ ID, Hwihyun Kim⁴ ID, Kathryn Kreckel¹ ID, Sharon E. Meidt⁵ ID,

Ralf S. Klessen⁶,⁷ ID, Simon C. O. Glover⁶ ID, Thomas G. Williams⁸,⁹ ID, Benjamin W. Keller¹⁰ ID, Adam K. Leroy¹¹,¹² ID,

Erik Rosolowsky¹³ ID, Janice C. Lee⁴,¹⁴ ID, Gagandeep S. Anand¹⁵ ID, Francesco Belfiore¹⁶ ID, Frank Bigiel² ID,

Guillermo A. Blanc¹⁷,¹⁸ ID, Médéric Boquien¹⁹ ID, Yixian Cao²⁰ ID, Rupali Chandar²¹ ID, Ness Mayker Chen¹¹,¹² ID,

Mélanie Chevance⁶,²² ID, Enrico Congiu¹⁸ ID, Daniel A. Dale³ ID, Sinan Deger²³ ID, Oleg V. Egorov¹ ID, Eric Emsellem²⁴,²⁵ ID,

Christopher M. Faesi²⁶ ID, Kathryn Grasha²⁷,²⁸ ID, Brent Groves²⁹ ID, Hamid Hassani¹³ ID, Jonathan D. Henshaw⁹,³⁰ ID,

Cinthya Herrera³¹, Annie Hughes³²,³³ ID, Sarah Jeffreson³⁴ ID, María J. Jiménez-Donaire³⁵,³⁶ ID, Eric W. Koch³⁴ ID,

J. M. Diederik Kruijssen²² ID, Kirsten L. Larson³⁷ ID, Daizhong Liu²⁰ ID, Laura A. Lopez¹¹,¹²,³⁸ ID, Ismael Pessa⁹,³⁹ ID,

Jérôme Pety⁴⁰,⁴¹ ID, Miguel Querejeta³⁵ ID, Toshiki Saito⁴² ID, Karin Sandstrom⁴³ ID, Fabian Scheuermann¹ ID, Eva Schinnerer⁹ ID,

Mattia C. Sormani⁶ ID, Sophia K. Stuber⁹ ID, David A. Thilker⁴⁴ ID, Antonio Usero³⁵ ID, and Bradley C. Whitmore⁴⁵ ID

¹ Astronomisches Rechen-Institut, Zentrum für Astronomie der Universität Heidelberg, Mönchhofstraße 12-14, D-69120 Heidelberg, Germany elizabeth.watkins@uni-heidelberg.de

² Argelander-Institut für Astronomie, Universität Bonn, Auf dem Hügel D-71, 53121, Bonn, Germany

³ Department of Physics and Astronomy, University of Wyoming, Laramie, WY 82071, USA

⁴ Gemini Observatory/NSF’s NOIRLab, 950 N. Cherry Avenue, Tucson, AZ 85726, USA

⁵ Sterrenkundig Observatorium, Universiteit Gent, Krijgslaan 281 S9, B-9000 Gent, Belgium

⁶ Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische Astrophysik, Albert-Ueberle-Straße 2, D-69120 Heidelberg, Germany

⁷ Universität Heidelberg, Interdisziplinäres Zentrum für Wissenschaftliches Rechnen, Im Neuenheimer Feld 205, D-69120 Heidelberg, Germany

⁸ Sub-department of Astrophysics, Department of Physics, University of Oxford, Keble Road, Oxford, OX1 3RH, UK

⁹ Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany

¹⁰ Department of Physics and Materials Science, University of Memphis, 3720 Alumni Avenue, Memphis, TN 38152, USA

¹¹ Department of Astronomy, The Ohio State University, 140 West 18th Avenue, Columbus, Ohio 43210, USA

¹² Center for Cosmology and Astroparticle Physics, 191 West Woodruff Avenue, Columbus, OH 43210, USA

¹³ Department of Physics, University of Alberta, Edmonton, AB, T6G 2E1, Canada

¹⁴ Steward Observatory, University of Arizona, 933 N Cherry Ave, Tucson, AZ 85721, USA

¹⁵ Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

¹⁶ INAF—Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50157, Firenze, Italy

¹⁷ The Observatories of the Carnegie Institution for Science, 813 Santa Barbara St., Pasadena, CA 91101, USA

¹⁸ Departamento de Astronomía, Universidad de Chile, Camino del Observatorio 1515, Las Condes, Santiago, Chile

¹⁹ Centro de Astronomía (CITEVA), Universidad de Antofagasta, Avenida Angamos 601, Antofagasta, Chile

²⁰ Max-Planck-Institut für Extraterrestrische Physik (MPE), Giessenbachstr. 1, D-85748 Garching, Germany

²¹ Ritter Astrophysical Research Center, University of Toledo, Toledo, OH 43606, USA

²² Cosmic Origins Of Life (COOL) Research DAO, coolresearch.io

²³ The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, AlbaNova, Stockholm, SE-106 91, Sweden

²⁴ European Southern Observatory, Karl-Schwarzschild-Straße 2, D-85748 Garching

²⁵ Univ Lyon, Univ Lyon1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230 Saint-Genis-Laval, France

²⁶ Department of Physics, University of Connecticut, 196A Auditorium Road, Unit 3046, Storrs, CT, 06269, USA

²⁷ Research School of Astronomy and Astrophysics, Australian National University, Canberra, ACT 2611, Australia

²⁸ ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Australia

²⁹ International Centre for Radio Astronomy Research, University of Western Australia, 7 Fairway, Crawley, 6009 WA, Australia

³⁰ Astrophysics Research Institute, Liverpool John Moores University, 146 Brownlow Hill, Liverpool, L3 5RF, UK

³² CNRS, IRAP, 9 Av. du Colonel Roche, BP 44346, F-31028 Toulouse cedex 4, France

³³ Université de Toulouse, UPS-OMP, IRAP, F-31028 Toulouse cedex 4, France

³⁴ Center for Astrophysics, Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138A, USA

³⁵ Observatorio Astronómico Nacional, Alfonso XII 3, E-28014, Madrid, Spain

³⁶ Centro de Desarrollos Tecnológicos, Observatorio de Yebes (IGN), E-19141 Yebes, Guadalajara, Spain

³⁷ AURA for the European Space Agency (ESA), Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

³⁸ Flatiron Institute, Center for Computational Astrophysics, NY 10010, USA

³⁹ Leibniz-Institut für Astrophysik Potsdam (AIP), An der Sternwarte 16, D-14482 Potsdam, Germany

⁴⁰ IRAM, 300 rue de la Piscine, F-38400 Saint Martin d'Hères, France

⁴¹ LERMA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, 75014, Paris, France

⁴² National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo, 181-8588, Japan

⁴³ Center for Astrophysics and Space Sciences, Department of Physics, University of California, San Diego, 9500 Gilman Drive, La Jolla, CA 92093, USA

⁴⁴ Department of Physics and Astronomy, The Johns Hopkins University, Baltimore, MD 21218, USA

⁴⁵ Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

Received 2022 October 20; revised 2022 November 25; accepted 2022 November 29; published 2023 February 16

Abstract

The first JWST observations of nearby galaxies have unveiled a rich population of bubbles that trace the stellar- feedback mechanisms responsible for their creation. Studying these bubbles therefore allows us to chart the

<figure><img src="image_2.png" /></figure>
interaction between stellar feedback and the interstellar medium, and the larger galactic flows needed to regulate star formation processes globally. We present the first catalog of bubbles in NGC 628, visually identified using Mid-Infrared Instrument F770W Physics at High Angular resolution in Nearby GalaxieS (PHANGS)-JWST observations, and use them to statistically evaluate bubble characteristics. We classify 1694 structures as bubbles with radii between 6 and 552 pc. Of these, 31% contain at least one smaller bubble at their edge, indicating that previous generations of star formation have a local impact on where new stars form. On large scales, most bubbles lie near a spiral arm, and their radii increase downstream compared to upstream. Furthermore, bubbles are elongated in a similar direction to the spiral-arm ridgeline. These azimuthal trends demonstrate that star formation is intimately connected to the spiral-arm passage. Finally, the bubble size distribution follows a power law of index $p = -2.2 \pm 0.1$, which is slightly shallower than the theoretical value by $1-3.5\sigma$ that did not include bubble mergers. The fraction of bubbles identified within the shells of larger bubbles suggests that bubble merging is a common process. Our analysis therefore allows us to quantify the number of star-forming regions that are influenced by an earlier generation, and the role feedback processes have in setting the global star formation rate. With the full PHANGS-JWST sample, we can do this for more galaxies.

*Unified Astronomy Thesaurus concepts: Superbubbles (1656); Stellar wind bubbles (1635); Infrared astronomy (786); H II regions (694); Stellar feedback (1602)*

# 1. Introduction

A ubiquitous feature of new mid-IR (MIR) JWST observations of nearby galaxies are the parsec- to kiloparsec-scale holes that riddle the interstellar medium (ISM; e.g., Figure 1). Many of these holes are carved out by stellar-feedback processes, and so trace HII regions, bubbles and superbubbles, although some holes are dynamically driven regions without a powering source or are remnants of older, sheared bubbles (Palous et al. 1990; Kim et al. 2002; Dobbs & Pringle 2013). Most of the smallest holes (hereafter bubbles) with sizes of ~6 pc are driven by pre-supernova (pre-SN) feedback and photoionization from individual high-mass stars and low-mass stellar OB associations (Lim et al. 2020), though some could be young (<10⁴ yr) supernova (SN) remnants. Larger bubbles (up to ~30 pc) are dominated by winds and SNe connected to OB associations (Ostriker & McKee 1988). The largest features, superbubbles, represent the combined impact of large or multiple OB associations that can drive bubbles to radii of ~1000 pc in 40 Myr (Nath et al. 2020).

Considering bubbles are such a dominant feature for a range of sizes spanning over two orders of magnitudes (see Figure 1), understanding what creates them, and their evolution, is crucial to understand how they impact the ISM, their role in driving the galactic-scale star formation cycle, and the processes that result in low star formation efficiencies compared to freefall timescales (Hopkins et al. 2014; Federrath 2015; Grudić et al. 2019; Keller et al. 2022). For example, the stellar populations that drive bubbles also contribute to the large-scale injection of turbulence that redistributes energy and matter, limiting the large-scale collapse of gas into stars (Ostriker et al. 2010; Faucher-Giguère et al. 2013; Krumholz & Burkhart 2016; Fisher et al. 2019) and leading to galactic-scale outflows (Fielding et al. 2018; Smith et al. 2021). The overall size distribution of bubbles, in particular, provides insight into how clustered star formation is within galaxies, and the pressure balance between different phases of the ISM (Nath et al. 2020).

In the Milky Way, we can reach the resolution needed to study and catalog bubbles at the earliest phases of their formation and evolution (Simpson et al. 2012; Anderson et al. 2014; Beaumont et al. 2014; Jayasinghe et al. 2019; Olivier et al. 2021). However, line-of-sight confusion within the Milky Way requires complex kinematic analysis (Reynolds & Ogden 1979; Ehlerova & Palous 1996; Ochsendorf et al. 2015; Tsujimoto et al. 2021; Zucker et al. 2022), whereas extragalactic observations provide the necessary context to easily connect their shells—and any cospatial

information such as powering sources—to larger scale features (Bagetakos et al. 2011; Watkins et al. 2022). Yet, replicating Milky Way studies in extragalactic environments has proven challenging due to the high resolution (~10 pc) required in atomic and molecular gas to resolve their edges, and in optical emission to resolve the hot ionized gas in their interiors.

With the high resolution (~10 pc) that JWST can reach tracing dust emission within nearby galaxies, we can finally bridge the scale gap between extragalactic and galactic studies. JWST observations are already revolutionizing our understanding of stellar feedback and the star formation cycle in nearby galaxies over a large range of evolutionary stages. MIR bands, such as the Mid-Infrared Instrument's (MIRI) F770W, trace hot dust heated by young stars and polycyclic aromatic hydrocarbons (PAHs), where PAHs are vibrationally excited in the presence of starlight (Sandstrom et al. 2023), especially when illuminated by UV photons. Therefore MIR observations allow us to identify new, young embedded clusters obscured at optical wavelengths (Rodríguez et al. 2023), large-scale filamentary structures containing dense, cold gas expected to host future star formation (Thilker et al. 2023), and hot dust emission shining in the presence of UV radiation emitted by OB stars (Leroy et al. 2023). Piecing these results together provides the observations needed to trace recent star formation histories within these galaxies (Kim et al. 2023). In this letter we provide a crucial piece needed to understand the star formation picture—connecting together small- and large-scale feedback processes—by cataloging the feedback-driven bubbles in NGC 628 using Physics at High Angular resolution in Nearby GalaxieS (PHANGS)-JWST.

# 2. Observations

## 2.1. NGC 628

NGC 628 (also known as M74 or “the Phantom Galaxy”) is a nearly face-on spiral galaxy at a distance of 9.84 ± 0.63 Mpc (McQuinn et al. 2017; Anand et al. 2021). Its star formation history and multiphase ISM have been investigated extensively in previous works (Sanchez et al. 2011; Grasha et al. 2015; Kreckel et al. 2018), and it has cospatial PHANGS observations—PHANGS-Atacama Large Millimeter/submillimeter Array (ALMA; Leroy et al. 2021), PHANGS-Multi Unit Spectroscopic Explorer (MUSE; Emsellem et al. 2022), and PHANGS-Hubble Space Telescope (HST; Lee et al. 2022)—providing a more comprehensive view of the star formation cycle within this galaxy.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Multi-scale bubble population in NGC 628 revealed with JWST observations. RGB illustration toward the center of NGC 628 highlighting bubble features and the bands used in their identification. Red traces MIRI F770W observations, green traces MUSE Hα observations, and blue traces B-band HST observations. White scalebar indicates physical size equal to 1 kpc.</figcaption></figure>

## 2.2. JWST Observations

NGC 628 is one of 19 galaxies being observed as part of the PHANGS$^{46}$–JWST Cycle 1 treasury program using the Near Infrared Camera’s (NIRCam) F200W, F300M, F335M, and

F360M, and MIRI’s F770W, F1000W, F1130W, and F2100W broadband filters (project-ID 02107; see Lee et al. 2023). It was observed on 2022 July 17 and the data have been reduced using the publicly available reduction pipeline along with additional reduction tools developed and outlined in Lee et al. (2023). These observations cover the central 3'8 times 2'2 (11 kpc × 6.3 kpc) of NGC 628, and for the F770W observations

<footnote>$^{46}$ http://www.phangs.org</footnote>
they have a sensitivity of 0.11 MJy sr⁻¹ and an angular resolution of ~0.′′25, which achieves a physical resolution of 12 pc (Lee et al. 2023).

## 2.3. Ancillary Data

NGC 628 was observed in five bands by the HST (WFC3/F275W, WFC3/F336W, ACS/F435W, ACS/F555W, and ACS/F814W) as part of the LEGUS survey (Calzetti et al. 2015), and reduced using the PHANGS-HST data pipeline (Lee et al. 2022). It was also observed in Hα (F658N; proposal 10402). In this Letter, we use the B-band observations (~0.′′08 ~ 4 pc) to visually identify younger sources to help separate the nested bubble population revealed in JWST observations.

Continuum-subtracted Hα imaging is taken from the PHANGS-MUSE survey (Kreckel et al. 2018; Emsellem et al. 2022), an optical integral field unit survey targeting 19 nearby galaxies (the same as PHANGS–JWST) at 0.′′9 (~50 pc) resolution. This survey traces local ionized gas properties and metallicities, which can be used to constrain the local star formation history and gas mixing scales. AstroSat (H. Hassani et al. 2022, in preparation) far-UV observations trace the ionized emission from older stellar populations at resolutions of ~1.′′5, and Spitzer observations of 8 µm emission (Dale et al. 2009) at ~2′′ allow us to compare the previous MIR observations to JWST. We also take into account HI data from the THINGS survey at resolutions of ~6′′ (Walter et al. 2008) and PHANGS-ALMA $^{12}$CO (J = 2–1) (Leroy et al. 2021) (at ~1′′) to investigate bubble shells in different bands.

# 3. A Panchromatic Picture of Bubble Features in NGC 628

In Figure 2, we present a panchromatic picture of NGC 628. Bubble features, such as shells, stellar sources, and hot emission enclosed by the bubbles, are detectable throughout the electromagnetic spectrum (from X-ray to radio), and especially pronounced at MIR bands, which is why MIR broadband filters have been used extensively to study bubbles in the Milky Way (Churchwell et al. 2006; Anderson et al. 2014; Jayasinghe et al. 2019).

For our panchromatic view, we chose bands that highlight the full range of different bubble features. For example, MIRI F770W and F1130W have similar emission features, but we only show F770W maps as they are at a higher resolution. Altogether, we consider AstroSat far-UV, HST *B*, *V*, and *I* bands, and narrowband Hα, MUSE Hα, NIRCam F335M, MIRI F770W and F2100W, Spitzer 8 µm, ALMA $^{12}$CO (J = 2–1), and THINGS HI moment 0 maps. For the NIRCam bands, we only show the F335M band since it contains both stellar sources and some weak PAH emission, potentially allowing us to see both the bubble shells and sources within them.

In F2100W emission, MUSE Hα, and AstroSat far-UV, we can see concentrated ionized emission, and some hints of ionized shells in MUSE Hα, but it is difficult to define their boundaries due to confusion with the more diffuse ionized gas (Belfiore et al. 2022). We also note that the vast majority are not resolved in the ~50 pc resolution MUSE observations (see, e.g., Barnes et al. 2021; Santoro et al. 2022).

In HST and JWST NIRCam imaging, we see bright clustered sources along the arms, though HST shows them more clearly

toward the center of the galaxy. While there is PAH emission in the NIRCam F335M map, it is faint compared to the MIRI F770W emission, which is expected given the relative intrinsic strength of the two features.

Molecular gas offers a clear constraint on bubble shells, particularly in the case of the larger (>50 pc) bubbles. The high spectral resolution achievable with molecular line observations also offers dynamical information (e.g., CO at 2.5 km s⁻¹; Leroy et al. 2021), providing kinematic confirmation of coherent bubbles structures, and constraining properties of the powering sources (Watkins et al. 2022). However, CO provides a limited sample of bubble structures. They must be large enough to be resolved by ALMA (~50 pc resolution), and they also need to be young enough (<8 Myr) to still contain molecular gas (Kim et al. 2021). In Watkins et al. (2022), we only find 12 unbroken bubbles in NGC 628 in CO due to these constraints. HI emission offers an alternative, as can be seen in the THINGS HI map (Bagetakos et al. 2011; Krause et al. 2015), but at 200 pc resolution it resolves only the largest bubble features (Barnes et al. 2023).

While lacking the kinematic confirmation, broadband dust (PAH) emission at MIR probes the three gas phases (Leroy et al. 2023). In addition, stellar populations found inside bubbles produce UV photons in abundance (Peeters et al. 2004). As the underlying dust and stellar continuum contributions are relatively weak at ~8 µm (Marble et al. 2010) and there is increased gas (and therefore dust) column densities in swept-up shells, bubble features stand out against the background. This makes them highly visible at a glance in MIRI F770W observations in Figures 1 and 2 at high (12 pc) resolution. While we can still see some larger bubble features toward the top of the image and in the lower left in the Spitzer 8 µm map, it lacks the resolution to probe smaller bubble features (Dale et al. 2009), limiting the number of bubbles we can see.

The high-resolution MIRI F770W shows that two physically different structures exist across the galaxy (also see Figures 1 and 3, and Figure 8 in Appendix A):

1. Spherical or elliptical “bubbles” that are closely linked downstream from a given arm (convex side) and are driven by feedback.

2. Very elongated elliptical “holes” that are situated between the arms, toward the upstream part of a given arm (concave side) and are old sheared bubbles, or dynamically created holes.

Dynamically, NGC 628 is rotating clockwise and the arms are trailing. Hence, inside corotation (4.4 ± 2.0 kpc; Williams et al. 2021), the fact that the disk material moves faster than the spiral puts recent star formation on the downstream side of the arms, with age increasing away from the arm (discussed further in Section 5). In this picture, the feedback from the newly formed star clusters preferentially flows into the bubbles that are already inflated by previous star formation associated with this arm.

The holes are therefore features formed by galactic dynamical processes (e.g., spurs or feathers; Kim et al. 2002; La Vigne et al. 2006; Dobbs & Pringle 2013; Williams et al. 2022) with morphologies shaped and elongated by galactic shear, or old superbubbles created by stellar feedback that are later sheared (e.g., Palous et al. 1990). The shells of older bubbles that have been dynamically sheared are likely the
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Multiwavelength observations covering the nearby star-forming galaxy NGC 628. In the top row of panels from left to right, we show the Hα map from PHANGS-MUSE (Emsellem et al. 2022), CO(2–1) peak intensity from PHANGS-ALMA, (Leroy et al. 2021), and the Very Large Array (VLA) Hi moment 0 map from the THINGS survey (Walter et al. 2008). In the middle row, left panel, we show an image produced from the *B*- (blue), *V*- (green) and *I*- (red) band filters from HST, and the continuum-subtracted HST Hα (658 nm, red; see Lee et al. 2022). In the middle center and right panels, we show the UV emission map from AstroSat (H. Hassani et al. 2022, in preparation) and NIRCam F335M emission from JWST (Lee et al. 2023). On the bottom row, from left to right, we show the 8 µm map from Spitzer (Dale et al. 2009), and the F770W and the F2100W maps from MIRI JWST (Lee et al. 2023). All panels have been overlaid with a dashed contour showing the coverage of the JWST observations.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Bubble locations identified in NGC 628 with multi-scale view zoom ins providing a more detailed view. Left: MIRI F770W band map of NGC 628 using a square root intensity stretch. Cyan ellipses show the location of the 1694 bubbles identified. Labeled lines indicate position of spiral arms derived by skeletonizing the arm environment masks from Querejeta et al. (2021). Solid circle markers on arm spines indicate the start position of the arm. White boxes labeled *a* and *b* indicate locations of close-up panels on right side of this figure. Right: same RGB as Figure 1 but focused on the regions indicated on the left panel. Bars at bottom left indicate the physical size scale.</figcaption></figure>

building blocks of dust lanes, as they get compressed upon entering the arm (Thilker et al. 2023). By contrast, bubbles are formed by recent stellar feedback.

To understand the link between star formation, feedback, and the larger environment, this Letter focuses only on bubbles currently being driven by feedback processes. While MIRI F770W provides the clearest view of holes, by itself it does not allow us to separate the feedback-driven holes from the dynamically created holes. Combining the MIRI F770W with HST imaging preserves the high-resolution view and helps delineate feedback-driven regions and dynamically created holes via the stellar populations within the shell. Furthermore, overlaying the HST Hα observations (shown as green and red in the three color image in Figures 1 and 2), we see that the spherical bubbles along the arms are usually traced by a shell of ionized gas—whereas the elongated holes typically have little ionized gas emission. If we instead combine the higher sensitivity MUSE Hα with F770W and HST B band in Figure 1 (the B band highlights the young OB stars that drive bubbles), we can separate feedback-driven bubbles from dynamically driven holes.

# 4. Identifying Bubbles

In this section, we outline how we identify bubbles using the multiwavelength data sets (MIRI F770W, MUSE Hα and HST B-band observations) outlined in the previous section. We note here that identifying shells in F770W MIRI observations is still the primary method we use, but for small, or highly elliptical features, we rely on HST sources and spatial MUSE Hα to help distinguish between old and dynamically driven holes, and those more likely to be feedback driven.

## 4.1. *Identification Method*

To find bubble shells, two approaches are possible: automated algorithms and manual methods. Automated algorithms are reproducible and have well-defined problems and biases (e.g., Thilker et al. 1998; Silburt et al. 2019; Van Oort et al. 2019; Collischon et al. 2021). However, they usually struggle in lower signal-to-noise ratio situations, or when presented with complex features such as broken or elliptical shells. Furthermore, automated methods usually require fine tuning or careful training sets to output reasonable catalogs.
Manual procedures are better able to utilize multiwavelength data when identifying bubbles, and human pattern recognition is superior to automated methods at identifying weak bubble features, but they are subjective and cannot be exactly reproduced.

JWST presents us with a brand new data set. Neither the exact features nor the associated scales are known. Automatic algorithms and citizen science approaches, such as machine learning (Van Oort et al. 2019) and Zooniverse projects (Jayasinghe et al. 2019), respectively, need to be trained with already tagged catalogs. We thus opt for manual bubble identification to deliver such a tagged catalog for future works. The aim of this catalog is not, therefore, to fully characterize all the feedback-driven bubbles present, but to identify the majority of them in an as unbiased way as possible. There are undoubtedly some bubbles that we have missed, and some which are false-positive detections. In Section 6, we discuss catalog completeness further.

We identify bubbles using an RGB image composite of MIRI F770W (red), MUSE Hα (green), and HST B-band (blue) observations using an asinh, log, and linear stretch, respectively, which we illustrate in Figure 1. We mainly focus on the shell-like features present in MIRI F770W. If a feature is easy to identify as a bubble, and has a complete shell in F770W, we usually identify it without other information. For highly eccentric features, or small, clustered features, we require HST sources or Hα emission within the bubble. For this, our definition of a shell includes any circular or elliptical feature with a radius greater than the physical resolution (6 pc), with no limit on the shell thickness in F770W. Partial features are included, such as continuous, but incomplete, circular arcs, or fragmented, clumpy circular shells. To find bubble features, we adjust the contrast for each band throughout the bubble-identification process. Finally, to identify bubble structures at different spatial scales, we use different zoom levels starting at a quarter of the image, down to a ~15" × 15" tile. For a more detailed discussion of how we identify bubbles, see Appendix A.

For each bubble, we use elliptical or circular apertures to trace the shell seen in MIRI observations (by eye), and we fit them using the shell ridgeline (as opposed to the inner or outer edge of the shell). We set no limits on the eccentricity or the position angles (PAs). To speed up the bubble identification, PAs are typically incremented in 10° steps.

We first perform this search using one member of the team (E.J.W.) who tried to identify most bubbles present. Two additional catalogs are generated by different team members (K.H., H.K.) using the same method to check the robustness of our results. Section 5 presents our results based on the first catalog, and Appendix B lists common findings and potential variations. We note here that all results pertaining to the bubble catalogs refer to the primary catalog (found by E.J.W.), unless otherwise stated.

## 4.2. Catalog Description

The bubble catalog contains 1694 bubbles, plotted in Figure 3; we tabulate the first 10 in Table 1.47 The properties listed include their ID, position (in R.A. and decl.), their semimajor and semiminor axes in parsec (we expect 10%

measurement uncertainties on their sizes via Watkins et al. 2022), the average radii in parsec, their PA (positive angles from north to east, with ±5° uncertainties), which spiral arm they are closest to (as defined in Section 5), the distance to this arm (where positive and negative distances are downstream and upstream, respectively), and their galactocentric distance. PAs are defined between -90° and 90°. While -90° and 90° are the same, we retain the minus sign to indicate which direction they are rounded from. We note here that while ellipses are generally used to find bubble candidates, we use circularized average radii when summarizing their sizes in the following sections.

# 5. Results

## 5.1. Global Catalog Properties

Bubble radii range between 6 and 552 pc with a mean and median of 34 and 27 pc, respectively, with a 16th–84th percentile spread of 17–45 pc. Since the average bubble diameter is 4.6–5.6 times the angular resolution of MIRI F770W, these are well resolved. The smallest bubbles also match our resolution limit, with a diameter of 12 pc, indicating that even smaller bubbles could be present for higher resolution observations.

On Figure 4, we plot the size distribution of all bubbles identified. The distribution should follow a power law of $N(R) \propto R^p$, where N is the number of bubbles, R are the bubble radii, and p is the power-law index of the distribution. At ~30 pc, the distribution visibly turns over and follows this power-law distribution. The turnover is not a physical limit to bubbles, but represents the size scale we can sample bubbles down to completely. If ambient pressures are low, some models predict the distribution can turn over at ~100 pc (Nath et al. 2020). However, we do not expect to see such low pressures toward the center of a main-sequence galaxy like NGC 628 (Sun et al. 2020), and the turnover point we measure is too low compared to the theoretical prediction.

To fit the power-law slope correctly, a precise measurement of the turnover point is needed. To find the optimal turnover point (and the corresponding p-value), we perform Pareto's maximum-likelihood estimate (MLE) using the package powerlaw and recalculate the power-law slope after removing the smallest bubble iteratively. For each solution, we perform a Kolmogorov–Smirnov (KS) test between the model and the data and choose the model that minimizes the KS value (Alstott et al. 2014). To quantify uncertainties on p and the turnover point, we bootstrap the bubble sizes for a 1694 bubble distribution, with replacement, 10,000 times and recalculate the optimal turnover point (i.e., the minimal KS distance), and the corresponding power-law index. We find $p = -2.2 \pm 0.1$ with a turnover of 29 ± 3 pc using the median bootstrapped solution (which is the same as the mean) and we quote the uncertainties using the 16th and 84th percentiles.

On smaller scales, we find bubbles are often nested within the shells of even larger bubbles with the major axis running parallel to the tangent to the shell of the larger parent bubbles. We quantify this nesting by counting the total number of bubbles that intersect with the boundary of larger bubbles, binned as a function of the larger bubble's radius, and we plot this distribution in the middle panel of Figure 4. We define the bubble boundary using an annular mask 0.9–1.1 times the radius for the inner and outer edge and cross-match them with

<footnote n="47">47 All tabulated properties can be found at https://www.canfar.net/storage/vault/list/phangs/RELEASES/Watkins_etal_2022.</footnote>
<table><caption>**Table 1**
Sample of 10 Bubbles Ordered by R.A. and their Defining Properties</caption><thead><tr><th>Bubble ID</th><th>R.A.</th><th>Decl.</th><th>Semimajor axis (pc)</th><th>Semiminor axis (pc)</th><th>Average Radius (pc)</th><th>PA (°)</th><th>Arm</th><th>Distance to arm (pc)</th><th>Galactocentric Radius (kpc)</th></tr></thead><tbody><tr><td>1</td><td>24°8′55″ 9</td><td>15°48′7″ 9</td><td>86</td><td>86</td><td>86</td><td>0</td><td>3</td><td>145</td><td>5.22</td></tr><tr><td>2</td><td>24°8′56″ 7</td><td>15°48′2″ 9</td><td>24</td><td>24</td><td>24</td><td>0</td><td>3</td><td>127</td><td>5.05</td></tr><tr><td>3</td><td>24°8′57″ 0</td><td>15°48′9″ 0</td><td>42</td><td>42</td><td>42</td><td>0</td><td>3</td><td>138</td><td>5.21</td></tr><tr><td>4</td><td>24°8′58″ 0</td><td>15°47′44″ 9</td><td>25</td><td>25</td><td>25</td><td>0</td><td>3</td><td>101</td><td>4.54</td></tr><tr><td>5</td><td>24°8′58″ 2</td><td>15°48′2″ 6</td><td>24</td><td>24</td><td>24</td><td>0</td><td>3</td><td>114</td><td>4.98</td></tr><tr><td>6</td><td>24°8′58″ 7</td><td>15°47′41″ 2</td><td>33</td><td>33</td><td>33</td><td>0</td><td>3</td><td>95</td><td>4.44</td></tr><tr><td>7</td><td>24°8′58″ 8</td><td>15°47′50″ 6</td><td>32</td><td>31</td><td>32</td><td>0</td><td>3</td><td>93</td><td>4.64</td></tr><tr><td>8</td><td>24°8′58″ 9</td><td>15°47′39″ 8</td><td>40</td><td>40</td><td>40</td><td>0</td><td>3</td><td>93</td><td>4.40</td></tr><tr><td>9</td><td>24°8′58″ 9</td><td>15°48′9″ 0</td><td>31</td><td>31</td><td>31</td><td>0</td><td>3</td><td>122</td><td>5.14</td></tr><tr><td>10</td><td>24°8′59″ 3</td><td>15°48′6″ 1</td><td>209</td><td>171</td><td>196</td><td>0</td><td>3</td><td>112</td><td>5.04</td></tr></tbody></table>

Note: Distance used to calculate physical lengths is 9.84 ± 0.63 Mpc. See Section 4.2 for more information about uncertainties.

the unaltered bubble masks. If any part of the bubble mask overlaps with the annular mask of a larger bubble, we counted it as a nested bubble. However, since such a comparison will result in a size bias (larger bubbles are more likely to overlap by chance due to their larger areas), we also calculate the same distribution after shuffling the positions of the bubbles 100 times, and compare the two distributions in the middle panel of Figure 4.

Altogether, we find 31% of bubbles overlap with at least one smaller bubble, and bubble nesting occurs 3.2 times more often than expected from random chance, in total, when comparing the two. In future work, we will explore the exact nature of this relationship in more detail (i.e., whether the nested bubbles represent regions of new star formation, or have been simply relocated or uncovered by the expansion of the larger bubble; see Barnes et al. 2023).

## 5.2. Azimuthal and Radial Trends

When considering the locations of bubbles, Figure 3 shows they are not distributed uniformly but follow the large-scale structure of the galaxy as seen by PAH emission. Specifically, bubbles are found closer to the bright emission associated with the spiral arms, with larger bubbles downstream from the arm. While an age gradient is not identified in the young stellar clusters (Shabani et al. 2018), we clearly see a systematic offset between star formation and molecular gas in Figure 1 in Kreckel et al. (2018).

To quantify what we see, we analyze the bubble properties in relation to the spiral arms. For the spiral-arm model, we define a spiral-arm ridge for each by skeletonizing the spiral-arm environment masks calculated in Querejeta et al. (2021), and show them in Figure 3. For each bubble outside of the galactic center (defined using the same environment masks), we identify which arm they are closest to using the difference between the bubble center and the closest approach to the spiral-arm ridge. We provide these assignments in Table 1, keeping track of which side of the arm (upstream, downstream) they are on. With these values, we measure the separate size distributions of bubbles upstream and downstream from their nearest arm on the bottom panel of Figure 4. While we find a similar number of bubbles downstream and upstream from the arms (794 versus 790, respectively), the bubbles downstream are somewhat larger on average. The mean and median radii are larger by 18% and 9%, respectively, downstream, and their distribution has a more pronounced high-end tail (the radii at

the 99th percentile is 220 pc for downstream bubbles, while the same for upstream bubbles is only 149 pc, a 48% difference).

To check how confident we can be that the two distributions are different—informing whether there is a statistical difference between bubbles downstream from an arm compared to upstream—we performed a Mann-Whitney test on them. The statistic indicates there is only a 0.01% chance the two distributions have the same underlying distribution, thus we reject this null hypothesis. As explained in Section 3, we expect to see larger bubbles downstream for stars forming inside corotation due to the flow of gas into and through spiral arms, causing the recently formed, young stars (and the bubbles they produce) to pile up downstream. These results also reinforce that almost all the bubbles found are feedback driven, rather than dynamical holes.

We also see that, in general, the PAs of bubbles correlate with the spiral arms. That is, the major axes of bubbles are somewhat parallel with the tangent angles of the spiral arms, though we also see some are perpendicular. If shear is responsible for noncircular bubbles, we expect PAs to correlate with the spiral arm tangent inside corotation (Palous et al. 1990). Bubbles offset perpendicular to the spiral arms instead represent where star formation was delayed and formed a bubble on the downstream side of the arm, which then blistered in the direction of the lower density medium downstream.

To quantify and confirm these trends, we calculate the difference ($\Delta\theta$) between the PAs of elliptical bubbles (i.e., with aspect ratios (ARs) >1) and the spiral-arm tangent angle for bubbles outside of the galactic center and plot their distribution on the top panel of Figure 5 in black. Six-hundred and fifty-four bubbles match these criteria. Positive $\Delta\theta$ angles are offset anticlockwise from the arm, while clockwise $\Delta\theta$ values are negative, resulting in a $\Delta\theta$ distribution between -90° and 90°.

We take the measurement uncertainty into consideration in Figure 5 by assuming the PA of each bubble has a Gaussian distribution with a standard deviation of 5° and summing their distributions together to create a kernel density estimation (KDE). We then quantify the uncertainty between the KDE and a binned histogram by calculating the mean difference between the data binned in 5° intervals and the KDE distribution generated using double the uncertainty (10°). This results in an uncertainty of $\pm 2.6$ bubbles.

Figure 5 reveals peaks at 0° and -30° and -90°. In addition, bubbles at 30° are underrepresented (and, in general, positive $\Delta\theta$ are lacking). The physical meaning behind the peak and trough at $\pm 30°$ is not immediately obvious. Due to our
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Size distributions of bubbles. Top: size distribution of the entire sample of bubbles. Solid gray and red lines are the mean and median of the distribution, respectively, shaded blue region is the 16th–84th sigma range (labeled as σ in the legend), and dashed cyan line is the power-law fit to the distribution equal to $p = -2.2$ derived using MLE. Inset shows the bootstrapped distribution of the power-law index in the bottom panel with a central value of −2.2 with a 16th–84th sigma range of −2.3–2.1. Middle: total number of nested bubbles within the boundary of a larger bubble. The x-axis corresponds to the radius of the larger bubble. White distribution shows the data while coral shows the same distribution if instead bubble locations were randomized, averaged over 100 different randomized realizations. Bottom: size distribution for bubbles upstream and downstream from a spiral arm, shown in orange and blue, respectively, with the same colored lines showing their 99th percentile values. A Mann–Whitney test comparing the two distributions indicates 0.01% chance that the two distributions are equal.</figcaption></figure>

selection criteria in Section 4, the most sheared bubbles are excluded, which, if included, should have positive $\Delta\theta$ inside corotation as they enter the arm upstream. The peak at −30° is likely tracing the bubbles downstream that lag behind the angle traced out by a spiral arm after it has passed.

To confirm that these peaks are statistically significant, we first estimate what Figure 5 would look like if the PAs were randomly distributed. We randomize the bubbles PAs 100 times, recalculate $\Delta\theta$, and average these distributions together. We find the randomization results in a uniform distribution, which we show using a dashed red line in Figure 5. This test confirms that the spiral arms do not have a preferred orientation, which would result in a bias in the random

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Offset angles between elliptical bubbles and the spiral-arm tangent angle. Top: solid black and blue lines with apertures show offset angles for bubbles with aspect ratios >1 and >1.15, respectively, for bubbles with an uncertainty in measured angle of $\pm 5^\circ$, resulting in a KDE-like distribution. The blue line has been renormalized to match the area of the black line (i.e., from 301 bubbles to 654 bubbles) to make comparisons easier. Translucent regions indicate the uncertainty. Dashed red line shows the distribution of offset angles if bubbles have random PAs. Bottom: same as blue line in top panel (i.e., AR > 1.15) but for bubbles >40 pc upstream (black) and downstream (purple) from an arm.</figcaption></figure>

distribution. Consequently, we can perform a Rayleigh test on $\Delta\theta$, which checks if a periodic distribution is nonuniform. The statistic indicates there is a 0.07% chance that the distribution of $\Delta\theta$ is uniform, therefore we reject the null hypothesis. Finally, to confirm any trends present are not driven by bubbles that might have less well-defined PAs (caused by having lower ARs), we repeat our analysis after excluding bubbles with an AR $\le 1.15$ and plot the result in blue in Figure 5 after renormalizing the area to improve the visual comparison. We see the same peaks and troughs, and the Rayleigh test has a $p$-value of 2.23%. Altogether, these tests confirm the peaks and troughs in Figure 5 are real and statistically robust.

To investigate if perpendicular $\Delta\theta$ represent a population of blistering bubbles, we remake Figure 5 on the bottom panel for bubbles upstream and downstream. If blistering is the cause, we expect more bubbles with perpendicular $\Delta\theta$ downstream than upstream. We focus on more elliptical bubbles (with AR > 1.15) since blistering should cause bubbles to reach higher eccentricities, and to ensure we are not dominated by smaller bubbles we limit the analysis to bubbles >40 pc. Again, we renormalize the KDEs to improve visual comparisons. Indeed, we find that bubbles downstream have perpendicular $\Delta\theta$, which are not present upstream. Moreover, the downstream distribution passes the Rayleigh test with a $p$-value of 2.99%, whereas the upstream distribution fails. However, in
<figure><img src="image_8.png" /><figcaption>**Figure 6.** KDE distributions of bubble sizes binned as a function of galactocentric distance in 1 kpc rings. Dashed and solid black lines show the median and mean values, while the solid pink, red, and cyan lines show the 84th, 95th, and the 99th percentiles, respectively. The color of the distributions shows the number of bubbles per unit surface area in annular rings at the distance indicated on the y-axis. The number labeling each distribution indicates the number of bubbles within the distribution.</figcaption></figure>

the independently generated catalogs presented in Appendix B we do not find a strong peak at ±90°. We expect the lower number of bubbles sampled in these catalogs is the cause. We therefore leave this as a tentative result and will explore whether bubbles align perpendicular to spiral arms in more detail in future work.

The last statistical trend we explore is how bubble radii change as a function of galactocentric radius. Typically, we expect bubble sizes to increase further away from the center due to higher ambient pressures in galaxy centers confining the bubble sizes (Bagetakos et al. 2011; Barnes et al. 2020). The trend can also be driven by the orbital time. Near the center, bubbles could be sheared before they have time to grow, resulting in smaller bubbles toward the center.

Using Kendall's τ correlation coefficient, we find there is a weak correlation of 0.10 (with a p-value of 1.0 × 10⁻⁸%) between bubble radii and the galactocentric radius. However, the trend is not obvious when plotting the two against each other. To view the trend, we split the catalog into 1 kpc rings as a function of galactocentric radius and plot their distributions in Figure 6. While weak, the average radii increase as a function of galactocentric distance, and a larger fraction of bubbles are found in the high-end tail at larger galactocentric distances. Therefore, the processes that impact bubble sizes radially have a secondary role in our observations. We expect the footprint of the MIRI F770W observations impacts our ability to measure increasing bubble size with galactocentric radii. Our data only cover the inner part of the galaxy, which might impede our ability to view this trend, and the shape of the footprint also results in a poor sampling of bubbles at larger galactic radii, which we indicate by coloring the distributions in Figure 6, by the number of bubbles per surface area.

# 6. Discussion

## 6.1. Understanding What Drives the Shallower Bubble Size Distribution

In Section 5.1, we find the bubble size distribution is $p = -2.2 \pm 0.1$. Numerical simulations predict a power law of −2.7 (Nath et al. 2020) for Milky Way-like pressures, and in nearby galaxies it has been shown that this power law can range between −2 and −4 (in HI; Bagetakos et al. 2011). A power law of −2.7 is significantly steeper than what we measure. But the size distribution of bubbles is expected to be proportional to the mechanical luminosity of OB associations (Oey & Clarke 1998), which is directly proportional to their $H\alpha$ luminosity function. For NGC 628, Santoro et al. (2022) recently measured that the HII region $H\alpha$ luminosity function follows a power law of $-1.7 \pm 0.1$. If we use this as a proxy for the OB luminosity function, it implies bubbles should follow a size distribution equal to $-2.4 \pm 0.2$ (see introduction in Nath et al. 2020). However, Nath et al. (2020) do not explore how different luminosity functions propagate in their models, so we present $-2.4 \pm 0.2$ as a lower limit on the theoretical value, −2.7. Therefore our result differs by 1–3.5σ, with our result leaning toward a shallower value compared to theoretical expectations (i.e., has more large bubbles and fewer small bubbles). This shows the size distribution potentially disagrees with theoretical expectations. Investigating the size distribution of more galaxies will help us confirm this conclusion.

If highly elliptical bubbles with aspect ratios >1.5 are instead two circular bubbles that we have misidentified, these bubbles could bias the size distribution to shallower values. To test the maximum impact this could have, we replace all bubbles with high aspect ratios (>1.5) with two circular bubbles with radii equal to half the semimajor axis and use MLE to find the power-law slope. We find 77 bubbles matching this criterion, and the new power-law slope is ~−2.3, meaning it could account for some of the difference between observations and theory. However, this represents an extreme scenario so it is very unlikely that it explains the discrepancy we potentially observe.

Counterintuitively, higher ambient pressures could explain the shallower slope we potentially observe (see Figure 14 of Nath et al. 2020). It is a consequence of the large range of evolutionary stages that—when averaged over time—create the size distribution in the first place. Higher ambient pressures decrease the time it takes for bubble expansion to stall, but, as a fraction of their expansion lifetime, bubbles powered by small stellar populations are more affected. Therefore the bubbles powered by larger stellar populations grow over a much longer timescale relative to the quickly stalled, smaller bubbles. As a result, at any given time, the growing bubbles occupy a larger fraction of the time-averaged distribution, which makes the size distribution top-heavy (Nath et al. 2020). Bubble merging can also flatten the slope as it decreases the number of smaller bubbles in favor of larger bubbles.

Given we find bubbles nest together within the shells of larger bubbles, with 31% of bubbles overlapping with at least one smaller bubble, if there is a discrepancy between theory and observations our results are consistent with bubble merging reducing the power-law slope.<fnref n="48" /> Bubble merging is expected in

<footnote n="48">⁴⁸ Although if we extrapolate the power-law index from Figure 14 of Nath et al. (2020), using the average dynamical pressure in NGC 628 of ~4.4 × 10⁻¹² dyne cm⁻² estimated in Barnes et al. (2021) as a proxy for the ambient ISM pressure, we recover that the power-law index would be ~−2.1, matching our results, though we concede there are many caveats to this estimate.</footnote>
galaxies, considering we see this process in the Milky Way (Krause et al. 2015), and without bubble merging the volume bubbles occupy can exceed the total volume of galaxies when modeling their porosity (Clarke & Oey 2002; Nath et al. 2020). Moreover, HI studies of nearby galaxies observe a bubble filling factor of only 10% (Bagetakos et al. 2011), meaning that bubbles must merge. If we assume the difference between the theoretical value is real and caused by merging, it predicts that bubble merging reduces the power law in NGC 628 by a minimum of ~0.2, providing some constraints on the impact of bubble merging for numerically derived theories. For example, the power-law index that was numerically generated in Nath et al. (2020) did not include the impact of bubble merging on the power-law index, even though they expect the power-law slope to become shallower if merging was included.

## 6.2. Expected Number of Bubbles

In the absence of merging, the number of bubbles, $N_b$, we expect to find in a galaxy is determined by three parameters: (i) the star formation rate (SFR) of the galaxy; (ii) the timescale, $t_{obs}$, over which we can observe bubbles, and (iii) the average cluster mass, $\bar{M}_*$, able to drive a bubble we can observe in this timescale (which itself depends on the form of the stellar initial mass function). Altogether, we can write

$$N_b = \frac{\text{SFR} \times t_{\text{obs}}}{\bar{M}_*}, \qquad (1)$$

as was first outlined in Clarke & Oey (2002).

While we do not know the precise value of $t_{obs}$, our identification method has an implicit time limit due to bubble shearing. Bubbles that are too sheared (and are missing an obvious powering stellar population) are not identified. Assuming that the deviation of the bubble axis ratio ($q$, which is the inverse of the aspect ratio) from unity originates entirely from the rotation curve, we can estimate the bubble expansion velocity, $v_{exp}$, as a function of shear and galactocentric radius, and hence the time over which bubbles are visibly identifiable. For the exact details of the assumptions and method used in this estimate, see Barnes et al. (2023), where it is presented in full.

We find the mean and median $q$ for elliptical bubbles (i.e., bubbles identified using elliptical apertures) are 0.8 and 0.9, respectively. Assuming the bubbles are driven by feedback, we employ a self-similar expansion model to relate the resulting expansion velocities, and bubble radii, $R_{bub}$, to bubble age, $t_{age}$, using

$$t_{\text{age}} = \eta \frac{R_{\text{bub}}}{v_{\text{exp}}}, \qquad (2)$$

and plot the resulting ages in Figure 7. The constant, $\eta$, is the self-similar scaling constant (Ostriker & McKee 1988) ranging between 2/3 and 1/4 for bubbles driven purely by adiabatic winds to SN-driven snowplow, respectively. Since we do not know the exact feedback mechanism driving the bubbles, we set $\eta$ to the model in the middle $\eta = 1/2$, which is the solution for radiative winds (Lancaster et al. 2021), and use 2/3 and 1/4 as uncertainties in Figure 7. We find the average bubble lies at 2.5 kpc from the galactic center, which, using the smallest and largest ages in Figure 7 at 2.5 kpc, predicts bubbles are visible between the ages of 7–42 Myr. We note that 42 Myr is nearly identical to the expected theoretical time limit for which the

largest bubbles grow. After 40 Myr, there are no OB stars left in the original stellar population that is driving the bubble.

The average cluster mass depends on the minimum cluster mass, $M_{min}$, that can drive bubbles above the completeness limit ($\ge 29$ pc; see Section 5.1) in 7 Myr, and the maximum number of stars within a stellar population, $N_{max}$, that drives bubbles,

$$\bar{M}_* = M_{\min} \ln(N_{\max}), \qquad (3)$$

shown in Section 2.4 of Clarke & Oey (2002). From studies of nearby galaxies, the maximum number of stars we expect within a single burst has a membership of ~14,000 stars (Larson et al. 2022). Using Figure 9 and Equation (11) in Nath et al. (2020), we predict the minimum cluster mass needed to drive a bubble $\ge 29$ pc is around $450 \pm 150 M_{\odot}$ and so $\bar{M}_* = 4300 \pm 1400 M_{\odot}$.

Finally, the SFR of NGC 628 is $1.8 M_{\odot} \text{yr}^{-1}$, and considering the footprint of our MIRI observations only cover around one-half of the total star formation (Leroy et al. 2021), we use a SFR of $0.9 M_{\odot} \text{yr}^{-1}$ for our estimate. With observable lifetimes of 7–42 Myr, NGC 628 should contain ($1000 \pm 300$)–($6200 \pm 2100$) bubbles $\ge 29$ pc, if we expect that 30% of bubbles that overlap merge together (Simpson et al. 2012; Krause et al. 2015). In our catalog we have 726 bubbles with average radii $\ge 29$ pc, which is close to the lower limit of what we expect to find. While there are many assumptions that go into such a calculation, it reinforces that the features we do find at $\ge 29$ pc are feedback-driven bubbles. It also shows that a more thorough exploration of bubbles with citizen science and machine-learning approaches could yield even more bubbles in NGC 628.

To quantify the observable timescale for feedback-driven bubbles we can instead ask: What timescale is needed to predict the number of bubbles present that match the number we find in a full bubble catalog (using citizen science, and machine-learning approaches in future work)? Such a timescale would constrain the average energy injection into bubble shells and would constrain gas recycling times, both of which are needed to quantify mixing processes between gas and the recently enriched material provided by the SN exploding in the bubbles.

### 6.2.1. Observable Bubbles in the Remaining 18 PHANGS Targets

Assuming the number of bubbles we find represents the complete number of observable bubbles in NGC 628, we can predict the total number of bubbles we should detect across the 19 PHANGS–JWST galaxies. If 787 bubbles (independently found by H.K. in catalog C; see Appendix B) represents a lower bound, and 1694 an upper bound, we find the number of bubbles found per SFR is ~900–1900 ($M_{\odot} \text{yr}^{-1}$)$^{-1}$. Using the SFRs given in Leroy et al. (2021), adjusted for the footprint of our JWST observations$^{49}$ the total SFR across the remaining 18 galaxies (see Lee et al. 2023 for details on the remaining 18 targets) is 42.1 $M_{\odot} \text{yr}^{-1}$. Thus, we predict a further 37,000–79,000 bubbles could be found using the 18 additional PHANGS–JWST galaxies. NGC 628 is one of the closest and most face-on galaxies in the PHANGS–JWST sample, and therefore has optimal parameters for identifying bubbles. Still, this represents an order of magnitude improvement on the number of bubbles we can use to investigate stellar-feedback

<footnote>$^{49}$ Footprints estimated using the WISE3 correction versus ALMA from Leroy et al. (2021), which are approximately the same size as the JWST footprints.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Bubble age as a function of galactocentric distance. The dashed blue and solid orange lines are the bubble ages for mean and median bubble ellipticity, respectively.</figcaption></figure>

mechanisms and its connection to the gas, even when compared to bubbles identified in the Milky Way (5106 MIR bubbles found in Simpson et al. 2012).

## 7. Conclusions

In this Letter, we present the first high-resolution (12 pc) view of feedback-driven bubbles using PAH emission in NGC 628. To find bubbles, we focus on JWST MIRI F770W observations, which contains the most shell-like features compared to the other JWST bands. We combine MIRI F770W with PHANGS-HST B-band and PHANGS-MUSE Hα observations—tracing young stellar sources and diffuse ionized emission, both of which are signposts of recent star formation and feedback—to help separate feedback-driven bubbles from dynamically driven features. We find that the three data sets together help to identify the complex, hierarchical nature of the ISM, with nested bubble structures being a key characteristic of the complex multiphase ISM.

With these three data sets, we visually identify a rich population of ~1700 bubbles with sizes spanning two orders of magnitude (6–552 pc). Locally, bubbles are highly nested, with 31% of bubbles containing at least one smaller bubble within their shell.

We find that the size distribution follows a power-law index of $p = -2.2 \pm 0.1$. While within the observational expectations outlined in HI studies, our power-law slope might be shallower than theoretical predictions derived using the HII Hα luminosity function of NGC 628. Owing to the significant number of bubbles within the shells of larger bubbles, we conclude that, if real, bubble merging best explains why the distribution is shallower compared to theoretical predictions, which do not include the impact of bubble mergers.

At large scales, we find azimuthal trends in the size distribution and the PAs of bubbles. Specifically, bubbles increase in size downstream from an arm compared to upstream. Since their size is related to age (i.e., their evolutionary state) bubbles are more evolved downstream, and therefore their evolution is linked to the spiral-arm passage. We find that the PAs are preferentially parallel or perpendicular to the spiral arms. Parallel PAs are induced when expanding bubbles shear, which stretches them parallel to the arm.

Potentially, the lower density gas downstream from the arm allows bubbles to blister, resulting in bubbles perpendicular to the arm, though further analysis is needed to confirm this.

To characterize how well we have sampled the population, we estimate the number of bubbles we would expect to identify. We find that bubble selection is primarily limited by shear as we are usually choosing bubbles that are less elliptical. We are also limited by the resolution and find we are complete for bubbles that have radii $\ge 29$ pc (where 29 pc represents our completeness limit). Given these selection criteria, we predict that we trace bubble lifetimes of 7–42 Myr and so we should find (1000 ± 300)–(6200 ± 2100) bubbles, respectively. With 726 bubbles identified with radii $\ge 29$ pc, we are just within this range.

Finally, by extrapolating our result to the remaining 18 PHANGS–JWST galaxies, we predict 37,000–79,000 more bubbles could be detected. Such a number represents an order of magnitude increase in the number of resolved bubble features in the nearby universe, which will form the basis for citizen science and machine-learning approaches. With comprehensive bubble catalogs, the population statistics explored in this Letter can provide quantitative constraints on energy injection into the ISM as well as stellar feedback and mixing timescales.

We would like to thank the referee for their constructive feedback that helped improve the quality of this work. This work was carried out as part of the PHANGS collaboration. Based on observations collected at the European Southern Observatory under ESO programmes 094.C-0623, 095.C-0473, and 098.C-0484. This work also makes use of the following ALMA data: ADS/JAO.ALMA#2012.1.00650.S, (N628/M74). It is also been carried out on observations made with the NASA/ESA/CSA JWST. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127. The observations are associated with JWST program 2107. The data presented in this paper obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute can be accessed via: 10.17909/9bdf-jn24. In addition, this research uses observations made with the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with program 15654 and can be accessed at: 10.17909/t9-r08f-dq31. It is also based on observations and archival data obtained with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. Finally, this publication uses the data from the AstroSat mission and the UVIT instrument of the Indian Space Research Organisation (ISRO), archived at the Indian Space Science Data Centre (ISSDC). This work is supported by a grant 19ASTROSA2 from the Canadian Space Agency. E.J.W. acknowledges the funding provided by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—Project-ID 138713538–SFB 881 (“The Milky Way System,” subproject P1). K.K., O.E., and F.S. gratefully acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) in the form of an Emmy Noether
Research Group (grant No. KR4598/2-1; PI: Kreckel). J.M.D.K. gratefully acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program via the ERC Starting Grant MUSTANG (grant agreement number 714907). COOL Research DAO is a Decentralized Autonomous Organization supporting research in astrophysics aimed at uncovering our cosmic origins. M.C. gratefully acknowledges funding from the DFG through an Emmy Noether Research Group (grant No. CH2137/1-1). T.G.W. acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 694343). G.A.B. acknowledges the support from ANID Basal project FB210003. M.B. acknowledges support from FONDECYT regular grant No. 1211000 and by the ANID BASAL project FB210003. E.C. acknowledges support from ANID Basal projects ACE210002 and FB210003. F.B. would like to acknowledge funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 726384/Empire) R.S.K., S.C.O.G., and M.C.S. acknowledge funding from the European Research Council via the ERC Synergy Grant "ECOGAL" (project-ID 855130), from the Deutsche Forschungsgemeinschaft (DFG) via the Collaborative Research Center "The Milky Way System" (SFB 881—funding ID 138713538—subprojects A1, B1, B2, and B8), and from the Heidelberg Cluster of Excellence (EXC 2181-390900948) "STRUCTURES," funded by the German Excellence Strategy. R.S.K. also thanks the German Ministry for Economic Affairs and Climate Action for funding in project "MAINN" (funding ID 50002206). M.Q. acknowledges support from the Spanish grant PID2019-106027GA-C44, funded by MCIN/AEI/10.13039/501100011033. E.R. acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC), funding reference number RGPIN-2022-03499. E.W.K. acknowledges support from the Smithsonian Institution as a Submillimeter Array (SMA) Fellow and the Natural Sciences and Engineering Research Council of Canada. K.G. is supported by the Australian Research Council through the Discovery Early Career Researcher Award (DECRA) Fellowship DE220100766 funded by the Australian Government. K.G. is also supported by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. A.K.L. and N.M.C. gratefully acknowledge support by grant Nos. 1653300 and 2205628 from the National Science Foundation, award JWST-GO-02107.009-A, award SOSP SOSPADA-010 from the NRAO, and by a Humboldt Research Award from the Alexander von Humboldt Foundation.

**Facilities:** HST, JWST, VLT.

**Software:** Python associated: Astropy (and its affiliated packages) (A. Collaboration et al. 2013; Astropy Collaboration et al. 2018; Astropy Collaboration 2022), GeoPandas (Jordahl et al. 2020), scikit-image (van der Walt et al. 2014), NumPy (van der Walt et al. 2011; Harris et al. 2020), Matplotlib (Hunter 2007), colorcet, descartes, Shapely, MultiColorFits (Cigan 2019), SciPy (Virtanen et al. 2020), powerlaw (Alstott et al. 2014). Other: SAOImage DS9 (Smithsonian Astrophysical Observatory 2000; Joye & Mandel 2003).

# Appendix A
## More Precise Rules for Identifying Bubbles

### A.1. Bubble Definition

We begin with our definition for a bubble. These are (note, this applies to the MIRI F770W data set):

1. Any full or partial circular feature. Partial features can be an arc of a bubble shell, or be a clumpy but round feature.

2. They can be any thickness, so both thinner and thicker bubbles are identified.

3. They can be any physical size, limited only by the resolution limit of 12 pc.

4. They can overlap and touch each other.

If we just used these criteria, we would identify empty holes, as well as bubbles, so to differentiate between them we loosely followed this set of rules in addition to the above:

1. We only identify stretched-out features as bubbles if they had obvious concentrated HST *B*-band sources within them, or required MUSE $H_α$ inside the hole or on the inner edge of the shell feature seen in MIRI (see Figure 8, where we show sheared features). A similar criterion was used when identifying large (≥300 pc) bubbles.

2. Small partial bubble features (usually <30 pc) needed either spatial HST *B* band or MUSE $H_α$ present within them to be identified.

Even with these rules, identifying all bubble features is a monumental task at 12 pc resolutions due to the hierarchical nature of the ISM. Bigger bubbles were often no longer visible when we looked at them more closely (breaking up into smaller bubbles, or were not empty in emission in MIRI F770W, making them hard to define), and we often found highly nested and complex bubble populations that were difficult to distinguish into individual bubbles (see left panel of Figure 8). We also found that different bubbles became visible after adjusting the contrast of the MIRI F770W or the HST *B*-band observations, further changing what bubbles we can detect. We note this to again emphasize that when generating our catalog, the goal was not to identify every single feature that perfectly matched these criteria but to find a representative set of bubbles at all scales bubbles are present at. A representative sample allows us to study bubble nesting at small scales, link the location and properties of bubbles to the larger scale environment, and investigate their evolution. We also note that if, for example, we added some bubbles at ~70 pc scales afterwards, we would also have to search for many more smaller bubble features to keep the catalog representative and scale unbiased.

### A.2. Detailed Method for Identifying Bubbles

Our identification flow followed this method:

1. Combine the three bands into an RGB image (red MIRI F770W using an asinh stretch to better show diffuse shell features; green MUSE $H_α$ using a log stretch since the emission spanned many orders of magnitude; blue HST *B* band using a linear stretch to focus on the point-like sources).

2. Look for bubbles, starting at a field of view that covered one-quarter of the MIRI footprint.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Close-up view of bubble features in NGC 628. From left to right, we show complex, nested bubbles where we might have mischaracterized the bubbles present; bubbles downstream from an arm with sheared features upstream; and a closer look at the sheared features shown on the previous panel. These RGB images have the same bands and scaling as in Figures 1 and 3.</figcaption></figure>

3. Fit ellipses (or circles if no strong eccentricity is visible) to the shell ridge of the bubbles present, incrementing the PA in steps of 10° to speed up the characterization process until it lined up with the shell feature. If 10° steps are insufficient, use finer increments until the ellipse traces the shell ridge.

4. After identifying bubbles and fitting their shells, increase and decrease the contrast of the MIRI F770W band, and adjust or add more bubbles if necessary after doing this.

5. Blink between the RGB image with and without HST B band to highlight any potentially missing bubbles and to check if the bubble still looked real without HST. Remove or add and adjust the bubbles that no longer match the definition of a bubble or whose shells now become clear, respectively. This step is subjective to the individual tolerance of what someone considers a bubble to be. We note that in Watkins et al. (2022), 325 superbubbles were identified both morphological and kinematically in the same set of galaxies that will be observed within PHANGS–JWST, and so our team already has expertise in identifying bubble features in lower resolution data sets.

6. Zoom in and repeat (down to a tile ~15" × 15"). If the larger scale bubbles no longer exist after zooming in, but instead break up into bubble complexes, fit the smaller features and remove the larger bubble if it is no longer needed to characterize the bubble features.

### A.3. Difficulty in Defining Bubble Features

When identifying bubbles using MIRI F770W, we exclude stretched features that we are unable to confirm as bubbles using the MUSE and HST data sets. However, we note that this is subjective and there are likely features that we missed, or that are even wrongly identified as a bubble. In Section 6 we show

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Close-up view of a bubble potentially missing from the catalog. Translucent black ellipse highlights a large (>500 pc) feature that could be a missing large bubble that is too sheared to confirm unambiguously. RGB image with the same bands and scaling as in Figures 1 and 3.</figcaption></figure>

that we are likely missing some bubbles with radii $\ge$29 pc, which we use to infer that it is unlikely that we have misidentified many bubbles $\ge$29 pc but are instead missing bubbles.

We therefore provide more close-up examples of bubbles we identify next to features that we deem to be too sheared to identify to help clarify our bubble finding process in Figure 8. In the same figure, we also show where features are ambiguous and are highly nested bubbles, and in Figure 9 we show a potentially missing large (>500 pc) bubble. These close-up views illustrate there is no simple solution to defining bubbles, which motivates using citizen science, or machine-learning techniques in the future, to build up a statistically motivated catalog.
## Appendix B
Generating Independent Catalogs

To confirm the results drawn in this letter in Section 5, two additional team members independently identified bubbles, shown in Figures 10(a) and 10(b), following the same method outlined in Section 4 and Appendix A. These catalogs are constructed using the full map area but with a less extensive search, and we use these catalogs to confirm that independent samples lead to the same conclusions (e.g., they generate the same size distribution). The purpose therefore was not to find as many bubbles as possible but to sample the bubble population present. K.H. and H.K. found 837 and 787 bubbles, respectively, and we label them as catalog B and C, with the original catalog labeled as A for the purposes of comparison in this section.

For catalog A, we found 759 out of 1694 bubbles overlap (treating each bubble as a mask) with catalog B, and from the perspective of catalog B 612 out of 837 bubbles overlapped with catalog A. The difference when switching perspective reflects where a bubble overlaps with multiple bubbles. Comparing catalogs A with C, we found 676 bubbles overlapped, and when comparing C with A we found 502 out of 787 bubbles overlapped. While there is good agreement on deciding which locations have bubbles present, there are significant differences in the size, shape, or number of bubbles identified in a specific region. For instance, in some locations catalog A characterizes the structures using smaller bubble apertures than catalogs B and C, while at the same locations catalogs B and C use as fewer but larger bubble apertures. The reverse (i.e., identifying fewer but larger bubbles in catalog A) also occurs.

To understand if our results are robust, we repeated the analysis and figures made in the results section for the two catalogs (B and C). We find they have mean and median bubble sizes of 35 and 46 pc, and 31 and 40 pc, respectively. The

higher average values compared to catalog A indicate that the peer-reviewed catalogs did not identify as many small bubble features, which might impact how reliably we can compare results drawn from the small-scale features.

We compare the overall size distributions between the three catalogs, and the power-law index derived from them in the top panels of Figures 11(a) and (b). Following the analysis in Section 5.1, we find all three catalogs turn over at ~30 pc (B: $30_{-1}^{+20}$ pc, C: $31_{-4}^{+6}$), confirming that all three samples have a similar completeness limit. When comparing their power laws, B and C have a power law equal to $-2.1_{-0.3}^{+0.2}$ and $-2.1_{-0.2}^{+0.1}$, respectively (matching catalog A), confirming that the power-law index we measure is robust. We note here that catalog B has a larger upper limit for its turnover point. When inspecting the cause, we find there are two solutions for the turnover point, one at ~30 pc and one at ~50 pc. This indicates there is a small break in the power law at ~50 pc, with the size distribution following a slightly steeper power law for values >50 pc. We do not think this is physical, but instead implies catalog B has a size bias.

In the middle panels of Figures 11(a) and (b), we find bubble nesting is 2.9 and 2.7 times more frequent than for a random distribution in catalogs B and C, strongly reinforcing this result. While the number of overlaps found is smaller than in catalog A (at 3.2), we again explain the difference by the lower number of small bubbles identified in B and C.

We next check the difference in bubble properties downstream and upstream from the arms in the bottom panels of Figures 11(a) and (b). Visually, the size distributions for B and C show a strong high-end tail and the lower number of small bubbles identified results in similar low-end tails between the upstream and downstream bubbles. As a result, the Mann-Whitney test returns a 4 and 40% chance (for B and C, respectively) that the two distributions sample the same underlying population. Catalog B still allows us to reject this

<figure><img src="image_12.png" /><figcaption>Figure 10. Left: same as Figure 3 for catalog B. Right: same as Figure 3 for catalog C.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 11. Left: same as Figure 3 for catalog B. Right: same as Figure 3 for catalog C.</figcaption></figure>

null hypothesis at low confidence, but at 40%, we are unable to do the same for C. Though, we note that both B and C find more bubbles in downstream (432 versus 362, 371 versus 316 downstream and upstream for B and C, respectively), which we would expect if more bubbles are smaller upstream (given that smaller bubbles are less represented in the two catalogs). Therefore we leave our initial conclusion unchanged from Section 5 given that we can explain why catalogs B and C show a less pronounced difference.

When comparing $\Delta\theta$ for catalog B in Figure 12(a), we find similar peaks and troughs with catalog A except catalog B does not peak at $\pm 90^\circ$ and it fails the Rayleigh test, with a $p$-value of 6.79%. Another difference is the number of bubbles with $AR > 1.15$ in catalog B, with 60 bubbles. In general catalog B uses more circular apertures and lower ARs to define elliptical bubbles. The poorer sampling likely causes the missing peak at $\pm 90^\circ$ and the failed Rayleigh test.

Catalog C (Figure 12(b)) also differs from Catalog A. We still find an excess of bubbles with $\Delta\theta < 0$, and it also passes the Rayleigh test with a $p$-value of 0.09%. However, the peak at $\pm 90^\circ$ is not present and the trough at $30^\circ$ is no longer a dominant feature. Instead, we see a peak at $45^\circ$. When inspecting the cause, we find the peak at $45^\circ$ comes from bubbles in the outer arms (the bottom-left and top-right corners near Arms 2 and 3, respectively, as defined in Figure 3). In fact,

catalog C finds most bubbles around Arm 3, whereas the other two catalogs find the most in Arm 4 (see Figures 3, 10(a) and (b)). If we remake Figure 12(b) excluding bubbles from Arms 2 and 3, it now matches catalog B. The outer bubbles are likely outside corotation defined at $4.4 \pm 2.0$ kpc (Williams et al. 2021), and therefore their $\Delta\theta$ might not correlate with the spiral arm passage.

In any case, the results from catalogs B and C inform us that we need much more complete bubble samples when analyzing trends in PAs. In future work, we will therefore undertake a more careful analysis when exploring these trends. We also note that since we do not find a peak at $\pm 90^\circ$, we do not compare perpendicular $\Delta\theta$ downstream versus upstream, and so we remove this panel from Figures 12(a) and (b).

Finally, we plot the radial trend with bubble size for catalogs B and C in Figures 13(a) and (b). Kendall's $\tau$ is 0.1 for both distributions, with $p$-values of 0.01% and 0.42% for B and C, respectively, matching A. These indicate a weak correlation exists. We see the high-end tails of the distributions weakly increase as a function of galactocentric radius, but the number statistics at distances >4 pc are too poor to show a definitive increase, similar to what we see in Figure 6. These results therefore reinforce the conclusions drawn previously (i.e., that the increase in bubble radius as a function of galactocentric
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Left: same as Figure 3 for catalog B but with the bottom panel removed. Right: same as Figure 3 for catalog C but again with the bottom panel removed.</figcaption></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 13.** Left: same as Figure 3 for catalog B. Right: same as Figure 3 for catalog C. We note here that the color of the distribution between 0 and 1 kpc is 31 kpc⁻² for catalog C.</figcaption></figure>

distance is weak, indicating that the process responsible for increasing the radii is not dominant).

## ORCID iDs

Elizabeth J. Watkins https://orcid.org/0000-0002-7365-5791

Ashley T. Barnes https://orcid.org/0000-0003-0410-4504

Kiana Henny https://orcid.org/0000-0001-7448-1749

Hwihyun Kim https://orcid.org/0000-0003-4770-688X

Kathryn Kreckel https://orcid.org/0000-0001-6551-3091

Sharon E. Meidt https://orcid.org/0000-0002-6118-4048

Ralf S. Klessen https://orcid.org/0000-0002-0560-3172

Simon C. O. Glover https://orcid.org/0000-0001-6708-1317

Thomas G. Williams https://orcid.org/0000-0002-0012-2142

Benjamin W. Keller https://orcid.org/0000-0002-9642-7193

Adam K. Leroy https://orcid.org/0000-0002-2545-1700

Erik Rosolowsky https://orcid.org/0000-0002-5204-2259

Janice C. Lee https://orcid.org/0000-0003-0946-6176

Gagandeep S. Anand https://orcid.org/0000-0002-5259-2314

Francesco Belfiore https://orcid.org/0000-0002-2545-5752

Frank Bigiel https://orcid.org/0000-0003-0166-9745

Guillermo A. Blanc https://orcid.org/0000-0003-4218-3944

Médéric Boquien https://orcid.org/0000-0003-0946-6176

Yixian Cao https://orcid.org/0000-0001-5301-1326

Rupali Chandar https://orcid.org/0000-0003-0085-4623

Ness Mayker Chen https://orcid.org/0000-0002-5993-6685

Mélanie Chevance https://orcid.org/0000-0002-5635-5180

Enrico Congiu https://orcid.org/0000-0002-8549-4083

Daniel A. Dale https://orcid.org/0000-0002-5782-9093

Sinan Deger https://orcid.org/0000-0003-1943-723X

Oleg V. Egorov https://orcid.org/0000-0002-4755-118X

Eric Emsellem https://orcid.org/0000-0002-6155-7166

Christopher M. Faesi https://orcid.org/0000-0001-5310-467X

Kathryn Grasha https://orcid.org/0000-0002-3247-5321

Brent Groves https://orcid.org/0000-0002-9768-0246

Hamid Hassani https://orcid.org/0000-0002-8806-6308
Jonathan D. Henshaw https://orcid.org/0000-0001-9656-7682

Annie Hughes https://orcid.org/0000-0002-9181-1161

Sarah Jeffreson https://orcid.org/0000-0002-4232-0200

María J. Jiménez-Donaire https://orcid.org/0000-0002-9165-8080

Eric W. Koch https://orcid.org/0000-0001-9605-780X

J. M. Diederik Kruijssen https://orcid.org/0000-0002-8804-0212

Kirsten L. Larson https://orcid.org/0000-0003-3917-6460

Daizhong Liu https://orcid.org/0000-0001-9773-7479

Laura A. Lopez https://orcid.org/0000-0002-1790-3148

Ismael Pessa https://orcid.org/0000-0002-0873-5744

Jérôme Pety https://orcid.org/0000-0003-3061-6546

Miguel Querejeta https://orcid.org/0000-0002-0472-1011

Toshiki Saito https://orcid.org/0000-0001-9016-2641

Karin Sandstrom https://orcid.org/0000-0002-4378-8534

Fabian Scheuermann https://orcid.org/0000-0003-2707-4678

Eva Schinnerer https://orcid.org/0000-0002-3933-7677

Mattia C. Sormani https://orcid.org/0000-0001-6113-6241

Sophia K. Stuber https://orcid.org/0000-0002-9333-387X

David A. Thilker https://orcid.org/0000-0002-8528-7340

Antonio Usero https://orcid.org/0000-0003-1242-505X

Bradley C. Whitmore https://orcid.org/0000-0002-3784-7032

## References

Alstott, J., Bullmore, E., Plenz, D., et al. 2014, PLoSO, 9, e85777

Anand, G. S., Lee, J. C., Van Dyk, S. D., et al. 2021, MNRAS, 501, 3621

Anderson, L. D., Bania, T. M., Balser, D. S., et al. 2014, ApJS, 212, 1

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167

Bagetakos, I., Brinks, E., Walter, F., et al. 2011, AJ, 141, 23

Barnes, A., Watkins, E. J., Meidt, S. E., et al. 2023, ApJL, 944, L22

Barnes, A. T., Longmore, S. N., Dale, J. E., et al. 2020, MNRAS, 498, 4906

Barnes, A. T., Glover, S. C. O., Kreckel, K., et al. 2021, MNRAS, 508, 5362

Beaumont, C. N., Goodman, A. A., Kendrew, S., Williams, J. P., & Simpson, R. 2014, ApJS, 214, 3

Belfiore, F., Santoro, F., Groves, B., et al. 2022, A&A, 659, A26

Calzetti, D., Lee, J. C., Sabbi, E., et al. 2015, AJ, 149, 51

Churchwell, E., Povich, M. S., Allen, D., et al. 2006, ApJ, 649, 759

Cigan, P. 2019, ascl soft, record ascl:1909.002 v2.1.0

Clarke, C., & Oey, M. S. 2002, MNRAS, 337, 1299

A. Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Collischon, C., Sasaki, M., Mecke, K., Points, S. D., & Klatt, M. A. 2021, A&A, 653, A16

Dale, D. A., Cohen, S. A., Johnson, L. C., et al. 2009, ApJ, 703, 517

Dobbs, C. L., & Pringle, J. E. 2013, MNRAS, 432, 653

Ehlerova, S., & Palous, J. 1996, A&A, 313, 478

Emsellem, E., Schinnerer, E., Santoro, F., et al. 2022, A&A, 659, A191

Faucher-Giguère, C.-A., Quataert, E., Hopkins, P. F., et al. 2013, MNRAS, 433, 1970

Federrath, C. 2015, MNRAS, 450, 4035

Fielding, D., Quataert, E., & Martizzi, D. 2018, MNRAS, 481, 3325

Fisher, D. B., Bolatto, A. D., White, H., et al. 2019, ApJ, 870, 46

Grasha, K., Calzetti, D., Adamo, A., et al. 2015, ApJ, 815, 93

Grudić, M. Y., Hopkins, P. F., Lee, E. J., et al. 2019, MNRAS, 488, 1501

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Hopkins, P. F., Kereš, D., Oñorbe, J., et al. 2014, MNRAS, 445, 581

Hunter, J. D. 2007, CSE, 9, 90

Jayasinghe, T., Dixon, D., Povich, M. S., et al. 2019, MNRAS, 488, 1141

Jordahl, K., den Bossche, J. V., Fleischmann, M., et al. 2020, Zenodo, Geopandas/Geopandas: v0.8.1, doi:10.5281/zenodo.3946761

Joye, W. A., & Mandel, E. 2003, in ASP Conference Series, Vol. 295, Astronomical Data Analysis Software and Systems XII, ed. H. E. Payne, R. I. Jedrzejewski, & R. N. Hook (San Francisco, CA: ASP), 489

Keller, B. W., Kruijssen, J. M. D., & Chevance, M. 2022, MNRAS, 514, 5355

Kim, J., Chevance, M., Kruijssen, J. M. D., et al. 2023, ApJL, 944, L20

Kim, J.-G., Ostriker, E. C., & Filippova, N. 2021, ApJ, 911, 128

Kim, W.-T., Ostriker, E. C., & Stone, J. M. 2002, ApJ, 581, 1080

Krause, M. G. H., Diehl, R., Bagetakos, Y., et al. 2015, A&A, 578, A113

Kreckel, K., Faesi, C., Kruijssen, J. M. D., et al. 2018, ApJ, 863, L21

Krumholz, M. R., & Burkhart, B. 2016, MNRAS, 458, 1671

Lancaster, L., Ostriker, E. C., Kim, J.-G., & Kim, C.-G. 2021, ApJ, 914, 89

Larson, K. L., Lee, J. C., Thilker, D. A., et al. 2022, arXiv:2212.11425

Lee, J., Sandstrom, K. M., Leroy, A. K., et al. 2023, ApJL, 944, L17

Lee, J. C., Whitmore, B. C., Thilker, D. A., et al. 2022, ApJS, 258, 10

Leroy, A., Bolatto, A. D., Sandstrom, K. M., et al. 2023, ApJL, 944, L10

Leroy, A. K., Schinnerer, E., Hughes, A., et al. 2021, ApJS, 257, 43

Lim, W., De Buizer, J. M., & Radomski, J. T. 2020, ApJ, 888, 98

Marble, A. R., Engelbracht, C. W., van Zee, L., et al. 2010, ApJ, 715, 506

McQuinn, K. W., Skillman, E. D., Dolphin, A. E., Berg, D., & Kennicutt, R. 2017, AJ, 154, 51

Nath, B. B., Das, P., Oey, M. S., et al. 2020, MNRAS, 493, 1034

Ochsendorf, B. B., Brown, A. G. A., Bally, J., & Tielens, A. G. G. M. 2015, ApJ, 808, 111

Oey, M. S., & Clarke, C. J. 1998, AJ, 115, 1543

Olivier, G. M., Lopez, L. A., Rosen, A. L., et al. 2021, ApJ, 908, 68

Ostriker, E. C., McKee, C. F., & Leroy, A. K. 2010, ApJ, 721, 975

Ostriker, J. P., & McKee, C. F. 1988, RvMP, 60, 1

Palous, J., Franco, J., & Tenorio-Tagle, G. 1990, A&A, 227, 175

Peeters, E., Spoon, H. W. W., & Tielens, A. G. G. M. 2004, ApJ, 613, 986

Querejeta, M., Schinnerer, E., Meidt, S., et al. 2021, A&A, 656, A133

Reynolds, R. J., & Ogden, P. M. 1979, ApJ, 229, 942

Rodriguez, J., Lee, J. C., Bradley, B. C., et al. 2023, ApJL, 944, L26

Sanchez, S. F., Rosales-Ortega, F. F., Kennicutt, R. C., et al. 2011, MNRAS, 410, 313

Sandstrom, K., Koch, E. W., Leroy, A. K., et al. 2023, ApJL, 944, L8

Santoro, F., Kreckel, K., Belfiore, F., et al. 2022, A&A, 658, A188

Shabani, F., Grebel, E. K., Pasquali, A., et al. 2018, MNRAS, 478, 3590

Silburt, A., Ali-Dib, M., Zhu, C., et al. 2019, Icar, 317, 27

Simpson, R. J., Povich, M. S., Kendrew, S., et al. 2012, MNRAS, 424, 2442

Smith, M. C., Bryan, G. L., Somerville, R. S., et al. 2021, MNRAS, 506, 3882

Smithsonian Astrophysical Observatory 2000, ascl soft, ascl:0003.002

Sun, J., Leroy, A. K., Schinnerer, E., et al. 2020, ApJ, 901, L8

Thilker, D., Lee, J. C., Deger, S., et al. 2023, ApJL, 944, L13

Thilker, D. A., Braun, R., & Walterbos, R. M. 1998, A&A, 332, 429

Tsujimoto, S., Oka, T., Takekawa, S., et al. 2021, ApJ, 910, 61

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, CSE, 13, 22

van der Walt, S., Schonberger, J. L., Nunez-Iglesias, J., et al. 2014, PeerJ, 2, e453

Van Oort, C. M., Xu, D., Offner, S. S. R., & Gutermuth, R. A. 2019, ApJ, 880, 83

La Vigne, M. A., Vogel, S. N., & Ostriker, E. C. 2006, ApJ, 650, 818

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Walter, F., Brinks, E., de Blok, W. J. G., et al. 2008, AJ, 136, 2563

Watkins, E. J., Kreckel, K., Groves, B., et al. 2022, submitted

Williams, T. G., Sun, J., Barnes, A. T., et al. 2022, ApJL, 941, L27

Williams, T. G., Schinnerer, E., Emsellem, E., et al. 2021, AJ, 161, 185

Zucker, C., Goodman, A. A., Alves, J., et al. 2022, Natur, 601, 334