Relative Alignment Between the Magnetic Field and Molecular Gas Structure in the Vela C Giant Molecular Cloud using Low and High Density Tracers

For the analysis in this work we utilize the magnetic field orientation inferred from linearly polarized dust emission measured by the BLASTPol balloon-borne polarimeter, during its last Antarctic science flight in December 2012 (Galitzki et al., 2014). BLASTPol observed Vela C in three sub-mm bands centered at 250, 350, and 500 $\mu$m, for a total of 54 hours. Due to a non-Gaussian telescope beam the maps required additional smoothing. In this paper we focus solely on the 2$\farcm$5-FWHM-resolution 500-$\mu$m maps previously presented in Fissel et al. (2016).¹¹1We note, however, that the inferred magnetic field orientation angles are largely consistent between the three BLAST bands, as discussed in Soler et al. (2017). This resolution corresponds to 0.7 pc at the distance of Vela C.

We assume that the orientation of $\mathbf{\langle\hat{B}_{\perp}\rangle}$, the magnetic field orientation projected on the plane of the sky, can be calculated from the Stokes parameters as

which corresponds to the polarization orientation $\mathbf{\hat{E}}$ derived from the BLASTPol 500 $\mu$m Stokes $Q$ and $U$ data rotated by $\pi$/2 radians.²²2In our coordinate system a polarization orientation angle of 0^∘ implies a Galactic North-South orientation, where the angle value increases with a counter-clockwise rotation towards Galactic East-West. Only BLASTPol measurements with an uncertainty in the polarization angle of less than 10^∘ are used in this analysis.

Fissel et al. (2016) discussed the several different methods for separating polarized emission due to diffuse ISM dust along the same sightlines as Vela C. This correction is important as the Vela C cloud is at a low Galactic latitude ($b$ = 0.5–2^∘). For our analysis, we use the “Intermediate” subtraction method from Fissel et al. (2016). In Appendix B.1, we show that the choice of diffuse emission subtraction method does not change our final results.

To study the density and velocity structure of Vela C we compare the BLASTPol data to results from a large-scale molecular line survey of Vela C made with the 22-m Mopra Telescope over the period from 2009 to 2013. The Mopra data presented here are the combination of two surveys: M401 (PI: Cunningham), which covered molecular lines at 3, 7, and 12 mm, and M635 (PI: Fissel), which mapped Vela C in the $J$ = 1 $\rightarrow$ 0 lines of ¹²CO and isotopologues ¹³CO and C¹⁸O. For the M401 observations the cloud was mapped in a series of square raster maps (5′, 10′, and 15′ respectively, for the 3-, 7- and 12-mm observations), while the M635 observations were taken using the Mopra fast-scanning mode, scanning the telescope in long rectangular strips of 6′ height in both the Galactic longitude and latitude directions.

For both surveys the UNSW-MOPS³³3The University of New South Wales Digital Filter Bank used for the observations with the Mopra Telescope was provided with support from the Australian Research Council. digital filterbank backend and the MMIC receiver were used, with multiple zoom bands covering 137.5 MHz each, with 4096 channels within the 8-GHz bandwidth. In this paper we present observations of the nine molecular rotational lines for which there is significant extended emission: the ¹²CO, ¹³CO, C¹⁸O, N₂H⁺, HNC, HCO⁺, HNC, and CS $J$ = 1 $\rightarrow$ 0 lines, as well as the NH₃(1,1) inversion line. Table 1 summarizes the observed lines including velocity resolution and beam FWHM $\theta_{\mathrm{beam}}$, which ranges from 33″ FWHM for the CO $J$ = 1 $\rightarrow$ 0 observations to 132″ FWHM for NH₃ (1,1). Our Mopra observations were bandpass corrected, using off source spectra with the livedata package, and gridded into FITS cubes using the gridzilla package⁴⁴4http://www.atnf.csiro.au/computing/software/livedata/index.html. Extra polynomial bandpass fitting was done with the miriad package,⁵⁵5http://www.atnf.csiro.au/computing/software/miriad/ and Hanning smoothing was carried out in velocity.

