Tracing Multi-Scale Magnetic Field Structure Using Multiple Chemical Tracers in Giant Molecular Clouds

TL;DR

This study demonstrates that the Velocity Gradient Technique, applied to multiple chemical tracers with different optical depths, effectively maps magnetic field structures in molecular clouds, aligning well with polarization observations.

Contribution

It validates the use of multiple chemical tracers in the Velocity Gradient Technique for detailed magnetic field mapping in giant molecular clouds.

Findings

Good correspondence between velocity gradient and polarization-based magnetic field maps.

Combining tracers improves magnetic field structure reconstruction.

Limited signatures of gravitational collapse detected in the cloud.

Abstract

Probing magnetic fields in Giant Molecular Clouds is often challenging. Fortunately, recently simulations show that analysis of velocity gradients (the Velocity Gradient Technique) can be used to map out the magnetic field morphology of different physical layers within molecular clouds when applied CO isotopologues with different optical depths. Here, we test the effectiveness of the Velocity Gradient Technique in reconstructing the magnetic field structure of the molecular cloud Vela C, employing seven chemical tracers that have different optical depths, i.e. 12CO, 13CO, C18O, CS, HNC, HCO+, and HCN. Our results show good correspondence between the magnetic field morphology inferred from velocity gradients using these different molecular tracers and the magnetic field morphology inferred from BLASTPol polarization observations. We also explore the possibility of using a combination of velocity gradients for multiple chemical tracers to explain the structure of the magnetic field in molecular clouds. We search for signatures of gravitational collapse in the alignment of the velocity gradients and magnetic field and conclude that collapsing regions constitute a small fraction of the cloud.

Tables

Table 1. Table 1: Observation parameters for seven chemical tracers used in our analysis, including the type of molecules, their respective transition levels, observing frequency in GHz, original telescope beam full width half maximum (FWHM) without any additional smoothing in arcsec, the smallest thickness of the velocity channel slice, and the approximate molecular hydrogen density (N H 2 subscript 𝐻 2 {}_{H_{2}} ) we expect each molecule to trace. Note that for all of molecules we have observed the ground state (J = 1 – 0) transitions. ξ 𝜉 \xi is the velocity channel width of each molecular line cube. Note that the cubes are Nyquist sampled so the true velocity resolution is 2 ξ 𝜉 \xi . δ ​ v 𝛿 𝑣 \delta v is the velocity dispersion and σ T R subscript 𝜎 subscript 𝑇 𝑅 \sigma_{T_{R}} is the per channel noise level of T R subscript 𝑇 𝑅 T_{R} for each data cube. Fissel et al. ( 2019 ) provides more details about the data reduction.

Molecular Line	Line of Transitions	Frequency [GHz]	FWHM [arcsec]	$\xi \left[km{s}^{-1}\right]$	$\delta v\left[km{s}^{-1}\right]$	Tracing Density	${\sigma }_{T_R}$[K]
12CO	J = 1 - 0	115.27	27.12	0.18	3.48	$\sim {10}^{2}c{m}^{-3}$	0.113
13CO	J = 1 - 0	110.20	28.37	0.18	4.75	$\sim {10}^{3}c{m}^{-3}$	0.053
C18O	J = 1 - 0	109.78	36.01	0.18	8.32	$\sim {10}^{4}c{m}^{-3}$	0.053
CS	J = 1 - 0	48.99	63.82	0.21	9.61	$\sim {10}^{4}c{m}^{-3}$	0.095
HNC	J = 1 - 0	90.66	34.48	0.22	4.97	$\sim {10}^5-{10}^{6}c{m}^{-3}$	0.039
HCO+	J = 1 - 0	89.19	35.05	0.23	4.79	$\sim {10}^5-{10}^{6}c{m}^{-3}$	0.018
HCN	J = 1 - 0	88.63	35.28	0.23	6.41	$\sim {10}^5-{10}^{6}c{m}^{-3}$	0.019

Table 2. Table 2: The expectation value μ 𝜇 \mu of the relative angle between rotated VCGs/VChGs and the magnetic field inferred from polarization, without sub-block averaging method. The uncertainty is given by the standard error of the mean, i.e. the standard deviation divided by the square root of the sample size, while the systematic uncertainty in polarization data is given in Appendix A .

Molecule Line	$\mu$ (VCGs)	$\mu$ (VChGs)	AM (VCGs)	AM (VChGs)
12CO	13.47${}_{}{}^{\circ }\pm {0.29}^{\circ }$	10.95${}_{}{}^{\circ }\pm {0.51}^{\circ }$	0.65$\pm 0.03$	0.58$\pm 0.02$
13CO	13.32${}_{}{}^{\circ }\pm {0.37}^{\circ }$	7.07${}_{}{}^{\circ }\pm {0.56}^{\circ }$	0.61$\pm 0.03$	0.66$\pm 0.03$
C18O	13.64${}_{}{}^{\circ }\pm {0.37}^{\circ }$	5.52${}_{}{}^{\circ }\pm {0.53}^{\circ }$	0.68$\pm 0.04$	0.70$\pm 0.03$
CS	7.24${}_{}{}^{\circ }\pm {0.22}^{\circ }$	5.68${}_{}{}^{\circ }\pm {0.27}^{\circ }$	0.63$\pm 0.02$	0.70$\pm 0.02$
HCN	13.84${}_{}{}^{\circ }\pm {0.44}^{\circ }$	5.35${}_{}{}^{\circ }\pm {0.40}^{\circ }$	0.66$\pm 0.04$	0.60$\pm 0.05$
HNC	7.42${}_{}{}^{\circ }\pm {0.36}^{\circ }$	5.07${}_{}{}^{\circ }\pm {0.38}^{\circ }$	0.63$\pm 0.04$	0.68$\pm 0.04$
HCO+	12.72${}_{}{}^{\circ }\pm {0.41}^{\circ }$	6.60${}_{}{}^{\circ }\pm {0.29}^{\circ }$	0.48$\pm 0.06$	0.62$\pm 0.04$

