The stellar population responsible for a kiloparsec size superbubble seen in the JWST "phantom" images of NGC628

TL;DR

This study investigates a large superbubble in NGC628, revealing its stellar population, gas properties, and star formation activity, highlighting the role of steady stellar feedback in maintaining the shell structure.

Contribution

It provides a detailed analysis of the stellar population and gas dynamics of a kiloparsec-sized superbubble using multi-band JWST and HST data, a novel combination for such a large structure.

Findings

The superbubble contains ~2x10^7 Msun of gas, mostly molecular.

A stellar population of 10^5 Msun formed over 50 Myr inside the bubble.

The shell is currently triggering new star formation.

Abstract

We here study the multi-band properties of a kiloparsec-size superbubble in the late-type spiral galaxy NGC628. The superbubble is the largest of many holes seen in the early release images using JWST/MIRI filters that trace the Polycyclic Aromatic Hydrocarbon (PAH) emissions. The superbubble is located in the interarm region ~3 kpc from the galactic center in the south-east direction. The shell surrounding the superbubble is detected in HI, CO, and Halpha with an expansion velocity of 12 km/s, and contains as much as 2x10^7 Msun of mass in gas that is mostly in molecular form. We find a clear excess of blue, bright stars inside the bubble as compared to the surrounding disk on the HST/ACS images. These excess blue, bright stars are part of a stellar population of 10^5 Msun mass that is formed over the last 50 Myr in different star formation episodes, as determined from an analysis of color-magnitude diagrams using a Bayesian technique. The mechanical power injected by the massive stars of these populations is sufficient to provide the energy necessary for the expansion of the shell gas. Slow and steady, rather than violent, injection of energy is probably the reason for the maintenance of the shell structure over the kiloparsec scale. The expanding shell is currently the site for triggered star formation as inferred from the JWST 21 micron (F2100W filter) and the Halpha images.

Tables

Table 1. Table 1: Description of imaging data used in this work

Telescope/	Filter	${\lambda }_c$	$\Delta \lambda$	Pixel	FWHM		ZP	mlim	tracer
instrument		$\mu$m	$\mu$m	″	″	pc	mag	mag
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
HST/acs	F435W	0.43	0.10	0.05	0.10	4.70	25.79	28.06	stars
HST/acs	F555W	0.53	0.12	0.05	0.10	4.70	25.73	27.33	stars
HST/acs	F814W	0.83	0.25	0.05	0.10	4.70	25.53	26.47	stars
JWST/nircam	F200W	1.99	0.22	0.031	0.06	3.01	26.38	25.21	stars
JWST/nircam	F300M	2.99	0.32	0.063	0.13	6.11	24.05	23.47	stars
JWST/nircam	F335M	3.36	0.35	0.063	0.13	6.11	23.82	23.37	stars
JWST/nircam	F360M	3.62	0.39	0.063	0.13	6.11	23.67	23.32	stars
JWST/miri	F770W	7.7	2.20	0.11	0.25	11.8	—	—	PAH
JWST/miri	F1000W	10.0	2.00	0.11	0.32	15.0	—	—	Silicate
JWST/miri	F1130W	11.3	0.70	0.11	0.36	16.9	—	—	PAH
JWST/miri	F2100W	21.0	5.00	0.11	0.67	31.5	—	—	hot dust

Table 2. Table 2: Datacubes used for kinematical and line maps analysis

Telescope	line	$\lambda$	$\Delta$V	Pixel	BEAM		tracer
/instrument			km/s	″	″	pc
(1)	(2)	(3)	(4)	(5)	(6)	(7)	8
VLT/muse	H$\alpha$	6563 Å	$\sim$15	0.2	0.6	28	H+
ALMA/phangs	CO	1.3 mm	2.5	0.2	2	94	H2
VLA/things	Hi	21 cm	2.6	1.5	6.9$\times 5.6$	300	H i

Table 3. Table 3: Variables entering the hierarchical model. Here dm, A V subscript 𝐴 𝑉 A_{V} and M 𝑀 M are the distance modulus, the visual extinction and the stellar mass, respectively.

Hyper-parameter	$𝒂=\{a_1,..,a_{N_I}\}$		Stellar fraction vector
Parameters	$f_j^k$	mag	predicted absolute magnitude
	${\hat{F}}_j^k=f_j^k+\mathrm{dm}+A_V$	mag	predicted apparent magnitude
Data ($𝑭$)	$F_j^k$	mag	Observed apparent magnitude
	$e_j^k$	mag	Apparent magnitude error
Fixed distributions	${ℱ}_i^k\left(M\right)$	mag	Isochrone absolute magnitude
	${\sigma }_i^k$	mag	Isochrone tolerance
	$\varphi \left(M\right)$		Initial mass function

Table 4. Table 4: Probability distribution function. Here Γ Γ \Gamma , ξ 𝜉 \xi are the Gamma function and the concentration parameter (See text).

Name	PDF	Eq.
Hyperprior	$p\left(𝒂\right)$	$\frac{\Gamma \left(\xi N_{iso}\right)}{\Gamma {\left(\xi \right)}^{N_{iso}}}{\prod }_{i=1}^{N_{iso}}{a}_i^{\xi -1}$
Likelihood	$P\left(F_j^k|{\hat{F}}_j^k,e_j^k\right)$	${\left(2\pi e_j^k\right)}^{-1/2}\exp {\left(\frac{F_j^k-{\hat{F}}_j^k}{e_j^k}\right)}^2$
Prior	$P\left(f_j^1{f}_j^2{f}_j^3\right|\text{a}$)	${\sum }_{i=1}^{N_{iso}}{a}_i{\int }_{M_{l,i}}^{M_{u,i}}\varphi \left(M\right){\prod }_{k=1}^3𝒩\left(f_j^k|{ℱ}_i^k\left(M\right),{\sigma }_i^k\right)dM$

Equations

$$V_{\mathrm{obs}}=V_{\mathrm{sys}}+V_{\exp }\sin (i)\sin (\theta )+V_{\mathrm{rot}}\sin (i)\cos ({\theta }^{\prime })$$

$$\tan (\mathrm{P}{\mathrm{A}}^{\prime })=-\frac{\mathrm{\Delta }{\alpha }_0-\mathrm{r}\sin \mathrm{PA}}{\mathrm{\Delta }{\delta }_0+\mathrm{r}\cos \mathrm{PA}}$$

$$p(\mathbf{a}|\mathbf{F})\propto p(\mathbf{a})\prod _{j=1}^{N_D}\prod _{k=1}^3\int \frac{S(F_j^k)\mathcal{L}(F_j^k|{\hat{F}}_j^k)p(f_j^1{f}_j^2{f}_j^3|\mathbf{a})}{\ell (\mathbf{a},S)}d{f}_j^k,$$

$$N_i=a_i\cdot N_s,$$

