BALLOON-BORNE SUBMILLIMETER POLARIMETRY OF THE VELA C MOLECULAR CLOUD: SYSTEMATIC DEPENDENCE OF POLARIZATION FRACTION ON COLUMN DENSITY AND LOCAL POLARIZATION-ANGLE DISPERSION

Standard techniques were applied to the bolometer time-ordered data (hereafter TOD) to remove detector spikes (mostly due to cosmic rays), deconvolve the bolometer time constant, and remove an elevation-dependent feature (see Matthews et al. 2014). The data were further preprocessed by fitting and removing an exponential function fit to each detector's TOD in the first 30 seconds after a HWP rotation or a telescope slew. A high-pass filter with power-law cutoff was used to whiten noise in the TOD below 5 mHz. Temporal gain variations were removed using the DC voltage level of each detector and periodic measurements of an internal calibration source. Pixel-to-pixel detector gain variations were corrected by frequent observations of the bright compact source IRAS 08470−4243.

Telescope attitude was reconstructed using pointing solutions generated from the BLASTPol optical star camera,²⁵ with payload rotational velocities from gyroscopes used to interpolate between pointing solutions (Pascale et al. 2008). Data having pointing uncertainties >5'' were discarded. The final on-sky pointing solution was calibrated to match the astrometry of publicly available Herschel SPIRE maps²⁶ at the same wavelength.

The BLASTPol 2012 beam differs from the beam predicted by our optics model. It has multiple elongated peaks, and the relative power in each peak varies from detector to detector. BLASTPol filters were designed for near-space conditions and the telescope far field is several kilometers away so it was not possible to map the far-field beam at sea level (Galitzki et al. 2014). Instead the beam had to be inferred from in-flight measurements of astronomical objects.

Our 2012 instrument beam model was defined in telescope coordinates and was informed by observations of two objects: IRAS 08470−4243, a warm compact dust source located in the Vela Molecular Ridge; and limited observations of the planet Saturn made on 27 December 2015. IRAS 08470−4243 was observed every 4–8 hr, with reasonable coverage for all detectors, but it is not a point source at BLASTPol resolution. BLAST observations of IRAS 08470−4243 in 2006 found a FWHM of ∼40'' (Netterfield et al. 2009). Saturn has a radius of 6 × 10⁴ km, which corresponds to an angular size of <20'', considerably smaller than the BLASTPol 2012 beam. Saturn was only observed early in the flight at telescope elevations of <25° and was only fully mapped by the bolometers near the center of the focal plane arrays. Three elliptical Gaussians were fit to the Stokes I maps of IRAS 08470−4243 and Saturn. The free parameters were the Gaussian amplitudes, centroids, widths, and position angles. Only pixels above an intensity threshold of 20%  of the peak intensity for IRAS 08470−4243 and 7.5% of the peak intensity for Saturn were used in the fits. The final 2012 instrument beam model used the centroid positions, amplitudes, and position angles from the fits to IRAS 08470−4243 and the Gaussian widths calculated from the fits to Saturn.

Next, an on-sky beam model was computed for the Vela C observations. The on-sky beam model is a time-weighted average of the instrument beam model rotated by the angle between the telescope vertical direction and Galactic north for each raster scan. The resulting average beam model for Vela C is shown in the left panel of Figure 2. A single elliptical Gaussian fit to this beam model gives FWHMs of 130'' by 64'', at position angle −51°. This beam is significantly larger than the expected diffraction limit of the telescope (FWHM = 60'').

Zoom In Zoom Out Reset image size

Lucy–Richardson (LR) deconvolution was previously used to correct a larger-than-expected beam from the BLAST 2005 flight (Roy et al. 2011). Here we used an iterative LR method and our beam model to deconvolve a simulated map consisting of a single Gaussian source with FWHM = 145''. The deconvolved map from this step when applied as a smoothing kernel to convolve the original beam model should restore the 145'' Gaussian. The success of this step can be judged from the right panel in Figure 2: it does produce a single-peaked source that is approximately round (FWHM = 144'' by 157''). This same smoothing kernel was used to convolve the I, Q, and U maps of Vela C, with a resulting resolution of approximately 25.

To determine the instrumental polarization (the polarization signal introduced by the instrument hereafter referred to as IP) we followed methods described in Matthews et al. (2014). In brief, the set of observations of Vela C was split into two bins based on parallactic angle, and maps were produced for each detector individually using the naivepol mapmaker (Moncelsi et al. 2014). The measured polarization is a superposition of one component fixed with respect to the sky and the IP, which is fixed with respect to the telescope. By comparing the polarization measurements at different parallactic angles the IP of each bolometer could be reconstructed. These IP terms were then removed during the mapmaking stage (Section 3.4). The 500 μm array has an average IP amplitude of 0.53%. To check the effectiveness of the IP correction the Vela C data were divided into two halves and IP estimates derived from the first half of the data were used to correct the second half of the data. By measuring the IP of the "corrected" second half of the Vela C data we estimate that the minimum value of the fractional polarization p measurable by BLASTPol at 500 μm is 0.1%.

Maps were made using TOAST (Time Ordered Astrophysics Scalable Tools),²⁷ a collection of serial and OpenMP/MPI parallel tools for simulation and map making. Specifically, the generalized least-squares solver was used, which iteratively inverts the map-maker equation using the preconditioned conjugate gradient method (see Cantalupo et al. 2010, a predecessor to TOAST). The map-maker's noise model was estimated using power spectra from observations of faint dust emission in the constellation Puppis (map centered at l = 2390, b = −17), with simulated astrophysical signal subtracted. The noise model is consistent with white noise plus ${(1/f)}^{\alpha }$ correlations that level out at low frequency due to data preprocessing. Correlations between detectors and non-stationarity of the noise were not required by the model. Instrumental polarization (Section 3.3) was removed as per Matthews et al. (2014). Per-pixel covariance matrices for Stokes $I,\;Q$, and U were estimated as the 3 × 3 diagonal block of the full pixel-pixel covariance matrix. Noise-only maps, both simulations and null tests (see Section 3.6), are consistent with the estimated covariances, up to a constant scaling factor due to unmodeled noise. A pixel size of 10'' was used for all signal and covariance maps.

To study the polarization properties and magnetic field morphology of Vela C it is necessary to isolate polarized dust emission originating in Vela C from the diffuse polarized emission associated with Galactic foreground and background dust as well as the Vela Molecular Ridge.²⁸ This separation requires extra care as previous studies show that diffuse sightlines, which may be used to estimate the foreground/background polarized emission, tend to have a higher average polarization fraction than dense cloud sightlines. In particular, Planck observations show that in the more diffuse clouds there is a range of p values with the maximum reaching 15%–20%, while such high values are not seen in the higher column density clouds (e.g., Orion, Ophiuchus, Taurus), where the average p is consistently lower (Planck Collaboration Int. XX 2015). Polarized emission from diffuse dust along dense cloud sightlines could therefore contribute significantly to the overall polarization measured. In this section we present two different diffuse emission subtraction methods, one conservative and one more aggressive with respect to diffuse emission removal.