Equations

$$\frac{l}{v_l}\sim \frac{l_{||}}{v_A}$$

$$l_{||}\sim l^{\frac{2}{3}}$$

$$C(x,y)=\frac{\int dv{T}_R(x,y,v)\cdot v}{\int dv{T}_R(x,y,v)}$$

$$\Delta v\le \sqrt{⟨\delta v^2⟩}$$

$$Ch(x,y)={\int }_{v_0-\Delta v/2}^{v_0+\Delta v/2}dv{T}_R(x,y,v)$$

$$AM=2(⟨co{s}^2{\theta }_r⟩-\frac{1}{2})$$

Full text

ISM:general—ISM:structure—magnetohydrodynamics(MHD)—radio lines:ISM—turbulence

1 Introduction

The magnetic field in the universe plays an essential role in multiple astrophysical processes, e.g., regulation of star formation (Burkhart et al., 2015; Chapman et al., 2013), guiding material and thermodynamic transfer between different media, and propagation and acceleration of cosmic rays (Fermi, 1949; Caprioli & Spitkovsky, 2014). One of the most critical roles of the magnetic field plays is to modify ubiquitous interstellar turbulence (Mac Low & Klessen, 2004; Ballesteros-Paredes et al., 2007) both in the diffuse interstellar medium (ISM), e.g. neutral hydrogen (H I), and in the molecular gas (Bell & Lin, 1994; Heitsch et al., 2007; Ostriker et al., 2001) which is distributed over an extensive range of density regimes in galaxies. While most of the sightlines towards Giant Molecular Clouds (GMCs) have low column density ($\leq 10^{22}\,cm^{-3}$), it is the properties of the dense gas components, such as magnetic field strength, turbulent energy, and angular momentum, that regulate the region where new stars are being born (Murray et al., 2010; Leroy et al., 2008; Chen et al., 2015; Li et al., 2005; Soler et al., 2017; Soler & Hennebelle, 2017). Thus, an exploration of their contribution on the dynamical evolution of molecular clouds is crucial to fully understand the process of star formation (Hopkins et al., 2012; Gao & Solomon, 2004; Solomon & Sage, 1988).

Many methods of tracing the magnetic fields have been proposed, although each method has limitations and biases. For instance, measuring polarization from dust grain alignment (Andersson et al., 2015; Lazarian, 2007; Planck Collaboration et al., 2015) requires high sensitivity either large scale far infrared or sub-mm polarization maps, which are extremely time-consuming from ground-based telescopes and often have low dust grain alignment efficiency (Lazarian, 2007). Surveys of magnetic fields from polarization due to the selective extinction of starlight passing through dust clouds are not highly effective for clouds with large dust columns. The synchrotron polarization method (Draine, 2011; Carilli & Taylor, 2002; Jansson & Farrar, 2012) is mostly used to trace magnetic field in warm and hot phases of ISM, while though Faraday rotation one can only measure the magnetic field along the line of sight towards regions where the ionization is significant (Hill et al., 2013).

Both numerical simulations and observations have shown that the Velocity Gradient Technique (VGT) is a promising method in studying magnetic fields (González-Casanova & Lazarian, 2017; Yuen & Lazarian, 2017a; Lazarian & Yuen, 2018a; Lazarian et al., 2018; Hsieh et al., 2019; Hu et al., 2019a, c, b; Zhang et al., 2019). The first suggestion to use velocity gradients to trace magnetic fields was made in González-Casanova & Lazarian (2017) (henceforth GL17) where velocity centroid gradients (VCGs) were used as proxies of velocity gradients. This technique was applied to a wide range of column densities from diffuse neutral hydrogen (H I) gas in Yuen & Lazarian (2017a) (henceforth YL17) and Hu et al. (2018), and later Lazarian & Yuen (2018a) proposed to use velocity channel gradients (VChGs) to trace magnetic fields. The application of the techniques to the diffuse H I data has proven the promise of this radically new way of magnetic field study (Lazarian & Yuen, 2018a). A recent exploration of the gradients technique in self-absorbing molecular gas data and self-gravitating media (Yuen & Lazarian, 2017b; González-Casanova et al., 2019; Hsieh et al., 2019; Hu et al., 2019a, c, b) have also demonstrated that VGT can be used as tracers of the magnetic field in regions with different physical conditions using 13CO as molecular tracer.

In this paper, we present an observational example for a low galactic latitude Giant Molecular Cloud Vela C using the recent advancement of the VGT made by Hsieh et al. (2019) in numerical simulations, which demonstrated that VGT can also be applied to CO isotopologues including 12CO and C18O. We use 7 molecular line maps for the analysis in Vela C, including the 3 isotopes of CO, CS, HNC, HCO+, and HCN. The latter four molecular tracers are high density tracers (number density of H2 around $10^{4}$ cm${}^{-3}-10^{6}$ cm*-3*) compared to the CO tracers 12CO, 13CO, and C18O, which typically trace number densities of H2 between $10^{2}$ cm${}^{-3}-10^{4}$ cm*-3* (Shirley, 2015; Fissel et al., 2019). To evaluate the success of the VGT we compare our inferred magnetic field orientation from the VGT to the magnetic field orientation inferred from a large scale 500 $\mu$m BLASTPol polarization map first presented in Fissel et al. (2016).