$$\varphi (M)=\{\begin{array}{l}{C}_1{M}^{-{\alpha }_1}\mathrm{if}{M}_l\le M\le M_c\\ C_2{M}^{-{\alpha }_2}\mathrm{if}{M}_c\le M\le M_{\lim }.\end{array}$$

$$N_{\lim }(i)={\int }_{M_l}^{M_{lim}(i)}\varphi (M)\cdot dM,$$

$$C_2$$

$$C_1$$

$$N_{\mathrm{corr}}(i)=N_i+N_{\lim }(i).$$

$$⟨M⟩\approx {\int }_{M_l}^{M_u}M\cdot \varphi (M)\cdot dM,$$

$$SF{R}_i=\frac{⟨M⟩N_{\mathrm{corr}}(i)}{\Delta t_i}=\frac{⟨M⟩N_{\mathrm{corr}}(i)}{2.30258t_i\Delta ({\log }_{10}{t}_i)},$$

$$R_{\mathrm{sb}}=A{\left(\frac{L_{\mathrm{mech}}}{{\rho }_{\mathrm{gas}}}\right)}^{1/5}{t}^{3/5},$$

$$v_{\exp }=\frac{3}{5}\frac{R_{\mathrm{sb}}}{t},$$

$$A={\left[\frac{375(\gamma -1)}{28(9\gamma -4)\pi }\right]}^{1/5},$$

$$L_{\mathrm{mech}}=\frac{125}{27}\frac{{\rho }_{\mathrm{gas}}{R}_{\mathrm{sb}}^2{v}_{\exp }^3}{A^5}.$$

$$L_{\mathrm{mech}}=\frac{7(9\gamma -4)}{27(\gamma -1)}\frac{M_{\mathrm{sh}}{v}_{\exp }^3}{R_{\mathrm{sb}}}.$$

$$f_j^k\sim \mathcal{N}({\mathcal{F}}_i^k,{\sigma }_i^k),$$

Full text

1 Introduction

Star clusters are fundamental blocks of star formation (Lada & Lada, 2003) and contain both low and high-mass stars whose relative numbers can be described by an Initial Mass Function (Kroupa, 2001). Massive stars in clusters continuously inject energy and momentum to the surrounding interstellar medium (ISM) through stellar winds and supernova explosions for the first 40 Myr (Leitherer et al., 1999, STARBURST99 henceforth). The injected energy can destroy the natal cloud (negative feedback), and/or trigger star formation in a neighbouring cloud (positive feedback), or under some circumstances modify the structure of the surrounding ISM into rings which are referred to as holes, cavities, bubbles, superbubbles, shells etc. (Tenorio-Tagle & Bodenheimer, 1988). We will refer to these structures as bubbles in the present work. The role of stellar feedback in the creation of bubbles has been numerically and analytically studied by several groups (e.g. Mac Low & McCray, 1988; De Young & Heckman, 1994; Mac Low & Ferrara, 1999; Silich et al., 2001; Koo & McKee, 1992; Steinwandel et al., 2020).

Bubbles are most easily traced as HI holes by the interferrometric HI observations (Brinks & Bajaja, 1986; Boomsma et al., 2008). They are also traced in the ionized gas as shells or filaments (e.g. Hunter & Gallagher, 1990). Late-type spirals and irregular galaxies are known to show a rich population of bubbles, with their radius varying from tens of parsecs to a kiloparsec (see e.g. Pokhrel et al., 2020). The number of bubbles in galaxies scales with the star formation rate (SFR) supporting the basic idea that most of the bubbles are created by the feedback from massive stars. Holes occupy as much as 15% of the surface area traced by H i in some late-type galaxies (Pokhrel et al., 2020), giving them a porous appearance.

Stellar feedback is found to provide enough energy to explain the observed sizes of bubbles (Vorobyov & Basu, 2005). However, systematic searches for the remnant cluster that created bubbles have failed to detect an obvious candidate (Warren et al., 2011). These authors conclude that the H i holes are likely formed from multiple generations of star formation and only under suitable interstellar medium conditions.

The MIRI camera mounted on the JWST provides a new tool to study the formation of bubbles (Bouchet et al., 2015; Wright et al., 2015; Wells et al., 2015). The F770W filter intercepts the Polycyclic Aromatic Hydrocarbon (PAH) features, which traces the Photo-Dissociated Region (PDR) phase and the diffuse medium associated with the neutral gas (e.g. Tielens, 2008; Montillaud et al., 2013; Hensley et al., 2022).

The PAH molecules are most commonly excited by soft ultraviolet radiation from stars (Tielens, 2008). The PAH molecules can be easily destroyed by hard radiation either from massive stars or an Active Galactic Nucleus (AGN) and high-velocity shocks from jets and or SNe explosions (Micelotta et al., 2010a; Zhang et al., 2022). Collision with electrons in pre-shock hot gas also destroys the PAH molecules (Micelotta et al., 2010b). However, PAH molecules are known to survive in the shells surrounding the expanding bubbles (see e.g. Pabst et al., 2019, 2021).

We here use the images of NGC 628 from Early Science data release of JWST principally in the F770W filter to detect a large bubble of $\sim$1 kpc diameter. The bubble is also traced in the F1000W, F1130W and F2100W filters that trace Silicate features, large PAHs and the hot dust continuum, respectively. The bubble studied in this work is the largest of numerous bubble-like structures in the F770W of this galaxy, which has given rise to the nickname “phantom galaxy” for galaxies that have porous appearance in this filter. NGC 628 is the first of such galaxies for which data are released. We complement the JWST data with the publicly available optical images from the HST/ACS filters and nebular line images from the MUSE spectral datacube. The detailed study of the largest bubble that we carryout here would help to address the correspondence between cluster feedback and bubble formation. NGC 628 is a late type spiral galaxy of solar metallicity (Kreckel et al., 2019) at a distance of 9.77 Mpc (Olivares et al., 2010, 47 pc arcsec*-1*). It is mildly star-forming with a global star formation rate of 1.7 M*⊙* yr*-1*  (Santoro et al., 2022).

Section 2 describes the data used in this study. Sec. 3 presents the morphology, kinematics and ionization state of the ionized gas associated with the bubble using the H$\alpha$ line. In Sec. 4, we present the results about the stellar population responsible for creating the bubble. We address the energetics, star formation scenario and the special circumstances that create kiloparsec bubbles in Sec. 5.

2 JWST observations of NGC 628 and the detection of the bubble

In Figure 1 we show the F770W image of NGC 628 from the JWST mission. The image is spectacular for the presence of numerous holes (or bubbles), giving rise to the use of a popular name “phantom galaxy” for galaxies with such an appearance. Most of these holes are physically coherent structures where the gas has been evacuated by stellar feedback processes. However, some large holes could be structures between spiral arm suprs (e.g. Williams et al., 2022). The largest of these bubbles is marked by a circle of 20 arcsec diameter, which corresponds to a physical size of 940 parsec at the distance of NGC 628. Interstellar structures of such large dimensions are normally reported in the H i images as H i holes (e.g. Pokhrel et al., 2020). Large bubbles can also be formed by the collective action of stellar winds and multiple supernovae (Mac Low & McCray, 1988), in which case they are referred to as superbubbles. The analysis carried out in this work clearly shows the presence of recently formed stars inside the hole that were part of a few dissolved clusters, and hence we refer to this hole as a superbubble.

In this work, we carryout a multiband analysis of this kiloparsec-scale superbubble using the dataset listed in Tables 1 and 2. The dataset includes stellar, as well as interstellar medium tracers, offered by the JWST and HST missions. In addition, we use the MUSE datacube observed within the PHANGS survey (Kreckel et al., 2018) to generate an H$\alpha$ map to trace the morphology of the bubble in the ionized gas, the CO (2-1) line map from ALMA-PHANGS (Leroy et al., 2021a, b) to trace the molecular gas, and HI 21 cm line map from VLA-THINGS (Walter et al., 2008) to trace the atomic gas. In Figure 2, we show an enlarged view of the bubble in four of the MIRI bands. The bubble is traced in all the four bands, but is best defined in the F770W and F1130W filters, which trace the PAH molecules. At the longest JWST wavelength (21 $\mu$m), compact sources, probably embedded stars, can be seen coinciding with the bubble boundary.

NGC 628 is one of sample galaxies in the PHANGS survey (Emsellem et al., 2022) and is observed by the JWST under the title “A JWST-HST-VLT/MUSE-ALMA Treasury of Star Formation in Nearby Galaxies” (Proposal ID: 2107; PI: Janice Lee). The JWST is used to carry out imaging observation of the galaxy in both NIRCAM and MIRI cameras. In Tables 1 and 2 we list the details of the publicly available images in the JWST as well as the HST, VLT, ALMA and VLA databases.

3 Physical characteristics of the bubble

3.1 Multi-phase morphology of the bubble

In Figure 3, we compare the bubble morphology seen in the F770W map with ionized (H$\alpha$), molecular (CO), and atomic (HI) gas tracers. We show two contours that delineate the inner and outer boundaries of the bubble in the F770W band and superpose these contours on the rest of the images to facilitate comparison of the multi-phase morphologies. These contours can be approximated by circles of inner and outer radii of 8.5 and 13.8 arcsec, which correspond to 400 and 650 pc, respectively. The morphology of the ionized gas very closely follows that seen in the F770W filter. The bubble is the brightest between PA=$-30^{\circ}$ to $150^{\circ}$ (north, north-west and west), with a well-defined circular-shaped boundary. On the other hand, the inner boundary of the bubble in the diagonally opposite direction (south-east) is corrugated and elongated radially. The western part of the bubble is brighter than the south-east segment in the CO emission also. At the relatively poor spatial-resolution of the H i image, the azimuthal emission is uniform.

The bubble is shown as an RGB image in Figure 4, where we superpose CO, F770W and H$\alpha$ images as red, green and blue components, respectively. The image facilitates to visualize spatial correspondences and offsets between the molecular, PAH and ionized gas emission. There is bright emission in the shell (delineated by the cyan contours) in all the three tracers. The brightest knots in the H$\alpha$, which are identified by letters A to E, are bright in the F770W image also. However, there are noticeable offsets in the emission in the three tracers at different azimuthal locations of the shell. For example, the CO emission in the bright north-western part is systematically offset outwards. There is a faint diffuse component in H$\alpha$ emission inside the bubble, with no such emission detected in any other ISM tracer. The inner edge of the shell (the thin region between the yellow and cyan contours) is filled with a narrow H$\alpha$ shell, with the PAH lying systematically outside this thin shell traced in H$\alpha$.

The observed non-circular morphology can be qualitatively understood as effects of an expanding bubble in a non-uniform medium. The bubble center lies in the inter-arm region, with the spiral arm present to the north-west (see Figure 1), and relatively diffuse gas on the southeast. Hence, the expanding bubble is bounded by the gas in the spiral arm, forming a semi-circular arc to the north-west, whereas it is expanding into relatively rarer medium in the opposite direction, which is the reason for the observed elongation in the direction perpendicular to the spiral arm. The bubble walls are intact even in this direction as evidenced by the almost uniform H i gas surrounding the yellow contours. We discuss the likely 3-D morphology of the bubble in Sec.5.1.

3.2 Radial and azimuthal structure of the bubble

In Figure 5, we show the azimuthally averaged radial intensity profiles of the bubble in the four chosen ISM tracers, as well as in the F2100W filter. The latter filter centered at 21 $\mu$m traces active star-forming locations through the continuum emission from hot dust. These profiles were obtained using concentric circular apertures on the disk background-subtracted images, where the average value in an annular zone between 1.1 and 1.4 kpc is taken as the disk background. It may be noted that the disk background is the major contribution to the observed flux in H i, F2100W and F700W images, whereas it is negligible in H$\alpha$ and CO images. In all the five plotted tracers, the profile has minimum and maximum intensity at the center of the bubble and the shell, respectively, with the peak emission occurring at a radius=0.55 kpc in H$\alpha$, H i, F770W and F2100W filters, and at 0.65 kpc in the CO line. The decrease of intensity outside the shell is the shallowest in H i, which is the typical characteristic of an H i hole. In comparison to H i, the three other tracers of neutral ISM (CO, F770W, F2100W) show a clear intensity enhancement at the shell, in addition to a dip at the bubble center with respect to the surrounding disk, which are clear characteristics of an expanding shell. The H$\alpha$ profile illustrates that the bubble center has detectable diffuse ionized gas (DIG), with a strength that is comparable in regions outside the shell. The normalized intensities in the two JWST filters inside the bubble are negative, which suggests that the PAH molecules are destroyed inside the bubble.

It can be inferred from images presented in Figure 3 that the shell shows considerable azimuthal structures. We show the azimuthal intensity profiles in the five chosen tracers in Figure 6 which allows us to compare the observed morphology in different phases of the ISM. All emission in an annular zone between radii 0.46 kpc to 0.93 kpc are summed to obtain these profiles. The profiles show narrow peaks superposed on a slowly varying extended emission. The extended emission comes from diffuse gas, which is expected to be predominantly the gas that is swept-up as the bubble expanded. On the other hand, the peaks correspond to compact knots present in the shell. The contrast between the extended emission and the compact knots is the highest in the H$\alpha$ profile, where we can identify five zones of emission, which we identify in the figure by letters A to E. Each of these five peaks in H$\alpha$ has a corresponding peak in F770W, and F2100W profiles. Examination of the image in these two JWST filters reveals that the compact sources dominate to the observed emission at these peaks rather than the diffuse component from PAH. The emission from compact sources at MIR wavelengths originates from the hot dust surrounding the recently formed stars. These knots are indeed new star-forming sites triggered by the positive feedback of the expanding bubble (Kuhn et al., 2019). We will address this issue later in the discussion section.

3.3 Kinematics of the gas in the shell

We used the velocity maps in CO, H i and H$\alpha$ to obtain kinematical information of the shell gas in molecular, atomic and ionized phases. We illustrate the results in Figure 7, where we plot the observed velocities against the position angle, PA. In the top panel, we show the velocities obtained in the three tracers. A systematic variation of velocity as a function of PA is seen in all the three tracers. A uniformly expanding ring is expected to follow sinusoidal curve in this plot. However, given the large size of the bubble, the contribution of the galactic rotation to the observed velocity profile cannot be ignored. We hence fitted the observed velocities $V_{\rm{obs}}$ with a sinusoidal model considering the expansion centered in the bubble, plus a rotation around the center of the galaxy, and obtained the expansion velocity $V_{\rm exp}$ using the following equation (van der Kruit & Allen, 1978):

[TABLE]

where $i$ is the inclination angle of the disk, $\theta$ is the azimuthal angle in the plane of the galaxy centered in the bubble, which is related to the measured position angle $\rm{PA}$ through the equation $\tan(\theta)=\tan(\rm{PA}-\rm{PA}_{0})/\cos(i)$, where $\rm{PA}_{0}$ is the position angle of the major axis. $\theta^{\prime}$ is the azimuthal angle in the plane of the galaxy centered in the galaxy since the rotation is around that center, and it is equally related to the position angle centered in the galaxy, PA*′. We can obtain the position angle centered in the galaxy, PA′*, in terms of the position angle centered in the bubble PA by solving:

[TABLE]

where $r=$650 pc is the radius of the bubble estimated from the azimuthal averaged CO intensity profile (black line in Fig. 6). $\rm{\Delta}\alpha_{0}$ and $\rm{\Delta}\delta_{0}$ are the sky coordinates difference between the center of the bubble, and the center of the galaxy: $(\rm{\Delta}\alpha_{0},\rm{\Delta}\delta_{0})=(\alpha_{\rm{bubble}}^{\rm{center}},\delta_{\rm{bubble}}^{\rm{center}})-(\alpha_{\rm{gal}}^{\rm{center}},\delta_{\rm{gal}}^{\rm{center}})$. We have fixed $V_{\rm{sys}}=653\rm{km/s}$ (Lu et al., 1993), $\rm{PA_{0}}=26.4^{\circ}$; $i=21.5^{\circ}$ (Daigle et al., 2006), $(\alpha_{\rm{gal}}^{\rm{center}},\delta_{\rm{gal}}^{\rm{center}})=(01\rm{h}36\rm{m}41.747\rm{s},15\rm{d}47\rm{m}01.18\rm{s})$, and $(\alpha_{\rm{bubble}}^{\rm{center}},\delta_{\rm{bubble}}^{\rm{center}})=(1\rm{h}36\rm{m}44.6172\rm{s},15\rm{d}46\rm{m}21.4\rm{s})$.

We have fitted the sinusoidal model from Eq. 1 in terms of PA. We fitted the observed velocities using the CO data since it offered us the best combination of spatial and kinematical resolutions among the available data. The observed data follows a sinusoidal curve in the PA range of $-60^{\circ}$ to $160^{\circ}$, where the maximum and minimum values of observed velocities are clearly visible. The observed velocities outside this PA range are almost approximately constant, which would require a model more complex than the model we are fitting. This part where velocity is almost constant corresponds to the western arc where the bubble grazes the spiral arm. The encounter of the swept-up material with the spiral arm is expected to decelerate the expansion of the bubble, which is probably the reason for the almost constant velocity in this part of the bubble. The results of the fits performed for PA$=-60^{\circ}$ to PA$=160^{\circ}$ (darker points) are shown in the bottom panel of Fig. 7. We show two fits, one where we force the expansion velocity, $V_{\rm{exp}}$, to be zero (dashed line) and the other where we let it free (solid line). It can be clearly seen the fit does not replicate the observed velocities when not considering expansion.

From these data, we determine $V_{\rm exp}=-$11.9$\pm$0.7 km s*-1*, and $V_{\rm rot}=108\pm$2 km s*-1*. The negative value of the expansion means that the approaching side of the expansion is in the positive range of PA, the eastern side of the bubble. This fact is consistent with the trailing spiral arms in this galaxy.

3.4 Mass of the gas in the shell

Most of the gas in the shell is in molecular form. Using a CO(2-1) to H2 conversion factor $\alpha_{\rm{CO}}=6.1\thinspace(\rm{K\thinspace km/s})^{-1}$ corresponding to the Galactic CO(1-0) value (Bolatto et al., 2013) and a CO(2-1)/CO(1-0) ratio of 0.7 (Leroy et al., 2009), we determine a total molecular hydrogen mass of 1.52$\times 10^{7}$ M*⊙* in the shell between the 400 to 700 pc annular ring. The atomic mass in the same zone is 5.4$\times 10^{6}$ M*⊙, giving a total gas mass of 2.06$\times 10^{7}$ M⊙. The CO and H i intensities per pixel is below the 3-$\sigma$ detection limit in pixels inside the bubble. We estimate an upper limit of $6\times 10^{6}$ M⊙* inside a radius of 460 pc.

4 Stellar populations inside the bubble

A superbubble is created due to the mechanical power deposited to the interstellar medium from winds from massive stars and multiple supernova explosions. Superbubbles are short-lived structures and hence the relatively lower-mass stars of the stellar population that contained these massive stars should still be present inside the bubble. Stars need to be more massive than 8 M*⊙* to explode as supernova, which implies that there should be stars younger than 40 Myr inside the bubble. On the other hand, the stellar population could be as young as 4 Myr if it contained stars as massive as 100 M*⊙* and the bubble got created in one of the first explosions of supernova in the cluster. In the STARBURST99 models, a Simple Stellar Population (SSP) produces its maximum mechanical power between 4 to 20 Myr. We used the HST/ACS and JWST/NIRCam images of NGC 628 to trace and characterize the stellar population that might be responsible for the creation of the superbubble.

Magnitude at the turn-off point of 8 M*⊙* star is F814W=24 mag, which is above the detection limit of the HST observations and hence post-main sequence stars of clusters as old as 40 Myr can be detected. For younger clusters, stars in the upper end of the main sequence can also be detected.

In Figure 8, we display the images in the F435W and F814W filters of HST and F200W, F300W, F335W and F360W of the JWST/NIRCam filters. The HST images needed a small astrometric correction which we carried out by visually identifying stars common in the F200W JWST image. We superpose the two F770W contours defined in Figure 3 in all the images which enables us to identify the location of stars with respect to the bubble. In the HST images, we clearly observe an overdensity of stars inside the bubble. In Figure 9 we show RGB images, where the HST/F435W and JWST/F200W images are chosen as the blue and red components, respectively. For the green component, we choose JWST/F770W image in the left and MUSE/H$\alpha$ image on the right. The figures illustrate that the blue stars are selectively confined inside the bubble, whereas the red stars are almost uniformly distributed inside and outside the bubble. The blue stars are in general brighter and are detected in both the HST and JWST filters, whereas the majority of the JWST-detected stars are mostly undetected in the HST/F435W filter. We carry out aperture photometry of all stars inside the bubble to determine the age of the blue population stars.

4.1 Detection and photometry of stars in the bubble

We used the SExtractor (Bertin & Arnouts, 1996) code to make a catalogue of stars detected in the HST and JWST/NIRCam images. Given the different spatial resolutions and pixel scales of the images used here (see Table 1), we ran SExtractor independently on each image and then obtained a list of common stars in any two images based on the RA and DEC coordinate matches. Stars in the two catalogues whose coordinates differ by less than 0.15 arcsec are considered counterparts of each other. For each catalogued source, we saved the FWHM, the AREA above the 4-$\sigma$ threshold, and magnitudes in concentric apertures of multiple radii. We used the zeropoints given in Table 1 to convert the image counts to apparent magnitudes in the VEGA system.

We defined a circular area of 450 pc radius centered on the bubble to create lists of bubble stars in each filter. This list contains stars physically inside the bubble as well as the disk stars in the line of sight. The latter contribution can be statistically determined and subtracted by analysing the color-magnitude properties of the stars in the immediate vicinity of the bubble. The disk brightness varies considerably around the bubble, and hence it is necessary to take into account this variation in order to infer a possible excess population of stars that are physically inside the bubble. We defined four regions around the bubble to sample the typical projected disk population at the location of the bubble. These regions are shown in Figure 10.

Given that the F200W has the best spatial resolution among all the images, we used the catalogue in this filter as the base catalogue and obtained the JWST colours of these stars. We found that the majority of the JWST stars are not detected in the F435W images and are only marginally detected in the F814W images. We hence made a separate catalogue of stars detected in both the F435W and F814W filters. We use the magnitudes and colours in aperture of 0.15 arcsec radius for constructing colour-magnitude diagrams. The 3-$\sigma$ limiting magnitudes in these filters are given in Table 1.

4.2 Colour-magnitude diagrams of the bubble and disk stars

In Figure 11, we compare the stellar populations in the bubble (left) with that of a typical disk (right) in the immediate vicinity in a CMD using the HST magnitudes in the Vega systems. In both the graphs, we plot isochrones at selected ages using the PAdova and TRieste Stellar Evolution Code (PARSEC; Bressan et al., 2012). A comparison of the locations of stars in these plots reveals that most of the disk population contains stars of 50 Myr or older. A subset of stars in the region bounded by $-0.5\leq F435W-F814W\leq 1.5$ and $20\leq F814W\leq 24$ in the bubble CMD clearly stand out. We suspect that these stars, most of which are younger than 30 Myr, are responsible for the creation of the expanding bubble. We carryout below a quantitative analysis of the stellar populations in the bubble and the disk in order to determine the ages of this relatively younger stellar population in the bubble.

4.3 Stellar population analysis

We used the statistical model from Alzate et al. (2021) to compute the ages of the stellar populations in the bubble and disk samples. This approach uses the Bayes theorem in its hierarchical mode to infer the star formation history (SFH) of any resolved stellar population from the photometry of its stars. In this framework, the algorithm computes the probability distribution of the position of stars in the color-magnitude diagram as a linear combination of theoretical isochrones of distinct ages, each with a weight $a_{i}$, whose numerical value is fixed by the fraction of total number stars fitting a given isochrone (Small et al., 2013). The statistical approach we have followed here has advantages over the traditional isochrone fitting method, as this novel method obtains SFH rather than the best-fit age of the single population.

Each photometric magnitude from SExtractor is modeled as a random sampling from a statistically independent normal distribution function (or the likelihood), where the mean ($\mu$) of the distribution is the predicted magnitude from stellar models and its standard deviation ($\sigma$) is the magnitude error from SExtractor. The prior Probability Distribution Function (PDF), from which photometric magnitudes and stellar counts are predicted, is obtained using the Initial Mass Function (IMF) and the stellar isochrones, and it is probabilistically conditioned to $a_{i}$ through the ${\rm P}(f_{j}^{1}f_{j}^{2}f_{j}^{3}|\bm{a})$ term in Eq. 3. The product of the likelihood and the prior, properly normalized and computed for the $j^{\rm th}-$star of the sample, gives the probability that the magnitudes of an observed star in different filters were sampled from the prior set of isochrones. Tables 3 and 4 (see Appendix) show the parameters and the probability distributions entering in the statistical model. Using all these components and selecting a hyper-prior for $a_{i}$, the posterior PDF of stellar fractions is computed through

[TABLE]

where subscript $j$ indicates the $j^{\rm th}$ star and the superscript $k$ indicates the photometric magnitudes of the star; $(f^{1},f^{2},f^{3})$=F435W, F555W, F814W. In Eq. 3, $p$ and $\mathcal{L}$ stands for PDF and likelihood function, respectively. Thereby, $p(\bm{a})$ and $p(f_{j}^{1}f_{j}^{2}f_{j}^{3}|\bm{a})$ are the hyper-prior and the prior of the theorem. The normalization constant ($\ell$) of the posterior PDF will depend on the completeness of the sample. With an inappropriate normalization, the posterior will predict incorrect stellar counts, and the inference of the SFH will be biased. We considered two cases to deal with the incompleteness of any sample, computing the posterior PDF for: i) a complete sample ($S(F_{j}^{k})=1$) and ii) a magnitude-limited sample (star counts are complete only for stars brighter than the limiting magnitude, $S(F_{j}^{k})=1$ if and only if F555W$\leq$F555Wlim. In the second case, counts fainter than the limit will have negligible probability and do not contribute information to the solution. Once the components of the Bayes theorem are defined and restricted to the completeness intervals, we sampled the marginal posterior distributions of the stellar fractions performing a Montecarlo sampling of $10^{5}$ iterations. The $a_{i}$ median estimates the contribution of the $i$-th isochrone, and the (10%, 90%) percentiles delimit the credibility intervals of $a_{i}$.

We transform these $a_{i}$ values per interval of age (or per isochrone) to star formation rates (SFRs), assuming a stellar IMF. The number of stars per interval of age is given by

[TABLE]

where $N_{s}$ is the total number stars of the sample. For a magnitude-limited sample, we obtain $N_{i}$ as the number of stars brighter than the $F555W_{\rm lim}$, and the number stars fainter than this limit are unknown. Therefore, we calculate the number of stars below $F555W_{\rm lim}$ by integrating the double power law IMF:

[TABLE]

as

[TABLE]

where $M_{l}=0.1$ M*⊙* is the lower mass boundary of the IMF and $M_{\rm lim}(i)=M_{i}(F555W_{\rm lim})$ is the mass of the $i^{\rm th}$-isochrone at the magnitude limit. The IMF in Eq. 6 is normalized to $N_{i}$, such that,

[TABLE]

We use $M_{c}=0.5$ M*⊙*, slopes $\alpha_{1}=1.3$ and $\alpha_{2}=2.3$ according to Kroupa (2001). Thus, the corrected number of stars is

[TABLE]

Evolutionary effects are not taken into account in Eq. 8, that is, $N_{\rm corr}$ does not include evolved high mass stars, $[M_{\rm lim}(i),M_{u}]$, that are evolved off the CMD. Here $M_{u}=350$ M*⊙*  is the maximum mass that a single star may have in the IMF. However, these stars are fewer and their contribution to the number and the mass of the whole population is negligible.

In this theoretical framework, the expected mean mass of $N_{\rm corr}$ stars is

[TABLE]

and the star formation rate at the $i^{\rm th}$ age interval is

[TABLE]

where $t_{i}=13.7{\rm Gyr}-{\rm Age}$.

4.3.1 Isochrones

We acquired stellar isochrones from version 1.2S of PARSEC (Bressan et al., 2012). This version of the code computes solar-scaled evolutionary tracks for stars with masses from 0.1 M*⊙to 350 M⊙*, for metallicities ranging from Z=0.0001 to Z=0.06 (Chen et al., 2015). We used 23 isochrones111http://stev.oapd.inaf.it/cgi-bin/cmd selected with the following prescription

• •

$Z=Z_{\odot}$,

• •

$\log({\rm Age}/{\rm yr})=[6,8.2]$ dex,

• •

$\Delta\log({\rm Age}/{\rm yr})=0.1$ dex,

• •

$M=[0.1,350]$ M*⊙*.

We used PARSEC magnitudes in the HST/ACS photometric system anchored to Vega zero-point reference. We selected $\log({\rm Age}/{\rm yr})=8.2$ dex as the upper limit of our model grid, motivated by the fact that the turn-off locus of the isochrone at this age is around the limiting magnitude of the bubble and the background samples (F555W=27 mag)

4.3.2 Age of the bubble population

We analyse with our Bayesian code the star samples of the bubble and the four disk regions described in Sec. 4.1. Instead of cleaning disk stars in the bubble CMD, we infer the SFH of each sample independently, and then compare the relative differences in SFH between the bubble and the backgrounds. We found that the differences in the numbers of stars in the bg1, bg2 and bg4 samples are less than the expected Poissonian error, whereas the bg3 sample contains distinctly less number of stars. In order to increase the statistical significance of the results, we determine the SFH of the disk for a sample obtained by combining the bg1, bg2 and bg4 disk samples, which we refer to as bg124 sample henceforth, and dividing the resulting SFH by a factor of three, which we denote by $<{\rm bg124}>$. We correct the bubble and the $<$bg124$>$ samples for the foreground Galactic extinction, using $A_{V}=0.19$ mag (Schlafly & Finkbeiner, 2011).

The SFHs of the bubble and the $<$bg124$>$ samples are shown in Fig. 12. There is a clear excess SF in the bubble over the 50 Myr. In order to demonstrate this excess star formation, we subtract the normalized SFR obtained for the combined disk sample from that for the bubble sample, which is shown in the top panel of Figure 13. These plots illustrate that the SF started inside the bubble around 50 Myr ago, with the SFR decreasing steadily thereafter. The decrease in SF is more drastic over the last 10 Myr with no new stars formed in the most recent 4 Myr. The surrounding regions in the disk are also forming stars in the last 4–50 Myr, but at a rate smaller than that inside the bubble.

4.3.3 Mass of the bubble population

In the bottom panel of Figure 13 we show the cumulative mass of all the populations that contain an excess SFR in the bubble with respect to that in the surrounding disk, which is obtained by integrating the excess SFR in the top panel. The total mass of all stars formed inside the bubble over the last 50 Myr is $10^{5}$ M*⊙*, of which 35% was formed in the first 10 Myr (between 40 and 50 Myr ago). In comparison, only $\sim$20% of stellar mass is formed in the most recent 20 Myr years. The remaining 45% of stars were formed 20–40 Myr ago.

5 Discussion

5.1 The energetics of the superbubble

The creation and subsequent evolution of the superbubble by the correlated supernovae (SNe) explosions in disk galaxies has been studied in the classical work of Mac Low & McCray (1988, see also () and references therein). In the case of a homogeneous ambient gas distribution and a constant star cluster mechanical power $L_{\rm mech}$, the wind-driven superbubble radius $R_{\rm sb}$ and velocity $v_{\rm exp}$ at a time $t$ after the onset of expansion are given by,

[TABLE]

[TABLE]

where $\rho_{\rm gas}$ is the ambient gas density. Coefficient $A$ in Eq. 11 is:

[TABLE]

where \gamma{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=\frac{5}{3}}} is the adiabatic index. Eqn. 11 can be rewritten to solve for $L_{\rm mech}$ in terms of observable parameters $R_{\rm sb}$, $v_{\rm exp}$ and $M_{\rm sh}$, the mass of the star cluster wind-driven shell. To do that, we first eliminate $t$ in Eqn. 11 by using Eqn. 12:

[TABLE]

In the case of a homogeneous ambient medium $\rho_{\rm gas}=3M_{\rm sh}/4\pi R_{\rm sb}^{3}$. One can obtain then the cluster mechanical power as a function of the observed shell mass, radius and velocity:

[TABLE]

The value of the star cluster mechanical luminosity predicted by Eq. 15 could be confronted with that obtained from the analysis of the JWST images.

In the case when $R_{\rm sb}$=450 pc, $M_{\rm sh}=2\times 10^{7}$ M*⊙* and $v_{\rm exp}$=12 km s*-1*, the bubble age calculated using Eq. 12 is 22 Myr, which is consistent with the age of the enclosed stellar populations. The mechanical power required to create the bubble using Eq. 15 is $L_{\rm mech}=2.2\times 10^{38}$ erg s*-1*. The mechanical power injected into the ISM by the winds of O and WR stars and SN explosions in a stellar population of $1\times 10^{5}$ M*⊙* at solar metallicity is $L=5.0\times 10^{39}$ erg s*-1* (Leitherer et al., 1999), which is $\sim$20 times more than the mechanical power required to create the observed bubble in a homogeneous ISM.

We cannot exclude, however, the possibility that the bubble tops were accelerated and perforated via Rayleigh-Taylor instabilities, and that some fraction of the deposited energy has escaped into the galaxy halo as the 450 pc inner radius of the shell exceeds the typical scale height of the gas disk (see Mac Low et al., 1989; Tenorio-Tagle et al., 1990; Silich & Tenorio-Tagle, 1998; Lopez et al., 2011, and references therein). In this case, the swept-up shell expansion is maintained by the injected momentum rather than energy. This agrees with the fact that Chandra and XMM-Newton observations did not reveal a significant amount of diffuse X-ray emission from the discussed region (Soria & Kong, 2002). A more comprehensive bubble model that includes the disk-like ambient gas distribution and non-coeval star formation and a more realistic spatial distribution of stars within the bubble will be presented elsewhere.

