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

We report the detection of the ordered alignment between the magnetic field and kpc-scale bubbles in the nearby spiral galaxy, NGC 628. Applying the Velocity Gradient Technique on CO spectroscopic data from the ALMA-PHANGS, the magnetic field of NGC 628 is measured at the scale of 191 pc (∼4″). The large-scale magnetic field is oriented parallel to the spiral arms and curves around the galactic bubble structures in the mid-infrared emission observed by the James Webb Space Telescope. A total of 21 bubble structures have been identified at the edges of spiral arms with scales over 300 pc, which includes two kpc-scale structures. These bubbles are caused by supernova remnants and prolonged star formation and are similar to the outflow chimneys found in neutral hydrogen in galactic disks. At the edge of the bubbles, the shocks traced by the O iii emission present a curved magnetic field that parallels the bubble’s shell. The magnetic field follows the bubble expansion and binds the gas in the shell to trigger further star formation. By analyzing the larger sample of 1694 bubbles, we found a distinct radial-size distribution of bubbles in NGC 628 indicating the star formation history in the galaxy.


INTRODUCTION
Magnetic fields are a primary component of the interstellar media (ISM) that play an important role in star formation, and the evolution of spiral galaxies (Mc-Kee & Ostriker 2007;Beck & Wielebinski 2013).Magnetic fields have been observed over the past decades at all spatial scales in the Galaxy and in extra-galactic sources using far-infrared and radio observations (Kronberg 1994;Beck & Gaensler 2004;Beck 2015;Lopez-Rodriguez et al. 2020;Jones et al. 2020).Shocks, su-Corresponding author: Mengke Zhao, Jianjun Zhou mkzhao628@gmail.com,zhoujj@xao.ac.cn pernovae, star formation activities, and other physical processes may affect and distort the local magnetic field structures in the ISM (eg: Bell 2004;Li et al. 2017;Arzoumanian et al. 2021 and the structure of spiral arms and spurs (Elmegreen 1980;Balbus & Cowie 1985;Lee & Shu 2012;Lee 2014).
The recent introduction of the Velocity Gradient Technique (VGT; González-Casanova & Lazarian 2017; Lazarian & Yuen 2018;Hu et al. 2018) provides a new way to measure the magnetic field using Doppler-shifted emission lines in nearby galaxies (Tram et al. 2023;Hu et al. 2022a;Liu et al. 2023).In a turbulent magnetized medium, the turbulent eddies will extend along the local B-field and have a velocity gradient perpendic-Zhao et al ular to their semi-major axis (Goldreich & Sridhar 1995;Lazarian & Vishniac 1999).These velocity gradients of elongated eddies are expected to be perpendicular to the local B-field orientation.The PHANGS-ALMA survey provides high-resolution CO isotopolog data where anisotropies of the turbulence can be extracted and used to trace the velocity gradients of eddies.The ability of the VGT technique to trace molecular-gas-associated magnetic fields has been tested in observations in nearby galaxies (Hu et al. 2022a;Liu et al. 2023).
Early Release Observations (ERO; Pontoppidan et al. 2022) of the spiral galaxy NGC 628 with the James Webb Space Telescope (JWST) reveal new features compared with previous observations in the form of 'bubble structures' along spiral arms that have no mid-infrared continuum emission.Such features resemble the chimney outflows found in neutral hydrogen related to star formation in edge-on spiral galaxies.Combined with the observation of the magnetic field, these bubbles provide new understanding of how the magnetic field may affect the structure of the spiral galaxy.
In this work, we aim to study the magnetic field in and around the bubbles structures of the spiral galaxy NGC 628 at the distance of 7.3 Mpc (Karachentsev et al. 2004), which exhibits ongoing star formation but is not undergoing a central starburst (Luisi et al. 2018).The details of the archival data are shown in Section 2. In Section 3, we introduce the details of the VGT and VDA (Velocity Decomposition Algorithm; Yuen et al. 2021) techniques.The magnetic field structure and galactic bubbles in NGC 628 are described and discussed in Sections 3 and 4. We discuss the nature of the bubbles and their magnetic field structure in Section 5. A summary has been provided in Section 6.

The spectral line emission
In this work, we select 12 CO (2-1) emission as the main tracer to measure the magnetic field of the spiral galaxy NGC628 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 data with the RMS noise level of ∼ 0.30 ± 0.13 K km s −1 .The spatial resolution for NGC 628 is ∼ 1.12 ′′ (FWHM ∼ 53.5 pc) and the velocity resolution is around 2.5 km s −1 (Leroy et al. 2021).The rms level of the PHANGS-ALMA Survey is around 0.3 ± 0.17 K km s −1 .
In addition, the OIII emission in NGC 628 has been selected to trace the shocks, which comes from the PHANGS-MUSE Survey (Emsellem et al. 2022) with a sub-arcsec beam scale size.

The polarization data
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 are estimated by: where the Q and U are the Stokes parameters.

The anisotropy of the MHD turbulence
The velocity gradient technique (VGT; González-Casanova & Lazarian 2017; Lazarian & Yuen 2018;Hu et al. 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 re-connection (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): 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): A . (3) Eq. 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); Liu et al. (2023).