In what follows, in Section 2 and Section 3, we describe the theoretical foundation and the Velocity Gradients Technique used in our work. In Section 4, we discuss the ability of VGT to trace magnetic fields over a large range of densities in GMCs and the contribution from the foreground and background. In Section 5, we show how to trace the magnetic field through the combination of different molecular tracers and estimate the fraction of the collapsing regions in Vela C. In Section 6 we discuss the possible application of gradient technique to other molecular clouds. In Section 7, we give our conclusions.

2 Theoretical Perspective on Self-Absorbing Gradient Technique

2.1 Basic MHD Turbulence Theory

The MHD turbulence has been explored both theoretically and numerically in decades (Shebalin et al., 1983; Higdon, 1984; Montgomery & Matthaeus, 1995). In Goldreich & Sridhar (1995) (henceforth GS95) they formulated the theory of incompressible MHD turbulence and also predicted predicted the turbulent anisotropy which later becomes the foundation of several magnetic fields tracing techniques. GS95 described that the scaling of turbulent eddies is approximate $v_{l}\sim l^{\frac{1}{3}}$, where $v_{l}$ is the turbulence velocity at scale $l$ and $l$ is the size of eddies perpendicular to the magnetic field. However, GS95 is done with respect to the mean magnetic field where the anisotropic relation is not expected to be observed.

In Lazarian & Vishniac (1999) they illustrate that the motion of turbulent eddies enables the magnetic fields to mix with minimal resistance from magnetic tension for eddies at all scales and thus a much faster rate of magnetic reconnection is allowed compared to the traditional Sweet-Parker model. As a result of the mixing motions, the Alfvénic turbulence tends to move along the magnetic field directions. To find how turbulent eddies evolve in the direction parallel to the magnetic field, it is also necessary to consider the mixing motions associated with magnetic eddies and Alfvén waves with the period equal to the period of an eddy:

[TABLE]

where $l_{||}$ is the parallel scale of the eddy and $v_{A}$ is Alfvén velocity. The correlation between the parallel and perpendicular scales of sub-Alfvénic turbulence, i.e. $v_{l}\leq v_{A}$, eddies can be obtained as:

[TABLE]

This correlation shows that turbulent eddies are elongated along the direction of the magnetic field and it holds for the eddies that are aligned with the local direction of the magnetic field that surrounds them. Incidentally, the concept of local magnetic field frame explains why, unlike the original GS95 treatment, the anisotropy of turbulence actually reflects the direction of magnetic fields that percolates turbulent eddies. This result was confirmed by numerical simulations (Cho & Vishniac, 2000; Cho & Lazarian, 2003; Maron & Goldreich, 2001). Due to this particular property of MHD turbulence, the velocities associated with a turbulent eddy are anisotropic so that the largest change of the velocity is in the direction perpendicular to the local direction of the magnetic field. Thus, it is essential that VGT is tracing the local magnetic field around eddies rather than the mean magnetic field.

It is worth to mentioning that several approaches have been proposed to trace the magnetic field based on the MHD anisotropy relation Eq. (2). The correlation function analysis (CFA) of the velocities and the Principal Component Analysis for Anisotropy (PCAA) were first proposed to study magnetic field morphology (Heyer et al., 2008), to estimate magnetization (Esquivel et al., 2015; Esquivel & Lazarian, 2011), and to determine the contribution of the fast, slow and Alfven modes in observed turbulence (Kandel et al., 2016, 2017a, 2017b) using the theoretical understanding of MHD turbulence discussed in Section 2.1. However, numerical study in Yuen et al. (2018) showed that comparing with VGT, CFA and PCAA face several issues in particular, the anisotropy may be distorted, multi-centered eddies or the contours are not closed. These significantly degrade the determination of the direction of anisotropy, thus the inferred magnetic field orientation through CFA and PCAA. Later Clark et al. (2014) proposed the Rolling Hough Transform (RHT) to study the magnetic field in the diffuse region based on the fact that so-called narrow ”H I fibers” are align parallel to the magnetic field orientation. However, the RHT requires linear structures in the ISM (Clark et al., 2014).

Moreover, VGT as a superior technique has been successfully tested to trace the local magnetic field from both diffuse region and absorbing media for the case of 13CO emission with different abundances and densities (González-Casanova et al., 2019), while Hsieh et al. (2019) showed in numerical simulations that the VGT can map out the magnetic field structures of different physical depths in a molecular cloud using CO isotopologues with different optical depths and in the presence of weak self-gravity.

2.2 Molecules as Probes of Gas in Different Density Regimes

Giant Molecular Clouds (GMCs) are the sites for most of the star formation to take place across the Milky Way and also in other galaxies. The star formation regions are often observationally unresolved in GMCs. The gas in GMCs have an average number density of the order of $10^{2}cm^{-3}$ and temperature of 10 - 30K (Heyer et al., 2009), while the star formation clumps may have number densities around to $10^{7}cm^{-3}$ and temperature as low as 10K or even lower (Williams et al., 1994). Molecular tracers like 12CO typically become optically thick, and therefore mostly trace the outer (low density) regions of molecular clouds. 13CO, which is typically a factor of $\sim$100 less abundant than 12CO and therefore usually optically thin, can be used to trace intermediate density regions, as can the even less abundant isotopologue C18O. Transitions of the molecules CS, HNC, HCO*+* and HCN have higher critical densities and are therefore used to trace intermediate or high-density molecular gas. These molecules tend to either trace regions where the effect of self-gravitation is very strong or gas that are just outside the self-gravitating regions. However, regions with high-density gas do not necessarily host star formation, that is, high density does not necessarily indicate that the self-gravitation is strong enough to trigger the gravitational collapse.

