PHANGS–JWST First Results: Multiwavelength View of Feedback-driven Bubbles (the Phantom Voids) across NGC 628

High-mass stars (>8 M_⊙) are fundamental in driving the evolution of galaxies, due to the large amounts of energy and momentum (i.e., stellar feedback) that they inject into the interstellar medium (ISM) during their short lifetimes (e.g., Krumholz et al. 2014). At early times (<5 Myr), feedback in the pre-supernova (pre-SN) stages of high-mass stars (i.e., within H ii regions) plays a critical role in disrupting molecular clouds and forming bubbles within the ISM (e.g., Dale et al. 2012, 2013; Raskutti et al. 2016; Gatto et al. 2017; Rahner et al. 2017; Grasha et al. 2018; Kim et al. 2018; Grasha et al. 2019; Hannon et al. 2019; Rahner et al. 2019; Kannan et al. 2020; Barrera-Ballesteros et al. 2021a, 2021b; Kim et al. 2021b; Barnes et al. 2021; Jeffreson et al. 2021; McLeod et al. 2021; Chevance et al. 2022a; Barnes et al. 2022; Hannon et al. 2022). At later times, these bubbles can merge (Clarke & Oey 2002; Simpson et al. 2012; Krause et al. 2015), and feedback from supernovae (SNe) can further act to drive expansion, forming so-called (super)bubble structures with scales of tens to thousands of parsecs over timescales of tens of millions of years (e.g., McKee & Ostriker 1977; Mac Low & McCray 1988; Oey & Clarke 1997; Weisz et al. 2009a, 2009b; Bagetakos et al. 2011; Keller et al. 2014, 2015, 2016; Kruijssen et al. 2019; Chevance et al. 2020; Nath et al. 2020; Kim et al. 2021a; Chevance et al. 2022b; Kim et al. 2022b; Orr et al. 2022; Zucker et al. 2022), even driving fountains out of the galactic plane (e.g., Veilleux et al. 2005; Fraternali 2017).

With the advent of the JWST, we can now, for the first time, achieve the resolution, sensitivity, and coverage for an unprecedented view of bubble populations within galaxies beyond the Milky Way and the Local Group. Specifically, the wavelength coverage of the Mid-Infrared Instrument (MIRI) is perfectly suited to this task, as it is sensitive to several polycyclic aromatic hydrocarbon (PAH) emission features (e.g., at 7.7 μm; see Draine & Li 2007; Smith et al. 2007; Li 2020). Emission from PAHs is particularly useful in tracing the shells of feedback-driven bubbles (e.g., Pineda et al. 2022), due to (a) the increased gas densities found in swept-up shells (PAHs are generally well mixed with the gas and illuminated by the average interstellar radiation field such that they trace the gas column very sensitively; e.g., Regan et al. 2006; Leroy et al. 2013; Chown et al. 2021; Gao et al. 2022; Leroy et al. 2023), (b) the high number of ionizing photons emitted by the OB association powering the bubbles, leading to PAHs being destroyed in the photoionized interiors of the bubbles (e.g., Galliano et al. 2018; Chastenet et al. 2023a, 2023b; Egorov et al. 2023), and (c) the low optical depth from the shell interior to the edge, which will allow far-UV photons to easily heat the small dust grains (e.g., Draine & Li 2007; Draine 2011; Hensley & Draine 2021). Together, these cause the edges of bubbles to appear with high contrast against their interior PAH emission (e.g., Churchwell et al. 2006; Watson et al. 2008).

In a companion Letter in this Issue, Watkins et al. (2023) use a combination of JWST and Hubble Space Telescope (HST) observations to study the bubble population across the nearby (9.84 ± 0.63 Mpc; Anand et al. 2021a, 2021b), star-forming (1.8 ± 0.5 M_⊙ yr⁻¹), face-on (i ∼9°; Lang et al. 2020, and also see Blanc et al. 2013), massive (M_* = 10 ^10.3 M_⊙, M_{H I} = 10 ^9.7 M_⊙, M_H2 = 10 ^9.4 M_⊙; Walter et al. 2008; Querejeta et al. 2015; Leroy et al. 2019, 2021a), spiral galaxy Messier 74 (also known as NGC 628 or the Phantom Galaxy). These authors manually identify bubbles using a combination of 7.7 μm JWST-MIRI observations (Lee et al. 2023), B-band (438 nm) Physics at High Angular resolution in Nearby GalaxieS (PHANGS)-HST observations (Lee et al. 2022), and Very Large Telescope (VLT)-Multi Unit Spectroscopic Explorer (MUSE) Hα observations (Emsellem et al. 2022). For each bubble, Watkins et al. (2023) fit circular or elliptical apertures to the bubble boundaries (i.e., the shells) seen within this multiwavelength data set, and in doing so identify ∼1700 bubbles with radii ranging between 6 and 500 pc (e.g., left panel of Figure 1). These structures are clearly pervasive across the galaxy, form a complex nested structure (where smaller bubbles are preferentially located at the edges of large bubbles), and are among the most striking features in the initial JWST images (see Watkins et al. 2023 for an in-depth discussion of the nesting of these structures and their size distributions).

