Magnetic Fields and Velocity Gradients in L1551: The Role of Stellar Feedback

Magnetic fields play a crucial role in star formation, yet tracing them becomes particularly challenging, especially in the presence of outflow feedback in protostellar systems. We targeted the star-forming region L1551, notable for its apparent outflows, to investigate the magnetic fields. These fields were probed using polarimetry observations from the Planck satellite at 353 GHZ/849 $\mu$m, the SOFIA/HAWC+ measurement at 214 $\mu$m, and the JCMT/SCUPOL 850 $\mu$m survey. Consistently, all three measurements show that the magnetic fields twist towards the protostar IRS 5. Additionally, we utilized the Velocity Gradients Technique (VGT) on the $^{12}$CO (J = 1-0) emission data to distinguish the magnetic fields directly associated with the protostellar outflows. These were then compared with the polarization results. Notably, in the outskirts of the region, these measurements generally align. However, as one approaches the center of IRS 5, the measurements tend to yield mostly perpendicular relative orientations. This suggests that the outflows might be dynamically significant from a scale of approximately $\sim0.2$~pc, causing the velocity gradient to change direction by 90 degrees. Furthermore, we discovered that the polarization fraction $p$ and the total intensity $I$ correlate as $p \propto I^{-\alpha}$. Specifically, $\alpha$ is approximately $1.044\pm0.06$ for SCUPOL and around $0.858\pm0.15$ for HAWC+. This indicates that the outflows could significantly impact the alignment of dust grains and magnetic fields in the L1551 region.

The exact mechanism of protostellar outflows is still unclear; however, magnetic fields are thought to be a key component, as they can accelerate the transfer of angular momentum during protostar formation (Bally 2016).Some studies suggest that outflow feedback might modify the magnetic field structure on scales of around 1000 AU (Hull et al. 2014;Hull & Zhang 2019;Xu et al. 2022); however, ★ E-mail: yue.hu@wisc.edu† E-mail: alazarian@facstaff.wisc.eduothers have indicated that the outflow geometry and the plane-ofthe-sky (POS) magnetic field may not be interconnected, displaying random relative orientation on larger scales (Pattle et al. 2022).This discrepancy has yet to be resolved.The magnetic field is primarily inferred from the polarized thermal emission of interstellar dust grains, which tend to align their semi-major axes perpendicular to local magnetic fields due to radiative torques (RATs; Lazarian & Hoang 2007;Andersson et al. 2015).Dust grains may retain this alignment in both diffuse parts in the ISM and also the dense regions with a volume density of up to 10 4 cm −3 (Hough et al. 2008;Seifried et al. 2019).Meanwhile, however, the diffuse molecular outflows are typically traced by diffuse 12 CO and 13 CO corresponding to the critical density of 10 2 and 10 3 cm −3 , respectively, when only collisional (with H 2 ) de-excitation and spontaneous radiative de-excitation at kinetic temperature ∼ 10 K are considered (Goldsmith & Langer 1978;Nagahama et al. 1998;Draine 2011).Such differences imply that the dust polarization measurements may reflect the superposition of magnetic fields in both the diffuse and dense regions, not necessarily those explicitly associated with the molecular outflows.
To directly trace the magnetic field associated with outflows, the Velocity Gradient Technique (VGT; González-Casanova & Lazarian 2017; Yuen & Lazarian 2017;Hu et al. 2023) offers a promising way.VGT is founded on the understanding that in MHD turbulence, eddies elongate along the magnetic field, exhibiting anisotropy (Goldreich & Sridhar 1995;Lazarian & Vishniac 1999).The velocity gradient of these turbulent eddies acts as a detector for this anisotropy and can, therefore, be used to trace the magnetic fields.The anisotropy, in turn, is imprinted on spectroscopic observations (Lazarian & Pogosyan 2000;Kandel et al. 2017;Hu et al. 2021b;Hu & Lazarian 2023).Thus, VGT can adeptly utilize molecular emission lines to trace magnetic fields.However, the presence of outflow feedback may alter turbulence properties (Hu et al. 2022a) and impose additional velocity gradients that are not associated with MHD turbulence.This effect is expected to be scale-dependent as outflows are typically weak at large scales, but this relationship has not yet been explored.
To study the magnetic field and velocity gradient properties amidst outflow, it is vital to obtain multi-scale and multi-frequency magnetic field measurements.These are now readily available through the use of advanced polarimetry on dust emission, such as the Planck satellite 353 GHz/849 m survey (Planck Collaboration et al. 2020b), the JCMT/SCUPOL survey at 850 m (Matthews et al. 2009), and the SOFIA/HAWC+ survey at 214 m (Harper et al. 2018).These instruments have been utilized for observing the molecular cloud L1551 (Snell et al. 1980;Snell 1981;Yoshida et al. 2010;Lin et al. 2016), a low-mass star-forming cloud (M ≈ 50 M ⊙ ) exhibiting substantial bipolar molecular outflow from young stellar objects L1551 IRS 5 (hereafter IRS 5) and L1551 NE (hereafter NE).Given that L1551 is located approximately 160 pc away, SCUPOL and SOFIA/HAWC+ offer high spatial resolutions of 10 ′′ and 18.2 ′′ respectively, capably resolving scales of 0.008 pc and 0.014 pc associated with stronger outflows.Importantly, the measurements from these two distinct wavelengths are expected to trace the magnetic field associated with different gas/dust temperatures.Additionally, we apply VGT to 12 CO (J = 1-0) emission lines, while its application to 13 CO (J = 1-0) has been investigated (Hu et al. 2019).As 12 CO is a common tracer for protostellar outflows (Snell et al. 1980;Stojimirović et al. 2006;Bally et al. 2007;Bally 2016), the magnetic fields inferred by VGT from the molecular gas emission are primarily associated with outflows in the star-forming region.Comparison of results from VGT and polarimetry will yield unique insights into the properties of velocity gradients in the presence of outflow feedback.
The paper is structured as follows: § 2 provides an overview of the observational data used in our study.In § 3, we dive into the VGT, discussing its theoretical foundation in terms of MHD turbulence and the software framework.Our principal findings concerning the magnetic fields in L1551, as derived from polarization measurements and VGT, are revealed in § 4. Subsequently, § 5 offers a discussion of the various magnetic field measurements.Our conclusions are drawn in § 6.

