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

Here we explore the impact of all major factors, such as the nonhomogeneous gas distribution, galactic rotation, and gravity, on the observational appearance of superbubbles in nearly face-on spiral galaxies. The results of our 3D numerical simulations are compared to the observed gas column density distribution in the largest southeast superbubble in the late-type spiral galaxy NGC 628. We make use of the star formation history inside the bubble derived from the resolved stellar population seen in Hubble Space Telescope images to obtain its energy and demonstrate that the results of numerical simulations are in good agreement with the observed gas surface density distribution. We also show that the observed gas column density distribution constrains the gaseous disk scale height and the midplane gas density if the energy input rate can be obtained from observations. This implies that observations of large holes in the interstellar gas distribution and their stellar populations have the potential power to solve the midplane gas density–gaseous disk scale height degeneracy problem in nearly face-on galaxies. The possible role of superbubbles in driving the secondary star formation in galaxies is also briefly discussed.


INTRODUCTION
Since the pioneer papers by Heiles (1979Heiles ( , 1984Heiles ( , 1987) ) and Brinks & Bajaja (1986), it became clear that shelllike structures and holes are characteristic features of the Interstellar gas distribution in star-forming galaxies.
The large kinetic energies (up to 10 53 erg) of the observed shells, their sizes, shapes and radial distributions led Bruhweiler et al. (1980); Ehlerova & Palous (1996) to concluded that they are a natural by-product of the stellar feedback on the interstellar medium (ISM) of their host galaxies.These results are in agreement with Castor et al. (1975), who suggested that the gas ejected by massive stars is heated up to large temperatures (10 6 K -10 7 K) within the bubble volume due to multiple stellar wind collisions and supernovae (SNe) explosions.This enhances the inner thermal pressure that drives a strong shock into the surrounding ISM and sweeps it up into a dense shell.
Indeed, diffuse X-ray emission has been detected around young star-forming regions in the Large Magellanic Cloud (LMC, e.g Chu & Mac Low 1990;Dunne et al. 2001;Oey et al. 2002).Bagetakos et al. (2011) detected more than 1000 HI holes in a sample of nearby galaxies studied within "The HI Nearby Galaxy Survey" (THINGS) project.The hole sizes range from ∼ 100 pc to about 2 kpc while the ages of the embedded stellar populations were estimated to be in the range of (3-150) Myr.McLeod et al. (2020) made an important step forward in systematic studies of the stellar feedback on the host galaxy ISM.They presented MUSE integral field unit (IFU) observations of the nearby (∼ 2 Mpc) dwarf spiral galaxy NGC 300.These observations, in combination with Hubble Space Telescope (HST) photometry, allowed them to obtain characteristics of individual massive stars in five large HII regions and study the impact of stellar feedback on the ambient ISM in this particular case.Nath et al. (2020) studied the superbubbles size distribution and showed by means of numerical simulations that the largest superbubbles are likely driven by multiple OB associations and can reach ∼ 1 kpc in tens of million years.JWST Mid-Infrared Instrument (MIRI) observations now allows us to identify bubbles driven by young stellar clusters by tracing emission of the Polycyclic Aromatic Hydrocarbon (PAHs) molecules associated to shells around the bubbles (Rodríguez et al. 2023;Watkins et al. 2023).
Relatively little attention has been paid to the impact of the host galaxy disk rotation and inclination on the study of superbubbles.It is probably because this requires time-consuming 3D hydrodynamical simulations.The effects of the disk rotation were first discussed by Palous et al. (1990) who found, by means of 2D calculations, that at later times shells are distorted by the differential galactic rotation and the swept-up mass is concentrated at the opposite tips of the wind-driven shell.3D simulations by Silich (1992) confirmed these findings and also showed that in the case of a plane-stratified interstellar gas distribution the majority of the swept-up mass is concentrated in the bubble belt -a thin zone of the expanding shell next to the midplane of the host galaxy.
NGC 628 is a late-type nearby face-on spiral galaxy at a distance of about 9.8 Mpc (Kreckel et al. 2019;Anand et al. 2021).It was extensively observed in the past (e.g.Sánchez et al. 2011;Grasha et al. 2015) and recently within the frame of the PHANGS (Physics at High Angular resolution in Nearby GalaxieS) survey (e.g.Leroy et al. 2021;Emsellem et al. 2022;Lee et al. 2022), THINGS (Walter et al. 2008) and the PHANGS -JWST Cycle I treasury program (Lee et al. 2023).It presents a rich population of bubbles and holes whose radii vary from tens to thousand parsecs (see also Pokhrel et al. 2020).Mayya et al. (2023) and Barnes et al. (2023) thoroughly discussed the multiband properties of the largest, kiloparsec-size, South-East superbubble of NGC 628 (see Fig. 1) using multiband JWST, HST, ALMA and VLT observations.The shell that surrounds this superbubble is traced in CO, PAHs, HI and Hα emission.Molecular gas dominates over the neutral and ionized gas components.The shell total mass and diameter are ∼ 2 × 10 7 M⊙ and ∼ 1 kpc, respectively.Mayya et al. (2023) determined the star formation history inside the bubble and in the surrounding stellar disk and concluded that the mechanical power of the enclosed stellar population exceeds that required to form a spherical shell of such mass and size.In Fig. 1 we show a JWST/MIRI image of NGC 628 in the F770W filter, which trace the bubble.
Another issue that requires careful analysis is the observed shell thickness.Mayya et al. (2023) showed that the observed shell is well resolved in the JWST/MIRI images and that it is rather thick (∼ 200 pc) when compared with the bubble size (∼ 1 kpc).This seems to be in conflict with the standard wind-driven bubble theory.Indeed, the gas density behind an adiabatic bubble- driven shock is about 4 × ρ ISM , where ρ ISM is the gas density of the ambient ISM.Hence, the shell thickness should be about ∆R = R/12 as the swept-up mass is located within the shell.In more realistic, radiative shock models, the shell should be much denser and thinner.Does this imply that the stellar feedback-driven bubble model is not an appropriate explanation for the observed structure?We demonstrate by means of 3D numerical simulations that this is not the case and that the observed thick column density distribution is a natural by-product of a stratified interstellar gas distribution.
Here we present the results of comprehensive 3D numerical simulations of the evolution of superbubbles in a non homogeneous, disk-like ISM.The simulations include the ISM differential rotation, galactic disk inclination and the star formation history obtained from the analysis of the stellar population detected in the HST/ACS and JWST/NIRCam images.It is shown that the energy input rate determined by the star formation history within the South-East superbubble is consistent with the observed bubble properties.The comparison of the model-predicted projected gas column density distribution with the observed one allowed us to determine the galactic gaseous disk properties: its midplane density and vertical scale.The comparison of the hydrodynamic model predictions with the observed bub-ble properties thus allows one to resolve the midplane gaseous disk density -disk height scale degeneracy in nearly face-on galaxies.
The paper is organized as follows.The model setup is presented in section 2. Specifically, the model adopted for the interstellar gas distribution, the adopted rotation curve and the gravitational field used in the simulations are discussed in Section 2.1.In Section 2.2 we make use of the star formation history obtained by Mayya et al. (2023) to derive the energy input rate responsible for the bubble expansion.In Section 2.3 the hydrodynamic scheme used for simulations is presented.The impact of the input parameters on the column density distribution is discussed in section 3.Here we also present the model that best fits the observations.The major results and conclusions are summarized in Section 4.

MODEL SETUP
The theory of wind-driven bubbles was developed by Castor et al. (1975); Weaver et al. (1977); Tomisaka & Ikeuchi (1986, 1988); Mac Low & McCray (1988); Koo & McKee (1992); Suchkov et al. (1994) and others (see for a review Tenorio- Tagle & Bodenheimer 1988;Bisnovatyi-Kogan & Silich 1995), who considered the thermalized kinetic energy of individual stellar winds and SNe to be the major driving mechanism responsible for the formation and evolution of interstellar bubbles.The impact of additional physical processes, such as radiation pressure and even cosmic rays on the bubble dynamics and observational appearance has been considered later (see reviews by Krumholz et al. 2019;Rosen & Krumholz 2020).
Here we make use of the thin-shell approximation (section 2.3) to follow the superbubble evolution and shape transformation in a non-homogeneous galactic interstellar medium (ISM) under the influence of the galactic gravity and the ambient gas rotation (section 2.1).We also estimate the superbubble driving energy from the star-formation history derived from HST observations (see section 2.2).

The adopted interstellar gas distribution and rotation curve
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: where ρ h = µn h is the gas volume density in the galactic halo, ρ mp = ρ 0 +ρ 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): where x, y, and z are Cartesian coordinates, G is the gravitational constant, V rot is the rotation velocity, R 2 = x 2 + y 2 is the distance along the galactic plane, H d and Σ D are the stellar disc scale-height and surface density.
The parameters n h , n 0 , 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 energy input rate
The star formation history within the South-East superbubble was carefully analysed 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 (Fig. 13 top in Mayya et al. 2023), as a sequence of N tot instantaneous starbursts separated by ∆t = t max /N tot time intervals, where t max ∼ 50 Myr is the time passed since the superbubble starburst onset till 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 ministarbursts: where t * = (i − j − 1) ∆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 Ṁ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 supernovae upon the assumption of a solar metallicity starburst and a canonical Kroupa initial mass function with lower and upper cutoffs of M low = 0.1M ⊙ and M up = 100M ⊙ , respectively.
The mass of each mini-starburst was calculated as: The resulting mechanical luminosity as a function of the bubble age is presented in Fig. 2. The energy input rate first grows fast and then remains almost constant until the SFR becomes negligible at the 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 Mayya et al. (2023) found that recent star formation is not confined to just the inside-bubble zone (see their Fig.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 where L d (t) and S b (t) are the stellar disk mechanical luminosity (red line in Fig. 12 from 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.

Main equations
All simulations were performed with a 3D code based on the thin shell numerical scheme presented in Palous (1992); Silich (1992); Ehlerová & Palouš (2018); 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 inside-bubble thermal pressure P th is uniform.The shell is split into a number of Lagrangian elements and the equations of mass, momentum and energy conservation are solved for each Lagrangian element: where r i , M i , and dΣ i are the i-th 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 i-th 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 here presented are performed with N shell ∼ 1600 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 inside bubble thermal pressure is: One can find a detailed description of the thin shell method in Bisnovatyi-Kogan & Silich (1995).

RESULTS AND DISCUSSION
In this section we confront our model with the observed gas column density distribution in the largest superbubble at the South -East of NGC 628 (see Fig. 1) whose multi-band properties were thoroughly discussed by Mayya et al. (2023) and Barnes et al. (2023).Following Mayya et al. (2023), we adopted a galactocentric radius R ∼ 2.8 kpc for the superbubble center.At this radius Σ D = 30 M ⊙ pc −2 and the disk rotation velocity V rot is approximately 165 km s −1 (Aniyan et al. 2018).We assume that the bubble was formed 50 Myr ago, when the star formation increased at the present bubble center.The SFR decreased thereafter as described by the star formation history (SFH) reported by Mayya et al. (2023).However, in the simulations we also took into account the mechanical luminosity of the older disk component that contributes to the bubble energy balance as each generation of massive stars within the superbubble interior supplies energy during ∼ 40 Myr, when the least massive (∼ 8 M ⊙ ) stars explode as supernovae (see section 2.2).
Estimates of the stellar disk scale-height in nearly face-on galaxies are difficult and rather uncertain.We adopt H d = 400 pc as a reference value for the NGC 628 stellar disk scale-height (see Aniyan et al. 2018).However, in section 3.4 we discuss the impact of the stellar disk scale-height and gravity on the projected gas column density distribution.Another input parameter used in our simulations is the galaxy inclination, which does not affect the hydrodynamical simulations but does affect the simulated column density maps.This is discussed below in section 3.2.We adopt for the NGC 628 disk inclination i = 9 • (Kamphuis & Briggs 1992;Dutta et al. 2008;Aniyan et al. 2018).
The remaining parameters, the gas disk density n 0 , the halo gas density n h , and the scale-height H z , are not determined by observations.We vary them to find the best fit to the observed gas column density distribution.

The reference model
Fig. 3 presents the output from a simulation that does not take into consideration the ambient gas rotation and gravity.The selected halo and disk gas densities and the gaseous disk scale-height are n h = 5×10 −2 cm −3 , n 0 = 2 cm −3 , and H z = 200 pc, respectively.The energy input rate L(t) was obtained from the star formation history as was explained in section 2.2.
The left-hand column in Fig. 3 shows the superbubble cross-section by the x-z plane at the bubble age 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 mushroomlike morphology that remains unchanged thereafter, despite the growth in its size.This strongly affects the column density distribution calculated along line of sights which are normal to the galactic plane (see the righthand column in Fig. 3).Indeed, the gas column density distribution is getting thicker with the bubble age, while its maximum moves further away from the bubble center.
It is interesting to note that the reference modelpredicted 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 Fig. 5 in Mayya et al. 2023).

Effects of the galactic disk inclination
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: where N mod and N con are the model-predicted and the convolved gas column densities.The parameter σ is determined by the beam FWHM: σ = FWHM/ √ 8 ln 2. We select 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 Fig. 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 ρ 0 , 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).
The solid and dashed lines in Fig. 4 correspond to the disk inclination angles i = 10 • , i = 30 • , respectively.This plot demonstrates the major effects of the host galaxy inclination.The peak in the column density distribution moves towards the bubble center, the maximum column density slightly increases and the column density distribution becomes wider in models with larger inclination angles.

