What Holes in the Gas Distribution of Nearly Face-on Galaxies Can Tell Us about the Host Disk Parameters: The Case of the NGC 628 Southeast Superbubble

Hereafter, we adopt a Gaussian interstellar gas distribution in the galactic disk and assume that the disk component is surrounded by a low-density galactic halo:

$$\begin{eqnarray}&&{\rho }_{\mathrm{gas}}(x,y,z)={\rho }_{0}\exp \left[-{z}^{2}/{H}_{z}^{2}\right]+{\rho }_{h},\end{eqnarray} \tag{ 1 }$$

where ρ_h = μ n_h is the gas volume density in the galactic halo, ρ_mp = ρ₀ + ρ_h is the midplane gas density, H_z is the gaseous disk scale height, and μ = 14/11m_H is the mean mass per particle in the neutral gas with 10 hydrogen atoms per each helium atom.

The components of the gravitational field are (e.g., Kuijken & Gilmore 1989; Bisnovatyi-Kogan & Silich 1995; Ehlerová & Palouš 2018)

$$\begin{eqnarray}&&{g}_{x}=-\displaystyle \frac{{V}_{{\mathrm{rot}}^{2}}}{{R}^{2}}x,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{g}_{y}=-\displaystyle \frac{{V}_{{\mathrm{rot}}^{2}}}{{R}^{2}}y,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{g}_{z}=-2\pi {Gz}\left[\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{D}}}}{\sqrt{{z}^{2}+{H}_{d}^{2}}}+2{\rho }_{h}\right],\end{eqnarray} \tag{ 4 }$$

where x, y, and z are Cartesian coordinates, G is the gravitational constant, V_rot is the rotation velocity, R² = x² + y² is the distance along the galactic plane, and H_d and Σ_D are the stellar disk scale height and surface density.

The parameters n_h , n₀, and H_z are not fixed by the available observations. They are determined by the best numerical fit to the observed gas column density distribution.

The SFH within the southeast superbubble has been carefully analyzed by Mayya et al. (2023) and Barnes et al. (2023). In order to determine the bubble-driven energy and the mass deposition rate, we approximate the excess of the star formation rate (SFR) within the bubble over that in the disk (top of Figure 13 in Mayya et al. 2023), as a sequence of N_tot instantaneous starbursts separated by ${\rm{\Delta }}t={t}_{\max }/{N}_{\mathrm{tot}}$ time intervals, where ${t}_{\max }\sim 50\,\mathrm{Myr}$ is the time passed since the superbubble starburst onset until now (see Silich et al. 2002).

One can obtain then the corresponding mechanical luminosity and mass deposition rate at any time t = iΔt, where i ≤ N_tot, by adding the mechanical luminosities and mass deposition rates from all previous mini-starbursts:

$$\begin{eqnarray}&&{L}_{b}(t)=\displaystyle \sum _{j=0}^{j\lt i}{L}_{j}({t}^{* }),\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\dot{M}}_{b}(t)=\displaystyle \sum _{j=0}^{j\lt i}{\dot{M}}_{j}({t}^{* }),\end{eqnarray} \tag{ 6 }$$

where ${t}^{* }=\left(i-j-1\right){\rm{\Delta }}t$ is the time interval between the evolutionary time t = iΔt and one of the previous mini-starbursts that occurred at time t_j = jΔt.

The mass ${\dot{M}}_{j}({t}^{* })$ and energy L_j (t^*) input rates of the mini-starbursts were obtained from the stellar population synthesis model STARBURST99 (Leitherer et al. 1999), which includes stellar winds and SNe based upon the assumption of a solar metallicity starburst and a canonical Kroupa initial mass function with lower and upper cutoffs of M_low = 0.1 M_⊙ and M_up = 100 M_⊙, respectively.

The mass of each mini-starburst was calculated as

$$\begin{eqnarray}&&{M}_{j}={\rm{SFR}}({t}_{j}){\rm{\Delta }}t.\end{eqnarray} \tag{ 7 }$$

The resulting mechanical luminosity as a function of the bubble age is presented in Figure 2. The energy input rate first grows fast and then remains almost constant until the SFR becomes negligible at the bubble age of about 50 Myr. This is a cumulative effect of the different age stellar generations detected within the superbubble volume. The total stellar mass formed during 50 Myr is ${M}_{\mathrm{tot}}={\sum }_{j=0}^{j={N}_{\mathrm{tot}}}{M}_{j}\sim {10}^{5}$ M_⊙.