In the conservative method for diffuse emission subtraction, we considered most of the emission surrounding Vela C as defined by Hill et al. (2011) to be associated with the cloud. The Hill et al. (2011) Vela C cloud subregions are overplotted (solid white lines) on a map of 500 μm total intensity in Figure 1. The zero-point for the intensity of the cloud emission was set by specifying a region with relatively low intensity near Vela C and assuming that emission in this reference region also contributes to sightlines on the cloud with spatial uniformity. This low flux region is labeled "C" in Figure 1. We calculated the average Stokes I, Q, and U in that region, and the appropriate mean flux was then subtracted from each of the maps. The result was a set of maps of Vela C emission isolated from the Galaxy, assuming a minimal, uniform contribution from foreground and background emission.

In the aggressive method we considered most of the diffuse emission surrounding Vela C to be unassociated with the cloud, and accordingly a higher level of flux needed to be removed to isolate the cloud. Furthermore, we noted that there is significantly more emission to the south of Vela C than to the north; thus it is reasonable to assume that the true map of the region consists of the Vela C cloud superimposed on a varying Galactic emission profile. In this method of referencing, we defined two regions closely surrounding the cloud (marked "A1" and "A2" in Figure 1) and performed two-dimensional linear-plane fits to the Stokes I, Q, and U maps excluding all map pixels except those located in regions A1 and A2. The three free parameters in these fits were the linear slopes of the plane in the directions tangent to l and b and a map offset. The equations for each of the resulting plane-fits to the I, Q, and U maps were used to specify the intensity to be subtracted from each pixel in the maps. Note that for regions far from Vela C, the linear approximation of the Galactic emission profile breaks down, leading to an inappropriate extrapolation. Therefore we defined an area within which the linear fit referencing method is valid (blue quadrilateral in Figure 1), bounded on the north and south by the reference regions A1 and A2, and on the east and west by the edges of the well-sampled portion of the map. This "validity" area roughly coincides with the four southernmost regions of Hill et al. (2011). We note that some of the emission in A1 and A2 might in fact be associated with Vela C, so this method is likely to over-subtract the diffuse dust emission.

The true I, Q, and U maps of Vela C probably exist somewhere between our most extreme physically reasonable assumptions corresponding to the conservative and aggressive diffuse emission subtraction methods. In this paper we present results for an "intermediate" diffuse emission subtraction method, derived by taking the arithmetic mean of the I, Q, and U maps from the aggressive and conservative methods. Most of our analyses are then repeated using the aggressive and conservative methods as a gauge of the uncertainties associated with the diffuse emission subtraction.

To characterize possible systematic errors in our data, we performed a series of null tests, which are described in detail in Appendix A. In these, we split the 250 μm observations²⁹ into two mutually exclusive sets. If the polarization parameters from the two independent data sets agree, we can conclude that the impact of systematics is small, and the uncertainties are properly characterized by Gaussian errors produced by TOAST. The four methods of splitting are: data from the left half of the array versus the right half; data from the top half of the array versus the bottom half; data from earlier in the flight versus later in the flight; and alternating every other scan set sequentially throughout the flight.

For each null test we made separate maps of I, Q, and U, which were then used to calculate residual maps of the polarized intensity P, the polarization fraction p, and the polarization-angle ψ as described in Appendix A. (The quantities P, p, and ψ are defined in Appendix B). If our data had no systematic errors we would expect to see uncorrelated noise in the residual maps. For P and p if the residuals were less than one-third of the signal in a given map pixel then that pixel was said to pass the null test. For ψ the residuals had to be less than 10° to pass the null test. We examined each of the four null tests listed above for each of the two methods of diffuse emission subtraction (Section 3.5) in P, p and ψ giving a total of 24 checks that our measured polarization signal is significantly above the systematic uncertainty level.

