Submillimeter Polarization Spectrum of the Carina Nebula

Dust grains in the interstellar medium (ISM) have long been known to linearly polarize background starlight at visible and near-infrared wavelengths (Hall 1949; Hiltner 1949). The thermal emission from these dust grains in the far-infrared and submillimeter portion of the spectrum is observed to be linearly polarized as well (Cudlip et al. 1982; Hildebrand et al. 1984). It is believed that the aspherical, spinning dust grains align with their long axes perpendicular to the direction of the local ISM magnetic field. The mechanism for the alignment of the dust grains is an area of active research, but recently the theory of radiative alignment torques (RATs) has gained some observational support. Under the RAT mechanism (Lazarian & Hoang 2007), an external radiation field is able to provide a net torque to spin up dust grains, which develop a net magnetic dipole moment that aligns with the direction of the local magnetic field. A detailed review is given in Andersson et al. (2015). In particular, they note that RAT alignment requires the presence of a radiation field with wavelengths less than the grain diameter. This condition predicts a lower alignment efficiency, and hence a lower polarization fraction, with increasing dust extinction.

The aligned dust grains preferentially absorb and emit light that is polarized in the direction parallel to their long axes. As a result, transmitted starlight is preferentially polarized in the direction parallel to the local magnetic field, whereas the dust thermal emission is polarized in the direction perpendicular to the magnetic field. Dust polarimetry therefore allows the directions of the plane-of-sky projection of the ISM magnetic field to be traced. This measurement can serve as a powerful tool for investigating the role played by magnetic fields in the early stages of star formation. This is true particularly when the polarization is measured in emission at longer wavelengths, where it can be traced into the interiors of the dense clouds of molecular gas where star formation takes place. Measurements of polarized dust emission in the diffuse ISM are also important, since this emission is a source of foreground contamination in studies of the polarization of the cosmic microwave background (CMB) radiation. Data that provide a better understanding of the variation of the dust polarization fraction with wavelength and with dust environment are important to both of these scientific applications.

The submillimeter polarization spectrum p(λ) is the linear polarization fraction of the thermal dust emission as a function of wavelength (polarization fraction being defined in Section 2.5). Typically, the spectrum is divided by the polarization fraction at a reference wavelength λ₀. This normalization removes the dependence on the inclination angle of the magnetic field relative to the line of sight and on any other unknown factors that would affect the observed polarization fraction across all bands. The relevant observable for the study of dust grain alignment is then the shape of this normalized spectrum, $p(\lambda )/p({\lambda }_{0})$.

A number of models attempting to predict the shape of the polarization spectrum over the submillimeter spectrum, in various column density regimes, have been developed. Draine & Fraisse (2009) investigated the polarization of dust in the diffuse ISM using models of aspherical silicate and graphite grains that were constrained to reproduce the observed dust extinction and polarization of starlight as a function of wavelength. They produced model polarization spectra that were increasing from 100 to 1000 μm. For a scenario in which the dust is diffuse enough that all grains are exposed to the same interstellar radiation field, it is not unreasonable to assume that the dust temperature will depend only on grain material and size, and not on physical location. Under this assumption, the rising spectrum can be explained, at least in part, in terms of an anticorrelation between dust temperature and grain alignment. Larger dust grains are known empirically to be better aligned (Kim & Martin 1995). However, larger dust grains also have higher emissivity and thus tend toward lower equilibrium temperatures. Therefore, it is the well-aligned lower temperature grains, the emission of which peaks at longer wavelengths, that contribute predominantly to the polarization fraction in this scenario.

In contrast to this diffuse ISM study, Bethell et al. (2007) modeled polarized dust emission in dense, clumpy molecular clouds and cores using simulations of magnetohydrodynamic turbulence. They found a polarization fraction that is largely flat over submillimeter wavelengths and only begins falling off toward the far-infrared, decreasing by a factor of 2 in percentage polarization between 250 and 100 μm. A possible explanation for this lack of variation is known as the extinction–temperature–alignment correlation (ETAC) effect. Under this effect, grains in the interior of dense molecular clouds are more shielded from the interstellar radiation field than grains on the surface of these clouds. If the RAT mechanism is correct, the shielded interior grains should be both colder and less well aligned, while the surface grains should be warmer and better aligned. This is the inverse temperature–alignment correlation from the one discussed above for the diffuse ISM. The net result is that the average temperature of the grains contributing to polarized emission is closer to the average temperature of all grains, leading to a flatter polarization spectrum. Additional discussion of the ETAC effect appears later in this section and in Section 4.

Previous ground-based observations of the submillimeter polarization spectrum in dense clouds and cores have found large ratios in the polarization fraction between different bands. These observations typically span <0.01 deg², around bright sources. As shown in Vaillancourt & Matthews (2012), the combination of ground-based measurements of different targets in different bands produces a V-shaped polarization spectrum that falls very steeply in the far-infrared, shows a pronounced minimum at 350 μm, and rises very steeply toward millimeter wavelengths (see Figure 9). The models described above are not able to account for the steepness of the observed slopes nor the overall magnitude of the variation.

More recently, high-sensitivity, wide-area mapping observations by the Balloon-borne Large Aperture Submillimeter Telescope for Polarimetry (BLASTPol) at 250, 350, and 500 μm have produced polarization measurements for molecular cloud targets in the Galactic plane. For example, Gandilo et al. (2016) presented the submillimeter polarization spectrum of the Vela C giant molecular cloud (GMC). Ashton et al. (2018) computed the first submillimeter polarization spectrum of a translucent molecular cloud near Vela C on the sky, but having approximately an order of magnitude lower column density. Both of these studies combined data from the three BLASTPol bands with data from the Planck High Frequency Instrument (HFI) 353 GHz (850 μm) band. Both analyses resulted in polarization spectra that were flat to within ± 15% in $p(\lambda )/p({\lambda }_{0})$ over these four bands.