Zoom In Zoom Out Reset image size

In the Watkins et al. (2023) catalog, one hole stands out due to its size (over ∼1 kpc in diameter), circular shape, and strong contrast with respect to its environment (e.g., right panel of Figure 1). We refer to this impressive feature as the Phantom Void, which is the focus of this Letter (Table 1 summarizes the main properties of the Phantom Void that have been determined within this work). To help understand this structure, we also identify a nearby more compact, but still very well-defined bubble, which we call the Precursor Phantom Void and analyze in parallel (also see Table 1). To conduct our analysis, we combine all the data sets available taken as part of the PHANGS ⁴⁴ -JWST survey (Lee et al. 2023) and existing PHANGS multiwavelength observations to assemble a complete panchromatic, multiphase picture of the gas, stars, and dust in these bubbles.

2.1. PHANGS–JWST Observations

2.2. Hubble Space Telescope Observations

2.3. Ancillary Observations

3.1. A Detailed Look at Feedback Bubbles

Figure 2 presents a detailed, multiwavelength view of the Phantom Void and the Precursor Phantom Void regions. In the Phantom Void, we see that the main bubble is relatively devoid of any JWST 7.7 μm emission (values of <1 MJy sr pix⁻¹) or HST Hα emission (values of <10⁻¹⁷ erg s⁻¹ cm⁻¹ arcsec⁻²) toward the central cavity, and emission is discretely distributed around the bubble shell. We note that we use this nomenclature throughout this Letter, where "bubble" refers to the inside cavity, while the "shell" refers to the perimeter with some thickness (as labeled on Figure 2). Comparing to the broadband NUV HST and Astrosat observations, we find that many of these emission peaks in the shell are connected to compact stellar emission sources (discussed further in Section 3.3). When comparing to the high-sensitivity, but lower resolution, PHANGS–MUSE observations, we see that there is a low-level component to the diffuse Hα emission inside of the bubble. This emission could be connected to the more diffuse emission seen in broadband HST and Astrosat observations, and is consistent with ionization by an older generation of stars (with respect to those associated with compact Hα emission in the shell).

Zoom In Zoom Out Reset image size

In contrast, the Precursor Phantom Void appears to be associated with a single star-forming complex (seen in Hα and NUV emission) located at the edge of the shell. Moreover, there is clear evidence from the HST observations that the ionizing radiation is propagating through the bubble and illuminating its far edge, which appears bright in Hα emission. Low-level Hα is also seen in the higher sensitivity MUSE observations within the bubble cavity. Additionally, it appears that there are many more compact and less evacuated secondary bubbles toward the south of the region, which are associated with bright extended Hα emission.

Together, the structures seen in the Phantom Void and the Precursor Phantom Void suggest that the bubbles we are seeing in the JWST imaging are driven by multiple generations of high-mass stars, the youngest of which are located at the bubble edges (see Section 4 for a discussion of, for example, trigger star formation). These stars are then injecting their ionizing radiation into the bubble, which must have a relatively low density of hydrogen such that the radiation can propagate through the bubble and cause the diffuse Hα emission at the dusty shell boundary. Specifically, for a characteristic ionizing photon energy of 20 eV, we expect the optical depth of a 100 pc bubble with mean atomic hydrogen density n_H to be approximately τ ∼ 500 (n_H/cm⁻³), implying that n_H cannot be larger than ∼0.01 cm⁻³. ⁴⁶ That said, the Phantom Void is even larger than this, so could have an even lower density. Indeed, the Phantom Void is so large that it is even visible as a hole within comparatively low-resolution H I observations (∼500 pc beam; Walter et al. 2008), confirming a low atomic hydrogen density within its evacuated center (see Figure 2). Given that the estimated H I disk scale height of NGC 628 is ∼800 pc (Dutta et al. 2008), roughly the size of the the Phantom Void, it has likely burst out of the galaxy. Bubbles on this scale are one mechanism that could be responsible for shaping galactic-scale chimneys (Heiles 1984) and establishing gas recycling via galactic fountains (e.g., Fraternali 2017). An interesting avenue for the future would be to investigate if this cavity is filled with outflowing hot X-ray gas (e.g., with Chandra).