The Velocity Gradient Technique
The input of the VGT technique is the high-resolution position-position-velocity (PPV) data cube of the 12 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 channel 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: where ∆v is the velocity channel width and δ(v 2 ) is the velocity dispersion.
The thin velocity channel map Ch i (x,y) is applied on 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 pixelised gradient map ψ i g (x, y): 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) twodimensional position will be displayed in the (RA,DEC) plane of the sky for pixels whose spectral line emission Zhao et al has a signal-to-noise ratio greater than 3, which is high enough for credibility.The orientation of the magnetic field is perpendicular to the velocity gradient direction in each sub-region.A sub-block averaging (Yuen & Lazarian 2017) method has been used to export the velocity gradients from the raw gradients within a sub-block of interest and then plotted in a corresponding histogram of the raw velocity gradient orientations ψ i g .The size of the sub-block is set as 20 × 20 pixels, which determines the size of the final magnetic field resolution.Using sub-block averaging, the total n v processed gradient maps ψ i gs (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: Ch i (x, y)sin(2ψ i gs (x, y)), ( 10) 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 three-dimensional ψ i gs data into a two-dimensional pseudo polarization angle ψ g .The pseudo magnetic field orientation ψ B is perpendicular to the pseudo polarization angle ψ g on POS: This pseudo magnetic field orientation is the B-field orientation measured with the VGT method from the spectral line data.Applying the VGT technique on the 12 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).

The Velocity Decomposition Algorithm
The VGT technique is based on the position-positionvelocity 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),.By 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: 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 −1 .Using these velocity properties p v to replace the values for Ch(x,y) in Eq. 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 signal-to-noise ratio, particularly those present in structures of high density and small scale.

VGT method in Shock Fronts
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 −1 .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 degrees (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 Eq. 12).This physical 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 1st column.The position is shown in the 2nd and 3rd 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).
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 pseudo-magnetic field orientation 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 OIII emission.The local magnetic field orientation can be discerned from the polarization at lower resolution, which encompasses the shock front.

The magnetic fields in NGC 628
The magnetic fields in the galaxy NGC 628 have been measured using the Velocity Gradient Technique using the high-resolution 12 CO (2-1) spectral line emission.Using the VGT technique, the magnetic field is derived at a scale of 191 pc (∼ 4 ′′ ) with the sub-block size set at 20×20 pixels and with a pixel size of the 12 CO emission of 0.2 ′′ .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).

Bubbles in extragalactic sources
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 a Full Width at Half Maximum (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 beam width 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 largescale 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 Fig. 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.

Comparison of the B-fields from VGT and from
Synchrotron Polarization 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): 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.

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).
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 28 cm 2 s −1 and a halo size 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 one or two 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 disc 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.

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 OIII emission from 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 noshock region, the velocity gradient remains perpendicular to the magnetic field.
Different from Bubble 1, a shock exists at 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 OIII emission.The OIII emission may be smoothed to trace the shock front at a scale of 18 ′′ , where the intensity is over 3×10 −16 cm −2 erg s −1 .Here the shape of the shock front is parallel to the magnetic field derived by polarization.The VGT orientation in this region is re-rotated by 90 degrees to achieve the Bfield orientation, as dictated by Eq. 14.By re-rotating 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 re-rotated 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 Fig. 3 shows, we compare the magnetic field measured with VGT-VDA and 3.1 GHz polarization.The AM value of Bubble 1 and Bubble 2 is 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.

