Magnetic Field of Molecular Gas Measured with the Velocity Gradient Technique. II. Curved Magnetic Field in kpc-scale Bubble of NGC 628

In this work, we select ¹²CO(2-1) emission as the main tracer to measure the magnetic field of the spiral galaxy NGC 628 with the VGT technique. The spectral cube data comes from the PHANGS-ALMA survey (Leroy et al. 2021), which includes 90 nearby galaxies (d ≤ 20 Mpc).

High signal-to-noise ratio (S/N) data with the rms noise level of ∼0.30 ± 0.13 K km s⁻¹. The spatial resolution for NGC 628 is ∼112 (FWHM ∼53.5 pc) and the velocity resolution is around 2.5 km s⁻¹ (Leroy et al. 2021). The rms level of the PHANGS-ALMA Survey is around 0.3 ± 0.17 K km s⁻¹.

In addition, the O iii emission in NGC 628 has been selected to trace the shocks, which comes from the PHANGS-MUSE Survey (Emsellem et al. 2022) with a subarcsecond beam scale size.

The continuum emissions of NGC 628 are observed by the James Webb Space Telescope (JWST) Early Release Observations (ERO; Pontoppidan et al. 2022). The filters F770W (7.7 μm), F1000W (10 μm), F1130W (11.3 μm), and F2100W (21 μm) from the JWST Mid-Infrared Instrument (MIRI) will be used in this work.

The synchrotron polarization of NGC 628 has been observed with the KarlG. Jansky Very Large Array (JVLA) at S-band (2.6–3.6 GHz effective bandwidth and spatial resolution ∼18''; Mulcahy et al. 2017). The orientation of the magnetic field θ_B derived from the JVLA polarization data is estimated by:

$$\begin{eqnarray}&&{\theta }_{B}=\displaystyle \frac{1}{2}\,\arctan (U,Q)\,+\displaystyle \frac{1}{2}\pi ,\end{eqnarray} \tag{ 1 }$$

where the Q and U are the Stokes parameters.

The velocity gradient technique (VGT; González-Casanova & Lazarian 2017; Hu et al. 2018; Lazarian & Yuen 2018) is the main analysis tool used in this work. This has been developed based on the current theories of MHD turbulence (Goldreich & Sridhar 1995) and fast turbulent reconnection (Lazarian & Vishniac 1999). These studies revealed that turbulent eddies are anisotropic, which means that the eddies are elongated along the local magnetic field lines (Goldreich & Sridhar 1995; Lazarian & Vishniac 1999). If the scale of the eddies is decomposed into parallel (l_∥) and perpendicular (l_⊥) components with respect to the magnetic field, the anisotropy suggests l_⊥ ≪ l_∥. This property has been confirmed with numerical simulations (Cho & Vishniac 2000; Maron & Goldreich 2001; Hu et al. 2021b) and with solar wind observations (Wang et al. 2016; Matteini et al. 2020; Duan et al. 2021). Another important property of MHD turbulence is that the velocity fluctuations are scale dependent such that (Lazarian & Vishniac 1999):

$$\begin{eqnarray}&&{v}_{l}\simeq {v}_{\mathrm{inj}}{\left(\displaystyle \frac{{l}_{\perp }}{{L}_{\mathrm{inj}}}\right)}^{\tfrac{1}{3}}{M}_{{\rm{A}}}^{1/3},\end{eqnarray} \tag{ 2 }$$

where v_l is the fluctuation at scale l, v_inj is the velocity at the injection scale L_inj, and M_A is Alfvén Mach number. Together with the anisotropic relation l_⊥ ≪ l_∥, the scaling of the velocity gradients can be easily obtained (Yuen & Lazarian 2020):