3.2. Bubble Expansion

We aim to measure the expansion speed of the bubbles with a number of independent methods using our multiwavelength data sets. These are summarized in the following subsections.

3.2.1. Molecular Gas

We make two estimates of the expansion velocity from the CO (2–1) data (Leroy et al. 2021a). First, we make an estimate using the line width measured inside of the bubble, which probes the motion out of the plane (i.e., a difference between faint emission at the front and back side of the bubble). Second, we make an estimate using the velocity of the emission in the shell of the bubble, which probes the motion in the direction of the galactic plane.

In Figure 3, we show the average CO (2–1) spectrum taken from inside of both bubbles, which show a weak, yet significant, signal. Here, we use a data cube that has had the local velocity field subtracted, including projected circular and noncircular motions, and, hence, the signal is centered on ∼0 km s⁻¹. We fit the emission with a Gaussian function, and measure a FWHM of 31.9 km s⁻¹ for the Phantom Void and 12.1 km s⁻¹ for the Precursor Phantom Void. We make the assumption that this line width is the result of the line-of-sight component of the approaching and receding side of the bubble shell. ⁴⁷ Assuming spherical expansion, then, ${v}_{\exp }=\mathrm{FWHM}/2\,\sim 16$ km s⁻¹ for the Phantom Void and ∼6 km s⁻¹ for the Precursor Phantom Void, where the expansion velocity here is the difference in velocity between the receding and approaching side of the bubble out of the plane.

Zoom In Zoom Out Reset image size

In Figure 4, we show the molecular gas velocity field that highlights local, noncircular flows in the galactic plane at the location of the Precursor Phantom Void. Focusing on the centroid velocity of the emission around the bubble, we find that residual velocities are negative on the left side of the shell and positive on the right side of the shell. Given the orientation, rotation direction and inclination of NGC 628, this would correspond to expansion outward from the center of the bubble (i.e., the left part of the shell is moving toward a larger galactocentric radius and the right part of the shell is moving toward a smaller galactocentric radius). Specifically, from Figure 4 we measure that the v_resid on the left part of the bubble is ∼−5 km s⁻¹, whereas on the right part of the bubble it is +3 km s⁻¹, with a few positions at +5 km s⁻¹. This would give ∣v_resid∣ ∼ 3−5 km s⁻¹, and hence we estimate that for the Phantom Void ${v}_{\exp }=| {v}_{\mathrm{resid}}| /\sin (i)$ = 20−35 km s⁻¹ (where i = 89; Lang et al. 2020; see also Blanc et al. 2013).

Zoom In Zoom Out Reset image size

Given the simplistic assumptions on the geometry of the system, we consider that these estimates of the expansion speed perpendicular and parallel to the plane of the galaxy are consistent for the Phantom Void. Hence the expansion velocity of the molecular gas is taken to be in the range ${v}_{\exp }\sim$ 5–35 km s⁻¹. As the Precursor Phantom Void is located within the spiral arm, and due to its smaller size relative to the the Phantom Void, there is more confusion in separating the bubble from the surrounding gas within the residual velocity map. Hence, the in-plane expansion speed of the Precursor Phantom Void cannot be estimated, and we take the expansion speed estimate from the out-of-plane method: ${v}_{\exp }\sim 6$ km s⁻¹.

3.2.2. Atomic Gas

To constrain the expansion speed in the neutral gas, we use the VLA H I (natural-weighted) data cube (Figure 2). As with the CO (2–1), we measure the averaged spectrum for the center of the Phantom Void bubble (Figure 3, right panel), and again fit the emission with a Gaussian function. We measure a FWHM of 36.3 km s⁻¹, which is very similar to the estimate obtained from the CO (32.4 km s⁻¹). This gives a v_exp =FWHM/2 ∼18 km s⁻¹.

For completeness, we also estimate the atomic gas expansion speed of the Precursor Phantom Void using the same method. We measure a FWHM = 26.4 km s⁻¹, giving v_exp ∼13 km s⁻¹. This estimate should, however, be taken with caution, as the beam size of the VLA i˝ is larger than the size of the Precursor Phantom Void bubble, and, hence, could include a significant contribution from gas kinematics outside of the bubble.

3.2.3. Ionized Gas