12 CO emission
The 12 CO (J = 1-0, 115.27120GHz) emission line of L1551 was observed by the Nobeyama 45m radio telescope from December 2007 to May 2008 (Yoshida et al. 2010).At 115 GHz, the telescope has a half-power beam width (HPBW) of 15 ′′ .The region was mapped into a 45 ′ ×45 ′ sample with a 7. ′′ 5 pixel grid size.The velocity resolution of the data is 37.8 kHz in frequency or 0.098 km s −1 in velocity, along with an RMS noise level of 1.23 K in  mb (Lin et al. 2016).We applied the VGT to the cube with a grid resolution of 7. ′′ 5.All velocity channels within the velocity range from -9.9 km s −1 to 14.9 km s −1 , in which the cloud's emission is concentrated (see Appendix C for the gas emission spectrum), are considered in VGT calculation.

Polarized dust emission
We adopt the following polarized dust emission data: Planck 353 GHz data from the 3rd Public Data Release (Planck Collaboration et al. 2020b), SOFIA/HAWC+ observation, and JCMT/SCUPOL survey (Matthews et al. 2009).
Planck: the Planck polarization has beam resolution ≈ 5 ′ .The magnetic field orientation  B =  + /2 is derived from the polarization angle , which is calculated by the Stokes  and : where the negative sign in − converts the angle originally in the HEALPix system to the IAU system.Note the  periodicity in radian denoted by the arc-tangent function.The polarization fraction  is defined as: where I is the intensity of dust emission.However, the polarization needs to be further debaised to correct for the instrumental error.We employ the common debiasing recipe as in Wardle & Kronberg (1974) and Pattle et al. (2019): and  are the uncertainty in the Stokes  and .SOFIA/HAWC+: The observation was carried out in band E (214 m) with an HPBW of ∼ 18.
′′ 2 in November 2021.Similar to Planck, the magnetic field orientation is derived from Stokes  and .We filter the data to keep only pixels that satisfy: i) /  > 3; ii) / > 4, where   and  are the uncertainty of polarization fraction and intensity, respectively.In addition, it is observationally rare to obtain a polarization fraction higher than 30% in diffuse ISM (Fanciullo et al. 2022).For example, Planck Collaboration et al. (2020a) reported that the polarization fraction across the full sky has a maximum of approximately 22%.Therefore, pixels with a polarization fraction larger than 25% are blanked out.
JCMT/SCUPOL: The data used in this work was retrieved from the SCUBA Polarimeter Legacy Catalogue 1 with an angular resolution of 10 ′′ ×10 ′′ at 850 m (Matthews et al. 2009).The data was processed by the SCUPOL team to achieve the following criteria: i) /  > 2; ii)   < 4% ; iii) I > 0.