We compare the observed molecular line emission to the total hydrogen column density map $N_{\mathrm{H}}$ (in units of hydrogen nucleons per cm^-2) first presented in Section 4 of Fissel et al. (2016).⁶⁶6Note that in this paper $N_{\mathrm{H}}$ and $n_{\mathrm{H}}$ refer respectively to the column density and number density of hydrogen nucleons, while $N_{\mathrm{H_{2}}}$ and $n_{\mathrm{H_{2}}}$ refer to the molecular hydrogen column and number density. Assuming all of the hydrogen is in molecular form at the densities probed in this work the conversion is $n_{\mathrm{H_{2}}}\thinspace=\thinspace n_{\mathrm{H}}/2$. These maps are also used in Section 4.3 and Appendix C to estimate the abundances of our observed molecules. The $N_{\mathrm{H}}$ maps are based on dust spectral fits to four far-IR/sub-mm dust emission maps: Herschel-SPIRE maps at 250, 350, and 500 $\mu$m; and a Herschel-PACS map at 160 $\mu$m. Each Herschel⁷⁷7Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. dust map was smoothed to match the BLASTPol 500-$\mu$m FWHM resolution of 2$\farcm$5 before spectral fitting.

To further explore the emission and line of sight velocity structure of Vela C we calculate the first three moment maps for the cloud. The zeroth-moment map is the integrated line-intensity:

where $T_{\mathrm{R}}$ is the radiation temperature in velocity channel $v$. $T_{\mathrm{R}}$ can be calculated from the measured antenna temperature $T_{\mathrm{A}}$ corrected by the main beam efficiency for extended structure $\eta_{\mathrm{xb}}$ values determined from previous Mopra observations and listed in Table 1 ($T_{\mathrm{R}}$ = $T_{\mathrm{A}}/\eta_{\mathrm{xb}}$). Here $v_{0}$ and $v_{1}$ are the minimum and maximum velocities over which the line data are integrated. These velocity integration limits are listed for each line in Table 1, and are generally within the 0 km s${}^{-1}<v_{\mathrm{LSR}}<$ 12 km s^-1 range where the molecular line emission is likely associated with Vela C. For the HCN and N₂H⁺ $J$ = 1 $\rightarrow$ 0 lines we integrate over a larger velocity range to include additional hyperfine spectral components and increase the signal-to-noise.

We can use higher order moments to study the velocity structure of each data cube. The first moment map gives the intensity weighted average line-of-sight velocity $\left<v\right>$:

Similarly where the signal-to-noise of the line data is high enough we can calculate the second moment, which gives the line of sight velocity dispersion $\Delta v$:⁸⁸8Note that the second moment is written as $\sigma_{v}$ in some publications, but in this paper we use $\Delta v$ to avoid confusion with the measurement uncertainties which are labeled with $\sigma$.

Note that in the case of a Gaussian line profile Equation 4 would return the Gaussian width ($\sigma$).

Before calculating the moment maps we first smooth each channel map with a 2-D Gaussian kernel so that the resulting cube has 120″ FWHM resolution.⁹⁹9The exception is the NH₃ cube, which has an intrinsic FWHM resolution of 132″. For this cube we smooth instead to 150″. This smoothing is needed both to increase the signal to noise ratio of the extended structure and to minimize any narrow spurious map features due to differences in $T_{\mathrm{sys}}$ levels within the map; we show in Appendix B.2 that the choice of smoothed resolution does not significantly change our final results. Table 1 lists the smoothed FWHM resolution ($\theta_{\mathrm{sm}}$) for each cube. The pixel size for the smoothed cubes is the same as the pixel size in the original data cubes (see the last column in Table 1).

To estimate the uncertainty in $T_{\mathrm{R}}$ we select velocity channels in the spectra that have no apparent signal, and find both the standard deviation of all the voxels and the standard deviation for each pixel over all the velocity channels that have no signal. We take as the per velocity channel uncertainty $\sigma_{T_{\mathrm{R}}}$ the maximum of these two standard deviations for each pixel. Uncertainties in the moment maps $\sigma_{I}$, $\sigma_{\left<v\right>}$, and $\sigma_{\Delta v}$, are then estimated through a Monte Carlo method by taking the data cube and adding to each voxel in the cube a random number selected from a normal distribution centered at 0 K, with a width of $\sigma_{T_{\mathrm{R}}}$. We recalculate the moment maps using this method 1000 times, and take the per-pixel standard deviation in the resulting moment maps to be our uncertainty.

For the analysis of the $I$ (zeroth-moment) maps we only use data points where $I/\sigma_{I}\thinspace>$ 8, except for the N₂H⁺ and NH₃ maps, which have relatively low signal-to-noise, where we relax the signal-to-noise requirement to 6 and 5 respectively. For the $\left<v\right>$ (first moment) maps we additionally require that the uncertainty $\sigma_{\left<v\right>}$ be less than 0.4 km s^-1. More strict criteria are applied for the calculation of the $\Delta v$ (second moment) maps, which are very sensitive to noise spikes. Here we only use spectral channels where $T_{\mathrm{R}}\thinspace\geq\thinspace 3\sigma_{T_{\mathrm{R}}}$ and require the integrated line strength to be above a threshold $I$ signal-to-noise level listed for each line in Table 1.