It should be noted that Ashton et al. (2018) attempted to model the ETAC effect analytically using two different populations of grains—bulk and surface—with different temperature–size distributions and alignment fractions. They found that in their observed range of column densities, this implementation of the ETAC effect was not strong enough to account for the observed flatness of the polarization spectrum of the translucent molecular cloud. In other words, the diffuse ISM models of Draine & Fraisse (2009) with no shielding should be applicable to the cloud observed by Ashton et al. (2018). However, these models are rising with wavelength, rather than flat, and can disagree with the translucent cloud data by up to 30% at 250 μm. Further modeling work, such as that of Guillet et al. (2018), attempts to produce flatter polarization spectra by varying the composition, porosity, and oblateness of the aligned grains.

To summarize some of the most recent developments in polarization spectrum analysis: combined BLASTPol/Planck results have produced flat polarization spectra in two different clouds of very different densities. These results disfavor some of the Draine & Fraisse (2009) models in the diffuse case and are in sharp contrast with previous ground-based measurements in the case of dense clouds. Thus, it has become even more important to study the polarization spectra of other targets, preferably having a range of different cloud environments, so as to further explore these discrepancies and to better inform grain-alignment models.

The Carina Nebula (NGC 3372), the largest and highest surface-brightness nebula in the southern sky, appears in visible light as a giant H ii region spanning several square degrees. Located at an estimated distance of 2.3 kpc (Allen & Hiller 1993; Smith 2006), the nebula and surrounding molecular cloud are part of a GMC complex spanning some 150 pc. In the context of the BLASTPol observations, the Carina GMC is perhaps the most active and evolved source observed, the other targets being relatively more quiescent molecular clouds. An overview of the the Carina molecular cloud complex, including the structure of the submillimeter emission, is given in Li et al. (2006). As they note, the central open clusters, Trumpler 14 and 16, contain an unusual concentration of massive stars, including η Carinae and 6 of the 17 O3-type stars in the Galaxy that were known at that time. In contrast, the most massive sources of excitation for the H ii region RCW 36 in Vela C have been measured to be two stars of type O9 and O9.5 (Ellerbroek et al. 2013). For this reason, comparisons of submillimeter polarimetric observations of Carina with other molecular clouds, such as Vela C, might be regarded as a way to probe the effects of radiative environment and internal heating on dust grain alignment, particularly in the context of the RAT mechanism.

In this paper, BLASTPol polarization data from the Carina Nebula at 250, 350, and 500 μm are presented along with Planck HFI 353 GHz (850 μm) data from the same region. A submillimeter polarization spectrum of Carina is produced over these bands following an analysis similar to, but independent from, that of Gandilo et al. (2016) for Vela C. Section 2 describes the BLASTPol instrument, the 2012 science flight, and the steps of the data analysis, including a detailed description of the polarimetric analysis. Section 3 contains the main results of the polarization spectrum analysis for Carina. The implications of these results are discussed in Section 4, and the overall findings of this paper are summarized in Section 5.

2.1. The BLASTPol Instrument and 2012 Flight

2.2. Map Making

BLASTPol maps of the I, Q, and U Stokes parameters (Figure 1) were produced using Time-Ordered Astrophysics Scalable Tools (TOAST),²⁶ a set of code for map making and simulation that can be used serially or with OpenMP/MPI parallelization. The TOAST generalized least-squares (GLS) solver was used, which iteratively inverts the map maker equation using the preconditioned conjugate gradient method. The map maker's input noise model came from per-bolometer TOD power spectral densities estimated from data obtained while observing a low-signal region of sky. It was not necessary for the noise model to include non-stationarity or detector-to-detector noise correlations. Using the input noise model, the map maker produces a 3 × 3 matrix of the (I, Q, U) covariances for each pixel. The TOAST maps were produced using a 10'' pixelization. Data from input Planck HFI 850 μm all-sky maps²⁷ were processed using coordinate information from TOAST to produce 850 μm maps of the Carina Nebula with the same pixelization, angular extent, and map projection as the BLASTPol maps. All of the TOAST Carina maps were produced in equatorial coordinates, which is therefore the coordinate system to which Q and U are referenced throughout this work (Section 2.5). For this analysis, the BLASTPol signal maps were smoothed to a beam size of 48 FWHM to match the resolution of the Planck data. This resolution corresponds to a physical scale of 3.2 pc at a source distance of 2.3 kpc. This is also well above the scale of irregularities that were observed in the BLASTPol beam shape (Fissel et al. 2016). The Lucy–Richardson iterative deconvolution (Richardson 1972; Lucy 1974) was used to deconvolve a model of the BLASTPol beam from a symmetric 48 FWHM Gaussian beam. The result of this deconvolution was then used as the smoothing kernel with which the BLASTPol signal maps were convolved. The covariance maps were smoothed with the square of this normalized kernel.

Zoom In Zoom Out Reset image size

2.3. Calibration

The TOAST BLASTPol maps are initially produced in the same units as the TOD: raw analog-to-digital converter (ADC) counts. They must be calibrated into physical units of MJy sr⁻¹. Calibration is achieved using a dust spectral energy distribution (SED) obtained from the Planck all-sky dust model described in Planck Collaboration XI (2014). This model is also discussed in more detail in Section 3.5. The model dust morphology is defined by the optical depth at a reference frequency. Using the model SED, the dust total intensity at this frequency is scaled to and integrated over top-hat approximations to each of the three BLASTPol bands. This procedure produces model maps for the dust Stokes I_λ at λ = 250, 350, and 500 μm. The pixel values in each of these model maps are fitted to the pixel values in the corresponding BLASTPol map using linear least-squares regression:

$$\begin{eqnarray}&&{I}_{\lambda }(\mathrm{BLASTPol})={G}_{\lambda }\left({I}_{\lambda }(\mathrm{model})+{D}_{\lambda }^{I}\right).\end{eqnarray} \tag{ 1 }$$