Essential of MHD turbulence and the Outflow Effect
Anisotropy of MHD turbulence: In decades, our understanding of MHD turbulence has rapidly changed.Many studies revealed that turbulence in the presence of magnetic fields is isotropic rather than anisotropic.Especially the cornerstone theoretical frame developed in Goldreich & Sridhar (1995) for trans-Alfvénic turbulence (i.e., the Alfvén Mach number M  ∼ 1) found the turbulent eddies are elongating along the magnetic fields, exhibiting the scaling relation: where  ∥ and  ⊥ represent the components of the wavevector parallel and perpendicular to the mean magnetic field, respectively.Also, the "critical balance" condition, i.e., the cascading time ( ⊥   ) −1 equals the wave periods ( ∥   ) −1 , and Kolmogorov-type turbulence, i.e.,   ∝  1/3 , were considered.Nonetheless, it should be noted that the "critical balance" condition considered in Goldreich & Sridhar (1995) was based on a global reference frame, where the direction of wavevectors is defined relative to the mean magnetic field.
For the purpose of tracing magnetic fields, it is important to find the anisotropy in the local reference frame defined relative to the magnetic field passing through an eddy at scale .This was later achieved by Lazarian & Vishniac (1999) in the study of turbulent reconnection.According to Lazarian & Vishniac (1999), turbulent reconnection of the magnetic field, takes place over just one eddy turnover time, naturally giving the "critical balance" in the local frame:  ,⊥  −1 ⊥ ≈    −1 ∥ , where  ⊥ and  ∥ represent the perpendicular and parallel scales of eddies with respect to the local magnetic field, respectively.Also, the motion of eddies perpendicular to the local magnetic field direction adheres to the Kolmogorov law (i.e.,  ,⊥ ∝  1/3 ⊥ ), since this is the direction in which the magnetic field offers minimal resistance.The scale-dependent anisotropy scaling is then given by: where  ⊥ and  ∥ represent the perpendicular and parallel scales of eddies with respect to the local magnetic field, respectively. inj denotes the turbulence injection scale and M  =  inj /  .On the other hand, the scaling of turbulent velocity and its corresponding gradient are: Here the velocity gradient is dominated by the perpendicular component, due to the anisotropy (i.e.,  ∥ ≫  ⊥ ).One can therefore derive the local magnetic field orientation from the velocity gradient.
Outflow Effects: Our earlier discussion presupposes the dominance of MHD turbulence.Nonetheless, within star-forming regions, outflow feedback might alter the properties of velocity gradients.Outflow feedback has been shown to substantially change fluid velocity statistics (Hu et al. 2022a) and is anticipated to mirror the effects of inflows on velocity gradients, i.e., causing the orientation of the velocity gradient to shift from being perpendicular to the magnetic fields to being aligned with them (Hu et al. 2020).Since 12 CO only resolves outflows in the L1551 region (Stojimirović et al. 2006), we focus our analysis on the shift of velocity gradients by outflows rather than inflows.It is well-documented that outflows are expelled from star-forming sites at high velocities (Snell et al. 1980).At a protostar's core, gas velocity peaks and diminishes towards its periphery.As a result, additional velocity gradients not attributed to turbulence might emerge, pointing from the center to the outskirts.These outflow-induced gradients are highly scale-dependent because of their correlation with the significance of outflows.On smaller scales proximate to protostars, outflows significantly overshadow turbulence.Magnetic fields in this case are expected to follow the outflows.Thus, the velocity gradient is primarily dominated by the potent outflows, being parallel to the magnetic fields.Conversely, on broader scales distant from protostars, weakened outflows allow turbulence-induced velocity gradients to take precedence, rendering the gradients perpendicular to the magnetic fields.A transition may exist wherein velocity gradients gradually shift from parallel to perpendicular alignment with the magnetic fields as one moves further from the protostar's center, contingent upon the comparative prominence of outflows versus turbulence.