Zoom In Zoom Out Reset image size

Mayya et al. (2023) found that recent star formation is not confined to just the inside of the bubble zone (see their Figure 12), instead any region of the disk in the vicinity of the bubble has also been forming stars over the last 50 Myr. Thus, the increasing inner bubble volume already contains some stars whose feedback should be considered. Therefore, we adopt as the bubble driving energy the sum of

$$\begin{eqnarray}&&L(t)={L}_{b}(t)+{L}_{d}(t)\times {S}_{b}(t),\end{eqnarray} \tag{ 8 }$$

where L_d (t) and S_b (t) are the stellar disk mechanical luminosity (red line in Figure 12 in Mayya et al. 2023) normalized to the unit surface area and the surface of the bubble cross section by the galactic plane at a given time t, respectively. The disk component contribution L_d (t) to the bubble mechanical luminosity is calculated in the same manner as L_b (t). However, the second term in Equation (8) also depends on the bubble cross section and therefore must be calculated at each time step during the simulations.

All simulations were performed with a 3D code based on the thin shell numerical scheme presented in Palous (1992), Silich (1992), Ehlerová & Palouš (2018), and Jiménez et al. (2021). There are two major assumptions in this approach. The first one is that all swept-up interstellar gas is accumulated in a thin shell behind the leading shock. The second one is that the thermal pressure P_th inside the bubble is uniform. The shell is split into several Lagrangian elements and the equations of mass, momentum, and energy conservation are solved for each Lagrangian element:

$$\begin{eqnarray}&&\displaystyle \frac{{{dM}}_{i}}{{dt}}={\rho }_{\mathrm{gas}}\left(x,y,z\right)\left({{\boldsymbol{U}}}_{i}-{\boldsymbol{V}}\right){{\boldsymbol{n}}}_{i}d{{\rm{\Sigma }}}_{i},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d\left({M}_{i}{{\boldsymbol{U}}}_{i}\right)}{{dt}}={P}_{\mathrm{th}}{{\boldsymbol{n}}}_{i}d{{\rm{\Sigma }}}_{i}+{\boldsymbol{V}}\displaystyle \frac{{{dM}}_{i}}{{dt}}+{M}_{i}{\boldsymbol{g}},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{dE}}_{\mathrm{th}}}{{dt}}={L}_{b}\left(t\right)-\displaystyle \sum _{i=1}^{{N}_{\mathrm{shell}}}{P}_{\mathrm{th}}{{\boldsymbol{U}}}_{i}{{\boldsymbol{n}}}_{i}d{\rm{\Sigma }},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{{\boldsymbol{r}}}_{i}}{{dt}}={{\boldsymbol{U}}}_{i},\end{eqnarray} \tag{ 12 }$$

where r _i , M_i , and dΣ_i are the ith Lagrangian element position vector, mass, and the surface area, ρ_gas is the ambient gas density, U _i , and V are the Lagrangian element and the ambient gas velocities in the rest frame, g is the gravitational acceleration, n _i is the unit vector that is normal to the ith Lagrangian element surface, Ω_b is the bubble volume, E_th is the bubble thermal energy, and γ = 5/3 is the ratio of the specific heats. N_shell is the total number of Lagrangian elements. All these variables are calculated for each Lagrangian element at each time step. The calculations presented here are performed with N_shell ∼1500 Lagrangian elements.

The projected gas column density distributions were evaluated by making use of each Lagrangian element position r _i , orientation (determined by the unit vector that is normal to the Lagrangian element surface), mass M_i , and surface area dΣ_i .

The thermal pressure inside the bubble is

$$\begin{eqnarray}&&{P}_{\mathrm{th}}=\left(\gamma -1\right)\displaystyle \frac{{E}_{\mathrm{th}}}{{{\rm{\Omega }}}_{b}}.\end{eqnarray} \tag{ 13 }$$

One can find a detailed description of the thin shell method in Bisnovatyi-Kogan & Silich (1995).

Figure 3 presents the output from a simulation that does not take into consideration the ambient gas rotation. The selected halo and disk gas densities and the gaseous disk scale height are n_h = 5 × 10⁻²cm⁻³, n₀ = 2 cm⁻³, and H_z = 200 pc, respectively. The energy input rate L(t) was obtained from the SFH as was explained in Section 2.2.

Zoom In Zoom Out Reset image size

The left-hand column in Figure 3 shows the superbubble cross section by the x-z plane at the bubble ages of 10 Myr (upper panel), 15 Myr (middle panel), and 35 Myr (bottom panel), respectively.