Impact of the midplane gas density and gaseous disk scale-height
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 ρ 0 and the gaseous disk scaleheight H z .The impact of the gas midplane density is shown in the upper panel of Fig. 5 where we compare the model-predicted column density distributions in cases with n 0 = 1 cm −3 (solid line) and n 0 = 5 cm −3 (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 towards the bubble center in the simulations with a larger midplane gas density.This occurs because the bubble expansion is slower in denser ambient media.
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 so rapidly along the zaxis 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 degen-  eracy problem.Note that Fig. 5 presents simulated column densities convolved with a FWHM = 100 pc beam.

Impact of gravity and disk rotation
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 the 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 −1 (see Aniyan et al. 2018).
The bubble midplane cross-section evolution is shown in Fig. 6.Here the solid, dotted and dashed lines display the bubble cross-section shape at the age of 10 Myr, 25 Myr, and 50 Myr, respectively.The other input parameters used in the simulations are: Σ D = 30 M ⊙ pc −2 , n 0 = 2 cm −3 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 semi-major axis rotates with time.
Fig. 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 Fig. 7  in Fig. 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 Fig. 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.We also present in Fig. 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.7kpc).

Best-fitted model
We now fix the inclination angle to (Kamphuis & Briggs 1992;Dutta et al. 2008;Aniyan et al. 2018), and vary the midplane gas density n 0 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 South-East superbubble.We confront the results of our simulations to the sum of the neutral and molecular gas column densities obtained from the observed, azimuthally averaged HI and CO radial intensity profiles (see Fig. 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 HI intensity we make use of the following relation from n hal = 10 −3 cm −3 n hal = 10 −2 cm −3 n hal = 10 −1 cm −3 Figure 8. Impact of the halo gas density on the model-predicted column density distribution.The solid, dashed and dotted lines correspond to n h = 10 −1 cm −3 , 10 −2 cm −3 , and 10 −3 cm −3 , respectively.Walter et al. (2008): The molecular gas column density is: where X CO = 3.3 × 10 20 (K km s −1 ) −1 and I HI and I CO are the HI and CO intensities expressed as the velocity integrated surface brightness temperatures in units of K km s −1 .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).Fig. 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 South-East 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 −1 cm −3 to about 330 pc in simulations with n h = 10 −3 cm −3 .However, in simulations with a large halo density (n h = 10 −1 cm −3 ) 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 0 = 2.6 cm −3 , 2.4 cm −3 and 2.5 cm −3 in models with n h = 10 −3 cm −3 , 10 −2 cm −3 , and 10 −1 cm −3 , respectively).
The model that best fits observations is shown in Fig. 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 a FWHM = 100 pc beam, which corresponds to the approximately 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 H2 + HI column densities in the inner-most parts of the cavity are within the 3-σ noise and hence the plotted points correspond to the column density upper limit.
The model is in excellent agreement with observations at distances ≥ 300pc 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 observe that the ellipticity q = b/a, where a and b are the semi-major and semi-minor axes of the hole, is ∼ 0.7 in our calculations, in good agreement with the observed value derived from the star-forming clumps distribution in the shell.Furthermore, the shell diameter along the semi-major 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 Fig. 9.
Another aspect of the differential galactic rotation is shown in Fig. 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 on Fig. 8 by Mayya et al. (2023).The position angle (PA) in Fig. 10 is measured counter-clockwise 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 interstellar matter accumulation in the opposite tips of the superbubble belt yields in 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 Fig. 8 in Mayya et al. 2023).It is important to highlight that the comparison of the column density distributions derived from the model and the observed profiles restricts the midplane gas densities and the gaseous disk scale-heights (the best model requires n 0 ≈ 2.3 cm −3 , n h ≈ 5 × 10 −2 cm −3 , and H z ≈ 250 pc, respectively) and thus may solve the gaseous disk midplane density -scale-height degeneracy.