In this model, G_λ and ${D}_{\lambda }^{I}$ are, respectively, a single calibration slope and DC offset obtained for each band. The BLASTPol maps Q_λ and U_λ in each band are then calibrated using the above gain as well (${G}_{\lambda }^{-1}$ in MJy sr⁻¹ count⁻¹). However, they also have to be divided by ${\gamma }_{\lambda }$, the measured instrumental polarization efficiencies in each waveband, which are reported in Galitzki et al. (2014a). Although this procedure corrects the map slopes, the polarization maps also have arbitrary DC offsets ${D}_{\lambda }^{Q}$, ${D}_{\lambda }^{U}$ that cannot be determined using Equation (1), because equivalent dust model maps ${Q}_{\lambda }(\mathrm{model})$ and ${U}_{\lambda }(\mathrm{model})$ do not exist. As discussed in the next section, polarization spectrum analysis normally proceeds by subtracting diffuse background emission from the maps before computing polarization quantities. However, a case with no background subtraction is presented here as well, for which the determination of the DC offsets is important. To determine the offsets, the Planck HFI 850 μm map²⁸ of dust emission is color-corrected by scaling it to the dust model map in each of the BLASTPol bands using a linear least-squares regression:

$$\begin{eqnarray}&&{I}_{\lambda }(\mathrm{model})={g}_{\lambda }\left({I}_{850}(\mathrm{HFI})+{d}_{\lambda }\right).\end{eqnarray} \tag{ 2 }$$

In principle, the linearity of this scaling requires that the dust be isothermal at temperature T_d and that it have an emissivity with a constant power-law spectral index β_d (Section 3.5). In practice, it is found that rejection of map pixels that are high-temperature outliers in the all-sky dust model is sufficient for good linearity. Model maps ${Q}_{\lambda }(\mathrm{model})$ and ${U}_{\lambda }(\mathrm{model})$ can then be obtained by color-correcting Q₈₅₀ and U₈₅₀ to the BLASTPol bands using the scale factor g_λ. The DC offsets of the calibrated BLASTPol polarization maps are then obtained by regressing pixel values of ${Q}_{\lambda }(\mathrm{BLASTPol})$ versus ${Q}_{\lambda }(\mathrm{model})$, and likewise for U_λ. For example, the DC offset for the Q maps would be determined using a model

$$\begin{eqnarray}&&{G}_{\lambda }^{-1}{\gamma }_{\lambda }^{-1}{Q}_{\lambda }(\mathrm{BLASTPol})={Q}_{\lambda }(\mathrm{model})+{D}_{\lambda }^{Q},\end{eqnarray} \tag{ 3 }$$

where ${Q}_{\lambda }(\mathrm{model})={g}_{\lambda }{Q}_{850}(\mathrm{HFI})$, while G_λ is the slope obtained from the I calibration in Equation (1), and ${D}_{\lambda }^{Q}$ is the DC offset parameter being fitted for here.

2.4. Diffuse Emission Subtraction

2.5. Polarimetry

The Stokes parameters are used to compute the net linear polarization of the dust emission,

$$\begin{eqnarray}&&P=\sqrt{{Q}^{2}+{U}^{2}},\end{eqnarray} \tag{ 4 }$$

as well as the fractional linear polarization,

$$\begin{eqnarray}&&p=\displaystyle \frac{P}{I}.\end{eqnarray} \tag{ 5 }$$

The angle ψ defining the direction of the linear polarization on the sky is given by

$$\begin{eqnarray}&&\psi =\displaystyle \frac{1}{2}\arctan (U,Q),\end{eqnarray} \tag{ 6 }$$

where the two-argument form of the arctan function is used in order to evaluate the angle quadrant properly. The IAU polarization angle convention is used; for a polarization pseudovector viewed on the sky in equatorial coordinates, ψ increases counterclockwise from 0° in the the north–south direction and ranges from −90° to +90°.

The TOAST pixel covariances are used to compute the variances in these quantities, ${\sigma }_{p}^{2}$ and ${\sigma }_{\psi }^{2}$, using error propagation (see Appendix B). Since P and p are restricted to positive values, any noise in Q and U positively biases these quantities. The polarization fraction is debiased approximately using a rudimentary method that is acceptable for high signal-to-noise in p (Wardle & Kronberg 1974; Montier et al. 2015):

$$\begin{eqnarray}&&{p}_{\mathrm{db}}=\sqrt{{p}^{2}-{\sigma }_{p}^{2}}.\end{eqnarray} \tag{ 7 }$$

All of the p values used in the polarization spectrum analysis (Section 3) have been debiased in this way. Note, however, that the map of P shown in the lower-right panel of Figure 1 is presented for visualization purposes only and has not been debiased.

Once maps of p and ψ have been produced, data cuts that have been customary for submillimeter polarization spectrum analyses in the literature are applied. The first is a signal-to-noise cut on map pixels using a threshold of p_db > 3σ_p. Only pixels for which this condition holds simultaneously in all four bands are kept. The second data cut is intended to mitigate the circumstance in which the different wavebands sample different cloud components along the line of sight, each with differing line-of-sight components of the magnetic field and hence differing polarization fractions. This situation would lead to artificial variation with wavelength in the measured polarization spectrum for the sightline in question, variation that is not intrinsic to any particular physical location within the cloud. Under the assumption that a sightline having a plane-of-sky component of B that is constant with wavelength implies (at least in a statistical sense) a constant line-of-sight component of B as well, the condition is imposed that the difference $| {\rm{\Delta }}\psi | \lt 15^\circ$ between any two of the four wavebands. The stringent maximum angle difference of 10° used in past analyses, including Gandilo et al. (2016), was relaxed for Carina, in order to include more sightlines in the analysis.