One can observe that the bubble shape deviates significantly from a spherical one at the age of 15 Myr. By the time it reaches 35 Myr, the bubble acquires a mushroom-like morphology that remains unchanged thereafter, despite the growth in its size. This strongly affects the column density distribution calculated along lines of sight, which are normal to the galactic plane (see the right-hand column in Figure 3). Indeed, the gas column density distribution becomes thicker with the bubble age, while its maximum moves further away from the bubble center.

It is interesting to note that the reference model-predicted column density distribution at the bubble age of 35 Myr looks fairly similar to the emission profiles observed in the direction of the NGC 628 largest bubble (see Figure 5 in Mayya et al. 2023).

Disk inclination and differential galactic rotation destroy the azimuthal symmetry in the simulated gas surface density maps. Therefore, we first calculate the azimuthally averaged gas column density distributions (see Section 3.2 in Mayya et al. 2023) and then convolve them with a Gaussian kernel to simulate the impact of the beam characteristics on the model-predicted column density distribution:

$$\begin{eqnarray}&&{N}_{\mathrm{con}}({R}_{i})=\displaystyle \frac{{\sum }_{j}{N}_{\mathrm{mod}}({R}_{j})\times \exp (-{\left({R}_{j}-{R}_{i}\right)}^{2})/(2{\sigma }^{2})}{{\sum }_{k}\exp (-{\left({R}_{k}-{R}_{i}\right)}^{2})/(2{\sigma }^{2})},\end{eqnarray} \tag{ 14 }$$

where N_mod and N_con are the model-predicted and the convolved gas column densities. The parameter σ is determined by the beam FWHM: $\sigma ={\rm{FWHM}}/\sqrt{8\mathrm{ln}2}$. We select an FWHM = 100 pc as this is approximately the FWHM of the CO observations (Leroy et al. 2021).

The convolution slightly reduces the peak column density value and makes the column density distribution broader than that obtained in the hydrodynamical simulations. However, its impact on the peak position is insignificant.

The impact of the host disk inclination on the column density distribution is shown in Figure 4, which presents the azimuthally averaged and convolved surface density profiles obtained upon different assumptions regarding the host galaxy inclination angle. The selected bubble age is 35 Myr. The gas density ρ₀, the disk scale height H_z , and the energy input rate L(t) used in the simulations are the same as those used in the reference model (see Section 3.1).

Zoom In Zoom Out Reset image size

The solid and dashed lines in Figure 4 correspond to the disk inclination angles i = 10° and i = 30°, respectively. This plot demonstrates the major effects of the host galaxy's inclination. The peak in the column density distribution moves toward the bubble center, the maximum column density slightly increases, and the column density distribution becomes wider in models with larger inclination angles.

We now consider the effects of two input parameters that for face-on galaxies cannot be determined directly from observations because of their degeneracy: these are the midplane gas density ρ₀ and the gaseous disk scale height H_z . The impact of the gas midplane density is shown in the upper panel of Figure 5 where we compare the model-predicted column density distributions in cases with n₀ = 1 cm⁻³ (solid line) and n₀ = 5 cm⁻³ (dotted line) keeping the other input parameters identical to those in Section 3.1. As one can note, the peak in the column density distribution is larger and moves toward the bubble center in the simulations with a larger midplane gas density. This occurs because the bubble expansion is slower in denser ambient media.

Zoom In Zoom Out Reset image size

The bottom panel in this figure shows how the value of the disk scale height affects the calculated gas column density distribution. The solid line in this panel shows the column density distribution in the case of a smaller scale height (H_z = 180 pc). In the case of the larger scale height (H_z = 300 pc, dotted line), the peak in the column density distribution is located further away from the bubble center, and the maximum value of the column density slightly increases. This is because, in this case, the bubble does not propagate very rapidly along the z-axis leading to a larger inner bubble pressure and larger cylindrical radii.

These results show that observations of large bubbles in nearly face-on galaxies together with appropriate numerical models have the potential power to solve the midplane gas density–gaseous disk scale height degeneracy problem. Note that Figure 5 presents simulated column densities convolved with an FWHM = 100 pc beam.

Finally, in this section, we discuss the impact of gravity and differential disk rotation on the bubble appearance. Two parameters were added to the set of reference model input parameters: the stellar disk scale height and the gas rotation velocity (see Equations (2)–(4) in Section 2.1 and Equation (10) in Section 2.3). We adopted H_d = 400 pc and V_rot = 165 km s⁻¹ (see Aniyan et al. 2018).