In summary, we expect that each molecular line will trace a different set of temperature, density, and excitation conditions. Therefore, by using multiple molecular lines it is possible to study the cloud velocity and density structure in low, intermediate and high-density molecular gas. Table 1 summarizes the properties of the molecular lines used in this study. Note that we only use the ground state ($J\,=\,$1 – 0) transition of each molecular line in our analysis.

3 Method

3.1 Velocity Gradients Technique

We use molecular line data from a survey with the 22-meter Mopra telescope covering the Vela C Giant Molecular Cloud, covering 264.62∘ - 266.56∘ in galactic longitude, and 0.54∘ - 1.92∘ in galactic latitude. The observations are is described in detail in Fissel et al. (2019). The beam full width half maximum (FWHM) in arcmin and the thickness of velocity slice in $kms^{-1}$ of different tracers are shown also in Table 1.

Vela C is a massive ($M\,\sim\,10^{5}\,M_{\sun}$), relatively nearby GMC (GAIA-DR2 distance $\sim$ 900 pc), that appears to be relatively young and unevolved (Yamaguchi et al., 1999; Netterfield et al., 2009; Hill et al., 2011). This molecular cloud is mostly cold ($T_{\mathrm{dust}}\,<\,$15 K) and in early stages of star formation, though it has formed a 1 Myr cluster of stars which powers a compact H II region RCW 36 (Ellerbroek et al., 2013). Since Vela C is dominated by a single velocity component in the range 0 - 12 km/s (Fissel et al., 2019), it therefore provides a good opportunity for testing the Velocity Gradient Technique through the study of alignment between VCGs/ VChGs and the magnetic field orientation inferred from polarization data. We infer the orientation of the projected magnetic field by rotating measurements of polarization angle made with the Balloon-borne Large-Aperture Sub-millimeter Telescope for Polarimetry (BLASTPol) at 500 $\mu$m by 90*∘* (Fissel et al., 2016).111 The data reduction and removal of the contribution to the polarized emission from dust in the diffuse ISM surrounding Vela C is described in Fissel et al. (2016). In our analysis, we use only data using the ”Intermediate” diffuse emission subtraction method from Fissel et al. (2016). The beam FWHM of the BLASTPol 500$\mu$m polarization data is 2.5′, which corresponds to $\sim$0.6 pc at the distance to Vela C.

Velocity centroid maps C(x,y) for all tracers are produced by integrating along the velocity axis of the PPV (Position-Position-Velocity) cube:

[TABLE]

where $T_{R}$ is the the radiation temperature (in units of Kelvin), and v is the line-of-sight velocity.

For the analysis of observational data, the theory describing the statistics of fluctuations in the Position-Position-Velocity (PPV) cube is crucial (Lazarian & Pogosyan, 2000). It predicts the anisotropy of the velocity channel maps (Lazarian et al., 2002).

Even when the velocity slice is thin, the channels can record more contribution from turbulent velocities (Lazarian & Pogosyan, 2000; Yuen et al., 2019). Therefore, the gradients of thin velocity channel maps (VChGs) are expected to trace the magnetic field orientation with high accuracy (Lazarian & Yuen, 2018a). Thus, we also investigate the correspondence between VChGs and the magnetic field. To construct the velocity channel maps, we create integrated maps over a narrow velocity range $\Delta v$ satisfying:

[TABLE]

where $\delta v^{2}$ is the line of sight velocity dispersion. This is the criterion proposed in Lazarian & Pogosyan (2000) so that the velocity contribution in the velocity channel map dominates over the density contribution (See Tab 1). Then the velocity channel map can be calculated also by integrating along the velocity axis:

[TABLE]

where $v_{0}$ is the velocity corresponding to the central peak of the velocity profile along the line-of-sight. We follow the gradient calculation algorithm of velocity centroid maps from Yuen & Lazarian (2017a), and velocity channel maps from Lazarian & Yuen (2018a), with the sub-block averaging method (Yuen & Lazarian, 2017a) applied (Note that the sub-block averaging method is not just a smoothing method for suppressing noise in a region, but used to increase the reliability of the important statistical measure.). In this work, the sub-block size is selected as 24 pixels ($\sim 1pc$, 12 arcsec per pixel for CO isotopologues, HCN, HNC, HCO*+*, while 24 arcsec per pixel for CS) which is enough for statistical accuracy. The polarization vector is averaged over the same sub-block size. In addition, Lazarian & Yuen (2018a) suggested a refinement of the velocity gradient technique by using the Moving Window approach which is using a continuous sub-block averaging, rather than dividing up the map into discrete sub-blocks. We apply the Moving Window to all our gradient maps unless specifically mentioned.

In practical observation, the noise is present in spectroscopic data. Lazarian & Yuen (2018a) find that the white noise would alter the alignment of gradients to the projected magnetic field. However, the alignment remains fairly good for data with intensity signal to noise of 3 or greater, and only significantly diverges from the true gradient alignment angle when the signal to noise of the data approaches 1. Therefore, in the case of high noise situation, a Gaussian filter is proposed in Lazarian et al. (2017) to reduce the effect of noise in velocity centroid maps and velocity channel maps. In addition, in order to match the gradient after the sub-block averaging, the polarization data is firstly smoothed to the same angular resolution as molecular data We average the polarization data over the area as the sub-block size in order to match the resolution of the gradient maps. With a sub-block size 24 pixels, the effective resolution of both gradients and polarization is 4.8’ (9.6’ for CS) which is larger than the minimum resolution of polarization 2.5’. Note that neither the polarization nor the spectroscopic data are filtering out scales smaller than 4.8’ (9.6’ for CS), and there could be small-scale magnetic field structures that cannot be resolved with our sub-block averaging method.