Model uncertainties and simplifications
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 Fig. 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 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 Fig. 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 North-West 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 into 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-feedbackdriven 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.

SUMMARY AND CONCLUSIONS
Here we discussed the evolution of superbubbles in galaxies with a disk-like interstellar gas distribution.The impact of different input parameters, such as the host disk inclination, differential galactic rotation and gravity, on the observational bubble appearance was thoroughly discussed.The results of 3D numerical simulations were confronted to the observed properties of the largest, ∼ 1 kpc in size, South-East superbubble in the nearly face-on spiral galaxy NGC 628.
We made use of the star formation history derived from HST/ACS and JWST/NIRCam observations by Mayya et al. (2023) to obtain the inside-bubble stellar population mechanical power and then performed numerical simulations of the superbubble evolution upon different assumptions regarding the interstellar gas properties.The simulations showed that the mechanical power of the inside-bubble stellar populations is sufficient to explain the observed, ∼ 1 kpc, hole size.
We then made use of multiple numerical calculations to find the model that best fits the observed radial column density distribution.For each set of input parameters the results of the simulations were convolved with a FWHM=100 pc beam and confronted to the observed column density distribution.The results show that a certain midplane gas density and a certain gaseous disk scale-height are required to fit the observations.This implies that the comparison of large holes in the interstellar gas distribution and their stellar populations with the results of numerical simulations has the potential power to solve the midplane gas density -gaseous disk scale-height degeneracy problem in nearly face-on galaxies.
Two observational parameters, the position of the peak in the column density distribution and the value of the maximum column density, must be fitted simultaneously to determine the gaseous disk scale-height and the midplane gas density.In the particular case of the NGC 628 South-East superbubble this method led us to conclude that the pre-superbubble midplane disk density and the gaseous disk scale-height are n mp ≈ 2.3 cm −3 and H z ≈ 250 pc, respectively.The model also predicts a non homogeneous angular column density distribution with the swept-up mass concentrated at the opposite tips of the bubble-driven shell.It is remarkable that the model-predicted enhancements in the azimuthal column density distribution correspond to the sites of the enhanced Hα emission which mark the sites where the recent star formation occurred.This favors the hypothesis that superbubbles may trigger a secondary star formation in those zones of the bubble-driven shells where a major fraction of the swept-up interstellar matter is accumulated.