5.2 Star formation scenario inside the bubble

Our analysis of the bubble stars in the CMD has indicated that the bubble contains stars that are formed over the last 50 Myr, which contribute to total stellar mass of $10^{5}$ M*⊙*. We here analyze the spatial distribution of stars of different ages in order to understand the star formation scenario inside the bubble since its formation.

In Figure 14, we show the distribution of the brightest of the stars inside the bubble in an RGB image formed by combining CO, H$\alpha$ and F435W, respectively as red, green and blue components. In addition, we show the compact sources detected in the F2100W images by the contours. There is one-to-one correspondence between the locations of these 21 $\mu$m sources and the compact H$\alpha$ sources. We identify five of the most prominent of these common sources in F2100W and H$\alpha$ images by letters A, B, C, D and E. The youngest stars are identified by circles of yellow ($<$10 Myr) and magenta (10–20 Myr) colours. We do not indicate stars older than 20 Myr in the image to avoid congestion of points.

If stars comprising each population that was inferred from the analysis of SFH belonged to star clusters, we would have expected agglomeration of stars of same age (i.e. circles of same colour in the figure). No such age-dependent spatial clustering of stars is easily discernible. Nevertheless, a closer inspection reveals that the 10–20 Myr old stars (magenta circles) are not distributed uniformly throughout the bubble, and instead majority of them are concentrated in three recognizable groups. Two of these groups are located at the southern edge of the cavity between the regions identified as C and D. Both these groups are surrounded by their own tiny H$\alpha$ bubbles of $\sim$100 pc radius. The third stellar subgroup is situated closer to the center of the bubble to the south of the region B. Additionally, there is a crowded distribution of stars older than 20 Myr towards the CO-bright arc to the north-west. The distribution suggests that the younger stars are preferentially located near the shell rather than at the bubble center.