3.2 Alignment Measure

The orientation of gradients from VGT is compared with the magnetic filed inferred from BLASTPol polarization. The relative orientations between the 90∘ rotated gradients and project magnetic field directions from polarization angles are measured by the Alignment Measure (AM) used:

[TABLE]

where $\theta_{r}$ is the angular difference between the gradient vector rotated by 90∘ and the magnetic field vector derived from polarization in a single sub-block. $\langle...\rangle$ indicates the average over all sub-blocks.

The range of AM is [-1,1]. There are three important cases for the value of AM:

AM = -1: in this case the gradients are perpendicular to the projected magnetic field; 2. 2.

AM = 0: the gradients are neither perpendicular nor parallel to the magnetic field; 3. 3.

AM = 1: the gradients are parallel to the projected magnetic field, which implies a perfect alignment.

The AM is a parameter used to quantify the relative orientation between two vectors. Note that since AM is insensitive to the sign of the velocity gradient, e.g. AM($\theta$) = AM($-\theta$) , it is advantageous to test the performance of the VGT in terms of the orientation of the plane-of-the-sky magnetic fields, i.e. the larger AM, the better alignment between the measurement from VGT and polarization.

4 Tracing the Magnetic Field using VGT over a larger range of densities

4.1 Low-density tracers: 12CO, 13CO & C18O

12CO, 13CO, and C18O typically trace gas with H2 number density between $10^{2}cm^{-3}-10^{4}cm^{-3}$, which is the typical density for a young self-gravitating molecular cloud (Richard M. Crutcher, 2012). Because their lines tend to have lower optical depths, 13CO and C18O can trace molecular gas over a larger range of densities, while optically thick 12CO traces the lower density outer cloud regions. In Yuen & Lazarian (2017a), we see that by combining the gradient vectors and polarization measurements from the same region, we can determine whether a region is optically thick or thin, and whether the self-gravity is dominant in a particular region or not.

Fig. 1 shows the normalized distribution of the relative orientation between rotated VCGs/VChGs and the inferred magnetic field from polarization using 12CO, 13CO, and C18O as tracers. The distribution is constructed by calculating gradients for each map pixel, without using sub-block averaging. We find that the set of rotated gradient alignment angles is roughly consistent with a Gaussian distribution. The expectation values of VCGs are located near 13.5o, while for VChGs the expectation values are decreasing with the increasing of the critical density of each tracer. Considering the systematic uncertainty in polarization data is less than 10∘ (see Appendix A), the results support our theoretical consideration that the velocity gradient rotated by 90∘ tends to align with the magnetic field.

Fig. 2 shows the orientation of VCGs, VChGs, and the magnetic field obtained from polarization with sub-block averaging method applied. Fig. 2 also shows that using C18O VChGs reflects the central structure of the cloud, while 12CO and 13CO determine the outskirts structure. Considering the systematic uncertainty in polarization data is approximately 2.07∘ (see Appendix A), these three tracers shows similar agreements with the magnetic field inferred from polarization data based on the estimation of AM.

4.2 High-density tracers: HCN, HNC, HCO+ & CS

Molecular CS is well known for tracing dense clumps in the molecular clouds in the Milky Way. Also, recent multi-line surveys of nearby spiral galaxies (Bayet et al., 2009) of this molecule in extragalactic environments have revealed low J (J $\leq$ 4) transitions of CS trace gas densities of the order $10^{5}cm^{-3}$, even the number is close to $10^{6}cm^{-3}$. Estimates of the effective excitation density of CS from Shirley (2015) show that cold CS J = 1 – 0 typically traces densities of $\sim$104 cm*-3*.

HCN and HNC are two species that have routinely been used as tracers of star formation regions in molecular clouds (Turner et al., 1997). In fact, HCN is used to trace the same approximate densities as the low-J CS molecule, and hence the gas that it traces is not necessarily in a self-gravitating star formation region. HNC has a similar critical density to that of HCN. Therefore, HNC and HCN should be good tracers for higher density star formation gas which might be self-gravitating. Pety et al. (2017) showed that HNC is a better molecular tracer than HCN when probing low-J lines, based on the analysis of visual extinction and line integrated intensities. We note that the HCN J = 1 – 0 has hyperfine structure which complicates the interpretation of the VCGs. Also the HCN emission has a lower integrated line emission compared to HNC (Fissel et al., 2019).

Compared to the other three tracers, HCO*+* is more often to be used as a tracer of ionized gas (e.g., see Cleeves et al. 2015). However, 80% of the HCO*+* and HCN emission originates in non-self-gravitating regions of molecular clouds (Vollmer et al., 2017). HCN and the HCO*+* emission lines can therefore both be used to trace the dense gas around the self-gravitating regions.

Fig. 3 shows the normalized distribution of the relative orientation between rotated VCGs/VChGs and the inferred magnetic field from polarization, using CS, HCN, HNC, and HCO*+. We do not apply sub-block averaging when constructing the distributions, but use the raw gradients for each pixel. We find that the distribution is consistent with a Gaussian distribution. However, the expectation values of VChGs are $\sim$5.5∘, lower than the expectation values for the 12CO, 13CO, and C18*O lines. CS and HNC show smaller expectation values of VCGs compared to HCN and HCO+.

