PHANGS–JWST First Results: Interstellar Medium Structure on the Turbulent Jeans Scale in Four Disk Galaxies Observed by JWST and the Atacama Large Millimeter/submillimeter Array

In this paper we study the first four targets observed to date as part of the PHANGS–JWST Cycle 1 treasury project (ID 02107): IC 5332, NGC 628, NGC 1365, and NGC 7496 (Lee et al. 2023). IC 5332 is a flocculent dwarf spiral, NGC 628 is a moderate-mass, grand design spiral, NGC 1365 is a massive, strongly barred spiral galaxy, and NGC 7496 is a moderate-mass, barred spiral (Lee et al. 2023).

Our focus for this pilot study is on MIRI F1130W imaging (036 FWHM resolution) where the sensitivity to low column densities is highest (Leroy et al. 2023b). Images show exceptional sensitivity to spatially extended filamentary features of interest in this work (see Figures 1 and 2). This band also shows a tight empirical correlation with the CO emission in these four targets which may indicate that the PAH emission in this band is a sensitive tracer of multiple gas phases (Leroy et al. 2023a). These data were processed together with the suite of Near Infrared Camera (F200W, F300M, F335M, and F360M) and MIRI (F770W, F1000W, F1130W, and F2100W) imaging obtained for each target (see Lee et al. 2023). The field of view in each case overlaps existing coverage from the Hubble Space Telescope (HST; Lee et al. 2022), Very Large Telescope-Multi Unit Spectroscopic Explorer (MUSE; Emsellem et al. 2022), and ALMA (Leroy et al. 2021a) assembled by the PHANGS survey.

Table 1 lists a number of relevant properties for the four target galaxies (see also Lee et al. 2023). We note here that, at the distances of our targets, 036 corresponds to 16–34 pc.

We use PHANGS–ALMA CO(2–1) (Leroy et al. 2021b, 2021a) as our primary tracer of the dynamics of the gas, its global rotation, turbulent velocity dispersion, and mass surface density. These are necessary inputs for estimating the Jeans length and the Toomre length (see Section 3.2) throughout the gas disk. Although it may be a promising area in the future, CO is not employed in this work as a tracer of the properties of individual dusty filaments.

In three of the four targets, the gas disk covered by the MIRI imaging is dominated by the molecular gas traced by CO. In NGC 628, the azimuthally averaged molecular gas mass fraction is >70% over the MIRI field of view using the atomic gas traced by resolved H i observations available for NGC 628 (see Sun et al. 2022). A similar fraction of the inner ISM is molecular in NGC 7496, and this raises to almost 90% across the targeted field of view in NGC 1365, adopting a representative H i surface density of 10 M_⊙ pc⁻² (Bigiel & Blitz 2012). For these three targets, we can thus directly leverage the gas properties measured by PHANGS–ALMA. In the case of the molecule-poor galaxy IC 5332 we estimate that the ISM is roughly 10% molecular, and thus we supplement CO constraints with additional information, as described more below.

For NGC 628, NGC 1365, and NGC 7496, we adopt the rotation curve measured from the CO velocity field by Lang et al. (2020) as our estimate of the circular velocity, V_c . For IC 5332, we use the semiempirical rotation-curve model implied by this galaxy's stellar mass (see Leroy et al. 2021a), following Meidt et al. (2018). Each rotation curve is interpolated onto a finer radial grid so that each pixel (and resolution element) in the JWST map has an estimate of the circular velocity, V_c , associated with it. This approach may lead to misestimation of V_c where strong noncircular motions are present and contribute to the measured rotation curve (within the barred regions of NGC 1365 and NGC 7496, for example). The result is that the Toomre length may be underestimated at these radii (see Section 3.2).

For the gas properties (velocity dispersion and surface density), where possible we adopt the values tabulated for the molecular gas in the PHANGS mega-tables recently presented by Sun et al. (2022). Our interest is in approximating the Jeans and Toomre lengths in the gas initially, prior to fragmentation. We therefore adopt the "large-scale" area-weighted average molecular gas surface density, measured in 1 kpc hexagonal apertures, for NGC 628, NGC 1365, and NGC 7496. In a given hex aperture, these values are calculated from the total flux traced in the 12m+7m+tp CO(2–1) 150 pc-scale integrated intensity (moment-0) map within the hex boundary, divided by the area of the aperture. (See Sun et al. 2022, for the adopted conversion factor α_CO and other details.) For IC 5332, the flux recovery is considerably higher in the 7m+tp CO(2–1) cube (compared to the 12m+7m+tp cube), and so we construct our "large-scale" area-weighted molecular surface densities from the 7m+tp moment-0 map convolved to 1.5 kpc resolution, sampling the convolved map at the center of each hex aperture. We then apply an α_CO appropriate for each aperture following the prescription of Sun et al. (2022).