In order to demonstrate an age-dependent spatial segregation of stars, we plot in Figure 15, the radial distance from the bubble center of each star associated to a particular isochrone. The stars associated with older isochrones are found at all radial distances, whereas stars younger than 12 Myr are preferentially found at larger radial distances. We binned the data into 12 radial bins, each containing 16 stars, and joined with a solid line the age corresponding to the 5th youngest star in each bin. The age of the youngest population is systematically lower at larger radii up to the ring radius of 450 pc. This suggests that the zone of active star formation is systematically moving away from the center as the bubble expands. It is most likely that the new stars are formed from the material that was swept up by the expanding bubble.

There is a clear evidence for such a triggered star formation further out in the shell, which marks the location of the swept-up gas, in the form of the sources identified by letters A to E. As described above, these sources are traced in H$\alpha$ and 21 $\mu$m, with the former emission suggesting the presence of an ionizing star at these locations, and the latter suggesting the presence of hot dust surrounding that star. Faint F435W sources are detected at the location of these knots. Unfortunately, only region B is captured within the area over which we analysed ages. We infer an age of $<$10 Myr for two of the HST stars coinciding with the peak of the F2100W emission. Thus, it seems that the F2100W sources are not heavily embedded, in which case no stars in the HST images, and the H$\alpha$ emission would have been detected at these locations. The newly formed stars are expected to get rid of the obscuring dust, and become undetectable in the F2100W filter in a few million years. All these characteristics suggest that the sources A to E trace the most recent star formation event and we assign them an upper age limit of 4 Myr. Correspondingly, we show these sources in the figure by downward-pointing triangles, each identified with their names, at fixed upper age limit of 4 Myr. These points continue the trend seen inside the bubble of star formation progressively shifting to larger radius with time.