VGT pipeline
Velocity caustics effect: Employing VGT requires velocity information to be extracted at each pixel of observed images.For such purpose, we treat the Doppler-shifted emission line data into maps with thin velocity channels.This is inspired by the velocity caustics effect (Lazarian & Pogosyan 2000;Kandel et al. 2017), which indicates that in observations, physical (density or intensity) structures are sampled into different channels with respect to their LOS velocities.In turbulent ISM, the intensity structures in each channel are distorted.Such distortion becomes more significant as the channel width gets narrower.Velocity fluctuations will eventually dominate over intensity fluctuations on Δ < √︁ ( 2 ): the channel width Δ is smaller than the velocity dispersion √︁ ( 2 ).In such thin channels, calculating the gradient of intensity amplitudes would in fact generate gradient information of the velocity fluctuations.More details about the thin channel as well as the Position-Position-Velocity (PPV) theory can be found in Lazarian & Pogosyan (2000); Kandel et al. (2017); Hu et al. (2023).
VGT: The software workflow for VGT used in this work follows the procedure introduced in Hu et al. (2022a) which is summarized below: each individual thin channel map Ch i (x, y)  is convolved with the 3×3 Sobel convolutional kernels 2 : the asterisks denote convolutions.∇  Ch i (, ) and ∇  Ch i (, ) are the gradient components in each thin channel maps along the x and y axis.They are used to calculate the overall pixelized map  i g (, ) of the gradients: From here, only pixels with an intensity less than three times the RMS noise are kept towards the following steps.The resulting  i g is then processed with the sub-block averaging method (Yuen & Lazarian 2017): i) the pixelized map is divided into rectangular sub-blocks with an optimal size of 20×20 pixels (verified both empirically and 2 numerically: Lazarian & Yuen 2018;Hu et al. 2021a); ii) a histogram is produced for the mean gradient angle within each sub-block; iii) the histogram is fitted with a Gaussian distribution, the peak value of which is then defined as the most probable orientation of the gradient in that sub-block.Now the averaged gradient angle map  i gs (, ) constructs the pseudo-Stokes parameters   and   : thus the POS magnetic field angle probed by the VGT  B can be inferred from the pseudo-polarization angle  g as  B =  g + /2.