Figures 3 and 4 show the calculated moment maps of nine different molecular lines, with contours of $N_{\mathrm{H}}$ (Section 2.3). In general the molecular lines appear to trace different density, chemical, and excitation conditions within the cloud. The ¹²CO $J$ = 1–0 $I$ map shows little correspondence to the column density structure of Vela C, which is consistent with the expectation that the emission is optically thick, such that only the lower density outer layers are probed by the line.

We expect ¹³CO to have a lower optical depth than ¹²CO. The $I$ map of ¹³CO shows similar structure to $N_{\mathrm{H}}$, but does not show the dense filamentary structure seen in the Herschel observations. The even rarer isotopologue C¹⁸O shows a very similar structure to ¹³CO, although with lower signal-to-noise ratio and more contrast towards the highest column density regions where ¹³CO might be optically thick.

The HNC, HCO⁺, HCN, and CS $J$ = 1 $\rightarrow$ 0 lines show significant $I$ detections only toward higher column density structures. We note that these intermediate number density tracers show weaker emission in the Centre-Ridge sub-region to the right of RCW 36 compared to the Herschel-derived $N_{\mathrm{H}}$ map (contours in Figures 3 and 4). This could imply that molecular abundance or excitation conditions are different in the Centre-Ridge compared to the rest of Vela C. We also include two tracers that are often used to probe higher density gas, NH₃ (1,1) and N₂H⁺ $J$ = 1 $\rightarrow$ 0. These lines tend to have low signal to noise ratios (Figure 2), but show emission near the highest column density cloud regions.

Throughout the paper we refer to ¹²CO and ¹³CO as “low density” tracers because these molecules are optically thick towards high $N_{\mathrm{H}}$ sightlines and have high enough abundance levels to be detected in the low density envelope of Vela C. We refer to N₂H⁺, HNC, HCO⁺, HCN, CS $J$ = 1 $\rightarrow$ 0, and NH₃ (1,1) as intermediate or high density tracers because these molecules trace mostly higher column density regions, are not generally detected in the cloud envelope, and tend to have higher estimated characteristic densities (see discussion in Section 4.3). C¹⁸O $J$ = 1 $\rightarrow$ 0 is also only detected toward higher column density structures, however radiative transfer modeling in Section 4.3.1 suggest the line typically traces lower densities than our intermediate or high density tracers.

The first moment or $\left<v\right>$ maps within the $N_{\mathrm{H}}$ contours show that the molecular gas of Vela C has on average a line of sight velocity 1–2 km s^-1 higher in the Centre-Ridge compared to the rest of the cloud. As discussed in Section 2.3, many cloud sightlines, particularly towards the South-Nest and Centre-Nest sub-regions, have multiple spectral peaks centered at different line of sight velocities. Some of the structure in the $\left<v\right>$ maps is therefore likely the result of variations in the relative intensity of the different spectral components that contribute to the total cloud sightline emission. In addition, the HCN, N₂H⁺, and NH₃ lines have hyperfine structure, and so $\left<v\right>$ maps calculated for these lines could be influenced by the optical depth of the different hyperfine components.

The second moment $\Delta v$ maps show large apparent velocity dispersions where there are two nearly equal strength spectral peaks at different line of sight velocities (for example location A in Figures 1 and 2). For the C¹⁸O and the intermediate to high density tracers HNC, HCO⁺, and CS, which do not have hyperfine line structure, we see that the two “nest-like” regions identified in Hill et al. (2011) show much larger average values of $\Delta v$ than the two “ridge-like” regions. This suggests that in addition to having filamentary structure with a variety of orientations, the South-Nest and Centre-Nest also have more complicated line of sight velocity structure than the South-Ridge and Centre-Ridge regions.

Similar to the methods described in Soler et al. (2013, 2017), the orientation of structure in the Mopra moment maps is calculated by computing the gradient vector field of each map. The moment map is convolved with a Gaussian gradient kernel of FWHM width $\theta_{\mathrm{gr}}$, where $\theta_{\mathrm{gr}}$ was chosen to be larger than three map pixels to avoid spurious measurements of the gradient orientation due to map pixelization (see Table 1).