Fig. 4 and Fig. 5 show the orientation of VCGs, VChGs using HCN, HNC, HCO*+*, and CS, and the magnetic field obtained from polarization with sub-block averaging method applied. Fig. 5 shows that CS is able to detect a clear structure of Vela C by using VChGs, and shows perhaps the best agreement with the polarization inferred magnetic field orientation, with an AM=0.70. Among four tracers, HNC and CS appear to best trace the magnetic field.

4.3 Contribution from the Foreground and Background

In a low galactic latitude cloud, such as Vela C, foreground and background material can change the observed cloud magnetic field orientation. In contrast, the VGT method for tracing magnetic fields is only sensitive to the molecule probed, and molecular line cubes can be studied solely over the velocity range associated with the molecular cloud of interest. Using the VGT it should therefore also be possible to trace the magnetic field in different layers of the cloud, by targeting lines that trace different ranges of densities. Thus, VGT using molecular tracers provides the information of the localized magnetic field corresponding to species in the cloud, while polarization accumulates the information along the line of sight, i.e. VGT and polarization are tracing different components of the magnetic field. If we want to fairly compare VGT and polarization, we should consider the contribution from the foreground and background222 Note that Fissel et al. (2019) shows that Vela C is dominated by a single velocity component in the range 0 - 12 km/s, while there might be multiple velocity components in a wider velocity range. For a fair comparison with polarization, we did not constraint the PPV cube to a particular velocity range, but keep its original velocity range -20 - 30 km/s for the presented results. We repeat our VCGs and VChGs analyses but limit the velocity range as 0 - 12km/s. The AM (with error bar $\sim\pm 0.03$) is generally 0.05 lower than the presented results. AM values for single/multiple component analyses still give similar values. We, therefore, expect the bulk of the signal in both data can be reasonably well compared and is dominated by signals from Vela C, although both polarization and spectroscopic data are not entirely constrained to Vela C.

In Fissel et al. (2016) the authors attempted to remove the contribution of the foreground and background dust to the BLASTPol 500$\,\mu$m polarization data by using two different methods to model the diffuse polarized dust emission, and subtracted their model from the data. However the systematic uncertainty associated with removing the diffuse polarized emission is expected to most affect polarization angles towards regions of low column density within Vela C. We note that several locations where the BLASTPol inferred magnetic field direction differs from the velocity gradient inferred field direction are towards the edge of the map in low column density regions of Vela C.

Tab.2 shows the expectation value $\mu$ of the relative angle between the rotated gradients of each molecular tracer and the magnetic field inferred from dust polarization. From Table.2, we find that low-density tracers, i.e. 12CO and 13CO, indeed show larger deviation between the rotated VChGs and magnetic fields inferred from the polarization than high-density tracers. As mentioned above, for low-column density regions we expect a larger systematic uncertainty in the polarization angle associated with background/foreground polarized emission subtraction, so it is not surprising that 12CO and 13CO, which show emission across the entire mapped area, show a larger deviation from the VChG and VCG inferred magnetic field orientation.

In addition, we find that VChGs show approximately $\mu\sim 6^{\circ}$ offsets for dense tracers, which are smaller than the offsets of VCGs. Noise is one possible factor contributing to the larger deviations of the VCGs compared to VChGs, for example, the map of VCGs is integrating over the whole velocity range of the PPV cube, while for VChGs, the map is integrating over a narrow velocity range $\Delta v$. There is also the second possibility for the worse alignment measures for the VCG maps. Although Vela C is dominated by a single velocity component, there are multiple velocity components in some parts of Vela C (see Fig. 1 in Fissel et al. (2019)). The VCG analysis implicitly assumes only one velocity component. We note also that the HCN J = 1 – 0 line has hyperfine structure, which makes calculating line velocity centroids more difficult. In addition, RCW 36 is a compact H II region on Vela C. However, the H II region is compact ($\sim$1pc in size) and as such only compromises a small portion of the map Fissel et al. (2016). The contribution from RCW 36 is therefore slim.

5 Tracing the Magnetic Field from low-density gas Region to High-density gas Region

5.1 Low-density Gas Region

Fig. 2 shows that the integrated velocity channel maps, i.e. the moment = 0 maps, are able to trace the cloud’s column density structure. VChGs of 12CO and 13CO appear to mostly trace the low to intermediate density material of the cloud, as they become optically thick towards high column density cloud sightlines.

Due to the lower optical depth of C18O, its emission lines traces deeper into the cloud, and therefore includes cloud structures over a larger range of densities. In this case, we can use 12CO, 13CO to get a fairly good idea of the magnetic field structure in the outer layers of the cloud and use C18O to better trace higher-density structure.

We therefore propose a method that uses multiple molecular lines as a combined tracer for magnetic field in GMCs, even if part of the cloud are strongly self-gravitating. We separate each moment = 0 map of the Vela C Cloud into the low-intensity regions and the high-intensity regions. The high-intensity region is defined by the percentile of the integrated line intensity: we take as the high-intensity region all pixels that are above the 95th percentile in integrated line intensity, while our definition of low-intensity gas is all the pixels that are below this same threshold.

Fig. 6 compares the alignment measure for the low-intensity regions and high-intensity regions ( high-intensity regions for each tracer are indicated with green contours in Fig .2, Fig. 4, and Fig. 5), traced by 12CO, 13CO, and C18O. We find that for VChGs, 12CO shows little alignment towards the high-intensity region while C18O shows better alignment in high-intensity regions than in low-intensity regions. The weak integrated line intensities of C18O towards low column density cloud regions likely contributes to the lower AM values observed towards low column density sight lines.