Molecular outflows
Fig. 1 reproduced the bipolar molecular outflows observed in the L1551 region.As detailed by Yoshida et al. (2010), the LSR velocities of the red-and blue-shifted components span from 7.9 km s −1 to 17.4 km s −1 and -3.7 km s −1 to 5.5 km s −1 , respectively.These outflows are almost symmetric around L1551 IRS 5, the binary protostar system located at the intersection of the two components.The blue-shifted flow stretches southwest, while its red-shifted counterpart spans northeastward, with an outlying structure that extends westward for approximately 20 ′ .These redshifted outflows might originate from a newly formed star near HH 102, a Herbig-Haro diffuse reflection nebula situated 5 ′ west of L1551 IRS 5 (Strom et al. 1976;Stojimirović et al. 2006).Nevertheless, other studies, like Moriarty-Schieven et al. (2006), suggest that the east-west outflow's driving source is one of the stars in the multi-star system L1551 NE (Reipurth et al. 2002;Yoshida et al. 2010).
On the contrary, in the southwest area, the magnetic fields are seen to cross the 12 CO intensity structure, which is associated with the positions of the blue-shifted outflows (see Fig. 1).In particular, at the center of L1551, i.e. around the IRS 5 protostar , the magnetic fields are significantly twisted.We observed that Planck's polarization fraction around the cloud is negative after debiasing, as seen in Fig. 3, indicating a considerable depolarization effect.This could be due to (1) turbulent magnetic fields on the POS and (2) significant variation of magnetic fields along the LOS.The low polarization fraction implies a very low signal-to-noise ratio.However, this could be reduced by increasing the resolution of the observation, i.e., using a smaller beam.Therefore, we further investigated the SCUPOL polarization, which was measured at a similar wavelength as Planck but has a much smaller beam width of ∼ 10 ′′ .As seen in Fig. 2, the magnetic fields (∼ 10 ′′ ) inferred from SCUPOL polarization are less regular.The polarization fraction is distributed in the range of 0 -5%, approximately, as seen in Figs. 3 and 4. To compare it with Planck polarization, we further smoothed it to ∼ 5 ′ , and the results can be found in the Appendix.C. Similar to Planck, the SCUPOL measurements reveal a twisted morphology of magnetic fields near IRS 5.If the Planck data was significantly contaminated with noises, the two independent observations are not expected to yield consistent results of magnetic fields in the same area.Therefore, we expect the Planck data to be still reliable in reflecting the magnetic field structure in the L1551 region and thus have kept the Planck data in our analysis.
We further synergized the HAWC+ data obtained most recently with Planck and SCUPOL to investigate the magnetic fields in L1551.Planck and SCUPOL polarization are measured at 850 m, which is usually associated with cold, dense gas.HAWC+ provides measurements at a shorter wavelength of 214 m, which is linked to gas with a higher temperature (Huang 1987).As seen in Fig. 2, the magnetic fields inferred from HAWC+ are still turbulent, but do not appear to be similar to the other two polarization measurements.While the polarization fraction values at most pixels are no larger than 5%, the maximum value could be up to more than 20%, as shown in Figs. 3  and 4. The uncertainties of the HAWC+ polarization are included in the Appendix A. This suggests that the magnetic field morphology may vary in different gas phases (or temperatures) in L1551.While polarization fraction is irrelevant to the following analysis with VGT, it is a significant intrinsic property of polarization measurements and can also shed light on the dust grain alignment within this cloud, which will be discussed in Sec. 5.