The evolution of the midplane cross section of the bubble is shown in Figure 6. Here the solid, dotted, and dashed lines display the bubble cross-section shape at the ages of 10, 25, and 50 Myr, respectively. The other input parameters used in the simulations are Σ_D = 30 M_⊙ pc⁻², n₀ = 2 cm⁻³, and H_z = 200 pc. The inclination angle is i = 9°. One can note that the differential galactic rotation distorts the cross-section shape significantly after about 20 Myr of the bubble expansion. After this time the cross section obtains an elliptical form and becomes progressively more elongated with the bubble age. In addition, the cross-section semimajor axis rotates with time.

Zoom In Zoom Out Reset image size

Figure 7 presents the projected, azimuthally averaged, and convolved (FWHM = 100 pc) radial column density at the age of t = 35 Myr. One can compare the dashed line in Figure 7 with the right-hand bottom panel in Figure 3 (which presents the results of the simulations with the same input parameters but without gravity) to realize how gravity and galactic rotation affect the results. In the calculations without gravity and rotation presented in Figure 3, the peak in the column density distribution is located at ∼0.65 kpc, while in the simulations with gravity and rotation, it is located closer to the bubble center, at ∼0.6 kpc. In addition, the column density distribution becomes broader in the simulations with gravity and the ambient gas rotation.

Zoom In Zoom Out Reset image size

We also present in Figure 7 two models with different stellar disk scale heights, H_d = 200 pc (solid line) and H_d = 600 pc (dashed line), respectively. One can note that the impact of the stellar disk scale height on the results is not significant. It slightly changes the bubble expansion velocity along the z-axis and the bubble shape but does not affect the column density distribution significantly (one can note a small difference only at large radii ≥0.7 kpc).

We now fix the inclination angle to i = 9°, H_d = 400 pc, Σ_D = 30 M_⊙ pc⁻², and V_rot = 165 km s⁻² (Kamphuis & Briggs 1992; Dutta et al. 2008; Aniyan et al. 2018), and vary the midplane gas density n₀ and the scale height H_z in models with different halo densities n_h looking for the best fit to the observed column density distribution in the NGC 628 southeast superbubble. We compare the results of our simulations to the sum of the neutral and molecular gas column densities obtained from the observed, azimuthally averaged H i and CO radial intensity profiles (see Figure 5 in Mayya et al. 2023) as the contribution of the ionized gas to the total mass is negligible.

To derive column densities from the observed H i intensity, we make use of the following relation from Walter et al. (2008):

$$\begin{eqnarray}&&{N}_{\mathrm{HI}}[{{\rm{cm}}}^{-2}]=1.823\times {10}^{18}{I}_{\mathrm{HI}}.\end{eqnarray} \tag{ 15 }$$

The molecular gas column density is

$$\begin{eqnarray}&&{N}_{{{\rm{H}}}_{2}}[{{\rm{cm}}}^{-2}]={X}_{\mathrm{CO}}{I}_{\mathrm{CO}},\end{eqnarray} \tag{ 16 }$$

where ${X}_{\mathrm{CO}}=3.3\times {10}^{20}{({\rm{K}}\,{\rm{km}}\,{{\rm{s}}}^{-1})}^{-1}$ and I_{H i} and I_CO are the H i and CO intensities expressed as the velocity integrated surface brightness temperatures in units of K km s⁻¹. The value of X_CO is that of the Milky Way CO(1–0) (Bolatto et al. 2013) taking into account a factor of CO(2–1)/CO(1–0) of 0.61 measured for this galaxy (den Brok et al. 2021).

Figure 8 compares the results of the simulations with different halo densities that reasonably reproduce the value of the maximum column density and the column density peak position observed around the NGC 628 southeast superbubble. Larger scale heights H_z are required in simulations with smaller halo gas densities to obtain maximum column densities similar to the observed value and accommodate it near the observed position. Indeed, the required H_z increases from about 200 pc in simulations with n_h = 10⁻¹ cm⁻³ to about 330 pc in simulations with n_h = 10⁻³ cm⁻³. However, in simulations with a large halo density (n_h = 10⁻¹ cm⁻³) the column density distribution looks too flat at large radii. It becomes narrower in simulations with smaller halo gas densities. In these cases, the model-predicted gas column density also drops significantly in the central (r ≤ 0.3 kpc) zone. In contrast, the value of the halo gas density almost does not affect the required midplane gas density (n₀ = 2.6, 2.4, and 2.5 cm⁻³ in models with n_h = 10⁻³, 10⁻², and 10⁻¹ cm⁻³, respectively).