$$\begin{eqnarray}&&{\rm{\nabla }}{v}_{l}\propto \displaystyle \frac{{v}_{l}}{{l}_{\perp }}\simeq \displaystyle \frac{{v}_{\mathrm{inj}}}{{L}_{\mathrm{inj}}}{\left(\displaystyle \frac{{l}_{\perp }}{{L}_{\mathrm{inj}}}\right)}^{-\tfrac{2}{3}}{M}_{{\rm{A}}}^{\tfrac{1}{3}}.\end{eqnarray} \tag{ 3 }$$

Equation (3) suggests that the gradient amplitude increases when the scale l_⊥ decreases, i.e., the small resolved eddies correspond to strong velocity gradients. Therefore, in extragalactic sources, where the shear velocity and the differential rotation might contribute to the total velocity gradient, the contribution from the MHD turbulence would dominate at small scales (<100 pc), as demonstrated in Hu et al. (2022a) and Liu et al. (2023).

The input of the VGT technique is the high-resolution position–position–velocity (PPV) data cube of the ¹²CO emission in this work. The extraction of the velocity channel Ch(x, y) information from the PPV cube of the spectral line data has been done via the velocity caustics effect. The concept of the velocity caustics effect has been defined as the effect of the density structure distorting the turbulence and shear velocities along the line of sight (Lazarian & Pogosyan 2000). The density structure changes significantly for the different velocity channels and the resulting velocity fluctuations are thought to be most prominent in thin channel maps (Lazarian & Pogosyan 2000; Hu et al. 2021a). The difference for thin and thick channels is as follows:

$$\begin{eqnarray}&&{\rm{\Delta }}v\,\gt \,\sqrt{\delta ({v}^{2})},\mathrm{thick}\,\mathrm{channel},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}v\,\lt \,\sqrt{\delta ({v}^{2})},\mathrm{thin}\,\mathrm{channel},\end{eqnarray} \tag{ 5 }$$

where Δv is the velocity channel width and $\sqrt{\delta ({v}^{2})}$ is the velocity dispersion.

The thin velocity channel map Ch_i (x, y) is applied to the VGT technique (Hu et al. 2022b) to extract the velocity information from the PPV cubes via the velocity caustics effect. Each thin velocity channel map Ch_i (x, y) is used to calculate the pixelized gradient map ${\psi }_{g}^{i}(x,y)$:

$$\begin{eqnarray}&&{ \triangledown }_{x}{\mathrm{Ch}}_{i}(x,y)={\mathrm{Ch}}_{i}(x,y)-{\mathrm{Ch}}_{i}(x-1,y),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{ \triangledown }_{y}{\mathrm{Ch}}_{i}(x,y)={\mathrm{Ch}}_{i}(x,y)-{\mathrm{Ch}}_{i}(x,y-1),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\psi }_{g}^{i}(x,y)={\tan }^{-1}(\displaystyle \frac{{ \triangledown }_{y}{\mathrm{Ch}}_{i}(x,y)}{{ \triangledown }_{x}{\mathrm{Ch}}_{i}(x,y)},\end{eqnarray} \tag{ 8 }$$

where ▽_x Ch_i (x, y) and ▽_y Ch_i (x, y) are the x and y components of the gradient, respectively. The (x, y) 2D position will be displayed in the (R.A., decl.) plane of the sky for pixels whose spectral line emission has an S/N greater than 3, which is high enough for credibility.

The orientation of the magnetic field is perpendicular to the velocity gradient direction in each subregion. A subblock averaging (Yuen & Lazarian 2017) method has been used to export the velocity gradients from the raw gradients within a subblock of interest and then plotted in a corresponding histogram of the raw velocity gradient orientations ${\psi }_{g}^{i}$. The size of the subblock is set as 20 × 20 pixels, which determines the size of the final magnetic field resolution. Using subblock averaging, the total n_v processed gradient maps ${\psi }_{{g}_{s}}^{i}(x,y)$ with i = 1, 2,..., n_v are taken to calculate the pseudo-Stokes-parameters Q_g and U_g . Then, these pseudo-Stokes-parameters Q_g and U_g of the gradient-inferred magnetic field would be constructed by:

$$\begin{eqnarray}&&{Q}_{g}(x,y)=\displaystyle \sum _{i=1}^{{n}_{v}}{\mathrm{Ch}}_{i}(x,y)\cos (2{\psi }_{{g}_{s}}^{i}(x,y)),\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{U}_{g}(x,y)=\displaystyle \sum _{i=1}^{{n}_{v}}{\mathrm{Ch}}_{i}(x,y)\sin (2{\psi }_{{g}_{s}}^{i}(x,y)),\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{\psi }_{g}=\displaystyle \frac{1}{2}{\tan }^{-1}\displaystyle \frac{{U}_{g}}{{Q}_{g}},\end{eqnarray} \tag{ 11 }$$

where I_i (x, y) is the two dimensions integrated intensity of spectra cubes and n_v is the number of velocity channels. The ψ_g pseudo polarization is derived from the pseudo-Stokes-parameters Q_g and U_g , which projects the 3D ${\psi }_{{g}_{s}}^{i}$ data into a 2D pseudo polarization angle ψ_g . The pseudomagnetic field orientation ψ_B is perpendicular to the pseudo polarization angle ψ_g on POS:

$$\begin{eqnarray}&&{\psi }_{B}={\psi }_{g}+\pi /2.\end{eqnarray} \tag{ 12 }$$

This pseudomagnetic field orientation is the B-field orientation measured with the VGT method from the spectral line data. Applying the VGT technique on the ¹²CO emission from the PHANGS-ALMA survey (Leroy et al. 2021), the magnetic field orientations can be measured and be presented as an LIC (Cabral & Leedom 1993) as shown in Figure 1. The resolution of the B-field measured with VGT can be approaching 4'' (∼190 pc).

Zoom In Zoom Out Reset image size

The VGT technique is based on the position–position–velocity statistics (Lazarian & Pogosyan 2000), where the PPV cube has two components: a velocity and a density contribution. While the VGT technique is more sensitive to the velocity contribution, the Velocity Decomposition Algorithm (VDA) is a new technique that separates the velocity and density contribution from the PPV cube (Yuen et al. 2021). With this method, the accuracy of VGT tracing of the magnetic field will be improved (Zhao et al. 2022).

In this work, the gas traced by the CO emission may come from star formation regions tracing the spiral arms. These star formation regions may have supersonic motions and the properties of the velocity component V(X, v, Δv) may be calculated as:

$$\begin{eqnarray}&&V(X,v,{\rm{\Delta }}v)=-{c}_{s}^{2}\displaystyle \frac{\partial \mathrm{Ch}(X,v,{\rm{\Delta }}v)}{\partial v},\end{eqnarray} \tag{ 13 }$$

where Ch(X, v, Δv) is the channel of the PPV cube, X means the position, v is the local velocity, and Δv is the velocity channel width. The velocity of sound c_s can be calculated by assuming a uniform temperature ∼10 K, which results in a value c_s ∼186 m s⁻¹. Using these velocity properties p_v to replace the values for Ch(x, y) in Equations (6) and (7), this VGT-VDA method can trace the magnetic field better. The VDA algorithm exhibits a sensitivity to spectral lines characterized by a high S/N, particularly those present in structures of high density and small scale.

On smaller scales, the shocks induce velocity jumps that are perpendicular to shock fronts, which can interfere with the turbulence-induced gradients. While the effect of shocks would reduce with a decreasing velocity channel width, it should be noticeable for the CO data having a velocity channel width large enough at ∼2.5 km s⁻¹. At a shock front, the velocity perpendicular to the shock will be reduced while the velocity component parallel to the shock front will be unchanged. For oblique shocks, this changes the orientation of the velocity vector and the velocity gradient at the shock front (Hu et al. 2019). Since the anisotropy of the MHD turbulence does not affect the shock, a magnetic field component in the same direction as the propagation direction of the shock front will still suggest a velocity gradient perpendicular to the magnetic field, and a propagation direction of the shock perpendicular to the shock front. For a magnetic field that is perpendicular to the propagation direction of the shock front or has some angle, the relative orientation of the velocity gradient and magnetic field is not maintained at 90° (Xu et al. 2019). For a shock front that is parallel to the magnetic field, the direction of the magnetic field and the velocity gradient tend to be parallel. In cases where the shock front is oriented perpendicularly to the local magnetic field, the magnetic field that is unaffected by the local shock front aligns perpendicularly to the orientation of the mean velocity gradient (see Equation (12)). This physical process mirrors that observed in the VGT pipeline, as detailed in Section 3.2. Conversely, when the shock front aligns paralleling the local magnetic field, the magnetic field influenced by the shock aligns paralleling the mean velocity gradient orientation and perpendicularly to the pseudomagnetic field orientation

$$\begin{eqnarray}&&\mathrm{if}\,B\,//\,\mathrm{shockfront},{\psi }_{{Bs}}={\psi }_{B}+90^\circ ,\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\mathrm{if}\,B\,\perp \,\mathrm{shockfront},{\psi }_{{Bs}}={\psi }_{B},\end{eqnarray} \tag{ 15 }$$

where the ψ_B is the VGT orientation and ψ_Bs is the magnetic field orientation at the shock front. The shock fronts can thus be traced by the velocity field of the O iii emission. The local magnetic field orientation can be discerned from the polarization at a lower resolution, which encompasses the shock front.

The magnetic fields in the galaxy NGC 628 have been measured using the VGT using the high-resolution ¹²CO (2-1) spectral line emission. Using the VGT technique, the magnetic field is derived at a scale of 191 pc (∼4'') with the subblock size set at 20 × 20 pixels and with a pixel size of the ¹²CO emission of 02. The magnetic field structure of NGC 628 is shown in Figure 1.

On larger scales, the magnetic field is distributed along the spiral arms and the B-field orientations are nearly parallel to the tangential direction of the spiral arms. On smaller scales, the magnetic field displays curved structures around the spiral arms (ring-like shapes).

The JWST Early Release Observations (Pontoppidan et al. 2022) present details of the spiral galaxy NGC 628. The high-resolution mid-infrared emission from JWST, characterized by an FHWM of approximately 10 pc, presents a unique opportunity for investigating bubbles within extragalactic environments. This advanced capability allows for the observation of distinctive bubble features such as shells, stellar sources, and the hot emissions encapsulated by these bubbles (Churchwell et al. 2006; Anderson et al. 2014; Jayasinghe et al. 2019). In a recent study, Watkins et al. 2023 identified numerous bubbles within the extragalactic system NGC 628, demonstrating a broad range of scales with diameters spanning from 12 to 1164 pc. From the bubble catalog of Watkins et al. 2023, we have identified 21 sizable bubbles exceeding twice the beamwidth of the magnetic field measured with VGT (approximately 390 pc). Of these, one particular bubble situated at the galactic center does not align with the focus of our work; it may be influenced more by the galactic center's black hole rather than stellar feedback and star formation activities. Our attention centers on the remaining 20 large-scale bubbles within our sample, which are potentially linked to local star formation activities, as detailed in Table 1. These large-scale bubbles display the hole-like shape around the ring shape. The ring structure contains large amounts of gas and shocks traced by CO and O iii emission (see Figure 7). The scale of our structures is all larger than 390 pc, where 390 pc is twice the beam for measuring the magnetic field structure accurately. Details of the bubble structures are shown in Table 1.

Table 1. The Position and the Supernova Remnant Numbers of the Bubbles in NGC 628

ID	R.A.	Decl.	Size	Estimated SNR Number
			(pc)
1	1h36m4464	15d46m207	1164	∼194
2	1h36m4084	15d48m096	808	∼135
3	1h36m3595	15d48m061	418	∼70
4	1h36m3629	15d48m111	408	∼68
5	1h36m3739	15d48m014	520	∼87
6	1h36m3751	15d47m12s	460	∼77
7	1h36m3847	15d46m079	392	∼65
8	1h36m3863	15d46m5185	408	∼68
9	1h36m3933	15d48m085	762	∼127
10	1h36m3943	15d46m078	496	∼83
11	1h36m395	15d45m532	490	∼82
12	1h36m4071	15d48m357	394	∼66
13	1h36m4098	15d46m278	500	∼83
14	1h36m4225	15d48m098	484	∼81
15	1h36m4241	15d45m527	592	∼99
16	1h36m4281	15d46m471	588	∼98
17	1h36m4325	15d48m081	476	∼79
18	1h36m4364	15d46m384	512	∼85
19	1h36m4439	15d45m327	614	∼102
20	1h36m4549	15d46m358	510	∼85

Note. These are 20 large-scale bubbles. The scale of Bubble 1 is at one kpc and others are hundreds-pc scale. The ID of large-scale bubbles is shown in the first column. The position is shown in the second and third columns. The position and size of bubbles are defined from the JWST MIRI continuum data. The SN number shows the estimated number of SNRs required to make a bubble of this size (see Section 5.2).

Download table as:  ASCIITypeset image

The magnetic field measured with VGT for CO emission has a higher resolution (∼4''), compared with previously determined B-field structures derived from synchrotron polarization (∼18'') with JVLA at 3.1 GHz (Mulcahy et al. 2017). We smooth the magnetic field derived by VGT to the same scale as that from JVLA resolution (∼18'') and compare the two versions of the magnetic field. The alignment between the B-field orientations measured with polarization θ_B and VGT Φ_B is quantified as the Alignment Measure (AM; González-Casanova & Lazarian 2017):

$$\begin{eqnarray}&&\mathrm{AM}=2\left(\left\langle {\cos }^{2}({\theta }_{B}-{{\rm{\Phi }}}_{B})\right\rangle -\displaystyle \frac{1}{2}\right).\end{eqnarray} \tag{ 16 }$$

The range of AM values should be from −1 to 1, where AM values close to 1 mean that ϕ_B is parallel to ψ_g and AM values close to −1 indicate that ϕ_B is perpendicular to ψ_g . The uncertainty in the AM value σ_AM may be given by a standard deviation divided by the square root of the sample size.

4.3.1. The B-field in the Galactic Disk

Figure 2 shows the magnetic field structure measured with VGT for CO emissions and from synchrotron polarization with a beam size of 18''. The synchrotron radiation could originate from the warm ISM in star formation regions (Mulcahy et al. 2017), while the CO emission originates from the cold ISM in the disk (Hu et al. 2022a; Liu et al. 2023).

Zoom In Zoom Out Reset image size

Star-forming regions may be actively generating cosmic rays (CRs) and the cooling times of these CRs emitting at 3.1 GHz could be estimated as t_synch ∼ 1.5(B/mG)^−1.5(hv/kev)^−0.5 ∼ 13 Myrs, where the total magnetic field strength is around 10 μG (Mulcahy et al. 2017). The average diffusion time of ∼1 GeV CR protons in the Milky Way is around 10 Myrs for an assumed diffusion coefficient of 10²⁸ cm² s⁻¹ and a halo size of around 1 kpc (Ginzburg & Ptuskin 1976; Berezinskii et al. 1990; Strong & Moskalenko 1998; Evoli et al. 2019). However, in star formation regions and molecular clouds, the CR diffusion is suppressed as a result of the CR steaming along confusing turbulence magnetic fields (Xu & Lazarian 2022) or mirror diffusion (Lazarian & Xu 2021; Xu 2021), which is also observed in Milky Way. Assuming that the disk height is hundreds of pc, the diffusion time could thus be 1 or 2 orders of magnitude longer. Since the cooling time of CR electrons is less than the diffusion time, electrons at 3.1 GHz would not diffuse into the halo. Therefore, CR electrons are still actively generated but are confined in the vicinity of the galactic disk and the origin of the CO emission could be similar to that of synchrotron emission.

Figure 2 shows that the mean AM value is approaching 0.53 and is ≥0.5, which suggests that the magnetic field orientations from VGT are nearly parallel to those inferred from synchrotron polarization. The magnetic field derived by VGT is basically consistent with the previous results obtained from synchrotron polarization. Because the magnetic field measured with VGT for CO emission has a higher resolution in the galactic disk, the results may reflect the magnetic field at the spiral arms but also the small-scale structures around it.

4.3.2. The B-field at the Bubbles

The two largest bubbles, Bubble 1 and Bubble 2 (elaborated in Table 1), exceed the spatial resolution of previous synchrotron polarization observations (FWHM ∼ 800 pc; Mulcahy et al. 2017). These two large-scale bubbles offer exemplary instances for assessing the congruence between magnetic field orientation and VGT within minute-scale structures. The O iii emission from the PHANGS-MUSE Survey is the tracer of shocks and Figure 2 shows that the shock in Bubble 1 is close but does not overlap the gas shell of the bubble. Because of the shock front and the gas not being in the same position, the shock could not change the velocity gradient of the cold neutral medium traced by the CO emission. By comparison between the magnetic field for VGT and that for polarization, the mean AM is 0.70 (±0.01). This means that the magnetic field measured with VGT could be similar to that from polarization. In the no-shock region, the velocity gradient remains perpendicular to the magnetic field.

Different from Bubble 1, a shock exists in the gas shell of Bubble 2. As Figure 2 shows, part of the bubble shell in Bubble 2 coincides with a shock front and the region has massive O iii emission. The O iii emission may be smoothed to trace the shock front at a scale of 18'', where the intensity is over 3 × 10⁻¹⁶ cm⁻² erg s⁻¹. Here the shape of the shock front is parallel to the magnetic field derived by polarization. The VGT orientation in this region is rerotated by 90° to achieve the B-field orientation, as dictated by Equation (14). By rerotating the VGT orientation at the shock front, the AM between two measures of magnetic fields is 0.54, and the two types of magnetic fields are parallel to one another. In the shock front, the accuracy of VGT tracing magnetic field has been tested and the velocity gradient is close to being parallel to the magnetic field. In the region of the shock front, the VGT orientation needs to be rerotated 90°.

For the two largest bubbles, we also use the VDA method to improve the accuracy of the VGT technique by removing the density contribution in the PPV data cube. As Figure 3 shows, we compare the magnetic field measured with VGT-VDA and 3.1 GHz polarization. The AM values of Bubble 1 and Bubble 2 are up to 0.74 and 0.64, which is larger than the AM of VGT (Bubble 1 = 0.7 and Bubble 2 = 0.54). This means that the VGT-VDA has higher accuracy in tracing magnetic fields than only the VGT method, but it also verifies that the magnetic fields measured with VGT for the CO emission. It also means that we can reliably use the VGT technique to study the small-scale structures in NGC 628.

Zoom In Zoom Out Reset image size

The bubble structures appear rather common around the inner spiral arms of NGC 628 as seen in Figure 1. To establish the nature of the structure of these large-scale bubbles and the role of the magnetic field, the VGT-VDA method has been used to study the effects of the high-resolution magnetic field on large-scale bubble structures. This method improves the accuracy of the VGT method as the application of VDA removes the effect of the density contribution to the velocity gradient, which is particularly important at the edges of the structures.

Supernova remnants (SNRs, Tenorio-Tagle & Bodenheimer 1988; Weiler & Sramek 1988) and prolonged star formation are the most likely cause for the creation of a cavity within the circular bubble structures. Supernova explosions could have enough energy to cause large-scale bubbles with scales up to hundreds of pc or even kpc. These bubbles are found to be surrounded by gas structures with much-triggered star formation activities as traced by 11.3 and 21 μm emissions. This makes the bubbles to be relics of past star formation that only now show up at mid-infrared frequencies (see Figure 1).

The kpc-scale Bubble 1 in Figure 3 may be used as a high-confidence example of a bubble, where the shape and magnetic field structure could be caused by SNRs. The O iii emission serves as a valuable indicator for tracing the impact of shocks, as shown in Figure 2. At the edge of Bubble 1, the shock front traced by the O iii emission is close to the bubble's shell but it does not overlap with the spiral arms traced by the CO emission. The magnetic field morphology in the curved structure of Bubble 1 is like a closed ring that is aligned around the shape of the bubble. Due to the shock front being close to the shell at the edge of the bubbles, the supernovae create the shock and compress the gas to cause the bubbles and the shell. As a result of the gas flow toward the shell, the magnetic field has been distorted and forced to be perpendicular to the direction of matter flow (see Figure 5). The curved magnetic field structure therefore represents the shocked boundaries of the large-scale bubble and serves as indirect evidence that the bubble structure and its curved magnetic field would be caused by the SNRs.

Alternative star formation processes, such as stellar feedback and the expansion of H ii regions, may exhibit a greater inclination toward instigating the smaller-scale structure of ISM. These formations might lack the requisite energy to give rise to structures on a scale of hundreds of pc. Similar phenomena may also manifest within subregions of the large-scale bubbles. As an illustration, the orientations of magnetic fields in the subregion of Bubble 2 (see Figure 4) are observed to be perpendicular to the gas structures in subregions with existing star formation as traced by 21 μm emissions. It is noteworthy that local magnetic fields may undergo distortion due to the influence of star formation activities in this specific region, a phenomenon akin to those observed within the Milky Way (Li et al. 2017; Zhao et al. 2024).

Zoom In Zoom Out Reset image size

In this paper, only the large-scale bubbles have been considered and their magnetic field structures have been determined (see Figure 5). A more complete sample of bubbles in NGC 628 includes 1694 bubbles with different sizes and is provided by Watkins et al. (2023). A histogram of the size distribution and the galactocentric distance of the bubbles in this extended sample is presented in Figure 5 and shows a peak in the extended distribution at a bubble size of 50 pc and at a galactocentric distance of 3 kpc. In the galactic distribution of bubbles, the number of bubbles increases nearly linearly up to the galactocentric distance of 3 kpc, and then they decrease sharply and nonlinearly. The large-scale bubbles considered are thus mainly distributed in the inner part of the galaxy within 3 kpc.

Zoom In Zoom Out Reset image size

The bubble number distribution along the galactocentric distance is similar to the star formation rate distribution in the Milky Way (Elia et al. 2022). The bubble scale of ∼50 pc could be typical for NGC 628 and is similar to the largest SNRs found in the Milky Way (Ostriker & McKee 1988; Becker et al. 2021). However, this typical bubble size in NGC 628 is much larger than the typical size of ∼6 pc for most SNRs in the Milky Way and it is unknown how many SNs actually contributed to this large size. The peak radial position of the bubble distribution exists at about 1/6 to 1/5 of the disk size and suggests a ring structure with an enhanced star formation history. The existence of this ring requires further investigation because of its implications for the radial density structure and the rotation curve in the disk, as well as for the location of the Inner Lindblad resonance in the galaxy.

The large-scale bubbles could be caused by the superposition of multiple supernova explosions for three reasons:

• 1.  
  The shock front traced by the O iii emissions is distributed around the shell of the bubbles (see Figure 2).
• 2.  
  There is triggered star formation at the edges of the bubbles, which is traced by the 11.3 and 21 μm emission (see Figures 1 and 3).
• 3.  
  The curved magnetic field structure is aligned with the structure of the bubbles.

The size of the large-scale bubbles (>300 pc) considered here clearly exceeds the size of even the largest SNR in the Milky Way (∼ 50 pc). Such large-scale bubble structures may only result from massive star formation and the contribution of multiple supernova explosions. The shock wave from an SN could compress the surrounding gas and merge with shells from many others (Watkins et al. 2023) to form a large-scale bubble with a surface close to spherical. Depending on the thickness of the disk, the bubble cavity will become cylindrical and form a chimney outflow perpendicular to the disk as seen in neutral hydrogen surveys. The diffuse gas in the disk would be compressed in a circular shell that is also constrained by the embedded local magnetic field and this enriched shell may lead to triggered star formation. A situation similar to that of NGC 628 is found, for instance, in the galaxy NGC 6946 (Heald 2012).

The radio emission from SNRs results from high-energy electrons accelerated during a supernova explosion and ends up spiraling along the magnetic field structure (Green 2019; Cotton et al. 2024). Mulcahy et al. (2017) provide a high-resolution synchrotron radiation map of NGC 628 at 3.1 GHz with a resolution of up to 8'' (FWHM ≈ 390 pc; twice of B-field beam in this work). As Figure 6 shows, the synchrotron emission is aligned with the spatial distribution of bubbles from Mulcahy et al. (2017), which includes both the large-scale bubbles (FWHM > 390 pc, the resolution of synchrotron radiation) and the small-scale bubbles (FWHM < 390 pc). These bubbles could be candidates for SNRs. The regions with richer synchrotron emissions have massive small-scale bubbles compared with the large-scale bubbles. The small-scale bubbles are located at the spiral arms with rich gas. In the spiral arms, despite impactful physical phenomena like supernova explosions, it proves challenging to propel massive amounts of gas and facilitate the formation of large-scale structures. On the other hand, the diffuse gas surrounding the large-scale bubbles could be pushed by supernovae over greater lengths and form the large-scale structure of bubbles.

Zoom In Zoom Out Reset image size

The diffuse gas and limited synchrotron emissions surrounding the expansive bubbles on a large scale imply that the number of supernovae contributing to the formation of such bubbles is comparatively modest when juxtaposed with equivalent areas in spiral arm regions. Assuming that the aggregate flux, I, of the large-scale bubbles and their scale size, R, follow a power-law relation:

$$\begin{eqnarray}&&I\propto {R}^{\alpha }\propto {\left(\displaystyle \frac{R}{{R}_{0}}\right)}^{\alpha },\end{eqnarray} \tag{ 17 }$$

where R₀ is the characteristic scale of one SNR, and α signifies the slope between R and I in log–log space. When the flux within the bubble is uniform, the value of α becomes 2 (since the area is proportional to the square of the radius). By fitting the total flux and scale of each large-scale bubble, the derived alpha converges around 1 (less than 2). The growth rate of the flux in the bubbles is below that of the growth in scale. As the bubble's scale expands, the synchrotron emissions from SNR may undergo dilution.

We adopt a maximum size of a single SNR in the Milky Way (∼6 pc) as the characteristic scale populating the bubbles. The SNR number of large-scale bubbles can be estimated by the scale R of the bubbles:

$$\begin{eqnarray}&&N\ (\mathrm{SNR})=k(R/{R}_{0})\propto (I/{I}_{0})\end{eqnarray} \tag{ 18 }$$

Because of the unknown characteristic flux I₀, we assume the factor k as 1. Then the largest bubble, Bubble 1 (scale ∼1164 pc) contains an estimated 194 SNRs. Details regarding the number of SNRs in other large-scale bubbles can be found in Table 1. The observed large-scale bubbles in NGC 628 provide a good picture of the succession of star formation and its history.