The relative orientation angle $\phi$ between the plane of the sky magnetic field $\mathbf{\langle\hat{B}_{\perp}\rangle}$ and a line tangent to the local iso-$I$ map contour is equivalent to the angle between the polarization direction $\mathbf{\hat{E}}$ and $\nabla\thinspace I$:

(Soler et al., 2017). With this convention, $\phi\thinspace=\thinspace$0^∘ indicates that the magnetic field and local $I$ structure orientations are parallel, while $\phi\thinspace=\thinspace$90^∘ indicates that $\mathbf{\langle\hat{B}_{\perp}\rangle}$ is perpendicular to the local $I$ structure. Because dust polarization can be used to measure only the orientation of the magnetic field, not the direction, the relative orientation angle $\phi$ is unique only within the range $\left[0^{\circ},\thinspace 90^{\circ}\right]$. That is, $\phi\thinspace=\thinspace$20^∘ is equivalent to both $\phi\thinspace=\thinspace-$20^∘ and $\phi\thinspace=\thinspace$160^∘.

We calculate the relative orientation angle $\phi$ for each Mopra molecular line $I$ map, sampling our data at the location of every Mopra map pixel (see Table 1 for pixel size information). Figure 5 shows the histograms of relative orientation (HROs). The black solid line shows the normalized histogram for all values of $\phi$ that have passed the $I$ map cuts described in Section 3.1 and have

where $\sigma_{\nabla I}$ is the measurement uncertainty of the gradient angle and $\sigma_{\mathbf{\hat{E}}}$ is the measurement uncertainty of the polarization angle.

The HROs for the nine observed molecular lines show different trends. For the ¹²CO HRO there are significantly more sightlines where the $I$ structure is parallel to the magnetic field than perpendicular. The other molecular lines show either slightly more sightlines parallel than perpendicular (¹³CO), a flat HRO indicating no preferred orientation with respect to the magnetic field (C¹⁸O and HCO⁺), or more sightlines perpendicular to the magnetic field than parallel (HCN, HNC, CS, N₂H⁺, and NH₃).

We test for changes in the shape of the HRO with $I$ by dividing our sightlines into seven bins based on their $I$ values, with the bins chosen such that each has the same total number of sightlines. Figure 5 shows no consistent trends in the shape of the HRO for different $I$ bins, in contrast with the Soler et al. (2017) application of HRO analysis to Herschel-derived column density maps, where there was a clear transition to a more perpendicular alignment with increasing column density. This could imply that our $I$ maps are not a direct proxy column density, or the difference could be due to the low resolution of our Mopra $I$ maps compared to the 36″ FWHM resolution $N_{\mathrm{H}}$ maps used in Soler et al. (2017). We discuss the change in relative orientation versus column density in Section 5.1.

As discussed in Jow et al. (2018) given a set $\{\theta_{i}\}$ of n independent angles distributed within the range [0, 2$\pi$], the Rayleigh statistic $Z$ can be used to test whether the angles are uniformly distributed

where $n_{\mathrm{ind}}$ is the number of independent data samples. This equation is equivalent to a random walk, with $Z$ characterizing the displacement from the origin if one were to take steps of unit length in the direction of each $\theta_{i}$. If the distribution of angles is uniformly random then the expectation value for $Z$ is zero.

To test for preferential parallel or perpendicular alignment we take $\theta\thinspace=\thinspace 2\phi$, where $\phi$ is the relative orientation angle calculated as described in Section 4.1. Here $\theta\thinspace=\thinspace 0$ corresponds to parallel alignment, while $\theta\thinspace=\thinspace\pi$ corresponds to perpendicular alignment. Jow et al. (2018) showed that the projected Rayleigh statistic (PRS) $Z_{x}$ can be used to test for a preference for perpendicular or parallel alignment:

$Z_{x}$ in Equation 8 represents the random walk component projected on the x-axis in a Cartesian coordinate system. If a measurement of $\mathbf{\langle\hat{B}_{\perp}\rangle}$ is parallel to the local iso-$I$ map contour then $\cos{\theta_{i}}$ = 1. If the two orientations are perpendicular then $\cos{\theta_{i}}$ = $-$1. Jow et al. (2018) used Monte-Carlo simulations to show that for uniformly distributed samples of {$\phi_{i}$} the expectation value of $Z_{x}$ converges to 0 with $\sigma_{Z_{x}}$ = 1. We also note that $Z_{x}$ in Equation 8 will increase proportionally to $n_{\mathrm{ind}}^{1/2}$. The PRS can therefore be thought of as quantifying the significance of a detection of relative orientation. Measurements of $Z_{x}$ $\gg$ 1 indicate a significant detection of parallel relative alignment, while measurements of $Z_{x}$ $\ll$ $-$1 indicate a strong detection of perpendicular relative alignment.