After the downsampling and data cuts, polarization data remained for 314, 285, and 261 sightlines, respectively, for the three cases of diffuse emission subtraction using the Far reference region, subtraction using the Near reference region, and no diffuse emission subtraction at all. Figure 2 shows maps of the polarization fraction ratios for those pixels surviving the data cuts in the Far reference region case. These are the ratios p_λ/p₃₅₀ for λ ∈ {250, 500, 850} μm. These discrete values sample the polarization spectrum over this wavelength range and are analyzed in detail in the next section.

Zoom In Zoom Out Reset image size

For all map pixels surviving the data cuts, the left panel of Figure 3 represents the linear polarization at 350 μm as pseudovectors of length p and direction ψ. The pseudovectors are overlaid on filled contours showing the 350 μm intensity over the cloud. The "polarization-hole effect" is evident here, in which the polarization fraction is lower near bright intensity peaks. This effect has been noted in past submillimeter polarimetry observations (Matthews et al. 2001). The right panel of Figure 3 shows the corresponding inferred directions (but not magnitudes) of the plane-of-sky component of the magnetic field, B_⊥. These directions are rotated by 90° relative to the ${\boldsymbol{E}}$-field polarization direction of the radiation.

Zoom In Zoom Out Reset image size

3.1. Median Polarization Ratios

Figure 4 shows histograms of the polarization ratios p_λ/p₃₅₀, specifically the distributions of these ratios over the cloud for the case of diffuse emission subtraction using the Far reference region. For consistency with Gandilo et al. (2016), the widths of the distributions were quantified using the median absolute deviation (MAD), defined as

$$\begin{eqnarray}&&\mathrm{MAD}\equiv \mathrm{median}\left(| {x}_{i}-\mathrm{median}({x}_{i})| \right),\end{eqnarray} \tag{ 8 }$$

where the quantities x_i are the measurements in question. Table 1 lists the median ratios and MADs for all three types of diffuse emission subtraction. Although the two cases with background subtraction show a slight minimum at 500 μm, none of the results are significantly different from a flat spectrum, with a polarization ratio of unity to within 15% in each band. This result is independent of the method of diffuse emission subtraction—a very similar outcome to that of Gandilo et al. (2016).

Zoom In Zoom Out Reset image size

Table 1.  Medians and MADs of Polarization Ratios (p_λ/p₃₅₀)

Diffuse Emission	250 μm	500 μm	850 μm
Subtraction Method
Far	1.00 ± 0.06	0.93 ± 0.10	1.01 ± 0.12
Near	1.02 ± 0.06	0.93 ± 0.08	0.99 ± 0.12
None	1.14 ± 0.08	0.96 ± 0.08	0.95 ± 0.12

Download table as:  ASCIITypeset image

3.2. Polarization Ratios from Scatter Plots of p_λ versus p₃₅₀

An alternative method for determining the polarization ratios, averaged over the cloud, is to produce linear fits to scatter plots of p_λ versus p₃₅₀. The polarization spectrum then consists of the best-fit linear slopes as a function of wavelength. For this fitting procedure, the least absolute deviation was used to optimize the fit parameters. This method is more robust to outliers than least-squares fitting. For each fit, an uncertainty on the slope was estimated using bootstrap resampling (Press et al. 1992). The fit was repeated for each of 10,000 random selections of the data points (with replacement), and the uncertainty was taken to be the standard deviation of this ensemble of fit parameter values. The linear fits are shown in Figure 5 for the case of diffuse emission subtraction using the Far reference region.

Zoom In Zoom Out Reset image size

Table 2 lists the slopes of the linear fits to each waveband, along with their uncertainties, for the cases of background subtraction using the Far reference region, using the Near reference region, and for the case of no background subtraction. For the Near and Far cases, the feature of a very slight minimum at 500 μm occurs using this method, just as it did for the median ratios. The polarization spectra obtained using the linear fitting are once again flat to within ±15%.

Table 2.  Slopes of Linear Fits to p_λ vs. p₃₅₀

Diffuse Emission	250 μm	500 μm	850 μm
Subtraction Method
Far	1.02 ± 0.01	0.98 ± 0.02	1.01 ± 0.02
Near	1.05 ± 0.01	0.95 ± 0.02	1.03 ± 0.03
None	1.14 ± 0.02	1.02 ± 0.02	0.90 ± 0.04

Download table as:  ASCIITypeset image

3.3. Fits to ${p}({\lambda })$

Two different functional forms for p(λ) were fitted to the per-pixel p measurements across the bands: a power law,

$$\begin{eqnarray}&&p(\lambda )={a}_{1}{\left(\displaystyle \frac{\lambda }{{\lambda }_{0}}\right)}^{{b}_{1}},\end{eqnarray} \tag{ 9 }$$

and a second-order polynomial,

$$\begin{eqnarray}&&p(\lambda )={a}_{2}\left[{b}_{2}{\left(\lambda -{\lambda }_{0}\right)}^{2}+{c}_{2}\left(\lambda -{\lambda }_{0}\right)+1\right].\end{eqnarray} \tag{ 10 }$$

Here, λ₀ = 350 μm, and for both models, a is an overall normalization constant (it is the fitted value of p₃₅₀). The models are intended to probe the shape of the per-pixel polarization spectra over the wavelength range of 250–850 μm. The power-law model investigates whether the spectra are increasing or decreasing. The quadratic model, in addition, allows for the polarization spectra to have minima or maxima somewhere within this wavelength range.

Figure 6 shows the results of the power-law and quadratic fitting for three example pixels. In these plots, the error bars in each band are derived from the TOAST covariances (Appendix B). After fitting to the per-pixel spectra, the distributions of the fit parameters were analyzed. The medians and MADs are listed in Table 3 for the fit parameters relevant to the spectral shape, i.e., b₁, b₂, and c₂. This table also lists the median ± MAD values of the scale factors a₁ and a₂ in units of percent polarization, as well as listing the ratios p₃₅₀/a. These ratios indicate the extent to which the model median value for the fractional polarization at 350 μm matches the measured median value. The distributions of the spectral-shape fit parameters for the Far reference region case are shown in Figure 7 for the power-law fit and in Figure 8 for the quadratic fit.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3.  Medians and MADs of the p(λ) Fit Parameters