As a best approximation for the initial velocity dispersion in the cold gas, where possible we select the average cloud-scale velocity dispersion in each hexagonal aperture (as opposed to a larger scale measure, which would likely reflect unresolved galactic motion). This assumes that principal fragmentation happens at/near the cloud scale rather than above it; any scale-dependent evolution in turbulent velocity dispersion is assumed to occur below the cloud scale, e.g., where it can lead to secondary fragmentation and eventually lead to star formation. An alternative choice would be to adopt a single value representative of the neutral gas as a whole, e.g., σ_gas ≈ 11 km s⁻¹, following Leroy et al. (2008) and others. We comment on the impact of this choice later in Section 4.

As our fiducial estimate of the molecular gas velocity dispersion in each aperture, we choose the average velocity dispersion measured specifically on 150 pc scales from the 12m+7m+tp CO(2–1) cubes for NGC 628, NGC 1365, and NGC 7496 (see Sun et al. 2022 for details). This allows us to consider the same physical scale in all galaxies at the highest resolution possible (given the distance of the most distant target studied here). For IC 5332, we expect the velocity dispersion measured from the 7m+tp moment-2 map at 1.5 kpc resolution to contain unresolved bulk motion. We therefore select the minimum velocity dispersion σ = 1.25 km s⁻¹ as a realistic upper bound on the velocity dispersion in the molecular gas at all locations in the map of IC 5332.

In all four targets, we use these measurements of the molecular gas surface density and velocity dispersion as our estimates of total gas surface density and velocity dispersion, but also explore the addition of an atomic gas component in the case of the molecule-poor galaxy IC 5332 (Leroy et al. 2023a), where it is most relevant. In this case, the total gas surface density is calculated by combining the molecular surface density with an atomic gas surface density that is assumed to exponentially decline from a maximum value of 10 M_⊙ pc⁻² with a scale length of 3.6 kpc. The integrated mass of atomic gas then equals the value estimated for this galaxy (Leroy et al. 2021a). When the atomic gas component is included, the gas velocity dispersion is assigned a value of 7 km s⁻¹ (e.g., Leroy et al. 2008).

Supplementing this set of measured gas properties, in each hexagonal 1 kpc mega-table aperture we also use the tabulated stellar mass surface density and the stellar scale height, z₀, assigned uniquely to each galaxy according to an empirical scaling relation between z₀ and disk scale length (see Sun et al. 2020, and references therein).

3.1.1. Identifying Filaments

3.1.2. Measuring Characteristic Spacings

We wish to compare the measured filament spacings with basic predictions for the characteristic fragmentation scale in rotating gas disks. The well-known Toomre length (Toomre 1964),

$$\begin{eqnarray}&&{\lambda }_{T}=\displaystyle \frac{4{\pi }^{2}G{\rm{\Sigma }}}{{\kappa }^{2}},\end{eqnarray} \tag{ 1 }$$

is the smallest scale stabilized by rotation in a two-dimensional (2D) rotating disk, and is typically hundreds of parsecs to a few kiloparsecs. We determine λ_T at the locations of the mega-table hexagonal apertures (see Section 2.2) for each galaxy using the observed gas surface density and an estimate for the radial epicyclic frequency,

$$\begin{eqnarray}&&\kappa ={\left(4{{\rm{\Omega }}}^{2}+R\displaystyle \frac{d{{\rm{\Omega }}}^{2}}{{dR}}\right)}^{1/2},\end{eqnarray} \tag{ 2 }$$

derived from an analytical fit to the observed rotation curve (Lang et al. 2020). Here, Ω = V_c /R, with V_c being the circular velocity (see Section 2.2).

We also compare with the most unstable 2D wavelength,

$$\begin{eqnarray}&&{\lambda }_{2{\rm{D}}}=\displaystyle \frac{2{\sigma }^{2}}{G{\rm{\Sigma }}}=2{\lambda }_{{\rm{J}},2{\rm{D}}},\end{eqnarray} \tag{ 3 }$$

which marks where the growth of instabilities is most rapid, positioned above the 2D Jeans length λ_J,2D (where growth is infinitely slow). In this expression, σ is the turbulent plus thermal velocity dispersion in the gas (following Chandrasekhar 1951).