Under the assumption that the uncertainty is dominated by the sample size, rather than by the measurement errors associated with the BLASTPol polarization angles or $I$ gradient angles, the variance of the $Z_{x}$ is

(Jow et al., 2018). For the null hypothesis of a uniform distribution of angles (no alignment), $\sigma_{Z_{x\thinspace\mathrm{uni}}}=1$, which is the standard against which $Z_{x}$ is tested. The behavior and convergence of the Rayleigh statistic and projected Rayleigh statistic are examined in detail in Jow et al. (2018).

In practice finding $Z_{x}$ for the set of relative orientation angles between the BLASTPol data and Mopra $I$ maps (as calculated in Equation 5) is complicated by the fact we measure $\theta_{i}$  for every map pixel, therefore our data is highly oversampled. In Table 2 we list the oversampled PRS $Z_{x}^{\prime}$ calculated for our measurements of {$\theta_{i}$} as

where $n_{pix}$ is the number of map pixels of size indicated in Table 1.

To correct for oversampling we calculate $Z_{x\thinspace WN}^{\prime}$ for a series of relative orientation angles $\{\phi_{\mathrm{WN}\thinspace i}\}$ where we replace $\nabla I$ in Equation 5 with $\nabla I_{\mathrm{WN}}$. $I_{\mathrm{WN}}$ is a white noise map smoothed to the same resolution as the Mopra $I$ maps. The gradient angles of $I_{\mathrm{WN}}$ should be random but will also have the same degree of oversampling as the Mopra $I$ maps. We calculate $Z_{x\thinspace\mathrm{WN}}^{\prime}$ for 1000 $I_{\mathrm{WN}}$ realizations and list the mean ($\langle Z_{x\thinspace WN}^{\prime}\rangle$) and standard deviation ($\sigma_{Z_{x\thinspace WN}^{\prime}}$) in Table 2. If every $\phi_{\mathrm{WN}}$ measurement was independent then $\sigma_{Z_{x\thinspace WN}^{\prime}}$ should approach 1. The value of $\sigma_{Z_{x\thinspace WN}^{\prime}}$ therefore gives an estimate for the factor by which the data is oversampled. We can therefore estimate the PRS corrected for oversampling by

while the number of independent data samples in the map is

Both quantities are listed in Table 2.

The statistical error bars for $Z_{x}$ listed in Table 2 are always $\simeq$ 1. However, these statistical error bars do not take into account potential systematic effects such as mapping artifacts associated with the Mopra telescope scanning strategy discussed in Section 2.2.

To quantify this we replaced $I$ in Equation 5 with $I$_noise, a “zeroth-moment” map made from velocity channels in the spectral data cube with no apparent molecular emission and recalculated $Z_{x}$. The map gradient angles should be random, and so we would expect these calculated $Z_{x\thinspace\mathrm{noise}}$ values to have a mean of 0 and a standard deviation of 1. The calculated values of $Z_{x}$ listed in Table 2 have a mean of $-$0.35 and a standard deviation of 1.18, which is consistent with our expectations.

In Appendix B we show our results are not sensitive to the resolution of the Mopra zeroth moment maps, the map sampling interval (provided the maps are sampled at least twice per smoothed Mopra beam FWHM $\theta_{\mathrm{sm}}$), or to the method used to remove the contribution of the diffuse ISM to the Vela C polarization maps.

Here we quantify the characteristic density traced by each of our observed molecular lines, in order to understand how the Vela C cloud structure is aligned with respect to the magnetic field over different number density regimes. We do this in three ways, by using a simple non-LTE radiative transfer model to calculate the $n_{\mathrm{H}_{2}}$ needed to reproduce our $I$ observations (Section 4.3.1), by calculating the critical density corrected for radiative trapping for the highly optically thick ¹²CO observations (Section 4.3.2), and by using the cloud width as a proxy for the depth in order to estimate characteristic number density from column density maps (Section 4.3.3).

Unlike the Herschel derived column density maps used in the analysis of Soler et al. (2017) and Jow et al. (2018), the $I$ maps in this work do not necessarily reflect the structure of the total gas column density. Instead the zeroth-moment maps shown in the left panels of Figures 3 and 4 are sensitive to the column density of the emitting molecules, the number density and average speed of particles colliding with the molecules (usually assumed to be H₂), the line optical depth, and the excitation temperature that characterizes the populations of the various rotational energy levels.