We use the MUSE Hα emission line kinematics to constrain the expansion speed in the ionized gas (Figure 5, left panels). Across the same region used to measure the molecular and atomic expansion speeds, we find that the average FWHM for the Phantom Void is 88.4 km s⁻¹ (${v}_{\exp }\,\sim 45$ km s⁻¹). ⁴⁸ Similarly, we find a FWHM for the Precursor Phantom Void of 85.1 km s⁻¹ (${v}_{\exp }\,\sim 45$ km s⁻¹). These values are significantly higher than what is measured in the neutral gas (∼6–20 km s⁻¹), even when accounting for the larger thermal contribution to the FWHM of the warmer ionized medium (e.g ∼20 km s⁻¹ at 10⁴ K), which could highlight the fact that the ionized medium is more directly connected to the source of feedback.

Zoom In Zoom Out Reset image size

Interestingly, we find systematically elevated line widths in the interior of both bubbles relative to their shells, and even measure FWHM values of more than 100 km s⁻¹ in the Phantom Void (though at these very low signal-to-noise ratios < 10 note there is a bias toward overestimating line widths). Such high values of the line width could be indicative of SN feedback contributing to increased turbulence within the center of the bubble (Egorov et al. 2023). Near the shell, we directly identify SN remnants via their broadened [S ii]6716 line emission and increased [O i]/Hα line ratios (Figure 5; J. Li et al, 2022, in preparation). The presence of SN explosions within the shell could be contributing to continued driving of the bubble expansion (Section 3.3).

3.2.4. Galactic Dynamical Constraints

In the presence of differential rotation or "rotation curve shear," expanding structures will become distorted over time (Palous et al. 1990). Here we use the ellipticity, , of shells of a given radius to place rough constraints on their expansion velocity. Assuming an isotropic bubble expansion velocity, v_exp, and bubble radius, ${R}_{{\rm{b}}}={\int }_{0}^{t}{v}_{\exp }\,{dt}$ (in the absence of shear), after bubble age t, the axis ratio (q = 1 − ) of shearing bubbles can be written

$$\begin{eqnarray}&&q={\left(1+\displaystyle \frac{{{dV}}_{{\rm{c}}}}{{dR}}{t}_{{\rm{b}}}\right)}^{-1}={\left(1+\displaystyle \frac{{{dV}}_{{\rm{c}}}}{{dR}}\displaystyle \frac{\langle {R}_{{\rm{b}}}\rangle }{\langle {v}_{\exp }\rangle }\right)}^{-1},\end{eqnarray} \tag{ 1 }$$

in terms of the rotation curve derivative dV_c/dR, the average (time-weighted) bubble radius,

$$\begin{eqnarray}&&\langle {R}_{{\rm{b}}}\rangle =\displaystyle \frac{{\int }_{0}^{t}{R}_{{\rm{b}}}{dt}}{t},\end{eqnarray} \tag{ 2 }$$

the average (time-weighted) bubble expansion velocity,

$$\begin{eqnarray}&&\langle {v}_{\exp }\rangle =\displaystyle \frac{{R}_{b}}{t}=\displaystyle \frac{{\int }_{0}^{t}{v}_{\exp }{dt}}{t},\end{eqnarray} \tag{ 3 }$$

(as would be estimated from the present-day size and age), and effective bubble age as estimated by t_b = 〈R_b〉/$\langle {v}_{\exp }\rangle$ =〈R_b〉/R_b t.

Using Equation (1), the v_exp required to generate bubbles with a given set of average sizes and ellipticities can be determined, provided that the rotation curve is known, that is,

$$\begin{eqnarray}&&\langle {v}_{\exp }\rangle =\displaystyle \frac{{{dV}}_{{\rm{c}}}}{{dR}}\langle {R}_{{\rm{b}}}\rangle {\left(\displaystyle \frac{1}{q}-1\right)}^{-1}.\end{eqnarray} \tag{ 4 }$$

Below we use this expression to obtain a basic estimate of ${v}_{\exp }$ under the assumption of negligible time evolution in the bubble expansion rate, such that the current bubble radius is a good approximation of <R_b >. Other assumptions for the time evolution of ${v}_{\exp }$ will yield different results, but this should be useful for illustration purposes. ⁴⁹

The Phantom Void, with a radius R_b ∼ 500 pc, sits at R_gal ∼ 2.5 kpc and has a low ellipticity corresponding to q ∼ 0.85. Figure 6 shows the expansion velocities that would be required for an expanding structure with these properties to retain such a high degree of circularity. Here we have adopted an analytical two-parameter model for this galaxy's rotation curve to estimate the rotation curve shear at all locations (Lang et al. 2020). To produce a structure like the Phantom Void, we estimate expansion velocities ${v}_{\exp }\sim$ 20 km s⁻¹ are required.

Zoom In Zoom Out Reset image size