Diffuse Emission	Power-law Fit			Polynomial Fit
Subtraction Method	b1	a1 [%]	p350/a1	b2 (×10−6)	c2 (×10−4)	a2 [%]	p350/a2
Far	−0.01 ± 0.10	4.8 ± 1.5	1.01 ± 0.03	0.8 ± 1.3	−3.1 ± 4.9	4.8 ± 1.5	1.03 ± 0.03
Near	−0.04 ± 0.10	5.5 ± 2.0	1.01 ± 0.03	0.9 ± 1.4	−4.2 ± 4.6	5.5 ± 2.0	1.03 ± 0.03
None	−0.14 ± 0.10	3.0 ± 1.0	0.96 ± 0.04	1.1 ± 1.1	−8.0 ± 4.6	3.0 ± 1.0	0.96 ± 0.03

Download table as:  ASCIITypeset image

3.4. Summary of Polarization Spectrum Measurements

The result of a flat spectrum for Carina computed using all of the previously described methods is shown in Figure 9 for the case of diffuse emission subtraction using the Far reference region. The median polarization ratios ± MADs (Section 3.1) are shown as red triangles. The polarization ratios from linear fits to scatter plots of p_λ versus p₃₅₀ are shown as red circles, with error bars based on bootstrap resampling (Section 3.2). Representative power-law and quadratic fits to per-pixel polarization spectra (Section 3.3) are also shown. For the power-law fit, the mean and dispersion among the per-pixel polarization spectra are demonstrated by plotting the power-law model corresponding to $\mathrm{median}({b}_{1})$ as a solid magenta line and plotting the models corresponding to $\mathrm{median}({b}_{1})\pm \mathrm{MAD}({b}_{1})$ as dashed magenta lines. Similarly, for the second-order polynomial fit, the parabola corresponding to the median values of b₂ and c₂ is plotted as a solid dark blue line. However, the fit parameters b₂ and c₂ are highly anti-correlated, with Pearson correlation coefficient ${\rho }_{{c}_{2},{b}_{2}}=-0.92$ for the Far reference region case. Therefore, it was not sufficient to simply plot extremal models using the median ± MAD values of each fit parameter individually. The 68% error ellipse of their joint distribution was constructed by diagonalizing the covariance matrix of b₂ and c₂ using eigenvalue decomposition. The major axis of the resulting ellipse had endpoints at $({c}_{2},{b}_{2})=(-1.8\times {10}^{-3},4.2\times {10}^{-6})$ and at $({c}_{2},{b}_{2})=(9.9\times {10}^{-4},-2.3\times {10}^{-6})$. The dashed dark blue lines in Figure 9 correspond to the parabolas having these extremal fit parameter values.

Zoom In Zoom Out Reset image size

3.5. Effect of Environment

An investigation was undertaken to probe whether the shape of the polarization spectrum over these four bands exhibits any dependence on the molecular cloud environment. Polarization spectrum parameters were correlated with two environmental parameters: dust temperature T_d and dust optical depth at 353 GHz τ₃₅₃. The latter is proportional to the dust column density. These parameters were obtained from the Planck all-sky dust model (Planck Collaboration XI 2014) first mentioned in Section 2.3. This model is generated by fitting a modified blackbody SED to the high-frequency dust I maps from HFI at 353, 545, and 857 GHz, along with the highest frequency map from IRAS 100 μm data. This SED is of the form

$$\begin{eqnarray}&&{I}_{\nu }(\nu )={\tau }_{{\nu }_{0}}{\left(\displaystyle \frac{\nu }{{\nu }_{0}}\right)}^{{\beta }_{d}}{B}_{\nu }(\nu ,{T}_{d}).\end{eqnarray} \tag{ 11 }$$

In Equation (11), β_d is the power-law spectral index of the frequency-dependent dust emissivity, B_ν is the Planck function, and T_d is the dust temperature. The parameter ${\tau }_{{\nu }_{0}}$ is the dust optical depth at a reference frequency of ν₀ = 353 GHz. In Planck Collaboration XI (2014), it is emphasized that these three parameters are only approximations to the true dust properties. A single-component model has been assumed, whereas in reality multiple temperature components could exist along any given line of sight. Therefore, the model parameters τ₃₅₃ and T_d are used here only to establish the relative ordering among sightlines in order to search for very obvious trends, which, if present, would lend themselves to a more detailed future investigation.

Figure 10 shows the results of plotting the power-law and polynomial fit parameters (b₁, b₂, and c₂ defined in Section 3.3) versus temperature and optical depth. In addition to the points for every individual sightline, the fit parameters are binned into seven evenly spaced bins in T_d or τ₃₅₃. The binned curves appear to be quite flat. For the most part, the mean values of the fit parameters within each bin (orange points) lie well within the median ± MAD range of the fit parameters over the whole cloud (bounded by the gray dashed lines). Exceptions to this are the first T_d bin for c₂ and b₂, and the last T_d bin for b₁. However, these bins only contain a handful of points each. Each panel of Figure 10 also shows the value of the Pearson correlation coefficient ρ for the scatter plot in question. Overall, strong evidence for trends in polarization spectrum shape with environment are not observed.

Zoom In Zoom Out Reset image size