The $I$ maps shown in Figures 3 and 4 exhibit noticeable differences in total sky area passing our signal-to-noise threshold requirements (described in Section 3.1). Emission from the lower density tracers ¹²CO and ¹³CO (which show on average a tendency to align parallel to the magnetic field) covers almost the entire map, while C¹⁸O and the intermediate or high density tracers, HCN, HNC, HCO⁺, and CS mostly show emission within the column density contour of $N_{\mathrm{H}}$ $=$ 1.2 $\times$ 10²² cm^-2 (this corresponds to the lowest $N_{\mathrm{H}}$ contour shown in Figures 3 and 4), and the weaker NH₃ and N₂H⁺ lines only show emission towards the highest $N_{\mathrm{H}}$ peaks.

Given the difference in map extent for each of our molecular line $I$ maps it is possible that the change in relative orientation between our molecular tracers is simply showing the same trend of $Z_{x}$ with $N_{\mathrm{H}}$ observed by Soler et al. (2017) and Jow et al. (2018) in Vela C. Soler et al. (2017) found that below $N_{\mathrm{H}}$ $\simeq\thinspace 1.2\thinspace\times\thinspace 10^{22}$ cm^-2 $\mathbf{\langle\hat{B}_{\perp}\rangle}$ is on average parallel to the $N_{\mathrm{H}}$ iso-contours. Since only ¹²CO and ¹³CO have significant emission at $N_{\mathrm{H}}$ $<\thinspace 1.2\thinspace\times\thinspace 10^{22}$ cm^-2 the differences in relative orientation between our observed lines could just be due to the difference in average $N_{\mathrm{H}}$ sampled by each line.

To test this hypothesis, in Figure 9 we recalculate $Z_{x}$ for ¹²CO and ¹³CO only for the sightlines where our intermediate and high density tracers were detected. We see that even when restricting ¹²CO and ¹³CO to the sightlines where higher density tracers are detected the behavior of $Z_{x}$ shows the same trends: structures in the ¹²CO $I$ map align preferentially parallel to $\mathbf{\langle\hat{B}_{\perp}\rangle}$; ¹³CO structures show a weak tendency to align parallel to $\mathbf{\langle\hat{B}_{\perp}\rangle}$, and intermediate to high density tracers show a weak preference to align perpendicular to $\mathbf{\langle\hat{B}_{\perp}\rangle}$. This suggests that the ¹²CO and ¹³CO preferentially trace lower density gas in outer cloud regions compared to the higher density molecular tracers. The only systematic difference in the $Z_{x}$ values for ¹²CO and ¹³CO shown in Figures 6 and 9 is that $Z_{x}$ is lower in Figure 9, which is expected as the intermediate and high density tracers have lower values of $n_{\mathrm{ind}}$ (see Table 2) and Equation 8 shows that $Z_{x}$ is proportional to $\sqrt{n_{\mathrm{ind}}}$.

We can also directly test for changes in $Z_{x}$ with column density by dividing our relative orientation angle $\phi$ data into seven groups binned by $N_{\mathrm{H}}$. The bins are chosen such that for a given molecular line each group has the same number of sightlines. We then calculate $Z_{x}$ for the sightlines in each group. Figure 10 shows the change in relative orientation $Z_{x}$ with increasing $N_{\mathrm{H}}$. Overall this figure gives the same impression as Figure 5, in that there is no consistent trend of relative orientation vs $N_{\mathrm{H}}$. The average $Z_{x}$ decreases with $N_{\mathrm{H}}$ for some tracers (e.g., ¹³CO and NH₃), but increases for other tracers (e.g., HCO⁺ and HNC).

In summary, our results are not consistent with a trend in relative orientation versus hydrogen column density, but suggestive of some relationship to volume density and/or excitation conditions. The magnetic field orientation probed by BLASTPol is always a sum along the line of sight weighted by the dust density, emissivity, and grain alignment efficiency within the volume probed by the telescope beam. For example, if the grain alignment efficiency and temperature were higher in low density cloud regions, the magnetic field orientation measured by BLASTPol could be more sensitive to the field direction in the low density rather than high density cloud regions within the sightline. This averaged $\mathbf{\langle\hat{B}_{\perp}\rangle}$ orientation measurement is what is compared to the orientation of the molecular structures, whether from a low density tracer or a high density tracer, wherever they happen to be along the line of sight. Thus it is important to keep in mind that the preference for intermediate and high density structures to appear aligned perpendicular to the magnetic field measured by BLASTPol does not imply that the magnetic field orientation is that of a field entirely within the volume highlighted by the molecules.