3.2.5. Mass and Energy of the Shell

We use the CO (2–1) ALMA observations to estimate that there is 3.8 × 10⁷ M_⊙ of molecular gas within the shell of the Phantom Void (1.3 × 10⁶ M_⊙ inside the bubble). Here we use the same shell mask shown in Figure 2, adopt a CO (2–1)/CO (1–0) ratio of 0.61 based on direct observations (den Brok et al. 2021), and a α_CO(1−0) = 3.93 M_⊙ pc⁻² (km s⁻¹)⁻¹ as predicted by empirical scaling relations (Sun et al. 2020). In addition, we estimate that there is 5.6 × 10⁶ M_⊙ of atomic gas within the shell of the Phantom Void (7.0 × 10⁵ M_⊙ inside the bubble), when using the VLA data and the conversion presented in Bigiel et al. (2010).

When taking a representative expansion speed determined in the previous section of ${v}_{\exp }\sim$ 20 km s⁻¹, we calculate that the energy required to drive a shell with this combined atomic and molecule mass is 1.7 × 10⁵³ erg. We find the total (H I + H₂) mass inside the shell of the Precursor Phantom Void is about 10 times less than the shell (4.2 × 10⁶ M_⊙), which—when assuming an estimated expansion speed of ${v}_{\exp }\sim$ 6 km s⁻¹—implies a kinetic energy of ∼10⁵¹ erg.

A plausible energy source for these bubbles is stellar feedback (as discussed further in the following section). For example, a typical value for the amount of kinetic energy released per SN is ∼10⁵¹ erg (see, e.g., Bethe 1990; Draine 2011), and, hence, approximately 100 SNe would be required to power the Phantom Void. On the other hand, the Precursor Phantom Void could be powered by a single SN. However, it is worth noting that this number could vary significantly depending on, for example, the contribution from pre-SN feedback (e.g., in the form of winds, which over the lifetime of a high-mass star have been suggested to provide an energy contribution similar to a SN; e.g., Chevance et al. 2022a), and coupling efficiency of these feedback mechanisms with the expanding bubble (e.g., values of the amount of kinetic energy retained in the surrounding medium, and not radiated away, ranges between a few to a few ten percent; e.g., Thornton et al. 1998; Tamburro et al. 2009; Sharma et al. 2014; Kim & Ostriker 2015; Martizzi et al. 2015; Gentry et al. 2017). We estimate the mass of a stellar association required to produce these number of SNe assuming ${f}_{* \to \mathrm{SN}}/\left\langle m\right\rangle \sim 0.01$, or that one SN is produced per 100 M_⊙ of stars that fully populate an initial mass function (see Tamburro et al. 2009). Hence, we estimate that the Phantom Void should contain >10⁴ M_⊙ of stars with ages less than a few 10 Myr (i.e., the timescale for the stars to explode as SNe after they are formed), and the Precursor Phantom Void >10² M_⊙. As we show in the following section, these criteria are fulfilled.

3.3. Stellar Populations Driving Bubbles

Within the Phantom Void, we find many compact and bright concentrations of NUV bright (young) stellar point sources within the shell (see HST F275W filter in Figure 2), which spatially correspond to regions of bright Hα and 7.7 μm emission. In addition, there appears to be an overabundance (with respect to a similar-size interarm region adjacent to the bubble) of fainter point sources distributed throughout the bubble cavity, which is also seen as diffuse emission in the Astrosat NUV map.

Figure 7 (upper panel) shows the flux distribution of the JWST 7.7 μm, HST Hα, HST NUV (F275W) filter, ALMA CO, and VLA H I emission maps measured radially from the center of the Phantom Void. We find a good correspondence between the 7.7 μm and Hα emission, which both peak within the shell radius at ∼600 pc (compare to dashed circles in Figure 2). We also see that the HST NUV (F275W) emission is fainter within the bubble. Interestingly, we find that the peak in the NUV emission is at a smaller (∼100 pc) radius to the 7.7 μm and Hα emission, hinting that these clusters may have already evacuated their immediate environment of dense atomic gas and ionized gas, respectively. We also see that the neutral gas tracers (CO and H I emission) are increased within the shell. The Precursor Phantom Void, on the other hand, shows a single compact cluster of point sources in the shell (again correspondent with bright Hα and 7.7 μm emission), and very little point-source emission within the bubble.

Zoom In Zoom Out Reset image size