Gandilo et al. (2016) carried out a very similar analysis for Vela C. The shape parameters from fits of the same models to the per-pixel polarization spectra were plotted against T_d and against column density N_H in order to search for a dependence of spectral shape on environment. No significant trends were found. That analysis used a dust SED model fitted to the Herschel SPIRE 160, 250, 350, and 500 μm maps of Vela C. This fit assumed β_d = 2 and obtained values for T_d and N_H. For the Carina Nebula, however, due to the limited spatial extent of the available Herschel maps, it was not possible to carry out such a fit using the exact same reference region for diffuse emission subtraction as was used for the BLASTPol maps. Thus, we opted to use the Planck all-sky dust model instead. As an additional benefit, using the Planck-derived τ₃₅₃ as the environmental density parameter, rather than computing N_H, also avoids making assumptions about the dust opacity, ${\sigma }_{353}={\tau }_{353}/{N}_{H}$, within our region.

As shown in Figure 9, the result of a polarization spectrum that is flat to within ±15% is inconsistent with the results of measurements of a number of molecular clouds obtained by ground-based submillimeter telescopes. Recall (Section 1) that the BLASTPol/Planck submillimeter polarization spectrum for the Vela C GMC presented by Gandilo et al. (2016) was also much flatter than the ground-based spectra. As a potential explanation for this discrepancy, Gandilo et al. (2016) posited that ground-based observations were limited to the densest sightlines within molecular clouds, in a very different regime of column density compared to the cloud-averaged case for BLASTPol. This can be quantified using the Planck 850 μm intensity as a proxy for column density. Although this proxy (unlike τ) is not completely free from dependence on dust temperature, it was the only method available for comparing the column density between the BLASTPol and ground-based measurements. Table 4 below shows a comparison of the median and interquartile range of the 850 μm intensity computed for Carina (this work), for Vela C (Gandilo et al. 2016), for a translucent molecular cloud (Ashton et al. 2018), and for ground-based measurements of 17 molecular clouds, which were calculated by Gandilo et al. (2016) using online data provided by Vaillancourt & Matthews (2012). The comparison shows that the median 850 μm intensity for the Carina Nebula is quite similar to that of the Vela C measurement, indicating a broadly similar regime in column density. In comparison, the median intensity of the ground-based measurements is nearly two orders of magnitude higher.

Another interesting comparison is to the analysis of the most recent Planck data release (PR3) in Planck Collaboration XI (2018). The information equivalent to the polarization spectrum is reported by Planck in terms of the difference in power-law spectral indices of the SEDs for total intensity and for polarization, ${\beta }_{d}^{I}$ and ${\beta }_{d}^{P}$. Under the simplistic assumption that the dust grains contributing to polarized and to overall emission are all isothermal at temperature T_d, these two SEDs can be written as $I(\nu )\propto {\left(\nu /{\nu }_{0}\right)}^{{\beta }_{d}^{I}}{B}_{\nu }(\nu ,{T}_{d})$ and $P(\nu )\propto {\left(\nu /{\nu }_{0}\right)}^{{\beta }_{d}^{P}}{B}_{\nu }(\nu ,{T}_{d})$. It can then be shown that the quantity computed in this work, the polarization ratio, is given by

$$\begin{eqnarray}&&\displaystyle \frac{p(\nu )}{p({\nu }_{0})}={\left(\displaystyle \frac{\nu }{{\nu }_{0}}\right)}^{{\beta }_{d}^{P}-{\beta }_{d}^{I}}.\end{eqnarray} \tag{ 12 }$$

Therefore, a measurement of ${\beta }_{d}^{P}\lt {\beta }_{d}^{I}$ would imply a polarization spectrum that was falling with frequency, or rising with wavelength. In the analysis of the PR3 maps, polarized emission was investigated for six nested sky patches at high Galactic latitude. A model was fitted to the angular cross-power spectra between Planck frequency maps. This model included a power-law dependence on the multipole moment (ℓ), a power-law frequency dependence for polarized synchrotron emission, and a modified blackbody SED with power-law emissivity for polarized dust emission at an assumed dust temperature T_d = 19.6 K. The result, averaged over all sky patches and multipole bins, was ${\beta }_{d}^{P}=1.53\pm 0.02$. No statistically significant difference from the spectral index for total intensity was found: ${\beta }_{d}^{P}-{\beta }_{d}^{I}\,=0.05\pm 0.03$ (Planck Collaboration XI 2018). This result implies a flat polarization spectrum, on average, at millimeter wavelengths.

In addition to this overall agreement with the high-latitude diffuse emission, there is now striking agreement among BLASTPol molecular cloud targets in the Galactic plane. The result of a flat spectrum in polarization fraction has now been reproduced for measurements of a cold, dense molecular cloud (Gandilo et al. 2016), a translucent molecular cloud (Ashton et al. 2018), and a more evolved, warmer molecular cloud being actively heated by many internal stellar sources (this work). A direct comparison of these three results is shown in Figure 11. Some care should be taken in the interpretation of this figure, since the polarization spectrum from Ashton et al. (2018) is normalized to the 850 μm band, whereas the other two polarization spectra are normalized to 350 μm. Even so, despite the differing radiative environments and densities of these targets, the result of a flat submillimeter polarization spectrum persists as a common property of these molecular clouds.

Zoom In Zoom Out Reset image size

This result of a lack of significant wavelength variation in the polarization spectrum for a cloud such as the Carina GMC, with significant internal heating from embedded sources, potentially has implications for theoretical grain-alignment mechanisms. For example, the RAT mechanism (Section 1) predicts a higher grain-alignment efficiency for dust grains that are less shielded from stellar radiation (Andersson et al. 2015). However, a detailed investigation of the implications of our result for grain-alignment theory is beyond the scope of this work. For the purpose of a baseline comparison to theory, we examine the model of a cold, dense molecular cloud with no internal sources presented in Bethell et al. (2007). Figure 12 shows this model over the wavelength range of interest. Overlaid are the previous ground-based polarization spectrum measurements first shown in Figure 9. In order to match the analysis method of Bethell et al. (2007), the BLASTPol/Planck data for Carina and for Vela C from Gandilo et al. (2016) are shown here as the "total polarized fraction," which is computed as $({\displaystyle \sum }_{j}{P}_{j})/({\sum }_{j}{I}_{j})$ using debiased values of P, where j indexes the pixels passing the data cuts. The results for these two BLASTPol targets are broadly consistent with the Bethell model, certainly showing a much closer correspondence than the previous ground-based studies.