We found that our map passed these tests for the majority of sightlines inside the "cloud" region shown in Figure 1 (blue solid line). The exceptions occurred in regions where the fractional polarization was small, so that a comparison of the scale of the polarization signal to the scale set by residuals in the null tests resulted in an apparent failure. The fact that we saw null test failures correlating with low p, but not with absolute difference in the null test I maps led us to the interpretation that the apparent low signal level compared to the null test residual is due to decreased signal and not increased systematic uncertainties. We did see significant structure in the null test residual maps of Q and U near the compact H ii region RCW 36, which coincided with null test residuals of one-fourth p, though the residuals in ψ were much smaller than 10°. These p measurements technically pass the null test criteria, but the systematic errors are larger than the statistical errors. We conclude that for the validity region shown in Figure 1 the null tests are passed, with the exception of a circular area centered on RCW 36 (l = 26515, b = 142 within a radius of 4').

Figure 3 shows Vela C 500 μm maps for the three Stokes parameters I, Q, and U. The maps have been smoothed to 25 resolution, as described in Section 3.2, and use the intermediate diffuse emission subtraction method (Section 3.5). Overlaid in gray are the outlines of the subregions of Vela C as defined in Hill et al. (2011) and labelled in Figure 8. The BLASTPol I map peaks at the location of RCW 36.

Zoom In Zoom Out Reset image size

Also included in Figure 3 is the derived map of the polarized intensity (P), which generally shows some signal where there is cloud emission. However, the correspondence is certainly not perfect, and varies considerably across the map. For example, along most of the Center-Ridge there is a corresponding peak in the P map along the main ridge. In the South-Ridge there are peaks at similar locations in the P and I maps. But in the South-Nest, prominent areas of polarized emission are only seen around the edge of the cloud structure seen in I. There are also some regions of significant P that stand out less in I, for example, along the north edge of the Center-Ridge.

Figure 4 shows the polarization fraction ($p=P/I$) for each of the three different diffuse emission subtraction choices discussed in Section 3.5. The conservative diffuse emission subtraction (top panel) results in p that is lower on average than from the aggressive diffuse emission subtraction (middle panel). This is expected as p is commonly observed to increase for regions of low dust emission. Thus, compared to the aggressive subtraction method that uses regions closer to the cloud with lower average p, the conservative method removes more P relative to I. The bottom panel shows the p map resulting from the intermediate diffuse emission subtraction method. Unless otherwise specified, for the remainder of the paper p, ψ, and Φ are calculated from Stokes parameter maps using the intermediate diffuse emission subtraction method.

Zoom In Zoom Out Reset image size

The mean value of p in our map is 6.0% with a median of 3.4%. For the map pixels within the dense cloud subregions defined by Hill et al. (2011) the mean polarization fraction is 3.5% with a median of 3.0% and a standard deviation of 2.4%. Previous submillimeter polarization maps having spatial coverage corresponding to the scales of entire clouds have yielded roughly similar values. Specifically, after subtracting the background/foreground emission, Planck Collaboration Int. XXXIII (2016) found mean 850 μm polarization fractions of 1.8%, 5.0%, and 6.1% for three nearby molecular clouds, while 450 μm polarization maps of four GMCs made by Li et al. (2006) with SPARO at the South Pole yield a mean polarization fraction of 2.0%. Our p map shows behavior that is broadly consistent with expectations from the P map. Values of p tend to decrease with increasing I, but there is not a one-to-one anticorrelation between p and I.

Figure 5 shows a detailed view of the magnetic field orientation projected onto the plane of the sky Φ, as inferred from the BLASTPol 500 μm data. This figure uses a "drapery" pattern produced using the line integral convolution method of Cabral & Leedom (1993) superimposed on the BLASTPol 500 μm I map.³⁰ Dotson (1996) showed that there is significant ambiguity in inferring the magnetic field lines from polarization data, particularly as polarization maps can sample multiple cloud structures along the line of sight, each potentially having a different magnetic field orientation. The drapery image is presented solely to show the range of orientations of Φ, and to give a sense of the range of spatial scales probed by BLASTPol. Figure 6 shows Φ as a series of line segments (approximately one line segment per 25 BLASTPol beam).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The projected cloud magnetic field direction appears to change across Vela C: at low Galactic latitudes the field is mostly perpendicular to the main cloud elongation direction, while at higher Galactic latitudes it bends to run mostly parallel to the cloud elongation direction. We also see some sharp changes in Φ, most noticeably in the South-Nest, and near the compact H ii region RCW 36.

Figure 7 shows that the dispersion in magnetic field orientation across the Vela C cloud is 28°. Novak et al. (2009) calculated dispersions on similar spatial scales by combining the large-scale GMC polarization maps of Li et al. (2006) with higher angular resolution submillimeter polarimetry data. They obtained 27°–28°, nearly the same result. In future publications we will present statistical studies of the correlations between magnetic field orientation, filamentary structure, and cloud velocity structure.

Zoom In Zoom Out Reset image size

To quantify the disorder of Φ in our Vela C maps at small scales we calculate the polarization-angle dispersion function S, implementing the formalism described in Section 3.3 of Planck Collaboration Int. XIX (2015). For each pixel in our map, S is defined as the rms deviation of the polarization-angle $\psi ({\boldsymbol{x}})$ for a series of points on an annulus of radius δ:

$$\begin{eqnarray}&&{S}^{2}({\boldsymbol{x}},\delta )=\displaystyle \frac{1}{N}\displaystyle \sum _{i=0}^{N}{S}_{{xi}}^{2},\end{eqnarray} \tag{ 3 }$$

where δ is the length scale of the dispersion, ${\boldsymbol{x}}$ is the position for which we evaluate the polarization-angle dispersion and

$$\begin{eqnarray}&&{S}_{{xi}}=\psi ({\boldsymbol{x}})-\psi ({\boldsymbol{x}}+{{\boldsymbol{\delta }}}_{i}).\end{eqnarray} \tag{ 4 }$$

Because S is always positive it is biased due to noise. We debias using the standard formula

$$\begin{eqnarray}&&{S}_{{db}}^{2}(\delta )={S}^{2}(\delta )\ -\ {\sigma }_{S}^{2},\end{eqnarray} \tag{ 5 }$$

where ${\sigma }_{S}^{2}$ is the variance of S.

Figure 6 shows S for δ = 25 (∼0.5 pc), the smallest scale that can be resolved with our smoothed beam. (Hereafter we refer to $S(2\buildrel{\,\prime}\over{.} 5)$ as S). The most striking features in the S map correspond to regions of sharp changes in Φ, which is indicated with line segments. These high dispersion regions sometimes occur near the locations of dense filaments (for example, the sharp bend in the South-Nest). More often they correspond to sightlines of lower than average p and do not appear to be coincident with any prominent cloud feature in I.

Our BLASTPol polarization maps were calculated from Stokes parameters smoothed to a resolution of ∼25 as discussed in Section 3.2. The resulting polarization maps were then sampled every 70'' to ensure at least Nyquist sampling. In total there are 4708 projected magnetic field sightlines over the validity region defined in Section 3.5.

In the following sections we attempt to model the polarization fraction p as a function of N, T, and S. For these detailed studies we restrict our analysis to sightlines that encompass the dense cloud regions as defined by Hill et al. (2011). These sightlines are better probes of the polarization structure in the cloud and are less sensitive to systematic uncertainties in our ability to separate the polarized emission emitted by diffuse dust foregrounds/backgrounds from the polarized emission emitted by dust grains in Vela C.

To ensure a robust sample, we use only p values that are large enough to be unaffected by uncertainties in instrumental polarization removal (p > 0.1%, see Section 3.3), and for which we have at least a 3σ detection of polarization ($p\;\gt \;3{\sigma }_{p}$) , which corresponds to an uncertainty in the polarization angle ${\sigma }_{\psi }\;\lt \;10^\circ$. To ensure that the polarization values are not dependent on our choice of diffuse emission subtraction method we require that ${p}_{{\rm{int}}}\;\gt \;3| {p}_{{\rm{int}}}-{p}_{{\rm{con}}}|$ and ${p}_{{\rm{int}}}\;\gt \;3| {p}_{{\rm{int}}}-{p}_{{\rm{agg}}}|$, where ${p}_{{\rm{con}}}$, ${p}_{{\rm{agg}}}$, and ${p}_{{\rm{int}}}$ are the polarization fraction values calculated using the conservative, aggressive, and intermediate diffuse emission subtraction methods respectively (see Section 3.5). Similarly, we require that $| {\psi }_{{\rm{int}}}-{\psi }_{{\rm{con}}}| \;\lt \;10^\circ$ and $| {\psi }_{{\rm{int}}}-{\psi }_{{\rm{agg}}}| \;\lt \;10^\circ$. We also exclude sightlines from a 4' radius region near RCW 36 as these show residual structure in our null tests (see Section 3.6). In total 2488 out of a 3056 possible Hill et al. (2011) sightlines meet these criteria. For our analysis of p versus N, T, and S in Sections 6 and 7 we also require at least 3-σ measurements of N, T, and S where the errors on N and T are derived from the SED fit covariance matrices (see Section 5). This results in a final sample of 2378 sightlines.

Figure 13 shows the median p (color map) for bins of S and N for ISRF-heated sightlines. There is a clear decrease of p with increasing N and S. Individual correlations are shown in Figures 12 and 14, and the derived associated power-law exponents are listed in Table 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Power-law Exponents of p vs. N and S

Diffuse Emission	${\alpha }_{N}$	${\alpha }_{S}$	${\alpha }_{{p}^{[S]}N}$	${\alpha }_{{p}^{[N]}S}$
Subtraction Method
Intermediate	−0.58 ± 0.02	−0.67 ± 0.02	−0.46 ± 0.01	−0.58 ± 0.01
Conservative	−0.53 ± 0.02	−0.67 ± 0.02	−0.39 ± 0.01	−0.56 ± 0.01
Aggressive	−0.66 ± 0.02	−0.66 ± 0.02	−0.57 ± 0.01	−0.59 ± 0.01

Note. The power-law exponents listed in this table are derived from linear fits of $\mathrm{log}p$ to $\mathrm{log}N$ and $\mathrm{log}S$ as described in Sections 6 and 7.1.

Download table as:  ASCIITypeset image

Decreasing p with increasing N has been seen in submillimeter polarization observations of many clouds and cores (e.g., Matthews et al. 2001; Li et al. 2006). The observed decrease in p is often attributed to either cancelation of polarization signal for high-N sightlines due to more disorder in the magnetic field, or to changes in grain alignment efficiency within the cloud. These possible explanations are discussed further in Section 8.

In Section 4.3, we showed that Vela C has high values of S in localized filament-like regions, where there are sharp changes in magnetic field direction. High S depends implicitly on spatial changes in the magnetic field locally in the map, and any related changes in the magnetic field direction within the volume sampled by the BLASTPol beam could lead to lower p. The top panel of Figure 14 shows p versus S. There is a clear anticorrelation between p and S (${\alpha }_{S}\;=$ −0.67 ± 0.02, the coefficient of determination R² = 0.47). We see no dependence of S on N (Figure 14, lower panel).

We showed in Section 6 that the dependence of p on N can be removed using Equation (8) to create p^[N]. Similarly we can normalize out the dependence of p on S by calculating

$$\begin{eqnarray}&&{p}_{i}^{[S]}={p}_{i}{\left(\displaystyle \frac{\bar{S}}{{S}_{i}}\right)}^{{\alpha }_{S}}.\end{eqnarray} \tag{ 9 }$$

The top panel of Figure 15 shows that by removing the dependence of S from p the degree of correlation of ${p}^{[S]}$ with N increases (R² = 0.35 compared to 0.30). Similarly, the bottom panel shows that the correlation of ${p}^{[N]}$ with S is better than the correlation of p with S (R² = 0.50 compared to 0.47). This indicates that both N and S contribute independently to the structure seen in p. The fitted power-law exponents tend to be systematically shallower for the decorrelated p^[S] and p^[N] than for trends with p (see the first row of Table 1), which might imply a weak underlying correlation between N and S (see Figure 14, bottom panel).

Zoom In Zoom Out Reset image size

As noted in Section 7.1, Figure 13 shows a color map of the median p binned two-dimensionally in S and N for ISRF-heated sightlines. The clear decrease of p with both increasing N and S is suggestive of a joint power-law relationship. Here we derive a function $p(N,S)$ that accounts for most of the structure seen in the p map. Specifically, we adopt the joint power-law form

$$\begin{eqnarray}&&\mathrm{log}p(N,S)={K}\;+\;{\alpha }_{N}\mathrm{log}N\;+\;{\alpha }_{S}\mathrm{log}S,\end{eqnarray} \tag{ 10 }$$

where ${K}$, α_N and α_S are the free parameters.

The exponents derived via a fit to Equation (10) are α_N = −0.45 ± 0.01 and  α_S = −0.60 ± 0.01. Just as in Sections 6 and 7.1, errors in fit parameters are derived via bootstrapping (Table 2). We note that, as expected, α_N and α_S derived from the two-variable power-law fit to N and S are identical within the error bars to ${\alpha }_{{p}^{[S]}N}$ (the power-law fit to ${p}^{[S]}$ as a function of N) and ${\alpha }_{{p}^{[N]}S}$ (the power-law fit to ${p}^{[N]}$ as a function of S), which were derived in Section 7.1 (also see Table 1).

Table 2.  Fit Parameters of $p(N,S)$ from Equation (10)

Diffuse Emission	${\alpha }_{N}$	${\alpha }_{S}$	${\boldsymbol{K}}$
Subtraction Method
Intermediate	−0.45 ± 0.01	−0.60 ± 0.01	8.42 ± 0.3
Conservative	−0.41 ± 0.01	−0.59 ± 0.01	6.92 ± 0.3
Aggressive	−0.58 ± 0.01	−0.60 ± 0.01	10.98 ± 0.3

Note. The power-law exponents (${\alpha }_{N}$ and ${\alpha }_{S}$) and fitted constant (K) listed in this table are calculated from a two-variable power-law fit to N and S as described in Section 7.2.

Download table as:  ASCIITypeset image

We can remove the dependence of p on N and S via

$$\begin{eqnarray}&&{p}_{i}^{[N,S]}=\displaystyle \frac{{p}_{i}\;\bar{p}}{p({N}_{i},{S}_{i})},\end{eqnarray} \tag{ 11 }$$

where ${p}_{i}^{[N,S]}$ is the decorrelated p for the ith data point, p_i is the measured polarization fraction for the ith data point, $p({N}_{i},{S}_{i})$ is the value of the two-variable power-law fit for the ith data point, and $\bar{p}=0.029$ is the median value of p. A comparison of the spatial distribution of p (middle panel) with ${p}^{[N,S]}$ (bottom panel) is shown in Figure 11. We discuss potential causes of residual structure in the ${p}^{[N,S]}$ map in Section 8.4.

Finally, we quantify the degree to which our two-variable power-law fit $p(N,S)$ can reproduce the observed dispersion in p. Figure 16 shows histograms of: our calculated $\mathrm{log}p$ (left panel); $\mathrm{log}p(N,S)$, the two-variable power-law fit values calculated for our N and S data points (center panel); and $\mathrm{log}{p}^{[N,S]}$, p with the derived dependence on N and S removed using Equation (11) (right panel). For each of the three cases the median p is 0.029. Histograms of $\mathrm{log}p$ rather than p are shown, because the fits are made in log space. The variances of $\mathrm{log}p$, $\mathrm{log}p(N,S)$, and $\mathrm{log}{p}^{[N,S]}$ are 0.068, 0.045, and 0.023, respectively. Our model $\mathrm{log}p(N,S)$ reproduces 66% of the variance in the $\mathrm{log}p$ map, which shows that our two-variable power-law fit model captures most of the physical effects that determine variations in fractional polarization in Vela C.

Zoom In Zoom Out Reset image size

The uncertainties of the exponents for the two-variable power-law fits ${\alpha }_{N}$ and ${\alpha }_{S}$ were estimated using the bootstrapping methods described in Section 6. However, as discussed in Section 3.5, the limiting uncertainty is our inability to precisely separate the contribution of background/foreground dust from the polarized emission of Vela C. We repeated our analysis for maps of p and S made with the "conservative" and "aggressive" diffuse emission subtraction methods, to gauge the systematic uncertainty of our derived power-law exponents. Table 1 gives power-law exponents derived from the individual correlations (Section 7.1) for the three different diffuse dust emission subtraction methods. Table 2 lists the exponents ${\alpha }_{N}$ and ${\alpha }_{S}$ of the two-variable power-law fits again for all three diffuse emission subtraction methods. Systematic uncertainties relating to the choice of subtraction method are seen to be ∼0.1 for ${\alpha }_{N}$ and ∼0.01 for ${\alpha }_{S}$.

In Section 6 we examined the dependence of the polarization fraction p on column density N and dust temperature T in Vela C. We divided our Vela C sightlines into two groups: those that show evidence of heating from the compact H ii region RCW 36, and those sightlines where the temperature decreases as ${e}^{-0.28N}$. For the latter sightlines we suggested that the dust temperature is primarily set by exposure to the interstellar radiation field (ISRF), with high N sightlines having on average more shielding and therefore receiving less heating per unit mass from the ISRF.

For the ISRF-heated sightlines, we find that p decreases with increasing N and also that p increases as T increases. Depolarization for higher column density sightlines has been seen in many studies (see Section 7.1). Vaillancourt & Matthews (2012) used the ratio of F(850 μm)/F(350 μm) as a proxy for dust temperature in two massive star forming clouds. They found that the polarization tended to decrease with increasing F(850 μm)/F(350 μm), implying that warmer dust grain populations tend to have a higher p. This agrees with our result, but Vaillancourt & Matthews (2012) caution that variations in F(850 μm)/F(350 μm) could be due to changes in dust spectral index, rather than just dust temperature. Our study is the first to fit p measurements within a molecular cloud as a function of both N and T. For the ISRF-heated sightlines, which are the majority of the sightlines, we find that the dependence of p on N is not separate from the dependence of p on T, since N and T are highly correlated.

There are two general classes of explanations for our observations of p versus N and T for the ISRF-headed sightlines. We may have greater magnetic field disorder along high N (and therefore lower T) dust columns, or we may have a decrease in the intrinsic polarization efficiency (see Section 1) for such sightlines. In the first explanation the increased field disorder could arise because of a higher field disorder at high particle densities n, or because high N sightlines pass through more cloud material and therefore may sample different field directions at different locations along the line of sight (Jones 1989). Regarding the second possibile explanation for our observed p versus N and T trends, note that in the RATs model of grain alignment, "alignment torques" from an anisotropic radiation field are responsible for aligning the dust grain spin-axes with the local magnetic field (Lazarian & Hoang 2007). Grains near the surfaces of molecular clouds (low N, high T) thus might be expected to show a higher average polarization fraction (Cho & Lazarian 2005). Alternatively, dust grain properties could change at high densities, e.g., grains could become rounder due to accretion of icy mantles (Whittet et al. 2008).

For our RCW 36 heated sightlines there is a significant anticorrelation between p and N (R² = 0.45). However, for these heated sightlines there is no correlation between N and T and no strong correlation between p and T (${R}^{2}=0.03$). This could indicate that the primary dependence of p is on N, rather than T  and that the correlation of p with T only appears when there is a strong correlation of T with N. However, we caution that we have relatively few sightlines near RCW 36 (143 Nyquist-sampled sightlines compared to 2235 ISRF-heated sightlines) so the lack of correlation between p and T could be caused by the angle of the magnetic field changing with respect to the line of sight, which would cause more spread in p. Indeed Figure 5 shows that near RCW 36 there are significant changes in the inferred magnetic field orientation projected onto the plane of the sky.

In Section 7.2 we fit a model that describes p as a function with a power-law dependence on two variables, hydrogen column density N, and S, the dispersion in the polarization-angle on scales corresponding to our beam FWHM (25, or 0.5 pc). The derived power-law exponents are α_N = −0.45 ± 0.10 for the dependence on N, and α_S = −0.60 ± 0.01 for the dependence on S (see Section 7.3). Our $p(N,S)$ fit is able to reproduce most of the structure seen in our $\mathrm{log}p$ maps.

The decrease in p with increasing S can be attributed to changes in the magnetic field direction within the volume probed by the BLASTPol beam. The mean magnetic field orientation, Φ, is an average over both the beam area (0.5 pc) and along the length of the cloud in the line of sight direction, and is weighted by the density and intrinsic polarization efficiency (Section 1). Large values of S indicate a substantial change in Φ on the scale of a beam, which implies a significant change in the orientation of polarization within the beam. This could be due to a sharp change in the magnetic field direction at some location within the cloud. Alternatively, it could indicate the overlap of two clouds, well separated along the line of sight, each with a different Φ. In either case we should see an overall decrease in the measured polarization fraction, since some of the polarization components cancel. Planck Collaboration Int. XX (2015) note a decrease of p with increasing S, both in their data and in corresponding MHD simulations (see Figure 19 of Planck Collaboration Int. XX 2015). However the Planck study sampled 5 × 10²⁰ cm⁻² < N < 10²² cm⁻² while our Vela C observations predominantly sample 10²² cm⁻² < N < 10²³ cm⁻². Also direct comparison with their derived power-law exponents is difficult, since they fit S versus p, thus minimizing the scatter in S, while we fit p versus S, which minimizes scatter in p. Nevertheless they do find a significant anti-correlation of p and S in their data that is reproduced in their MHD simulations. In these simulations there is by contrast only a weak correlation of S with N (F. Levrier 2016, private communication), just as we found in our data (Figure 14, lower panel).

In Section 8.1, we discussed two classes of explanations for the observed p versus N trend. The first class involves magnetic field disorder. An example is the work of Falceta-Gonçalves et al. (2008). These authors were able to reproduce a decrease in p with increasing N via synthetic polarization maps made from supersonic, sub-Alfvénic MHD molecular cloud simulations, assuming uniform intrinsic polarization efficiency. Their power-law exponents ${\alpha }_{N}$ ranged from 0 to −0.5, with models where the mean magnetic field was in the plane of the sky (γ = 0°) having the steepest slope and models where the mean field was parallel to the line of sight (γ = 90°) having no dependence of p on N. In this theoretical study, the decrease in polarization for higher column density regions is due to an increase in the dispersion of the magnetic field direction for high density regions. An analytic model by Jones (1989), is similarly able to reproduce a falling p versus N for a medium having uniform intrinsic polarization efficiency.

Our analysis shows only a weak correlation (or perhaps no correlation) between S and N (see Section 7). Thus on 0.5 pc scales, we find no significant increase in the dispersion of Φ for sightlines of increasing column density. Such an increase might be expected if disorder in the magnetic field direction increased in high density regions (for example due to accretion-driven turbulence as in Hennebelle and André 2013), or if the magnetic field were affected by large-scale gas motions near self gravitating filaments. In the above-mentioned theoretical study by Falceta-Gonçalves et al. (2008), the authors showed that rarefied cloud regions show little variation in polarization direction whereas significant fluctuations in direction do occur within dense condensations. In this case, one might expect a positive correlation between S and N, which we do not see in our observations.

The second class of explanations for the observed decrease in p with increasing N involves changes in intrinsic polarization efficiency. This idea derives support from the observations of Whittet et al. (2008). These authors measured the near-IR polarization of background stars in four nearby molecular clouds. For studies of polarization of background starlight the quantity that is analogous to fractional polarization of dust emission is referred to as the "polarization efficiency," defined as the fractional polarization of the starlight divided by the extinction optical depth at the same wavelength ${p}_{\lambda }/{\tau }_{\lambda }$. Whittet et al. (2008) found that the polarization efficiency in their clouds was consistent with a power-law dependence, ${p}_{\lambda }/{\tau }_{\lambda }\;\propto \;{A}_{V}^{-0.52}$. Because the inferred magnetic field direction is mostly uniform across the region studied, they attributed all of the decrease in polarization efficiency with increasing N to changes in the intrinsic polarization efficiency. It is interesting to note that our power-law exponent (α_N = −0.45 ± 0.10) is similar to that found by Whittet et al. (2008). Other starlight polarization studies have found power-law exponent values ranging from −0.34 to −1.0 (Gerakines et al. 1995; Goodman et al. 1995; Arce et al. 1998; Chapman et al. 2011; Alves et al. 2014; Cashman & Clemens 2014; Jones et al. 2015). Ground-based studies of polarized thermal dust emission yield similar results. For example Matthews et al. (2002) examined p versus I₈₅₀ for three clouds in Orion B south and found $p\;\propto {({I}_{850\mu {\rm{m}}})}^{\alpha }$, with α ranging from −0.58 to −0.95.

Which of the two general classes of explanations for the observed p versus N trend best explains our observations of Vela C? Naively, the absence of a correlation between S and N would suggest that magnetic field disorder does not increase toward high N sightlines, which would imply that variation in intrinsic polarization efficiency is the more likely explanation. However, if the increased disorder in the field occurs on much smaller scales than 0.5 pc, the scale probed by S, then S is not sensitive to the random component of the field and so we would not expect a correlation between N and S. We emphasize that detailed statistical comparisons with simulations of magnetized clouds that include variations in intrinsic polarization efficiency are needed to fully understand the origin of the p versus N anticorrelation (e.g., Soler et al. 2013). Such comparisons are beyond the scope of the present paper.

In Section 8.2 we advanced various explanations for the anticorrelation between p and N in Vela C. Here we consider an extreme case where all of the dependence of p on N is due to reduced intrinsic polarization efficiency in shielded regions. Our goal is to quantify the ability of our measurements to trace magnetic fields deep inside the Vela C cloud under this pessimistic assumption. If most of the polarized emission comes from the outer diffuse layers of the cloud then our derived magnetic field orientations will not be sensitive to changes in the magnetic field direction within dense structures embedded deep in the cloud.

We model the efficiency of the dust along a given cloud sightline in emitting polarized radiation with , where epsilon is normalized such that $\epsilon =\xi \;{\mathrm{cos}}^{2}(\gamma )I/{A}_{{\rm{V}}}$, where ξ is the intrinsic polarization intensity as defined in Section 1, ${A}_{{\rm{V}}}$ is the total dust extinction in the V band for that sightline, and γ is the angle of the magnetic field with respect to the plane of the sky. For these models $\epsilon =\epsilon (\chi )$, where χ is the parameterized depth into the cloud, which is equal to the integrated visual extinction to the nearest cloud surface as indicated in Figure 17. Note that these models make a number of assumptions: (a) that the cloud is isothermal; (b) that the dust emissivity does not change within the cloud; (c) that the magnetic field direction is uniform; and (d) that the geometry of the cloud is slab-like.

Zoom In Zoom Out Reset image size

As the p versus N and p versus S trends appear to be independent (see Section 7) we compare predictions from our models with ${p}^{[S]}$, the polarization fraction decorrelated from S (Figure 15, upper panel). Figure 18 shows the predicted p from three models of $\epsilon (\chi )$ compared to our measurements of ${p}^{[S]}$ versus A_V, where here N has been converted to A_V assuming N(H)/A_V = 2 × 10²¹ cm⁻² (Bohlin et al. 1978 with R_V = 3.1). Because of the normalization of ${p}^{[S]}$ the overall level of polarization is somewhat arbitrary, but of the right order.

Zoom In Zoom Out Reset image size

Here we describe the three plausible models of $\epsilon (\chi )$ shown in Figure 18:

Constant  Model: if the intrinsic polarization efficiency were constant throughout the cloud we would expect no dependence of p on N. The best fit to this model is shown as a dashed–dotted line in Figure 18.

Skin Depth Model: alternatively, we can consider a model where the intrinsic polarization efficiency is constant up to an extinction depth ${\chi }_{\mathrm{crit}}$ and zero thereafter; i.e., a diffuse layer near the cloud surface is responsible for all of the polarized emission and the dust at cloud depths above ${\chi }_{\mathrm{crit}}$ does not contribute to the polarized emission. For a sightline of total extinction A_V the maximum value of χ is ${A}_{{\rm{V}}}/2$. We express  as

$$\begin{eqnarray}\epsilon (\chi )=\left\{\begin{array}{ll}{\epsilon }_{0}, & \mathrm{for}\ \chi \leqslant {\chi }_{\mathrm{crit}};\\ 0, & \mathrm{for}\ \chi \gt {\chi }_{\mathrm{crit}},\end{array}\right.\end{eqnarray} \tag{ 12 }$$

where ${\epsilon }_{0}$ is a constant, and from this we can calculate the total polarized intensity:

$$\begin{eqnarray}P\; & = & \;2{\displaystyle \int }_{0}^{{A}_{{\rm{V}}}/2}\epsilon (\chi )d\chi \;\\ & = & \left\{\begin{array}{ll}{A}_{{\rm{V}}}{\epsilon }_{0}, & \mathrm{for}\ {A}_{{\rm{V}}}\;\leqslant \;{A}_{{\rm{V}}\mathrm{crit}};\\ {A}_{\mathrm{Vcrit}}{\epsilon }_{0}, & \mathrm{for}\ {A}_{{\rm{V}}}\;\gt \;{A}_{{\rm{V}}\mathrm{crit}},\end{array}\right.\end{eqnarray} \tag{ 13 }$$

where we define ${A}_{\mathrm{Vcrit}}=2{\chi }_{\mathrm{crit}}$. The percentage polarization for a given sightline is then

$$\begin{eqnarray}p=\displaystyle \frac{P}{I}=\left\{\begin{array}{ll}{p}_{0}, & \mathrm{for}\ {A}_{{\rm{V}}}\;\leqslant \;{A}_{{\rm{V}}\mathrm{crit}};\\ {p}_{0}\;\tfrac{{A}_{{\rm{V}}\mathrm{crit}}}{{A}_{{\rm{V}}}}, & \mathrm{for}\ {A}_{{\rm{V}}}\;\gt \;{A}_{{\rm{V}}\mathrm{crit}}.\end{array}\right.\end{eqnarray} \tag{ 14 }$$

In the "skin depth" model, p is constant for sightlines with ${A}_{{\rm{V}}}\;\leqslant \;{A}_{\mathrm{Vcrit}}$ and decreases with a power-law slope of −1.0 for sightlines with ${A}_{{\rm{V}}}\;\gt \;{A}_{\mathrm{Vcrit}}$. The dashed line in Figure 18 shows a fit to the skin depth model.

Power-law Model: finally we consider a model where the polarization efficiency is constant up to ${\chi }_{\mathrm{crit}}$ and thereafter decreases as a power-law with coefficient η:

$$\begin{eqnarray}\epsilon (\chi )=\left\{\begin{array}{ll}{\epsilon }_{0}, & \mathrm{for}\ \chi \leqslant {\chi }_{\mathrm{crit}};\\ {\epsilon }_{0}{\left(\displaystyle \frac{\chi }{{\chi }_{\mathrm{crit}}}\right)}^{\eta }, & \mathrm{for}\ \chi \gt {\chi }_{\mathrm{crit}}.\end{array}\right.\end{eqnarray} \tag{ 15 }$$

This model simulates a constant  for the diffuse outer cloud layers and a decreasing  at greater cloud depths. The polarized intensity for a given sightline described by the power-law model is:

$$\begin{eqnarray}P\;=\left\{\begin{array}{ll}{\epsilon }_{0}\;{A}_{\mathrm{Vcrit}}\;a, & \mathrm{for}\ a\;\leqslant \;1;\\ {\epsilon }_{0}\;{A}_{\mathrm{Vcrit}}\left(1\;+\;\tfrac{1}{\zeta }[{a}^{\zeta }-1]\right), & \mathrm{for}\ a\;\gt \;1,\end{array}\right.\end{eqnarray} \tag{ 16 }$$

where $\zeta =\eta \;+\;1$ and $a\;\equiv \;{A}_{{\rm{V}}}/{A}_{\mathrm{Vcrit}}$. The corresponding fractional polarization is then

$$\begin{eqnarray}p\;=\left\{\begin{array}{ll}{p}_{0}, & \mathrm{for}\ a\;\leqslant \;1;\\ {p}_{0}\;{a}^{-1}\left(1\;+\;\tfrac{1}{\zeta }[{a}^{\zeta }-1]\right), & \mathrm{for}\ a\;\gt \;1.\end{array}\right.\end{eqnarray} \tag{ 17 }$$

The power-law  model best fit parameters are ${p}_{0}=0.038$, ${A}_{\mathrm{Vcrit}}$ = 4.5 mag and $\eta =-1.21$ (Figure 18 solid line). Our power-law model fit would imply that at cloud depths of about two magnitudes or greater of visual extinction the intrinsic polarization efficiency decreases with depth into the cloud as $\sim {\chi }^{-1.21}$.

It can be seen that both the skin-depth and power-law model capture the negative slope of the ${p}^{[S]}$ versus N curve for high N. For the purposes of quantifying our ability to trace magnetic fields deep within the cloud, we will use the power-law model as it seems to more closely follow the data points in Figure 18. We also caution that these are all simple models, so the fits should be taken merely as indicative of the trends of  with χ.

Using Equation (16) for a sightline of total dust extinction A_V we can calculate the fractional contribution to the polarized intensity from cloud material at depths of $\lt \chi ^{\prime}$

$$\begin{eqnarray}&&f(\chi ^{\prime} )=\displaystyle \frac{P(\chi ^{\prime} )}{P({\chi }_{\mathrm{max}})}\end{eqnarray} \tag{ 18 }$$

where by definition ${\chi }_{\mathrm{max}}={A}_{{\rm{V}}}/2$. Figure 19 shows $f(\chi ^{\prime} )$ for a sightline of total dust extinction ${A}_{{\rm{V}}}=40\;\mathrm{mag}$ (about the largest found in Vela C) over the range of $\chi ^{\prime} =1$ to $\chi ^{\prime} ={\chi }_{\mathrm{max}}=20\;\mathrm{mag}$. Dashed and solid lines represent different assumptions for ${A}_{\mathrm{Vcrit}}$, and line colors represent different power-law slopes η in Equation (16). The solid black line is derived using the best-fit parameters. For comparison we also show the expected behavior for the constant  model (red dotted line). For our best fit parameters 27% of the polarized emission comes from the outer 2.2 mag of extinction, or the outer 2.2/20 = 11% of the cloud. A further 47%  of the total polarized emission comes from $2.2\;\mathrm{mag}\;\lt \;\chi \;\lt \;10\;\mathrm{mag}$, which accounts for 39% of the dust column. The most deeply embedded regions of the cloud ($10\;\mathrm{mag}\;\lt \;\chi \;\lt \;20\;\mathrm{mag}$) contribute 50% of the total dust column but only 27% of the total polarized emission. Figure 19 also shows that steeper power-law slopes and lower ${A}_{\mathrm{Vcrit}}$ values would imply that more of the total polarized intensity measured comes from the outer diffuse cloud layers.

Zoom In Zoom Out Reset image size

To estimate the fraction of the cloud that, from the perspective of contributing to polarization, is "hidden" we first calculate the  weighted mean cloud depth $\langle \chi \rangle$:

$$\begin{eqnarray}&&\langle \chi \rangle =\displaystyle \frac{{\int }_{0}^{{A}_{V}/2}\;\chi \;\epsilon \;d\chi }{{\int }_{0}^{{A}_{V}/2}\;\epsilon \;d\chi }=\displaystyle \frac{{{\boldsymbol{A}}}_{\mathrm{Vcrit}}}{4}\;\displaystyle \frac{h(\chi )}{g(\chi )},\end{eqnarray} \tag{ 19 }$$

where $g(\chi )=P/({\epsilon }_{0}{A}_{\mathrm{Vcrit}})$ (see Equation (16)) and $h(\chi )$ is given by

$$\begin{eqnarray}h(\chi )=\left\{\begin{array}{ll}a, & {\rm{for}}\ a\;\leqslant \;1;\\ 1\;+\;\tfrac{2}{\zeta \;+\;1}({a}^{\zeta +1}\;-\;1), & {\rm{for}}\ a\;\gt \;1.\end{array}\right.\end{eqnarray} \tag{ 20 }$$

If  were constant throughout the sightline then $\langle \chi \rangle$ would equal half of the maximum value of χ giving $\langle \chi \rangle ={A}_{{\rm{V}}}/4$. The fraction of the cloud that is hidden can then be roughly estimated as

$$\begin{eqnarray}&&{f}_{\mathrm{hidden}}=1.0-\displaystyle \frac{\langle \chi \rangle }{{A}_{{\rm{V}}}/4},\end{eqnarray} \tag{ 21 }$$

which is shown in Figure 20. For A_V = 10 mag (assuming our best fit parameters) only 16%  of the cloud is hidden. For a sightline of A_V = 40 mag about 48% of the cloud is hidden.

Zoom In Zoom Out Reset image size

In summary, for "moderate" dust column sightlines (${A}_{{\rm{V}}}\;\lt$ 10 mag) our polarization measurements sample most of the cloud (${f}_{\mathrm{hidden}}\;\lt \;16 \%$). So for sightlines with dust columns of ${A}_{{\rm{V}}}\;\sim \;10$ mag or less, the BLASTPol 500 μm measurement of the magnetic field orientation should be representative of the density-weighted average magnetic field orientation along the sightline. For higher dust column sightlines, the fraction of the cloud that is not well sampled by our polarization measurements increases (${f}_{\mathrm{hidden}}\;\sim \;34 \%$ for ${A}_{{\rm{V}}}=20$ mag) and for our highest column sightlines (${A}_{{\rm{V}}}\;\sim \;40$ mag) about half of the dust contributes little to the polarization measured by BLASTPol. For these latter sightlines BLASTPol would not be sensitive to changes in magnetic field direction in the most deeply embedded cloud material. Recall that our model assumes that all of the decrease in p with N is due to lower intrinsic polarization efficiency of material deep within the cloud. If some of the decrease in p with N is due to increased field disorder along high N sightlines then the  drop-off with χ would be shallower, which would decrease the fraction of the cloud that is hidden.

As noted earlier, our model has many implicit assumptions (slab-geometry, uniform dust temperature, power-law dependence of $\epsilon (\chi )$). In particular, the assumption of isothermality is clearly incorrect. Figure 10 shows that for the ISRF-heated sightlines included in this analysis the average temperature decreases with increasing column density (and thus increasing A_V). For the temperature extremes of 11 K and 15 K of the ISRF sightlines in Figure 10, we calculate that for the colder highest column density sightlines the dust on average emits half as much radiation per unit mass at 500 μm compared with dust on the warmer more diffuse cloud sightlines. It is quite likely that the more deeply embedded dust grains in Vela C are colder, which implies they will contribute less than warmer grains near the cloud surfaces to both the total intensity and the measured polarized intensity. We therefore expect that the average magnetic field orientation inferred from polarization data will be weighted more toward the orientation in the warmer regions of the cloud. This will increase ${f}_{\mathrm{hidden}}$: even assuming uniform intrinsic polarization efficiency if the outer half of the dust grains had T = 15 K and the inner half of the dust grains had T = 11 K then we find that ${f}_{\mathrm{hidden}}=20 \%$.

To some degree this problem can be reduced by measuring polarization at millimeter wavelengths where the intensity of thermal dust emission is less sensitive to temperature. For detailed statistical comparisons of submillimeter polarization data with synthetic observations of molecular clouds derived from numerical simulations it will be important to not only model the effects of grain alignment in simulation postprocessing but also include a realistic cloud temperature structure. Due to the aforementioned uncertainties that are related to the assumptions in our model, our values of ${f}_{\mathrm{hidden}}$ should be taken only as crude estimates.

Despite these uncertainties, we note that, our results are consistent with the findings of Cho & Lazarian (2005), who showed that dust grains can be aligned efficiently by RATs at cloud depths χ of up to 10 mag in visual extinction. Bethell et al. (2007) found that the exact depth to which grains are aligned depends on grain size and on the degree of anisotropy of the local radiation field. Our model is also consistent with recent observations by Alves et al. (2014) who argue that their observations of submillimeter polarization of a starless core suggest loss of grain alignment at column densities higher than ${A}_{{\rm{V}}}=30$ mag. If  changes appreciably with χ in the cloud then this might be revealed in the frequency dependence of p. Thus studying p at higher frequencies, as can be done using BLASTPol data, might provide further constraints.

In Section 7.2 we showed that we can account for most of the variations in p that we observe in Vela C with a simple two-variable power-law fit $p(N,S)$. Here, we consider a number of factors besides N and S which could contribute to the variance in p. The dispersion in the logarithm of the decorrelated fractional polarization $\mathrm{log}({p}^{[N,S]})$ is 0.15, which corresponds to a variance in ${p}^{[N,S]}$ of 1.0 × 10⁻⁴.

If the variance in p were entirely due to the measurement uncertainty, then we would expect the variance in p to be described by:

$$\begin{eqnarray}&&{\sigma }_{{\rm{p}}\;\mathrm{stat}}^{2}=\sqrt{\displaystyle \frac{{\displaystyle \sum }_{i}^{n}{\sigma }_{{p}_{i}}^{2}}{n}}=1.7\;\times \;{10}^{-6},\end{eqnarray} \tag{ 22 }$$

where ${\sigma }_{{p}_{i}}^{2}$ is the variance for each individual p_i and n is the total number of data points. This value for ${\sigma }_{p\;{stat}}^{2}$ is much smaller than the measured variance in ${p}^{[N,S]}$. Measurement uncertainties thus play a minor part in the observed variance of ${p}^{[N,S]}$.

A more likely possibility is that the variance in ${p}^{[N,S]}$ seen in Figure 16 and the bottom panel of Figure 11 is the result of changes in the direction of the magnetic field with respect to the line of sight. The observed polarization of a population of dust grains aligned with respect to a uniform magnetic field depends on γ, the angle between the magnetic field direction and the plane of the sky:

$$\begin{eqnarray}&&p={p}_{\mathrm{max}}{\mathrm{cos}}^{2}\gamma ,\end{eqnarray} \tag{ 23 }$$

where p_max is the polarization that would be observed if the magnetic field were parallel to the plane of the sky (γ = 0°). Our inferred magnetic field maps (Figure 5) clearly show several large scale changes in magnetic field direction Φ. Corresponding large scale changes in γ would add width to the $\mathrm{log}p$ distribution. In theory a detailed statistical comparison of S and ${p}^{[N,S]}$ on different angular scales could be used to gain insight into the three-dimensional structure of the magnetic field. However, such a treatment is beyond the scope of the present paper.