In Figure 2 (upper-right panel), we also overlay as colored contours the ages of the stellar associations (see Section 2.2; Deger et al. 2020; Larson et al. 2023). These indicate that the stellar associations throughout the bubble and shell are generally young (<20 Myr). We estimate the total stellar mass of these young stellar associations within the the Phantom Void region to be 10^6.0 M_⊙ (80% of this stellar mass is in the shell and 20% is in the bubble). These associations could be viable powering sources for these bubbles, either currently or, in the case of those inside the bubble of the Phantom Void, they may represent the past generation of stars that created the bubbles (Section 4).

Interestingly, several of the stellar associations are very young (<5 Myr). Both in the Phantom Void and the Precursor Phantom Void, these associations reside in the shells (as inferred from the Hα). Indeed, Figure 7 (lower panel) shows the radial distribution of ages from the association catalog across the Phantom Void, where the size of each point indicates the stellar mass. We see evidence that there is a higher number of younger and more massive associations situated within the shell, as opposed to inside of the bubble. We find that there are 10 associations with ages less than 5 Myr within the shell of the Phantom Void, which have a maximum mass of 10⁵ M_⊙ (total mass of 3.2 × 10⁵ M_⊙), and, therefore, are massive enough to harbor OB stars (i.e., strong sources of stellar feedback).

We consider whether the stellar population ages align with the previous estimates of the bubble expansion speed. To test this, we make a calculation of the expansion velocity based on the radius of the bubble (R ∼ 500 pc), and a mean age of the associations within the bubble t ∼ 10 Myr. According to the classical Weaver et al. (1977) model of expanding bubbles, ${v}_{\exp }\,\simeq \,0.6R/t$, assuming they are in an adiabatic stage of their evolution. This implies ${v}_{\exp }$ = 30 km s⁻¹, which is very close to the values estimated from our kinematic analysis (Section 3.2).

The large size, high expansion speed, and large energy required to power the the Phantom Void bubble (Section 3.2) suggests that the driving mechanism is stellar feedback, sustained by a combination of multiple mechanisms. We propose that early feedback first cleared a bubble, as we are now currently witnessing for the Precursor Phantom Void (and to some extent continue to see in the shell of the Phantom Void). Since then, SNe have been exploding within the cavity and have accelerated the shell (e.g., Keller et al. 2022). This sustained injection of energy and momentum is what is required to reach high (50 km s⁻¹) expansion speeds (i.e., seen in the ionized gas), because the continuous energy input keeps the pressure in the bubble high, whereas having only a single SN would cause the pressure to drop as the bubble expands (e.g., Mac Low & McCray 1988). Finally, the star formation and SNe that are observed in the shell could further contribute to the current and future expansion of the shell.

We have found that the Phantom Void contains a significant amount of ongoing star formation within its shell (as seen in the Hα emission and the young, massive stellar associations), and also that this shell contains a large reservoir of molecular gas (see Table 1). This could then imply several physical scenarios as to why star formation is preferentially located toward the shell, of which we provide an nonexhaustive list of examples below:

• (i)  
  Gravitationally unbound gas is collected in the shell of the bubble as it expands, which then becomes dense and then star formation is triggered (i.e., gas that would not have otherwise formed stars; e.g., Elmegreen & Lada 1977; Whitworth et al. 1994).
• (ii)  
  Gravitationally bound gas either from the parent cloud or surrounding the initial star-forming region is moved by the expansion of the bubble into the shell, which then forms stars (i.e., gas that would have otherwise formed stars but in a different location).
• (iii)  
  If star formation proceeds sequentially in a cloud/region in an inside-out fashion, then we would naturally get the oldest stars in the middle of the cloud, with younger stars near the edge.

Differentiating between these cause-and-effect scenarios for the formation of the feedback-driven bubbles and future episodes of star formation is nontrivial, and has been a very long-standing problem in the galactic star formation community. Nonetheless, below we outline methods to test this in future work.

The different cases make different predictions regarding the ages of the various young stellar clusters/associations. In scenario (i), where star formation is triggered by shell collapse, we would expect all of the young regions in the shell to have ages less than the age of the shell (and hence also less than the age of the central associations). However, this need not be the case in scenario (ii). If bound clouds are being swept up that are already capable of forming stars, then it is likely that some of them will already have stars when they get swept up, whereas others will not, and we would therefore expect to see a broader distribution of ages than in scenario (i).

Another interesting test would be to look at whether there is a population of young stellar clusters/associations just inside the shell. This is what one might expect if the expanding shell sweeps away molecular clouds that have already started forming stars, as the shock will act on the gas but not the stars, which should therefore get left behind. This would be much harder to bring about in scenario (i), as in this case we would expect that the radial outward velocity of the stars at their birth should be the same as that of their parent giant molecular clouds (i.e., they should move along with the shell, or even move out faster than the shell in the case where it is significantly decelerating). Interestingly, as seen in Figure 7, we do find that the young stars (seen in the HST NUV and in the cluster association catalog) are preferentially located on the inner edge of the shell (offset from the peak in 7.7 μm emission by 50 to 100 pc).