5.2 High-density Gas Region

As is described above, the CS, HNC, HCN, and HCO*+* are used as dense gas tracers. Although HCN and HCO*+* are usually generated in non-self-gravitating regions, they can be used to trace gas that is near the self-gravitating regions. Lower density tracers 12CO 13CO and intermediate density tracer C18O can also be applied to relatively diffuse regions that are near the strong self-gravitating molecular core.

From Fig. 5 and Fig. 4, it is obvious that both CS and HNC show excellent alignment between gradients and the magnetic field. Fig. 7 shows the cumulative plot of the alignment measure, when using VChGs to trace the magnetic field by CS, HCN, HCO+, and HNC. CS and HNC show a very strong increase in the cumulative fraction over the AM range [0.5,1.0].

Therefore, HNC and CS give a more robust performance than others. We therefore propose to use HNC or CS to trace the magnetic field in the dense gas regions (i.e. H2 density $\sim 10^{4}-10^{6}cm^{-3}$) of GMCs.

Furthermore, we propose including HNC to the inferred magnetic field orientation map made from a combination of multiple molecular line tracers. 12CO, 13CO, and C18O provide the structure of the magnetic field towards low and intermediate column density gas. Taking the dense region traced by HNC which traces higher densities on average than CS into account, we can, therefore, infer the plane of sky component of the magnetic field orientation over a wide range of cloud densities in GMCs.

Fig. 8 shows the magnetic field structure of Vela C obtained from VChGs using the combination of multi-tracers 12CO, 13CO, C18O and HNC. The map is produced as follows: using velocity channel maps, we compare the structure contours333 The structure contour is a curve connecting points which have the same particular value. We use the structure contour to highlight the main structure of Vela C for each tracer, as shown in Fig. 2. The value selected for 12CO is 50K$\cdot$km/s and above, 15K$\cdot$km/s for 13CO, 3K$\cdot$km/s for C18O, 5K$\cdot$km/s for CS, 4K$\cdot$km/s for HNC, while 2K$\cdot$km/s for HCN and HCO+. as discussed above from different molecular tracers and remove the gradients calculated from the tracer has lower critical density in the overlap region but keep the gradients calculated from the tracer has higher critical density. The combination of multiple-tracers shows a better alignment measure over all than any single molecular line tracer.

5.3 Re-rotation Test

In extremely high density cloud regions, self-gravity can be the main force affecting the dynamics of the gas in GMCs. We expect the gradients to be aligned parallel to the magnetic field with the presence of gravitational collapse, i.e., the direction of the gradient shoud flip 90∘ (Yuen & Lazarian, 2017b). In this section we search for the signature of gravitational collapse in Vela C using VCGs and VChGs.

We explore this by re-rotating the gradients by 90∘ again at the high-intensity cloud regions. This is equivalent to not rotating the gradients by 90∘ towards high intensity regions, but only low intensity regions. For 12CO, 13CO, and C18O, we vary the threshold for the definition of high-intensity region from the 75th percentile in integrated line intensity, for which AM is close to 0, to the 100th percentile which means that there is no re-rotation. For CS, HCN, HNC, and HCO*+, we change the range of high-intensity region from 90th* percentile, since with the threshold at the 90th percentile the AM is already close to 0.

Fig. 9 shows the AM variation with different moment = 0 map percentile thresholds defining the re-rotated region. It shows that the peak value of AM among all tracers is achieved at 100th percentile which means no re-rotation is required for both VCGs and VChGs to get good alignment with polarization. Thus the high-intensity region contributes positively to the whole AM value, among all tracers. Therefore we conclude that strongly self-gravitating, collapsing regions constitute only a small fraction of the cloud area within Vela C. This, however, does not prevent molecular gas in small scales regions below the resolution limit of the Mopra data to be collapsing and forming stars (Hu et al., 2019a).

6 Discussion

6.1 Extracting Magnetic Field Orientation for Gas in Different Density Regimes

Molecular tracer maps with different optical depth provide the spectroscopic information of gas dynamics up to certain line-of-sight depth. The concept of ”gradient tomography” was first discussed in Lazarian & Yuen (2018) by considering the effective accumulation line-of-sight deepness of synchrotron polarization data with the different wavelength. Both the synchrotron polarization data with the presence of strong Faraday Rotation effect and the gas spectroscopic data with the presence of optically thick radiative transfer effect share the same concept that the contribution of gas dynamics with line-of-sight deepness larger than some certain physical boundary would be effectively noise. Lazarian & Yuen (2018) showed that, by stacking the gradient maps from the polarized synchrotron intensities measured from different frequencies, one can create the 3D tomography information of the magnetic field. The number of layers in the gradient tomography completely depends on how many individual frequency measurements one has taken for the synchrotron data. Hsieh et al. (2019) implemented analogous idea in the case of multiple molecular tracer maps using SPARX radiative transfer code. For our analysis of the Vela C data each molecular tracer samples a different range of densities. Assuming the excitation temperature of Vela C is approximately 10K and CMB brightness temperature 2.725K, Fissel et al. (2019) then calculated the optical depth $\tau_{18}$ for C18O typically ranges from 0.015 to 0.18, with a median value of 0.026. Assuming a [13CO/C18O] ratio of 10 and a [12CO/ C18O] ratio of 400, this implies a typical $\tau_{12}$= [12CO/ C18O]$\tau_{18}$ in the range of 6 to 72, and $\tau_{13}$ in the range of 0.15 to 1.8. We can therefore expect that the velocity gradients tell us about the plane of sky component of the magnetic field over different density ranges, as shown Fig. 10.