Zoom In Zoom Out Reset image size

The negative slope of the ground-based polarization spectra in the far-infrared has been explained heuristically in terms of a very strong ETAC effect (Section 1), perhaps due to internal sources, in which hotter, less shielded dust grains are preferentially more aligned, leading to a higher polarization fraction for the dust emission at shorter wavelengths (Hildebrand et al. 1999; Zeng et al. 2013). The right-hand side of the V is more problematic, however: slopes of this steepness are not reproduced in the BLASTPol data. Furthermore, no theoretical models that can reproduce the very large polarization ratios that are seen, in general, in the ground-based spectra have been constructed. Regardless, we can conclude that the cloud-averaged submillimeter polarization spectra of two very different molecular clouds agree with the flat Bethell et al. (2007) prediction for externally illuminated, dense molecular clouds with no internal radiation sources. Thus, it appears, at least from the BLASTPol measurements, that the internal sources do not significantly affect the cloud-averaged submillimeter polarization spectrum.

In light of the discrepancy between BLASTPol and ground-based submillimeter polarization spectrum results, a logical step for future work would be to repeat the measurements of the specific ground-based targets using more advanced ground-based polarimeters of higher sensitivity and mapping speed that are slated to come online in the near future. One example is the TolTEC camera (Bryan et al. 2018), which will probe wavelengths of 1 mm and longer. Another interesting avenue for future research is a more detailed investigation of the spatial variation of the polarization spectrum from point to point within a cloud, and the correlation of the spectral shape with environment. The BLASTPol data show some evidence of spatial variation (Figure 2), but no strong trends with environment. Ultimately, this investigation is limited by the angular resolution of the data. The BLAST-TNG experiment (Galitzki et al. 2014b), which is currently scheduled for an Antarctic balloon flight in the Austral summer of 2018–2019, will greatly aid this effort. This polarimeter offers an order of magnitude more detectors than BLASTPol and will observe in the same bands with a much higher angular resolution of 31''–59'' FWHM. The combination of high resolution and full cloud-scale coverage offered by this experiment may also help test the hypothesis stated above that the discrepancy between the BLASTPol/Planck and ground-based measurements is due to the fact that the latter were only sensitive to the densest clumps within GMCs.

Measurements of the linear polarization along 314 sightlines were made by BLASTPol in the Carina Nebula in the 250, 350, and 500 μm wavebands. These data were combined with Planck 353 GHz (850 μm) data from the same region in order to produce submillimeter polarization spectra. These spectra were calculated using several methods, including the median polarization ratios, slopes from linear fits to scatter plots of p_λ versus p₃₅₀, and by fitting quadratic and power-law models p(λ) to the per-pixel polarization spectra. No strong evidence was found for variation of the fitted parameters of these models as a function of cloud environment, as quantified by the dust temperature T_d and dust optical depth τ₃₅₃, which were obtained from the Planck all-sky dust model. The cloud-averaged polarization spectrum of the Carina Nebula appears flat to within ±15% in the polarization ratio quantity p_λ/p₃₅₀, where p_λ is the fractional linear polarization in a given waveband. This is at odds with previous ground-based measurements of the polarization spectrum of other molecular clouds, which showed a V-shaped spectrum, with a negative slope in the far-infrared, a positive slope toward millimeter wavelengths, and a pronounced minimum near 350 μm. The flatness of the polarization spectrum in Carina is, however, in remarkably close agreement with BLASTPol/Planck measurements in other molecular clouds, including the measurement of the Vela C GMC (Gandilo et al. 2016), and of a translucent molecular cloud near Vela C on the sky (Ashton et al. 2018). The shapes of the Vela C and Carina polarization spectra are both in relatively good agreement with the Bethell et al. (2007) theoretical prediction for an externally illuminated, dense molecular cloud with no internal radiation sources.

The BLASTPol collaboration acknowledges support from NASA (through grant Nos. NAG5-12785, NAG5-13301, NNGO-6GI11G, and NNX0-9AB98G, and the Illinois Space Grant Consortium), the Canadian Space Agency, the Leverhulme Trust through the Research Project Grant F/00 407/BN, Canada's Natural Sciences and Engineering Research Council, the Canada Foundation for Innovation, the Ontario Innovation Trust, and the US National Science Foundation Office of Polar Programs. This work was based in part on observations obtained with Planck (http://www.esa.int/Planck), an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada. C.B.N. also acknowledges support from the Canadian Institute for Advanced Research. F.P.S. is supported by CAPES grant 2397/13-7. F.P. acknowledges the European Union's Horizon 2020 research and innovation programme under grant agreement number 687312 (RADIOFOREGROUNDS). P.A. is supported through Reach for the Stars, a GK-12 program supported by the National Science Foundation under grant DGE-0948017. L.M.F. is a Jansky Fellow of the National Radio Astronomy Observatory (NRAO). NRAO is a facility of the National Science Foundation (NSF), operated under cooperative agreement by Associated Universities, Inc. Finally, we thank the Columbia Scientific Balloon Facility staff for their outstanding work.

During the initial polarization spectrum analysis of Carina, anomalously high values were found for the polarization ratio p₈₅₀/p₃₅₀ for sightlines in the northern part of the cloud. The problem was determined to be in the BLASTPol data: receiver 1/f noise led to the BLASTPol I maps having a large-scale gradient in the north–south direction, which is perpendicular to the predominant scan direction during Carina observations. Systematically high values of I₂₅₀, I₃₅₀, and I₅₀₀ led to systematically low values of p₂₅₀, p₃₅₀, and p₅₀₀ in the north of the cloud. This bias in turn inflated the Planck-to-BLASTPol polarization ratio. Since the striping is predominantly in I, it could be corrected by differencing calibrated BLASTPol maps with publicly available Herschel SPIRE maps taken of Carina in the same bands.