Magnetic fields inferred from VGT
Fig. 5 shows the magnetic fields inferred from VGT using the 12 CO emission line.The VGT-12 CO approach uniquely permits the separation of magnetic fields correlated with distinct cloud structures: (1) total: the full 12 CO emission, integrated from -15 km s −1 to 25 km s −1 ; (2) redshifted: exclusively for the redshifted component in the gas emission, spanning 7.9 km s −1 to 17.4 km s −1 ; (3) blueshifted: exclusively for the blueshifted gas, ranging from -3.5 km s −1 to 5.5 km s −1 .The effective resolution of these VGT-12 CO measurements is 20×20 pixels.
In the fully integrated VGT-12 CO measurements, the magnetic fields broadly align with the 12 CO intensity configurations in both the northeastern and southwestern parts.Remarkably, the magnetic fields in VGT-12 CO are not twisted near the IRS 5 protostar.This stands as different from the Planck polarization.As for the red-shifted and blue-shifted ones, the magnetic fields are less regular but still generally follow the gas structures, except for the red-shifted west tail.The magnetic fields are also twisted towards the center around IRS 5 in both red-shifted and blue-shifted cases.
To quantify the agreement between VGT-12 CO measurements and polarization, we utilize the Alignment Measure (AM; González-Casanova & Lazarian 2017), expressed as: Here,  r = | B − B |.An AM value of 1 implies parallel alignment of  B and  B , while -1 indicates perpendicularity.
Fig. 6 shows the spatial distribution of the AM between the fully integrated VGT-12 CO measurement and the magnetic field orienta-  To resolve the inconsistency in the spatial resolutions of the polarimetry surveys, a Gaussian filter was applied to each of the polarimetry data so that every resulting map has the same 20×20 pixel resolution.The graphs are different in shape with the vectors in Fig. 3, which are averaged over several pixels for visualization purposes, while each map here was produced by calculating AM at every pixel.tion inferred from Planck polarization.In the vicinity of the IRS 5 protostar, anti-alignments (AM = -1) are observed.This is also seen in SCUPOL and HAWC+.We attribute this to the strong molecular outflows characteristic of this area.As delineated in Sec. 3, velocity gradients stemming from turbulent gas motion are typically perpen-dicular to magnetic fields and are used in VGT to trace the magnetic fields.This type of gradient is denoted as turbulence-asscociated velocity gradients.However, vigorous outflow dynamics can create pronounced velocity differences, or non-turbulence-associated velocity gradients, which we anticipate to align parallel to the magnetic fields in the presence of substantial outflows.
When extracting gradients from observational data, we inherently deal with the combined velocity field (turbulence plus nonturbulence).Consequently, the VGT-12 CO -which rotates the velocity gradients by 90 degrees to estimate magnetic field orientations -may result in a perpendicular anti-alignment with polarization when outflows predominate.Moving outward from the center, where outflow influence wanes and turbulence takes precedence, the VGT-12 CO measurements begin to align with magnetic field orientations inferred from polarization, indicating the relative increase in the significance of turbulent motions.
Furthermore, Fig. 7 shows the histogram of the VGT-Planck AM for the red-shifted component, the blue-shifted component and the total gas emission.An intriguing feature is that in both redshifted gas and total emission, the VGT-Planck alignments (AM = 1) are statistically more significant, while in the blue-shifted component, AM concentrates more on ∼ 0.5.We expect this is due to the Planckinferred magnetic fields do not follow the blueshifted outflows, as seen in Fig. 2.
The AM map between the fully integrated VGT-12 CO measurement and both SCUPOL and HAWC + results has been studied, as seen in Fig. 8. Since the polarimetry surveys observed the protostar systems in L1551 at close-in scales, the magnetic fields they trace are likely related to the protostellar outflows.However, because of the limited sampling ranges of SCUPOL and HAWC+, their measurements are not compared with the magnetic fields associated with individual outflow components.Similarly to the VGT-Planck AM map, anti-alignments of magnetic field measurements inferred by the VGT and the two polarizations are mostly detected near the IRS 5 protostar, suggesting that strong outflows exist in the area, resulting in a difference between VGT and the polarization measurements.Moreover, the anti-alignment is observed in the three polarization measurements with different beam resolutions of ∼ 5 ′ , 10 ′′ , and 18.2 ′ .Since the trend of negative AM was observed in both comparisons with the low-resolution Planck data and the high-resolution SCUPOL and HAWC+, we infer that this may imply the outflow motions dominate the gas dynamics of the L1551 star-forming region starting from ∼ 5 ′ , corresponding to a physical scale of ∼ 0.2 pc.Nonetheless, other observations are needed for cross-checking to confirm this speculation.

Dust grain alignment in the presence of outflows
The polarization fraction is usually inversely related to the Stokes  in a power-law fashion:  ∝  −  .Index  is usually used to assess the effectiveness of dust grain alignment in star-forming regions.It is believed that a value of  = 0 implies perfect alignment between the grains and the magnetic field, while an index of  = 1 suggests a total loss of alignment (Whittet et al. 2008;Alves et al. 2014;Jones et al. 2015).In Fig. 9, we present the correlation between  and  for SCUPOL and HAWC+ data.We find the index is ∼ 1.044 ± 0.06 and ∼ 0.858 ± 0.15 for SCUPOL and HAWC+, respectively, indicating a lack of grain alignment in the L1551 region.This, along with the negative polarization fractions in Planck results, induced us to provide an additional analysis, in which we observed similarities between the SCUPOL and Planck measurements in vital statistics.As shown in Appendix C, Planck, and SCUPOL, which are nearly at identical wavelengths (∼ 850 m) with different spatial resolutions indicate similar magnetic field orientations.While the magnetic field structure inferred from HAWC+ at 214 m shows differences with the former two, all three polarization measurements have negative AM with VGT in the central star-forming region.Should dust grain alignment fail to give reliable indications of magnetic fields, the polarization measurements are not expected to show similar magnetic field orientation as well as AM values.Since three independent polarization measurements give similar results, we expect the uncertainties in polarization to be insignificant and that the comparison with VGT stands valid.