Figure 1 .
Figure 1.The JWST/MIRI F770W image showing the largest, kpc-size hole in the gas column density distribution in the South-East zone of NGC 628.The inner and outer edges of the region with an enhanced gas column density around the hole are displayed by the white ellipses.The radial vector joining the hole center to the center of the galaxy is also shown.

Figure 3 .
Figure 3.The reference superbubble evolution.The left-hand column presents the bubble shape at the age 10 Myr (upper panel), 15 Myr (middle panel), and 35 Myr (bottom panel).The right-hand column displays the model-predicted column density distributions at the same bubble ages.The input parameters selected for these simulations are n0 = 2 cm −3 , Hz = 200 pc, i = 0 • , and the energy input rate derived from the star formation history in the NGC 628 South-East superbubble (see section 2.2).

Figure 4 .
Figure 4.The azimuthally-averaged gas column density distribution after the convolution with a FWHM = 100 pc beam at the age of 35 Myr.The solid and dashed lines correspond to the disk inclination angles i = 10 • and i = 30 • , respectively.All other input parameters are the same as in the case of the reference model (see section 3.1).

Figure 5 .
Figure5.The azimuthally-averaged gas column density distribution in models with different values of n0 and Hz and fixed inclination angle i = 9 • for the host galaxy.The bubble age is 35 Myr.The upper panel displays the projected column density distributions in the case when the gaseous disk scale-height is fixed to Hz = 200 pc, but densities n0 are different: n0 = 1 cm −3 (solid line) and n0 = 5 cm −3 (dotted line).The lower panel presents the projected column density distributions in the case when the value of the gas density is fixed to n0 = 2 cm −3 , but the disk scale-heights are different: Hz = 180 pc (solid line) and Hz = 300 pc (dotted line), respectively.

Figure 7 .
Figure 7.The azimuthally-averaged radial column density distribution in models with different stellar disk scale-heights at the bubble age of 35 Myr.The solid, dashed and dotted lines correspond to H d = 200 pc, H d = 400 pc, and H d = 600 pc, respectively.The other input parameters are the same as in the model presented in Fig. 6.

Figure 9 .
Figure9.The best-fitted bubble model.The left panel presents the projected gas column density map at full numerical resolution in the case when the NGC 628 disk inclination angle is i = 9 • .The corresponding azimuthally-averaged radial gas column density distribution after convolving with a beam of 100 pc is shown in the right panel.The model input parameters are n0 = 2.3 cm −3 , n h = 5 × 10 −2 cm −3 , Hz = 250 pc, H d = 400 pc, ΣD = 30 M⊙ pc −2 , and Vrot = 165 km s −1 .The results of the simulations are presented by the solid line while triangle symbols display the sum of the observed neutral and molecular gas column densities.

Figure 10 .
Figure 10.The best-fitted model column density angular distribution.The position angle (PA) is measured counterclockwise from the positive x-axis.