A total of $\sim 10^{5}$ M*⊙* of stars have been formed inside the bubble over the last 50 Myr in several star formation episodes. Each episode is expected to be associated with one or more star clusters of size not exceeding $\sim$10 pc, the typical size of open clusters. The stellar groups discussed above are the only identifiable agglomerations, with most of the stars inside the bubble distributed over almost a kiloparsec scale. The logical explanation for this is that the stars were formed in clusters but the cluster got disrupted quickly. Clusters are prone to disruption when the gas in the cluster volume gets expelled due to the mechanical power input from the massive stars in the form of stellar winds and supernova explosions (e.g Krause et al., 2016). This starts to happen at $\sim$4 Myr. Cluster is expected to expand and survive if the expelled gas mass forms only a small fraction of the total mass within the cluster volume. The fraction of gas mass in a cluster depends on the star formation efficiency (SFE), with less residual gas for more efficient SF. According to hydrodynamical calculations of Krause et al. (2016), a $10^{5}$ M*⊙* cluster produces enough mechanical power to expel all its gas if it is not very compact (e.g. $r_{h}$=10 pc), even for large residual gas masses (or low SFE). It is likely that the progenitor cluster that was present at the center of the present-day bubble was not very compact and the binding energy of the cluster was not strong enough to counter the pressure exerted by its stars, following the expulsion of residual gas. If the cluster was in Virial equilibrium before it got dissolved, cluster stars at the most would have had a dispersion velocity $v_{\rm disp}$=10 km s*-1*.