VGT
The magnetic field mapped by VGT is subject to two sources of uncertainty: (1) systematic errors in the observational cube and (2) the uncertainty inherent in the VGT algorithm itself.The latter uncertainty arises from the subblock-averaging approach utilized in the algorithm.Specifically, the VGT fits a Gaussian histogram to the orientation of the gradient within a subregion and outputs the angle corresponding to the peak value of the histogram.The associated uncertainty can be quantified as the error   s (, ) from the Gaussian fitting algorithm within a 95% confidence level.The uncertainty maps of the VGT-12 CO measurement are given in the Appendix B, indicating a globally low level of error.

Velocity gradients in the presence of outflows
When assessing to what extent VGT-12 CO aligns with the polarization measurements, a common trend is revealed that in the outskirts of the L1551 region, these measurements tend to align with each other reasonably well (AM → 1); yet approaching the center of the star-forming region, especially near the protostellar system IRS 5 which is believed to be the major powering source of outflows in L1551 (Fridlund & Liseau 1998;Itoh et al. 2000;Rodríguez et al. 2003), these measurements seem to produce mostly perpendicular orientations at the same position (AM → -1).As detailed in Sec. 4, this could be attributed to the presence of strong protostellar outflows in such regions.Energetic outflows ejected from the protostar cores generate an outbound acceleration that induces significant velocity gradients.At positions close to the protostars, such non-turbulenceassociated gradients generated by the outflow gas motions are the greatest and thus dominate over the gradients induced by turbulent motions in the gas, causing the VGT-inferred magnetic fields to appear perpendicular to the actual magnetic field orientations, which eventually are reflected by the anti-alignments with polarization measurements.Such dominance of non-turbulence-associated gradients stays profound starting from ∼ 0.20 pc.At the outer edge away from the protostellar systems, the influence of outflows decreases and turbulence-induced gradients become dominant again.Thus, a high degree of alignment in the VGT and polarization is observed in the outskirts of L1551.This emphasizes the fact that VGT and polarimetry measurements respond differently to the presence of outflows.While polarization could be altered by rotational disruptions of dust grains by both RATs and mechanical torques (METs; Hoang et al. 2022), velocity gradients are closely associated with the kinematic properties of the cloud and thus are mainly affected by the dynamical process.More efforts should be dedicated to investigating in detail what effects those factors, as well as other vital physical properties such as radiations and shocks, have on polarization and gradients, as that can not only better determine the scenario where each technique becomes the most effective, but also promote our insights on the fundamentals of star formation.

Prospects of VGT
In this work, magnetic fields directly associated with the protostellar outflows in the L1551 star-forming region have been investigated with VGT.This novel technique utilizes velocity information from spectroscopic observations to access magnetic fields, hence is capable of separating magnetic fields associated with red-shifted and blue-shifted gas with specific velocities.This way of applying VGT was inspired by Liu et al. (2022) and Hu et al. (2022b), where the technique was used to decompose the magnetic field associated with different velocities components in the supernova remnant W44 and the Central Molecular Zone, respectively.This capacity gives VGT unprecedented advantages in studying astrophysical systems that have complicated velocity components.For instance, as proposed in Liu et al. (2023), different parts of a galaxy, such as spiral arms, typically rotate with different velocities.With VGT, it is possible to isolate them in a velocity space and investigate their magnetic properties individually.