The Phantom Void is one of the most striking features within the new JWST data, which was not as pronounced in our existing PHANGS multiwavelength data sets. Only with the resolution and sensitivity afforded by JWST do such coherent, large-scale features become clearly apparent (Figure 1). It is then interesting to consider if phantom voids are a common feature of other galaxies.

Watkins et al. (2023) estimates that there should be ∼1900 bubbles identifiable in the PHANGS–JWST observations per 1 M_⊙ yr⁻¹ of star formation rate (Section 3.1). Hence, using their size distribution, we expect to find ∼1 with a radius ∼500 pc per 1 M_⊙ yr⁻¹ of star formation rate. This is more or less in line with the observations, given that the star formation rate contained within the mapped region is ∼1 M_⊙ yr⁻¹. Across the remaining galaxies to be mapped as part of the PHANGS–JWST treasury program (total 42 M_⊙ yr⁻¹), we could then expect to uncover many phantom voids (albeit these may be harder to identify in galaxies with higher inclinations than NGC 628, because of greater foreground confusion).

We offer a brief comparison between the new JWST observations and theoretical predictions/simulations. Figure 8 compares our 7.7 μm emission map and snapshots from two sets of galaxy simulations (with matched physical scales).

Zoom In Zoom Out Reset image size

The first is the isolated (noninteracting) galaxy simulation from Tress et al. (2020, 2021). This simulation represents a generic isolated galaxy and is not tuned to reproduce the particular properties of NGC 628. The aim of the simulation was to study molecular clouds and their associated star formation in a galactic environment. The simulation includes a live N-body stellar and dark matter component. The ISM is modeled using a time-dependent chemical network that keeps track of hydrogen and carbon chemistry, a physically motivated model for the formation of new stars using sink particles, and Type Ia and Type II SN feedback. The simulation reaches subparsec resolution in the dense regions (see Figure 3 in Tress et al. 2020) and self-consistently follows the formation of individual molecular clouds from the large-scale flow and their embedded star formation. More details about the simulation can be found in Tress et al. (2020).

The second simulation is a dedicated hydro-dynamical simulation of an isolated main-sequence star-forming galaxy, using known properties of NGC 628 to generate initial conditions (stellar and gas mass profiles, velocity profile, viewing angles). The simulation was run using Ramses (Teyssier 2002), an adaptive-mesh refinement code with a maximum sampling of 3.6 pc for the gas cells. It includes live particles for dark matter, old and new stars, atomic and molecular cooling, constant efficiency per freefall time star formation above a given gas density, with recipes for UV background, feedback from Type Ia and Type II SNe, radiative pressure, and stellar winds.

The similarity in morphology between simulations and observations is striking. We show as ellipses in the bottom two panels of Figure 8 several large-scale (∼1 kpc-sized) voids, similar to the Phantom Void, that have been identified in the simulations. ⁵⁰ In addition, we see many smaller bubbles in the simulations that are similar in size to the the Precursor Phantom Void. Indeed, simulations have been predicting a bubble-rich ISM morphology for decades (e.g., Wada 2008; Grisdale et al. 2017), and H I observations have hinted at a similar morphology (see, e.g., Figure 1 in Boomsma et al. 2008). However, H I observations have limited resolution outside of the Milky Way, and it is only with the advent of JWST that we can for the first time confirm the predictions of simulations so spectacularly. In simulations which include perturbations from satellites, such as the M51 analog shown in Pettitt et al. (2017), large-scale voids are also clear in the gas column density. The voids seen in Pettitt et al. (2017) are driven by SN feedback, and also track the edges of the spiral arms, as we see in M74. The strong similarity between observations and simulations suggests that the physical processes included in the simulations (gas self-gravity, SN feedback, and differential rotation) are at the origin of the observed morphology. Further analysis, which is out of the scope of this paper, will be necessary to quantify in more detail the statistical properties of simulations versus observations (see, e.g., Sandstrom et al.2023).