The currently observed spatial extent of the dissolved cluster allows us to check the above hypothesis. After dissolution, cluster stars are expected to move radially outward with a mean velocity which should be equivalent to the dispersion velocity of the progenitor cluster. At the inferred $v_{\rm disp}$=10 km s*-1*, stars would reach 500 pc in 50 Myr. Thus those formed in the first episode had enough time to occupy the entire area of the present bubble. Stars that are formed at subsequent epochs formed at some location between the bubble center and the shell, in such a way that these stars also would have had time to reach the boundaries of the bubble. Thus, the large spatial extent of the recently-formed stars inside the bubble can be reconciled with the scenario of disruption driven by gas expulsion and subsequent movement of the stars with the velocities that they had when they were bound to the cluster.

While the scenario described above explains the reason for the distribution of stars over the kiloparsec scale, it does not satisfactorily explain the observed complete absence of stars outside the bubble, and the strict confinement of stars inside the bubble. We here propose an alternative scenario, known as the gravitational feedback, that consistently explains the cluster disruption, and the confinement of stars inside the bubble. The gravitational feedback is the process in which the gravitational potential well of the shell gas is deeper than that of the cluster, which not only destabilizes the cluster, but also accelerates the stars towards the shell. Zamora-Avilés et al. (2019) showed that this happens in expanding bubbles when the swept-up mass is orders of magnitude more than the cluster mass. They found that the “Hubble-flow” nature of the observed velocity fields of stars in some of the Galactic open clusters (see e.g. Kounkel et al., 2017; Román-Zúñiga et al., 2019) can be explained by the process. The gravitational force exerted by the mass collected in the expanding ring not only pulls apart the stars in an initially bound cluster, but also accelerates the stars towards the shell.