Zoom In Zoom Out Reset image size

The model that best fits observations is shown in Figure 9, where the left-hand panel displays the simulated column density map at the bubble age of 50 Myr and the right-hand panel shows the theoretical and the observed azimuthally averaged radial column density distributions. Here the solid line displays the simulated radial column density distribution after the model results were convolved with an FWHM = 100 pc beam, which corresponds to the approximately an FWHM = 2'' beam of the CO observations at the distance of NGC 628. The observed column densities are shown by the triangle symbols. The observed H ii + H i column densities in the innermost parts of the cavity are within the 3σ noise and hence the plotted points correspond to the column density upper limit.

Zoom In Zoom Out Reset image size

The model is in excellent agreement with observations at distances ≥300 pc from the bubble center, where most of the swept-up gas is located. At smaller radii, the model-predicted column densities also agree with observations as all the model points fall below the upper limits obtained in Mayya et al.(2023).

It is interesting to note that the model-predicted ellipticity q = b/a, where a and b are the semimajor and semiminor axes of the hole, is ∼0.7, in good agreement with the observed value derived from the star-forming clumps distribution in the shell. Furthermore, the shell diameter along the semimajor axis is d = 2a ∼ 1.2 kpc, also close to the observed values.

We also performed numerical simulations with the best-fitted model input parameters upon the assumption of an exponential vertical gas distribution and did not find significant differences with the results presented in Figure 9.

Another aspect of the differential galactic rotation is shown in Figure 10. Here we present the normalized angular column density distribution within the ring 460 pc ≤r ≤ 930 pc, to compare it with the normalized flux azimuthal distribution presented in Figure 8 by Mayya et al. (2023). The position angle (PA) in Figure 10 is measured counterclockwise from the positive x-axis. The figure clearly demonstrates that the majority of the swept-up matter is accumulated at the opposite tips of the major axis of the oval-shaped bubble-driven shell, within the shell sectors around PA ∼ 30° and PA ∼ 200°. It is likely that the accumulation of interstellar matter in the opposite tips of the superbubble belt yields a secondary star formation in these regions. This suggestion is supported by the fact that the positions of the enhancements in the model-predicted surface density distribution are in remarkably good agreement with peaks in the observed Hα emission, which trace the sites of recent star formation (see Figure 8 in Mayya et al. 2023).

Zoom In Zoom Out Reset image size

It is worth noting that the comparison of the model-predicted column density distributions with the observed profiles restricts the midplane gas densities and the gaseous disk scale heights (the best model requires n₀ ≈ 2.3 cm⁻³, n_h ≈ 5 × 10⁻² cm⁻³, and H_z ≈ 250 pc, respectively) and thus may solve the gaseous disk midplane density–scale height degeneracy.

It was assumed in our simulations that the bubble expands into a smooth interstellar gas distribution, while the NGC 628 galactic disk has a very complex density structure, as evident in the 770W JWST image (see Figure 1). However, it is unlikely that small-scale (in comparison with the superbubble size) inhomogeneities affect our conclusions significantly. It is expected that the bubble-driven shell just overtakes sufficiently smaller bubbles. Certainly, this should result in a more rippled, less smooth shell structure, but it should not affect the shell dynamics significantly. The situation becomes more complicated if the superbubble was not driven initially from a single center, but instead formed via the merging of several bubbles comparable in size. We cannot exclude this scenario, but it is difficult to believe that in this case, the resulting shell would have an almost perfect elliptical shape (see Figure 1). Collisions with spiral arms can distort the superbubble shape, and in addition, induce star formation at the sites of collisions. This seems to be happening at the northwest side of the superbubble where the recent star formation is concentrated, but we still do not see a significant shell distortion there, probably because the shell reached the spiral arm only recently. The collisions with spiral arms could be included in the model, but this requires more information regarding the position and density distribution in the spiral arms.

It is worth noting, however, that despite the simplified assumptions, the model predictions align well with observations. This leads us to assert that the gas density stratification, and the known galaxy disk rotation velocity and gravity, are the main factors required for modeling the observational appearance of stellar feedback-driven bubbles. Furthermore, together with strong constraints on the energy input rate obtained from observations, they allow one to fit the observed properties of the stellar feedback-driven bubbles and obtain the host disk parameters in nearly face-on galaxies.