The Herschel maps were first regridded to match the BLASTPol pixelization. A least-squares fit of a fifth-order polynomial was applied to the (BLASTPol – Herschel) intensity values (for pixels within the left Far reference region rectangle) as a function of the pixel y-coordinate. The fitted intensity as a function of the map y-coordinate was extended horizontally across the entire map to produce a model for the vertical gradient. This model was then subtracted from the BLASTPol I maps in each band.

Only pixels within the left Far reference region rectangle were used in the fit, because there was a slight beam mismatch between the smoothed BLASTPol and smoothed Herschel intensity maps, leading to leakage of cloud structure (including intensity peaks) into the (BLASTPol – Herschel) difference map. The left Far reference region rectangle happens to be free of such structure, meaning that most of the flux variation within it is due to the gradient alone. As the analysis progressed, more refined data cuts excluded some of the more extreme northern sightlines from the final data set. As a result, repeating the analysis of this paper with diffuse emission subtraction using the Far reference region, and without the Herschel correction, results in relatively small changes to the median ratios of ${\rm{\Delta }}\left({p}_{\lambda }/{p}_{350}\right)=-0.007,+0.019$, and −0.048 at 250, 500, and 850 μm, respectively. Since the systematic does not change the final result of a polarization spectrum that is flat to within ±15% over all four bands, it was not deemed to be necessary to repeat the Herschel correction more carefully using maps of matching resolution.

The per-pixel variances of the polarization quantities, P, p, and ψ, are computed from the per-pixel covariance matrix of the Stokes parameters on the sky (as estimated by the TOAST map maker) using the following procedure. Each of the above three quantities can be expressed as a function $f(I,Q,U)$. If this function is approximated by its first-order Taylor series expansion about the mean of the pixel Stokes parameter values, then the variance in f is given by

$$\begin{eqnarray}{\sigma }_{f}^{2}=\left[\begin{array}{ccc}\tfrac{\partial f}{\partial I} & \tfrac{\partial f}{\partial Q} & \tfrac{\partial f}{\partial U}\end{array}\right]\left[\begin{array}{ccc}{\sigma }_{I}^{2} & {\sigma }_{{IQ}} & {\sigma }_{{IU}}\\ {\sigma }_{{IQ}} & {\sigma }_{Q}^{2} & {\sigma }_{{QU}}\\ {\sigma }_{{IU}} & {\sigma }_{{QU}} & {\sigma }_{U}^{2}\end{array}\right]\left[\begin{array}{c}\tfrac{\partial f}{\partial I}\\ \tfrac{\partial f}{\partial Q}\\ \tfrac{\partial f}{\partial U}\end{array}\right].\end{eqnarray} \tag{ 13 }$$

Applying Equation (13) to each of Equations (4)–(6) above yields the following results for the variances of the polarization quantities:

$$\begin{eqnarray}&&{\sigma }_{P}^{2}=\displaystyle \frac{1}{{P}^{2}}\left[\left({Q}^{2}\right){\sigma }_{Q}^{2}+\left({U}^{2}\right){\sigma }_{U}^{2}+\left(2{QU}\right){\sigma }_{{QU}}\right]\end{eqnarray} \tag{ 14 }$$

and

$$\begin{eqnarray}&&{\sigma }_{\psi }^{2}={\left(\displaystyle \frac{180^\circ }{\pi }\right)}^{2}\displaystyle \frac{1}{4{P}^{4}}\left[\left({U}^{2}\right){\sigma }_{Q}^{2}+\left({Q}^{2}\right){\sigma }_{U}^{2}-\left(2{QU}\right){\sigma }_{{QU}}\right].\end{eqnarray} \tag{ 15 }$$

For notational convenience, we can define normalized Stokes parameters for linear polarization, given by $q\equiv Q/I$ and $u\equiv U/I$. Applying Equation (13) to these, we have

$$\begin{eqnarray}&&{\sigma }_{q}^{2}={q}^{2}\left(\displaystyle \frac{{\sigma }_{I}^{2}}{{I}^{2}}+\displaystyle \frac{{\sigma }_{Q}^{2}}{{Q}^{2}}-2\displaystyle \frac{{\sigma }_{{IQ}}}{{IQ}}\right),\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{\sigma }_{u}^{2}={u}^{2}\left(\displaystyle \frac{{\sigma }_{I}^{2}}{{I}^{2}}+\displaystyle \frac{{\sigma }_{U}^{2}}{{U}^{2}}-2\displaystyle \frac{{\sigma }_{{IU}}}{{IU}}\right),\end{eqnarray} \tag{ 17 }$$

and

$$\begin{eqnarray}&&{\sigma }_{{qu}}={qu}\left(\displaystyle \frac{{\sigma }_{I}^{2}}{{I}^{2}}+\displaystyle \frac{{\sigma }_{{QU}}}{{QU}}-\displaystyle \frac{{\sigma }_{{IQ}}}{{IQ}}-\displaystyle \frac{{\sigma }_{{IU}}}{{IU}}\right).\end{eqnarray} \tag{ 18 }$$

The variance in p can therefore be expressed compactly as

$$\begin{eqnarray}&&{\sigma }_{p}^{2}=\displaystyle \frac{1}{{p}^{2}}\left[\left({q}^{2}\right){\sigma }_{q}^{2}+\left({u}^{2}\right){\sigma }_{u}^{2}+\left(2{qu}\right){\sigma }_{{qu}}\right].\end{eqnarray} \tag{ 19 }$$