Fragmentation occurs at λ_2D, subject to the condition that the Toomre parameter

$$\begin{eqnarray}&&Q=\displaystyle \frac{\sigma \kappa }{\pi G{\rm{\Sigma }}}={\left(\displaystyle \frac{{\lambda }_{2{\rm{D}}}}{2{\lambda }_{T}}\right)}^{1/2}\leqslant 1.\end{eqnarray} \tag{ 4 }$$

When the disk is stable and Q ≳ 1, estimates for λ_2D exceed λ_T . Since rotation entirely prevents growth at and above λ_T , this scenario underlines the 2D stability predicted by Q ≳ 1.

Accounting for the three-dimensional (3D) nature of gas disks, the turbulent Jeans length (Jeans 1902),

$$\begin{eqnarray}&&{\lambda }_{J}=\displaystyle \frac{\sigma {\pi }^{1/2}}{{\left(G\rho \right)}^{1/2}}={\lambda }_{2{\rm{D}}}\displaystyle \frac{{f}_{g}^{1/2}}{\pi },\end{eqnarray} \tag{ 5 }$$

is predicted to be the more natural scale for fragmentation at the midplane in rotating gas disks (Meidt 2022). The relation to λ_2D in this expression is written in terms of the gas fraction f_g = ρ/(ρ + ρ_b ), where ρ_b is the background density. Three-dimensional instability can occur on this scale even when rotation is predicted to lead to 2D stability with Q > 1. The presence of rotation slightly lengthens gravitationally unstable fragments (Meidt 2022), but the Jeans length serves as a good lower bound on the 3D fragmentation scale.

In the molecular disks of nearby galaxies, the turbulent Jeans length is typically on the order of tens to hundreds of parsecs, systematically smaller than λ_2D. In some scenarios, λ_J approaches or exceeds λ_T , i.e., when the Toomre criterion is not satisfied (Q ≳ 1; see Equations (3) and (5)), signaling that only 3D instability may be possible.

With constraints from PHANGS–ALMA for the gas surface density and turbulent velocity dispersion, estimating λ_J throughout each galaxy requires a determination of the gas volume density and thus the vertical scale height. Where possible we do this in two ways: either assuming a uniform 100 pc scale height (i.e., the height of the MW disk; Heyer & Dame 2015) or solving for the scale height required for hydrostatic equilibrium in the presence of self-gravity and a background potential, i.e.,

$$\begin{eqnarray}&&\displaystyle \frac{{dP}}{{dz}}=-\rho \int {{\rm{\nabla }}}_{z}^{2}{\rm{\Phi }}{dz},\end{eqnarray} \tag{ 6 }$$

in terms of the gas pressure P = ρ σ² and total gravitational potential Φ that obeys Poisson's equation ∇²Φ = 4π G(ρ + ρ_b ).

Assuming that the gas is isothermal and follows a Gaussian vertical distribution $\rho \propto \exp (-{z}^{2}/(2{h}^{2}))$ with scale height h (such that ${\rm{\Sigma }}=\rho \sqrt{2\pi }h$), and making the typical assumption that the vertical terms dominate Poisson's equation, this can be written as the approximation

$$\begin{eqnarray}&&\displaystyle \frac{{\sigma }^{2}}{{h}^{2}}z=\int \left(\sqrt{8\pi }G\displaystyle \frac{{\rm{\Sigma }}}{h}+4\pi G{\rho }_{b}\right){dz},\end{eqnarray} \tag{ 7 }$$

yielding the following quadratic equation to solve for h:

$$\begin{eqnarray}&&\sqrt{8\pi }G{\rm{\Sigma }}h+4\pi G{\rho }_{b}{h}^{2}-{\sigma }^{2}=0.\end{eqnarray} \tag{ 8 }$$

For this work, the background is assumed to be dominated by the underlying stellar distribution, omitting the contribution from atomic gas and the dark-matter halo (unless otherwise specified). This may lead to overestimation of the gas scale height and Jeans length determinations, but should be suitable for the preliminary estimates presented in this work. Following Sun et al. (2022), the stellar volume density is estimated from the observed stellar mass surface density together with the galaxy's empirically assigned stellar scale height. The gas scale heights estimated for our four targets lie typically in the range 30 to 100 pc and exhibit a modest amount of radial variation. Replacing the nominal velocity dispersion with the conservative σ_gas ∼ 11 km s⁻¹ (following Leroy et al. 2008), brings this to a typical value of 100 pc on average. Using either scale height, the midplane gas volume density is calculated as $\rho ={\rm{\Sigma }}/(\sqrt{2\pi }h)$, assuming a Gaussian vertical distribution.