Using the density estimates presented in Sections 4.3.1 to 4.3.3 we can probe the characteristic number density at which the relative orientation of the cloud structure changes with respect to the magnetic field, as traced by the $I$ maps. Figure 11 shows $Z_{x}$ versus $n_{\mathrm{H}_{2}}$ for our two number density estimation techniques. The top panel shows $Z_{x}$ vs. $n_{\mathrm{H}_{2}\thinspace\mathrm{rad}}$, which was derived from the RADEX models (or in the case of ¹²CO $J$ = 1 $\rightarrow$ 0 the critical density corrected for radiative trapping). The bottom panel shows $n_{\mathrm{H}_{2}\thinspace\mathrm{x-sect}}$, where we use the molecular column density cross-sections shown in Figure 8 to estimate the cloud depth and give a rough estimate of the average molecular hydrogen density along the cross-section. In both panels the values for $Z_{x}$ are the same as those discussed in Section 4.2, listed in Table 2, and shown in Figure 6.

Figure 11 shows a transition from a clear detection of preferentially parallel alignment to $\mathbf{\langle\hat{B}_{\perp}\rangle}$($Z_{x}$ $\gg$ 0) for ¹²CO to no preferred orientation or a weakly perpendicular alignment ($Z_{x}$ $<$ 0) for intermediate and high density tracers. As discussed in Section 4.2 while the intermediate and high density tracers with characteristic densities $n_{\mathrm{H}_{2}}\geq\thinspace$10³ cm^-3 tend to individually have low significance values of $Z_{x}$, this is partially explained by the lower number of independent relative orientation angle measurements $n_{\mathrm{ind}}$ compared to ¹²CO and ¹³CO. When we calculate the averaged projected Rayleigh statistic accounting for the correlations in $I$ map structure between  N₂H⁺, HCO⁺, HCN, HNC, CS, and NH₃ we obtain an average $Z_{x}$ $=\thinspace-4.03$, showing that on average intermediate and high density gas structures do preferentially align perpendicular to $\mathbf{\langle\hat{B}_{\perp}\rangle}$.

Our results show that the change in $Z_{x}$ from cloud structures aligned parallel to structures aligned perpendicular to the magnetic field takes place at molecular gas densities between between those traced by ¹³CO and C¹⁸O. For both number density estimation methods, this transition number density $n_{\mathrm{H}_{2}\mathrm{tr}}\sim\thinspace$10³ cm^-3, though with the spread in density estimates the uncertainty in the value of $n_{\mathrm{H}_{2}\mathrm{tr}}$ could be up to a factor 10.

Above $n_{\mathrm{H}_{2}}\sim\thinspace$10³ cm^-3, there are significant inconsistencies between the characteristic density estimated for the same molecules in the two panels of Figure 11. For the molecules N₂H⁺, HCO⁺, HNC, HCN, and CS $n_{\mathrm{H}_{2}\thinspace\mathrm{x-sect}}$ is at least an order of magnitude lower than $n_{\mathrm{H}_{2}\thinspace\mathrm{rad}}$. For HNC and HCN $n_{\mathrm{H}_{2}\thinspace\mathrm{x-sect}}$ is more than a factor of 100 lower than $n_{\mathrm{H}_{2}\thinspace\mathrm{rad}}$. This discrepancy may in part be due to the estimated density being averaged over the width of the cross-section, and also partly because we assume a molecular gas volume filling factor of unity. If the molecular gas filling factor is less than unity then $n_{\mathrm{H}_{2}\thinspace\mathrm{x-sect}}$ will be less than the true characteristic number density probed by the molecular line. In the astrochemical models of a molecular cloud simulation presented in Gaches et al. (2015), the volume filling factor for these molecules ranges from 0.005 (N₂H⁺ $J$ = 1 $\rightarrow$ 0) to 0.40 (HCN $J$ = 1 $\rightarrow$ 0).