Large-scale Bubbles as Supernova Remnants
The bubble structures appear rather common around the inner spiral arms of NGC 628 as seen in Fig. 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 µm 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 Fig. 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 OIII emission serves as a valuable indicator for tracing the impact of shocks, as shown in Fig. 2. At the edge of Bubble 1, the shock front traced by the OIII 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 towards the shell, the magnetic field has been distorted and forced to be perpendicular to the direction of matter flow (see Fig. 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 HII 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 sub-regions of the large-scale bubbles.As an illustration, the orientations of magnetic fields in the sub-region of Bubble 2 (see Fig. 4) are observed to be perpendicular to the gas structures in sub-regions 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).

Bubble Size and Galactic Distribution
In this paper, only the large-scale bubbles have been considered and their magnetic field structures have been determined (see Fig. 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 galacto-centric distance of the bubbles in this extended sample is presented in Fig. 5 and shows a peak in the extended distribution at a bubble size of 50 pc and at a galacto-centric distance of 3 kpc.In the galactic distribution of bubbles, the number of bubbles increases nearly linearly up to the galacto-centric distance as 3 kpc, and then they decrease sharply and non-linearly.The large-scale bubbles considered are thus mainly distributed in the inner part of the galaxy within 3 kpc.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 supernova remnants 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 origin of large-scale bubbles
The large-scale bubbles could be caused by the superposition of multiple supernova explosions for three reasons: 1.The shock front traced by the OIII emissions is distributed around the shell of the bubbles (see Fig. 2).2. There is triggered star formation at the edges of the bubbles, which is traced by the 11.3µm and 21µm emission (see Fig. 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 highenergy electrons accelerated during a supernova explosion and end up spiraling along the magnetic field structure (Green 2019;Cotton et al. 2024).Mulcahy et al. 2017 provides a high-resolution synchrotron radiation map of NGC 628 at 3.1 GHz with resolution 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 of SuperNova Remnants.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 largescale 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.
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: where R 0 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: Because of the unknown characteristic flux I 0 , 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 largescale bubbles can be found in Tab. 1.The observed large-scale bubbles in NGC 628 provide a good picture of the succession of star formation and its history.

SUMMARY
The magnetic fields of NGC 628 has been measured with the VGT technique using the 12 CO (2-1) emission at the scale of 4 ′′ (FWHM ∼ 191 pc).The magnetic field is in alignment with the structure of the bubbles in NGC 628 and forms a surrounding ring-like structure leaving a central cavity at the Mid-infrared wavelengths.
1.The high-resolution magnetic field measured with VGT for the CO emission is ordered and distributed along the spiral arms at the large scale.The result is Zhao et al in agreement with the preview observations of the synchrotron polarization at the low resolution of 18 ′′ .For smaller-scale bubble sources, the VGT-VDA method has been used to trace the magnetic field and improve accuracy, and the magnetic fields display curved structures.
2. 20 large-scale bubbles (390 pc -1164 pc) in NGC 628 display a hole-like shape surrounded by a ring structure, which includes rich gas reservoirs and shock wave trails.These large-scale bubbles are all distributed on the edge of spiral arms and are the result of multiple supernova activities.
3. The distribution of bubble size versus galactocentric distance of a large sample of 1694 bubbles displays a peak for a size around 50 pc at a galactocentric distance of 3 kpc.The number of bubbles in the galactic distribution increases linearly up to 3 kpc and rapidly decreases beyond this distance.The distribution of bubbles in NGC 628 forms a ring structure that is densely filled with bubbles and reveals the intense star formation history in the central region of NGC 628.The distribution of bubbles is similar to the star formation rate distribution in Milky Way.The bubbles observed in NGC 628 have been observed to be associated with outflow chimneys in neutral hydrogen in other galaxies.
4. The apparent ring of 50 pc bubble structures at about 3 kpc that surrounds regions of intense star formation in the inner part of the disk in NGC 628 has implications for the radial density structure the rotation curve in the disk, and the location of the Inner Lindblad resonance in the galaxy.
5. The bubble structure appears to be caused by repeated supernova explosions and prolonged star formation, which squeezes gas towards the edge of the bubbles and triggers new star formation in these compressed regions.The organization of the magnetic field in the shock boundary surrounding the bubble structures will be perpendicular to the propagation direction of the SNR shocks and parallel to the shell of OIII gas.While these bubbles are relics of past star formation, they continue to trigger new star formation.
6.The large-scale bubbles could be caused by the contribution of massive/multiple supernova explosions.
The multiple supernovae drive the diffuse gas over a large length and cause the large-scale bubbles.Scaling of the power needed to form such large-scale bubbles suggests contributions from large numbers (up to hundreds) of single SNR as observed in the Milky Way.The number and size of the bubbles observed in NGC 628 reveal the long-term star formation history in the galaxy, which is a common evolutionary process in galaxies.
Further plan are to observe the magnetic field using dust polarization at sub-arcsec resolution in order to study the small-scale magnetic field structure of the large-scale bubbles.Analysis of a face-on galaxy will promote understanding of similar physical processes in the Milky Way.

Figure 1 .
Figure1.The curved magnetic field aligns with the bubble structures in NGC 628.This map displays the magnetic field structure measured with VGT for the CO emission of the nearby spiral galaxy NGC 628, where the magnetic field orientation is indicated as a line-integral-convolution (LIC) map.The background is the RGB image at the mid-infrared with red: F1130W; green: F1000W; blue: F770W.The cyan circle is the one-kpc-scale bubble.The white circles show the position of 20 bubbles with scale sizes > 390 pc.

Figure 2 .
Figure 2. The magnetic field distribution measured with VGT and taken from the JVLA synchrotron polarization.The red and black vectors in the left panel display the magnetic field orientation measured with VGT and the JVLA polarization at 3.1 GHz (Mulcahy et al. 2017), respectively.The background is the 12 CO (2-1) emission intensity map.The mean value for the Alignment Measure AM is 0.53 with an uncertainty of 0.01.The middle panel shows the distribution of the magnetic field at the two kpc-scale bubbles (see Fig. 1) associated with the white boxes in the left panel.The cyan contours show the flux of OIII emissions at [4, 5, 6]×10 −16 cm −2 erg s −1 as determined from the PHANGS-MUSE survey data (Emsellem et al. 2022).The OIII lines overlap the large region of Bubble 2 and are located in a small part of Bubble 1 The orange contour shows the shock front identified by smoothing the OIII emissions to 18 ′′ .The right panels display the comparison of magnetic field orientations between the VGT-VDA and JVLA polarization at 3.1 GHz.The yellow vectors show the magnetic field orientation measured by VGT-VDA.

Figure 3 .
Figure 3.The magnetic field structure measured with VGT-VDA for kpc-scale Bubble in NGC 628.The panels display the curved magnetic field at bubble 1.The beam size is about 4 ′′ and is displayed as a cyan circle in the left bottom corner.The image background is a three color map in the mid-infrared with red: F770W, green: F1130W, and red: F2100W.The yellow vectors represent the magnetic field orientations measured by the VGT-VDA technique for 12 CO.The white contours are the contour levels of the 12 CO integrated intensity at 5, 10, and 15σ.The cyan line shows the bubble structure.

Figure 5 .
Figure 5.The distribution of the size of the bubbles versus galactocentric distance to the galaxy center using a probability density function and a 2D histogram.

Figure 6 .
Figure 6.Distribution of bubbles aligned with synchrotron radiation at 3.1 Ghz.The main panel shows the distribution of bubbles in NGC 628(Watkins et al. 2023).The black rings display the large-scale bubbles (FWHM > 390 pc) and the red rings show the small-scale bubble(FWHM < 390 pc) in this work.The background is mapping the 3.1 GHz emissions at the scale of 8 ′′(Mulcahy et al. 2017).The white dash line displays the intensity map of 12 CO (2-1) emissions.The sub-panel shows the distribution between scale R and flux I of large-scale bubbles.The red line is the fitting Eq. 17.

Table 1 .
The position and the supernova remnant numbers of the bubbles in NGC 628