Traditional theoretical arguments for the development of gravitational instability in rotating disks (e.g., Toomre 1964; Goldreich & Lynden-Bell 1965) predict that rotation slows the growth of structures forming on scales λ > λ_J,2D (with slowest growth on the smallest scales) and entirely suppresses fragmentation above the Toomre length (Equation (1)). The result is a preference for fragmentation on the intermediate-scale λ_2D coinciding with the Toomre threshold Q = σ κ/(π GΣ) = 1.

In the disks of the present target sample, λ_2D lies near or above λ_T , which occurs when Q ≳ 1 (see Equation (4)). Indeed, the gas disks here (like most gas disks in nearby galaxies) already appear to be generally Toomre stable (see Figure 4; Leroy et al. 2008; Williams et al. 2022). This makes 2D fragmentation near λ_2D or any scale above λ_T highly unlikely.

In this case, the structure on scales below λ_2D but above or near to λ_T found in this work would suggest that fragmentation may be occurring in a mode in which factors are active that reduce the influence of rotation in favor of gravity. For example, in magneto-Jeans instability magnetic forces can act to counter the stabilizing influence of rotation (Elmegreen 1987; Kim & Ostriker 2001). However, this is most effective in diffuse ISM phases where the Alfvén speed is high compared to turbulent plus thermal velocity dispersion. In cold neutral gas with characteristic high densities and equilibrium pressures (Sun et al. 2020), the Alfvén speed is low and magnetic forces will mostly act as a source of stability (Kim & Ostriker 2001), disfavoring structures near λ_J,2D and lowering the stability threshold. Given the Q values estimated for these targets (Figure 4) magnetic forces would even more strongly disfavor 2D fragmentation near λ_2D.

Turbulent dissipation (Elmegreen 2011) has been proposed as an avenue for shifting (and possibly removing) the Toomre condition that can also lead to growth below λ_2D. Taking into account thickness, the stability threshold in 2D rotating gas disks with dissipation shifts to Q ∼ 2, allowing for weak instability above Q = 1. The instability in this dissipation-assisted regime is characterized by a fastest growing length that shifts to smaller scales as Q increases, though with growth rates that are typically an order of magnitude lower than when Q < 1.

Another way that gravity has been shown to gain a strong hold over rotation is in a fully 3D scenario, including the additional degrees of freedom associated with the third (vertical) direction (Meidt 2022). Accounting for the perturbed vertical pressure and gravitational forces absent from 2D treatments, 3D fragmentation can be triggered even where 2D fragmentation is suppressed, subject to a modified stability threshold:

$$\begin{eqnarray}&&{Q}_{M}=\displaystyle \frac{{\kappa }^{2}}{4\pi G\rho }=\sqrt{\displaystyle \frac{\pi }{2}}\left(\displaystyle \frac{{\lambda }_{J}}{{\lambda }_{T}}\right){f}_{g}^{1/2}\equiv 1,\end{eqnarray} \tag{ 9 }$$

or roughly Q = 2, again writing f_g = ρ/(ρ + ρ_b ). The 3D structures that form in this scenario are able to grow rapidly, as fast as Toomre instabilities, down to scales λ_J ≲ λ_2D (taking into account the smaller scale height in the presence of a background potential with density ρ_b , as in Equation (5); Meidt 2022). Even where the disk is stabilized in 2D by rotation at and above λ_T , 3D instability can occur near λ_J up to the threshold Q_M = 1, where λ_J ≈ λ_T (see Equation (9)).

From this perspective, the observed gas structures on hundreds of parsec scales suggests that disk fragmentation can indeed occur in a mode in which shear plays little role compared to pressure and gravity, as predicted (Elmegreen 1987, 2011; Kim & Ostriker 2001; Meidt 2022). Rotation could nevertheless be important for regulating the formation of structure on larger scales.

