PHANGS–JWST First Results: A Statistical View on Bubble Evolution in 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.

NGC 628 is one of 19 galaxies being observed as part of the PHANGS ⁴⁶ –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 38 times 22 (11 kpc × 6.3 kpc) of NGC 628, and for the F770W observations they have a sensitivity of 0.11 MJy sr⁻¹ and an angular resolution of ∼025, which achieves a physical resolution of 12 pc (Lee et al. 2023).

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 (∼008 ∼ 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 09 (∼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 ∼15, 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 ¹²CO (J = 2−1) (Leroy et al. 2021) (at ∼1'') to investigate bubble shells in different bands.

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.

The bubble catalog contains 1694 bubbles, plotted in Figure 3; we tabulate the first 10 in Table 1. ⁴⁷ 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.

Table 1. Sample of 10 Bubbles Ordered by R.A. and their Defining Properties

Bubble	R.A.	Decl.	Semimajor	Semiminor	Average	PA	Arm	Distance to	Galactocentric
ID			axis (pc)	axis (pc)	Radius (pc)	(°)		arm (pc)	Radius (kpc)
1	$24^\circ 8^{\prime} 55\buildrel{\prime\prime}\over{.} 9$	$15^\circ 48^{\prime} 7\buildrel{\prime\prime}\over{.} 9$	86	86	86	0	3	145	5.22
2	$24^\circ 8^{\prime} 56\buildrel{\prime\prime}\over{.} 7$	$15^\circ 48^{\prime} 2\buildrel{\prime\prime}\over{.} 9$	24	24	24	0	3	127	5.05
3	$24^\circ 8^{\prime} 57\buildrel{\prime\prime}\over{.} 0$	$15^\circ 48^{\prime} 9\buildrel{\prime\prime}\over{.} 0$	42	42	42	0	3	138	5.21
4	$24^\circ 8^{\prime} 58\buildrel{\prime\prime}\over{.} 0$	$15^\circ 47^{\prime} 44\buildrel{\prime\prime}\over{.} 9$	25	25	25	0	3	101	4.54
5	$24^\circ 8^{\prime} 58\buildrel{\prime\prime}\over{.} 2$	$15^\circ 48^{\prime} 2\buildrel{\prime\prime}\over{.} 6$	24	24	24	0	3	114	4.98
6	$24^\circ 8^{\prime} 58\buildrel{\prime\prime}\over{.} 7$	$15^\circ 47^{\prime} 41\buildrel{\prime\prime}\over{.} 2$	33	33	33	0	3	95	4.44
7	$24^\circ 8^{\prime} 58\buildrel{\prime\prime}\over{.} 8$	$15^\circ 47^{\prime} 50\buildrel{\prime\prime}\over{.} 6$	32	31	32	0	3	93	4.64
8	$24^\circ 8^{\prime} 58\buildrel{\prime\prime}\over{.} 9$	$15^\circ 47^{\prime} 39\buildrel{\prime\prime}\over{.} 8$	40	40	40	0	3	93	4.40
9	$24^\circ 8^{\prime} 58\buildrel{\prime\prime}\over{.} 9$	$15^\circ 48^{\prime} 9\buildrel{\prime\prime}\over{.} 0$	31	31	31	0	3	122	5.14
10	$24^\circ 8^{\prime} 59\buildrel{\prime\prime}\over{.} 3$	$15^\circ 48^{\prime} 6\buildrel{\prime\prime}\over{.} 1$	209	171	196	0	3	112	5.04

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

Download table as:  ASCIITypeset image

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) ∝ 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.

Zoom In Zoom Out Reset image size

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 ± 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 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).

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 (Δθ) 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 Δθ angles are offset anticlockwise from the arm, while clockwise Δθ values are negative, resulting in a Δθ distribution between −90° and 90°.

Zoom In Zoom Out Reset image size

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 ±2.6 bubbles.

Figure 5 reveals peaks at 0° and −30° and −90°. In addition, bubbles at 30° are underrepresented (and, in general, positive Δθ are lacking). The physical meaning behind the peak and trough at ±30° is not immediately obvious. Due to our selection criteria in Section 4, the most sheared bubbles are excluded, which, if included, should have positive Δθ 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 Δθ, and average thse 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 distribution. Consequently, we can perform a Rayleigh test on Δθ, which checks if a periodic distribution is nonuniform. The statistic indicates there is a 0.07% chance that the distribution of Δθ 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 ≤ 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 Δθ 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 Δθ 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 Δθ, 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 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.

Zoom In Zoom Out Reset image size

In Section 5.1, we find the bubble size distribution is p = −2.2 ± 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α luminosity function. For NGC 628, Santoro et al. (2022) recently measured that the Hii region Hα luminosity function follows a power law of −1.7 ± 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 ± 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 ± 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. ⁴⁸ Bubble merging is expected in 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.

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, ${\overline{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

$$\begin{eqnarray}&&{N}_{{\rm{b}}}=\displaystyle \frac{\mathrm{SFR}\times {t}_{\mathrm{obs}}}{{\overline{M}}_{* }},\end{eqnarray} \tag{ 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

$$\begin{eqnarray}&&{t}_{\mathrm{age}}=\eta \displaystyle \frac{{R}_{\mathrm{bub}}}{{v}_{\exp }},\end{eqnarray} \tag{ 2 }$$

and plot the resulting ages in Figure 7. The constant, η, 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 η to the model in the middle η = 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.

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&{\overline{M}}_{* }={M}_{\min }\mathrm{ln}({N}_{\max }),\end{eqnarray} \tag{ 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 ≥29 pc is around 450 ± 150 M_⊙ and so ${\overline{M}}_{* }=4300\pm 1400$ M_⊙.

Finally, the SFR of NGC 628 is 1.8 M_⊙ yr⁻¹, 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_⊙ yr⁻¹ for our estimate. With observable lifetimes of 7–42 Myr, NGC 628 should contain (1000 ± 300)–(6200 ± 2100) bubbles ≥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 ≥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 ≥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

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 cospatial HST B band or MUSE Hα present within them to be identified.

Zoom In Zoom Out Reset image size

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.

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.
• 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.

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 that we are likely missing some bubbles with radii ≥29 pc, which we use to infer that it is unlikely that we have misidentified many bubbles ≥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.

Zoom In Zoom Out Reset image size