Molecular abundance variations with density are not accounted for in either technique for estimating the characteristic density. For example, CO, the primary reservoir of carbon with molecular clouds, is expected to “freeze-out” onto dust grains at intermediate densities. In pre-stellar cores Bacmann et al. (2002) estimate that freeze-out becomes important above $n_{\mathrm{H}}\thinspace\sim\thinspace$10⁴ cm^-3 (corresponding to $n_{\mathrm{H_{2}}}\sim$ 5 $\times$ 10³ cm^-3). Lower levels of carbon in the molecular phase can then reduce the abundance of other carbon-bearing molecules such as CS, HCN, HNC, and HCO⁺ (Bergin & Tafalla, 2007). In contrast nitrogen-bearing molecules such as N₂H⁺ and NH₃ are not expected to freeze-out onto dust grains, and because these molecules tend to be destroyed in interactions with CO and HCO⁺, their abundance can increase towards high densities where CO is depleted (Aikawa et al., 2001; Tafalla et al., 2002; Jørgensen et al., 2004). These abundance variations most likely are not important at our estimated transition density $n_{\mathrm{H}_{2}\mathrm{tr}}\thinspace\sim\thinspace 10^{3}$ cm^-3. However, studies of whether the observed trend of decreasing $Z_{x}$ continues with increasing $n_{\mathrm{H_{2}}}$ continues beyond $n_{\mathrm{H}_{2}\mathrm{tr}}$ will need to consider the possibility of molecular abundance variations with density.

In the PRS analysis presented in this work we have shown for the first time a clear change in the average orientation of gas structures of different characteristic number density with respect to the magnetic field. Previous comparisons with synthetic observations of magnetized cloud formation show that this change of relative orientation has implications for the magnetization of Vela C. This was first shown by Soler et al. (2013), who analyzed three RAMSES-MHD adaptive mesh refinement simulations with self gravity for low, intermediate, and high magnetization cases (specifically, initial thermal to magnetic pressure ratio $\beta$ = ($c_{\mathrm{s}}/v_{\mathrm{A}}$)² = 100.0, 1.0, and 0.1). After beginning the simulation and allowing turbulence to decay they found that only the highest magnetization simulation (initially sub-Alfvénic) showed a change in relative orientation from parallel to perpendicular with increasing density/column density. The intermediate magnetization simulation, where the turbulence was initially close to equipartition with the magnetic field, showed the alignment changing from preferentially parallel at low values of $n$ or $N$, to showing no preferred orientation at high densities.

Our PRS results thus imply that the cloud-scale magnetic field in Vela C is at least trans-Alfvénic in strength, and therefore strong enough to have played an important role in the formation of global cloud structure. This same conclusion was also reached in the studies of Soler et al. (2017) and Jow et al. (2018), which revealed a change in relative orientation of column density iso-contours and magnetic field orientation with increasing column density (see Section 1).

Does the observation of a change in the project Rayleigh statistic $Z_{x}$ for gas tracers of different densities give us any additional information about the cloud magnetic field structure compared to the studies of $Z_{x}$ vs. $N_{\mathrm{H}}$ presented in Soler et al. (2017)? One advantage of studying the change in relative orientation with density rather than column density is that the observed column density distribution will change for different cloud viewing angles. This is shown in Figure 10 of Soler et al. (2013), where different viewing angles resulted in different transition column densities $N_{\mathrm{t}r}$, even though in both cases the magnetic field is parallel to the plane of the sky. Studies of $Z_{x}$ versus $n_{\mathrm{H}_{2}}$ remove this projection effect; however, this method is still sensitive to yet another projection effect, because the polarization data is only sensitive to $\mathbf{\langle\hat{B}_{\perp}\rangle}$, the orientation of the magnetic field projected on the plane of the sky. If the mean direction of the cloud magnetic field is exactly parallel to the line of sight then $\mathbf{\langle\hat{B}_{\perp}\rangle}$ will only measure the disordered components of $\mathbf{B}$ and no average correlation of the $\mathbf{\langle\hat{B}_{\perp}\rangle}$ direction with cloud structure is expected.

Comparisons of the probability distribution functions of the fractional polarization $p$, and the dispersion of polarization angle $S$ on 0.7-pc scales with those from synthetic observations of cloud-forming simulations suggest either that the magnetic field in Vela C is highly turbulent and disordered, or that the mean-field direction is highly inclined with respect to the plane of the sky (King et al., 2018). The first explanation of a disordered (i.e., relatively weak) magnetic field is in conflict with the PRS observations presented in this work and Soler et al. (2017). The latter explanation of a highly inclined magnetic field is therefore more likely and might explain why the $Z_{x}$ versus $N_{\mathrm{H}}$ trend in Vela C appears to be shallower than the same curves for many of the clouds discussed in Planck Collaboration Int. XXXV (2016). However, we note that the simulations considered in King et al. (2018) are highly idealized and did not cover a wide range of cloud physical parameters. A more comprehensive parameter study is being conducted and will be published in a separate paper.