This letter demonstrates the potential for new PHANGS–JWST observations to be used to identify (morphologically) stellar-feedback-powered bubbles, and highlights the need for multiwavelength data sets to study and understand their physical properties. Here we provide a detailed case study of two regions of interest, one of which contains the most prominent bubble in the galaxy (the Phantom Void, over 1 kpc in diameter), and the other being a smaller region that may be the precursor of the larger bubble (the Precursor Phantom Void). We see that the compact Hα typically sits on the edge of the shells, highlighting that the youngest stars are forming at the boundaries of the bubbles, coincident with the youngest (∼1 Myr) and most massive (∼10⁵ M_⊙) HST stellar associations that are also found toward these regions. We also find an older generation (∼20 Myr) of stellar associations that is present within the bubble of the Phantom Void. From our kinematic analysis of the H I, H₂ (CO), and H ii gas across the Phantom Void, we infer a high expansion speed of around 15 to 50 km s⁻¹. The large size and high expansion speed of the Phantom Void suggest that the driving mechanism is sustained stellar feedback due to multiple mechanisms. We propose a scenario where early feedback first cleared a bubble (as we observe now in the Precursor Phantom Void), and since then SNe have been exploding within this cavity and have accelerated the shell (which can also be enhanced by shear). Finally, comparison to simulations shows a striking resemblance to our JWST observations, and suggest that such massive stellar-feedback-driven bubbles should be common with other galaxies.

We would like to thank the referee for their constructive feedback that helped improve the quality of this paper.

This work was carried out as part of the PHANGS collaboration. This work is based on observations made with the NASA/ESA/CSA JWST and NASA/ESA Hubble Space Telescopes.

Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) 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 for JWST and NASA contract NAS 5-26555 for HST. The JWST observations are associated with program 2107, and those from HST with program 15454. The specific observations analyzed can be accessed via 10.17909/9bdf-jn24 and 10.17909/t9-r08f-dq31.

This work also makes use of observations collected at the European Southern Observatory under ESO programmes 094.C-0623 (PI: Kreckel), 095.C-0473, 098.C-0484 (PI: Blanc), 1100.B-0651 (PHANGS–MUSE; PI: Schinnerer), as well as 094.B-0321 (MAGNUM; PI: Marconi), 099.B-0242, 0100.B-0116, 098.B-0551 (MAD; PI: Carollo) and 097.B-0640 (TIMER; PI: Gadotti). 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 center (ISSDC). This work is supported by a grant 19ASTROSA2 from the Canadian Space Agency.

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2012.1.00650.S, ADS/JAO.ALMA#2017.1.00886.L. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.

A.T.B. and 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). 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). J.P.e. acknowledges support by the DAO iSM grant No. ANR-21-CE31-0010 and by the Programme National "Physique et Chimie du Milieu Interstellaire" (PCMI) of CNRS/INSU with INC/INP, cofunded by CEA and CNES. 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. 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 No. 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). R.S.K. acknowledges 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 the project "MAINN" (funding ID 50OO2206). T.G.W. acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant No. 694343). K.K. and O.E. 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). G.A.B. acknowledges the support from ANID Basal project FB210003. S.J. is supported by Harvard University through the ITC. M.B. acknowledges support from FONDECYT regular grant 1211000 and by the ANID BASAL project FB210003. E.R. acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC), funding reference number RGPIN-2022-03499. This research was supported by the Excellence Cluster ORIGINS which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy—EXC-2094-390783311. Some of the simulations in this paper have been carried out on the computing facilities of the Computational Center for Particle and Astrophysics (C2PAP). E.E. would like to thank Alexey Krukau and Margarita Petkova for their support through C2PAP. 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 supported by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. M.Q. acknowledges support from the Spanish grant No. PID2019-106027GA-C44, funded by MCIN/AEI/10.13039/501100011033. F.R. acknowledges support from the Knut and Alice Wallenberg Foundation. C.E. acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG) Sachbeihilfe, grant No. BI1546/3-1 A.K.L. gratefully acknowledges support by grant Nos. 1653300 and 2205628 from the National Science Foundation, by award JWST-GO-02107.009-A, and by a Humboldt Research Award from the Alexander von Humboldt Foundation. J.S. acknowledges support by the Natural Sciences and Engineering Research Council of Canada (NSERC) through a Canadian Institute for Theoretical Astrophysics (CITA) National Fellowship.

Facilities: HST (Hubble Space Telescope) - , JWST - , ALMA (Atacama Large Millimeter/submillimeter Array) - , VLT-MUSE (Very Large Telescope-Multi Unit Spectroscopic Explorer) - , VLA (The Karl G. Jansky Very Large Array) - , Spitzer Space Telescope - , Astrosat (Astronomy satellite). -

Software: Astropy (Astropy Collaboration et al. 2013,2018, 2022), SAO i mageDS9 (Joye & Mandel 2003), CARTA (Comrie et al. 2021), APLpy (Robitaille & Bressert 2012; Robitaille 2019).