Previously, Jeans-scale fragmentation has been proposed for describing the spacing of clumps in low-shear spiral-arm filaments (Elmegreen et al. 2018). For the filaments that form a disk-wide network and originate with disk fragmentation, the next levels of fragmentation might be scale dependent and ending at the thermal Jeans length, according to dissipation and the scale dependence of turbulence (as observed in MW filaments; e.g., Mattern et al. 2018) and also sensitive to stability in the new filamentary geometry (e.g., Nagasawa 1987; Freundlich et al. 2014). The internal kinematics of these structures should be studied in order to probe how subsequent fragmentation can proceed, eventually leading to star formation. In this context we note that, although the formation of small-scale filamentary structures through disk fragmentation may not be sensitive to shear, they are not necessarily insensitive to galactic orbital motions: the growth rates of 3D fragments still depend on shear (Meidt 2022) and the turbulent velocity dispersion responsible for the gas pressure that counters self-gravity partially reflects (orbital) motion in the galactic potential in equilibrium (Meidt et al. 2018; Sun et al. 2020). Structures formed through 3D disk fragmentation may thus not be kinematically decoupled from their host galaxy.

As alluded to already in Section 4, the filamentary structure evident in the MIRI images is spaced with a characteristic scale that is 1–1.5 orders of magnitude larger than the disk scale height. Often spatially coincident with these filaments, there is an abundance of rich structure evident on smaller scales, with shell-like and clumpy morphologies (Watkins et al. 2023), that can be linked to directly to star formation (Barnes et al. 2023; Watkins et al. 2023). Indeed, the MIRI-traced filamentary web harbors 75%–80% of H ii regions and 60% of star clusters younger than 5 Myr in NGC 628 (Thilker et al. 2023).

The planned detailed studies of star formation feedback leveraging the full multiwavelength PHANGS data set will be key for understanding how feedback-driven structures form and evolve and interact with structures that are primarily dynamical in origin. We propose that the scale height is a useful reference for distinguishing gas structures that are dynamical in origin from those that are feedback driven. Theory and numerical simulations predict that many bubbles driven by supernovae feedback can stall at or below h (whether in an energy- or momentum-driven regime; Mac Low & McCray 1988; Kim et al. 2017; Fielding et al. 2018; Orr et al. 2022) and even the bubbles that do break out of the plane may expand only slightly beyond h (considering that wind momentum goes to driving a vertical outflow rather than to in-plane expansion; see Orr et al. 2022). Consistent with this expectation, the catalog of bubble candidates in NGC 628 morphologically identified in the MIRI images and visually matched with a stellar counterpart (Watkins et al. 2023) have average radii 20–30 pc, just at the scale height estimated here (see Figure 3), although they can reach as large as 550 pc. In total, 60% of the bubbles identified in NGC 628 have radii below 1h and ∼90% are within 2h.

With some exceptions (like the "Phantom Void"; Barnes et al. 2023), many of the feedback-driven shells traced by MIRI appear to be stalling at or near the scale height, supporting the idea that structures present on larger scales (with separations >100 pc) are mostly dynamical in origin.

In this light, we speculate that the structure originating with dynamical factors like disk fragmentation and shear acts as a scaffold that organizes star formation and feedback, which then shapes the large-scale structure and influences the process of disk fragmentation. Feedback acts not only to break up fragments (e.g., most clearly in NGC 7946 and NGC 628; see also Barnes et al. 2023; Watkins et al. 2023), but serves as a coordinated source of perturbations that stimulate further fragmentation. From the visual appearance of the four targets studied here, and taking the present distribution of star formation as representative of the time-averaged activity in the recent past, star formation seems to prompt fragmentation with both radial and azimuthal components (e.g., NGC 628 and NGC 7496). In contrast, ring-like fragmentation may be sustainable in the absence of star formation or other sources of azimuthal perturbations (e.g., within the bar in NGC 1365). Studies in a greater number of targets will be important for recognizing how certain types of perturbations (local feedback versus underlying large-scale patterns) act to trigger disk fragmentation with different radial and azimuthal properties.

Spatial power spectra (e.g., Elmegreen & Scalo 2004; Grisdale et al. 2017; Koch et al. 2020) can be expected to provide an indispensable view of the rich multiscale web of gas structures created by gravitational instability and feedback visible in the MIRI images. However, the large-scale structure of the galaxy and the resolution of the observations can masquerade as characteristic features in power spectra (Koch et al. 2020). The high spatial resolution of the MIRI data are crucial for testing whether there are genuine features in the power spectra associated with the scale height as previously proposed (Lazarian & Pogosyan 2000; Elmegreen et al. 2001; Dutta et al. 2009; Combes et al. 2012). Characteristic Jeans-scale filamentary structure may also leave an imprint on the shape of a power-law power spectrum, but given that the structure follows a preferred radial–azimuthal organization, this signature is perhaps best recovered with a non-Cartesian (non-sky-plane) coordinate system. Spatial variation in the Jeans length and scale height also calls for radial and/or azimuthal segmentation of the input maps (e.g., Elmegreen et al. 2001).