CONCLUSIONS
Magnetic fields are integral to the process of star formation, though tracing them in regions with outflow feedback presents complexities.In this work, we use polarization observations from the Planck at 849 m, the SOFIA/HAWC+ measurement at 214 m, and the JCMT/SCUPOL at 850 m, as well as VGT with the 12 CO (J = 1-0) emission data, to trace the magnetic fields in the star-forming region L1551.Our major discoveries are: (i) Upon targeting the star-forming region L1551, known for its dis-cernible outflows, polarization consistently revealed the magnetic fields to be twisted in the direction of the protostar IRS 5. (ii) By applying VGT to the 12 CO (J = 1-0) emission in specific velocity ranges, we were able to separate magnetic fields that are directly tied to the protostellar outflows (iii) While VGT and polarization showed a general alignment on the outskirts of the L1551 region, a trend of mostly perpendicular relative orientations was observed closer to the center of IRS 5. (iv) The observed perpendicular orientations imply that outflows could be dynamically significant from a scale of around ∼ 0.2 pc.This dynamic significance may be causing a change in the direction of velocity gradients by 90 degrees.(v) A distinct correlation  ∝  −  was observed between the polarization fraction  and the total intensity .Values for  were found to be 1.044 ± 0.06 for SCUPOL and 0.858 ± 0.15 for HAWC+, suggesting that outflows may play a crucial role in determining the alignment between dust grains and magnetic fields within the L1551 region.revealing that in a larger scale, the magnetic fields present a similar trend as in the Planck map that the orientations bend dramatically near IRS 5.

Figure 2 .
Figure 2. Magnetic field maps inferred from Planck (left), SCUPOL (middle), and HAWC+ (right).The magnetic field orientations are represented by colored line segments, overlaid upon the integrated intensity map of the 12 CO emission.SCUPOL and HAWC+ only surveyed a small portion of the L1551 region, thus the maps are zoomed in to the main intensity structure for better illustration, and so does every figure that contains results from the two surveys in the following sections.

Figure 6 .
Figure 6.Spatial distributions of AM between VGT measurements on total 12 CO emission and Planck (left), SCUPOL (middle), and HAWC+ (right).The positions of L1551 NE (white circle) and L1551 IRS 5 (white square) have been highlighted on each map.To resolve the inconsistency in the spatial resolutions of the polarimetry surveys, a Gaussian filter was applied to each of the polarimetry data so that every resulting map has the same 20×20 pixel resolution.The graphs are different in shape with the vectors in Fig.3, which are averaged over several pixels for visualization purposes, while each map here was produced by calculating AM at every pixel.

Figure 7 .
Figure 7. Histogram of AM between Planck and VGT measurements on total 12 CO emission (black line), redshifted emission (red line), and blueshifted emission (blue line).

Figure 8 .
Figure 8. Histogram of AM between VGT measurements on total 12 CO emission and SCUPOL (orange line) and HAWC+ (purple line).

Figure 9 .
Figure 9. Polarization fraction  as a function of intensity (Stokes ) of the SCUPOL (left) and HAWC+ (right) data in the log-log scale.Gray dots represent original data points, and colored dotted lines mark the best-fit trend lines with slopes indicated, which correspond to the  index that reflects the efficiency of dust grain alignment with magnetic fields.

Figure A1 .
Figure A1.Uncertainty in polarization fraction of the Planck data.

Figure A2 .
Figure A2.Uncertainty in polarization fraction of the HAWC+ data.

Figure A3 .
Figure A3.Uncertainty in polarization angle of the Planck data.

Figure A4 .
Figure A4.Uncertainty in polarization angle of the HAWC+ data.

Figure B2 .
Figure B2.Uncertainty in the VGT measurements on the redshifted 12 CO emission.

Figure B3 .
Figure B3.Uncertainty in the VGT measurements on the blueshifted 12 CO emission.

Figure C1 .
Figure C1.Spectral line of 12 CO emission.Red and blue dotted lines mark the redshifted and blueshifted regions, respectively.

Figure C2 .
Figure C2.Smoothed magnetic field orientation map inferred from the SCUPOL data.
Fig. C4 is the histogram of polarization angle (east from the north) in the SCUPOL and HAWC+ data.Fig. C5 shows the I-p correlation of the debiased Planck data of the L1551 region.Its low polarization fraction along with the  index close to 0 imply that it has a questionable accuracy for tracing the dense part of the cloud.

Figure C3 .
Figure C3.Smoothed magnetic field orientation map inferred from the HAWC+ data.

Figure C4 .
Figure C4.Histogram of polarization angle in the SCUPOL (orange line) and HAWC+ (purple line) data.Note that the angle is defined in the IAU convention.

Figure C5 .
Figure C5.Polarization fraction  as a function of intensity (Stokes ) of the Planck data in the log-log scale.