5.3 The conditions for the formation of kiloparsec-size bubbles

In the analysis carried out in this work, we comprehensively showed that the $\sim 10^{5}$ M*⊙* stars formed in multiple episodes over the last 50 Myr, and the feedback from these stars, are responsible for the observed characteristics of the bubble. Such star-forming events are common in galaxies unlike the presence of a kiloparsec size bubble. So, what is the most crucial ingredient for the creation of such large bubbles?

Vorobyov & Basu (2005) found that the bubbles created by stellar feedback from multiple star formation events sustain a kiloparsec-diameter organized bubble as compared to the one created by a massive cluster or exotic events like a hypernova explosion associated with a gamma ray burst (GRB; Efremov et al., 1998), and the vertical impact of a high velocity cloud (HVC; Tenorio-Tagle et al., 1987), all of the same total energy output. Their computations suggest that the location and mass of the first-formed cluster in the NGC628 bubble was the most crucial ingredient to create the observed large bubble. The SFR in the last bin in Figure 13 suggests a mass of at the most $3.5\times 10^{4}$ M*⊙* of stellar mass in the first cluster that was formed. The location of this cluster in the interarm region provided conditions for a uniform expansion of the bubble in its initial stages, which then formed new generations of stars in the shell that was formed from the swept-up material. Given the lower ISM density in the interarm region as compared to a spiral arm, the cluster formed in the swept-up gas mass was also of relatively smaller mass, thus allowing for a uniform, rather than a violent, expansion of the bubble. This process of steady expansion, collect and collapse, continued in each subsequent star formation event, each time forming moderate mass clusters. The bubble at the current time has reached the location of the spiral arm in the north-west side, where the bubble expansion seems to have been stalled. In the diagonally opposite direction, the bubble is growing asymmetrically and is less defined. We believe the bubble has reached the maximum size over which it can maintain a coherent structure.

6 Conclusions

We analysed the multiband data for a kiloparsec size bubble in the late type spiral galaxy NGC 628, to understand the creation and maintainance of large ISM bubbles. The bubble is the largest of many such structures seen in the JWST F770W images of this and other late-type galaxies, which has given them the popular name of “phantom galaxies”. The bubble seen in the F770W filter corresponds to the PAH molecules in the diffuse ISM, and is detected in other tracers of multi-phase ISM, such as in molecular gas by the CO emission, neutral atomic gas by the H i line and ionized gas by the H$\alpha$ emission. A kinematical analysis of the CO velocity field reveals that the bubble is currently expanding at a velocity of $\sim$12 km s*-1*. The mass of the gas in the shell around the bubble is $2\times 10^{7}$ M*⊙, which is mostly in molecular form. We find a clear evidence for the presence of a spatially resolved stellar population that is confined to the inner boundaries of bubble in the HST images. An isochrone-fitting analysis of these stars in the colour-magnitude diagram using a Bayesian technique, reveals that they are formed in multiple star formation events over the last 50 Myr. We find evidence for age-dependent spatial segregation of star formation with recent star formation happening closer to the inner shell of the bubble, rather than in its center, and on-going star formation at the shell itself. We propose a scenario wherein star formation started at the center of the present bubble with the formation of a star cluster of $10^{4}$ M⊙*  which got dissolved when the gas in the cluster volume got expelled due to the combined energy output of stellar winds and supernovae as the massive stars in the cluster ended their lives. The mechanical power created by these massive stars gave rise to an expanding bubble. The swept-up mass in the shell enclosing the expanding bubble got dense enough to form new generations of stars. This cycle of cluster formation, its destruction, subsequent creation of an expanding shell, and secondary star formation in the shell gas continued until the present day, in the process forming a large bubble. We propose that the location of the seed cluster in the interarm region, and its relatively low mass are the critical ingredients required to set in motion the creation of a kiloparsec size interstellar bubble.

Acknowledgements

We thank the PHANGS team for making the ALMA processed data publicly available. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2012.1.00650.S. We thank the referee for the comments that have contributed to improvements of the manuscript. We also thank CONACyT for the research grant CB-A1-S-25070(YDM), A1-S-28458 (SS), CB-A1-S-22784 (DRG) and CF2022-320152 (DFA). JAA acknowledges support from the IRyA-UNAM computer department and for granting him the required CPU time in the computers acquired through CONACyT grant CB2015-252364. LLN thanks Fundação de Amparoà Pesquisa do Estado do Rio de Janeiro (FAPERJ) for granting the postdoctoral research fellowship E-40/2021(280692).

7 Data availability

The paper is based on publicly available archival data. The processed fits files and data tables will be shared on reasonable request to the first author.

Appendix A Definition of parameters used in the Bayesian method of age analysis

In Table 3, we give the mathematical definitions of random variables entering the hierarchical model. Subscripts $i=1,2,...,N_{I}$ and $j=1,2,...,N_{D}$ indicate that the computation is made using the $i^{\rm th}$ isochrone and the $j^{\rm th}$ star, where $N_{I}$ and $N_{D}$ are the number of isochrones and observed stars, respectively. Superscript $k=1,2,3$ specifies the HST photometric magnitudes: $(f^{1},f^{2},f^{3})$=F435W, F555W, F814W. The model establishes that the pdf of the predicted absolute magnitude $f^{k}$ is a function of the isochrone magnitudes $\mathcal{F}^{k}(M)$ as follows

[TABLE]

taking the isochrone tolerance as constant, independent of $i$ and $k$. Distance modulus and visual extinction, dm and $A_{V}$ respectively, are known a priory and they are not variables.

Table 4, we give the equations for the PDFs. The hyper-prior $p(\bm{a})$ is the symmetric Dirichlet distribution, with concentration parameter $\xi$. When $\xi=1$, the hyper-prior is uniform, and no $a_{i}$ is preferred. See Sec. 3.4 of Gelman et al. (2013) for more details of the Dirichlet distribution, and Sec. 3.2 of Walmswell et al. (2013) for an example of its application within an astronomical context. Prior PDF is built as a mixture of model (isochrones) and its assemble the mass PDF and the $f^{k}$ PDF.