If the velocity gradient can be used to determine the magnetic field orientation over different cloud density regimes, this also has implications for the efficiency of dust grain alignment as a function of density. Dust grains in molecular clouds are thought to be aligned by radiative alignment torques (RATs) (see Lazarian & Hoang 2007), and in deeply embedded cloud regions the alignment efficiency may be lower as the photons of the interstellar radiation field with the wavelength comparable with the grain size would be selectively extincted. However, some observational studies and numerical simulation indicate that grains can be effectively aligned for moderately extincted dust sightlines (Bethell et al., 2007; Alves et al., 2014; Fissel et al., 2016). Moreover, the radiation of the embedded stars can play an important role for the RAT alignment (Whittet et al., 2008). The fact that the BLASTPol-inferred magnetic field orientation shows a higher degree of alignment with the VGT of intermediate or high density tracers towards high column density sightlines would seem to indicate that dust grains are efficiently aligned in molecular gas with number density $\approx\,10^{4}\,$cm -3 or greater, and therefore that dust polarization can trace the magnetic field of intermediate and high density gas. While in numerical simulation it is confirmed that the grain align efficiency is still high fith number density $\geq 10^{3}$ (Seifried et al., 2019), this issue requires further studies with higher resolution molecular line observations and polarization data.

6.2 Applicability of the velocity gradient technique to other molecular clouds

The gradient technique, which shows very good alignment with the magnetic field orientation obtained from polarization data, is a promising method to infer the magnetic field morphology over a wide range of physical conditions. Low resolution polarization data is now universally available from the Planck all-sky data, while higher resolution polarimetry data has been obtained from stratospheric and ground-based instruments.

However, dust polarization may not be sensitive to the magnetic field within deeply embedded high density regions, where as the efficiency of grain alignment due to radiative torques is expected to be less efficient (Lazarian & Hoang, 2007). In addition, it is very difficult to determine the magnetic field structure from polarization data when there are multiple dust clouds along the line of sight. However, gradient technique can make use of readily available large scale molecular line surveys, such as CHAMP (Stolle et al., 2006) and ThRuMMs (Nguyen et al., 2015) survey, or (HI4PI Collaboration et al., 2016) neutral hydrogen atom distribution survey and the COMPLETE survey (Ridge et al., 2006), to provide measurements of the magnetic field orientations, which can also be cross-checked by using different tracers towards the same cloud.

7 Conclusion

As one of the most versatile methods for probing magnetic field in the molecular clouds, the gradient method has been proposed to trace magnetic field orientations in multiple scenarios, including diffuse media, shocks, and self-gravitating regions. By examining the gradients of seven chemical tracers, and judging from the distribution of gradient vectors that are rotated by 90∘, we present here an example making uses of the recent advancement of gradient techniques from numerical and observational studies. Moreover, this work also suggests the promise of using multiple tracers with different tracing abilities to probe magnetic fields in different scales according to the density of the molecular clouds. To summarize:

The Velocity Gradients Technique, which is for the first time applied to multiple tracers, is able to trace the variations in the plane-of-sky component of the magnetic field within different density regimes of the Vela C cloud. 2. 2.

The Velocity Gradients Technique opens new a way of exploring the localized magnetic field in Giant Molecular Clouds without the systematic uncertainty inherent in polarization observation associated with subtracting out foreground and background dust emission. 3. 3.

We see no column density threshold above which re-rotation of the velocity gradients is required to match the BLASTPol inferred magnetic field orientation. This implies that at the resolution of the Mopra data, we do not see the expected signature of the velocity gradients aligning parallel to the magnetic field in regions of gravitational collapse. We therefore infer that the collapsing regions constitute a small fraction of the Vela C cloud.

Acknowledgements. A.L. acknowledges the support of the NSF grant AST 1715754, and 1816234 as well as a distinguished visitor PVE/CAPES appointment at the Physics Graduate Program of the Federal University of Rio Grande do Norte, the INCT INEspao and Physics Graduate Program/UFRN. K.H.Y. acknowledges the support of the NSF grant AST 1816234. Y.H. acknowledges the support of the Future Generation Fellowship by Hu’s family. L.M.F. is a Jansky Fellow of the National Radio Astronomy Observatory (NRAO). NRAO is a facility of the National Science Foundation (NSF operated under cooperative agreement by Associated Universities, Inc.). The Mopra radio telescope is part of the Australia Telescope National Facility, which is funded by the Australian Government for operation as a National Facility managed by CSIRO.

Appendix A Statistical error in polarization data

We show the statistical angle uncertainty in polarization data in Fig. 11. The errors are not constant across the BLASTPol map and the error is higher when it gets close to the boundary of Vela C. However, the mean error for the polarization data used in comparison with the CO isotopologues data is $\sim 2.07^{\circ}$, while $\sim 1.44^{\circ}$ for the one used for CS, and $\sim 1.37^{\circ}$ for the one used for HCN, HNC, and HCO+. The statistical polarization angle errors depend on the signal to noise of the polarization fraction. For Gaussian polarization errors, the error $\sigma_{\phi}\sim 28.6^{\circ}\times\sigma_{p}/p$, where $\phi$ is the polarization angle and p is the polarization fraction (Fissel et al., 2016). Thus, for our polarization data, which is a three-sigma polarization detection, the statistical error is less than 10 ∘ (Fissel et al., 2019).