Magnetic fields in the Horsehead Nebula

Polarized dust emission from the Horsehead Nebula was observed using the POL-2 polarimeter inserted in front of the SCUBA-2 bolometer camera (Holland et al., 2013) on the JCMT at 450 and 850 $\mu$m simultaneously, although we consider only the 850 $\mu$m data in this work. The field was mapped 10 times in 32-minute observing blocks with SCUBA-2/POL-2 on 2018 November 21. The total on-source time was about 5 hours 20 minutes. The observations were conducted in JCMT Weather Band 1, in which the atmospheric opacity at 225 GHz ($\tau_{225}$) $\leq$ 0.05. The effective beam size of the JCMT is 14$\farcs$1 at 850 $\mu$m (Dempsey et al., 2013).

We reduced the data using the Submillimetre User Reduction Facility (SMURF) package (Chapin et al., 2013) from the Starlink software suite (Currie et al., 2014). The $pol2map$⁴⁴4http://starlink.eao.hawaii.edu/docs/sun258.htx/sun258ss73.html command is used to reduce POL-2 data. In the first step of the data reduction process, $pol2map$ is used to convert raw bolometer timestreams for each observation into separate Stokes Q, U, and I timestreams, and to create initial Stokes $I$ maps from the Stokes $I$ timestreams using the SMURF routine $makemap$ (Chapin et al., 2013). The coadd of the initial Stokes $I$ maps is used to make masks defining areas of astrophysical emission. In the second step of the data reduction process, $pol2map$ is used to create final Stokes $I$, $Q$ and $U$ maps and a polarization segment (or half-vector) catalogue. In this second step, we used the $skyloop$ and $mapvar$ parameters in $pol2map$. $skyloop$ improves signal-to-noise ratio and image fidelity by reducing all 10 observations concurrently rather than consecutively, while $mapvar$ calculates the variances of the Stokes $I$, $Q$ and $U$ maps from the spread of pixel data values between the individual observations. We used the August 2019⁵⁵5https://www.eaobservatory.org/jcmt/2019/08/new-ip-models-for-pol2-data/ instrumental polarization model to correct the Stokes $Q$ and $U$ maps. A detailed description of the the POL-2 data reduction process is given in, e.g., Pattle et al. (2021); Könyves et al. (2021). In the second step of the data reduction process, we used a pixel size of 8^′′ instead of the default 4^′′ pixel size to increase the signal-to-noise ratio (SNR) of the final Stokes $I$, $Q$ and $U$ maps, but used masks defined using the 4^′′-pixel-size Stokes $I$ map from step 1, in order to provide sufficient constraints on the mapmaker. A detailed investigation into the effect of pixel size on POL-2 data reduction and the choice to use 8^′′ pixels is presented by Karoly et al. (in prep.).

The native units of the reduced Stokes $I$, $Q$ and $U$ maps are pW. We converted the maps to Jy beam^-1 using a Flux Conversion Factor (FCF) of 495 Jy beam^-1 pW^-1 (Mairs et al., 2021), multiplied by a factor of 1.35 to account for additional flux losses in POL-2 (Friberg et al., 2016). As we used 8^′′ rather than 4^′′ pixels, we further multiplied the maps by a factor of 1.12, determined from SCUBA-2 calibration data (Karoly et al. in prep.). Figure 2 shows an image of the 850 $\mu$m total intensity (Stokes $I$) map of the Horsehead Nebula in units of Jy beam^-1, on which polarization segments are overplotted.

We binned the polarization segments obtained using the data reduction processes described above to a pixel size of 12^′′, which is the JCMT 850 $\mu$m primary beam size, and is similar to the JCMT 850 $\mu$m effective beam size of $14\farcs 1$. The polarization fractions $p$ listed in the polarization half-vector catalogues are debiased, and are given by

where $\delta Q$ and $\delta U$ are the square roots of the measured variances on $Q$ and $U$ respectively. The polarization angles $\theta$ are given by

Uncertainties on $p$ and $\theta$ are given by

and

respectively. We rotated these polarization angles by 90 degrees to show magnetic field orientations in the Horsehead Nebula in Figure 2 and Figure 2, and in subsequent analysis of magnetic field direction. The maps and half-vector catalogs used in this work are available at [DOI to be inserted in proof].

Three CO isotopologues, ¹²CO ($J=3-2$), ¹³CO ($J=3-2$) and C¹⁸O ($J=3-2$), were observed in the Horsehead Nebula under project code M18BD002 using HARP. The ¹³CO and C¹⁸O observations were made simulatenously on 2018 December 30. All observations of the CO isotopologues were carried out in weather band 2 (0.05 $<$ $\tau_{225}$ $\leq$ 0.08). In this work, we use only the C¹⁸O data, which like the other isotopologues has a critical density of $\sim 10^{4}$ cm^-3, and which, with the lowest fractional abundance of the three isotopologues, is able to trace gas density to the highest optical depth (e.g. Rigby et al., 2016). With a rest frequency of 329.331 GHz $\approx 911$ $\mu$m, the C¹⁸O data has comparable resolution to the SCUBA-2 850 $\mu$m data.

To estimate velocity dispersion values in the Horsehead Nebula, we reduced the C¹⁸O data using the ORAC Data Reduction (ORAC-DR) pipeline and the Kernel Application Package (KAPPA; Currie et al. 2008) in the Starlink software suite (Jenness et al., 2013). We obtained a reduced C¹⁸O data cube in which the pixel size, beam size and spectral resolution are 7^′′, 15.3^′′ and 0.05 km s^-1, respectively. To improve the SNR, we smoothed the data cube to a spectral resolution of 0.15 km s^-1. Figure 3 shows the velocity dispersion of the C¹⁸O data, which was obtained by fitting the C¹⁸O spectra with a single Gaussian profile. The contours in the figure show total intensities at 850 $\mu$m of 9, 100 and 150 mJy beam^-1. We analyzed the two regions within the lowest contour level, in which velocity dispersions of C¹⁸O are similar to each other.

Both SMM1 and SMM2 have ordered magnetic fields, but the fields in the two regions have very different magnetic field morphologies. We select those polarization segments for which $p$/$\delta p\geq 2$ and $p<20\%$ for analysis. Figure 4 shows magnetic field orientations in the Horsehead Nebula as a polar bar chart. The left-hand side of the plot shows magnetic field orientations in SMM1 and SMM2 obtained by the JCMT, while the right-hand side shows magnetic field orientations in the vicinity of the Horsehead Nebula obtained by Palomar and Planck, such that diametrically opposite angles on the plot agree. The magnetic field orientations in SMM2 have a peak at 115 $\pm$ 30 degrees, approximately perpendicular to the major axis of the source, and are broadly aligned with those measured with Planck, and similar to those measured at the Palomar Observatory. The magnetic field orientations in SMM1 are roughly perpendicular to the curved structure of SMM1 in the plane of the sky. The magnetic field orientations in SMM1 thus occupy a wide range of angles in Figure 4.

Pattle et al. (2018) presented magnetic field orientations within the photoionized pillars known as the ‘Pillars of Creation’ in M16 at a distance of $\sim 1.8$ kpc. They found that magnetic fields are aligned along the length of the pillars, but were unable to resolve the magnetic fields in the PDRs at the pillars’ tips. Our observations of the Horsehead nebula are a factor $\sim 5$ higher in linear resolution than the M16 observations, and the PDR of the Horsehead appears to be more extended in the plane of the sky than those of the pillars in M16, allowing us to resolve the magnetic field within the PDR. The magnetic field in SMM1 is perpendicular to the major axis of the PDR structure. This is quite dissimilar to POL-2 observations of the magnetic field in the Orion Bar in OMC-1, at a comparable distance to the Horsehead, in which the magnetic field runs along the length of the PDR (Ward-Thompson et al., 2017; Pattle et al., 2017). A possible cause of the morphology which we observe in SMM1 is that we are seeing a magnetic field in the PDR which is folded back on itself along the line of sight. A detailed schematic view of the magnetic field morphology of SMM1 is presented in Section 4.2.

We calculated the column density of molecular hydrogen in SMM1 and SMM2 using (Hildebrand, 1983)

where $I_{\nu}$ is the total intensity at frequency $\nu$, $\mu$ is the mean molecular weight, $m_{\text{H}}$ is the mass of a hydrogen atom, $\kappa(\nu)$ is the dust opacity and $B_{\nu}(T)$ is the Planck function for a dust temperature $T$. We took $\mu$=2.8 (Kauffmann et al., 2008). The dust opacity is estimated using $\kappa(\nu)=\kappa_{\nu_{0}}(\nu/\nu_{0})^{\beta}$ where $\kappa_{\nu_{0}}$ is the dust opacity at the reference frequency $\nu_{0}$, and $\beta$ is the dust opacity spectral index. We took $\kappa_{\nu_{0}}$ = 0.1 cm² g^-1 at $\nu_{0}$ = 1000 GHz, and $\beta$ = 2, thus assuming a dust-to-gas mass ratio of 1:100 (Beckwith et al. 1990; Motte & André 2001; André et al. 2010).

The dust temperatures in SMM1 and SMM2 are 22 and 15 K, respectively (Ward-Thompson et al., 2006). We estimated mean column densities for SMM1 and SMM2 by substituting the mean Stokes I intensities for SMM1 and SMM2 of 71 and 69 mJy beam^-1 at 850 $\mu$m, as measured from our observations, into equation (5). The mean Stokes $I$ intensities are estimated by averaging intensities within the orange contours in Figure 2, in which the left contour outlines SMM2 and the right contour outlines SMM1. The mean column densities estimated in SMM1 and SMM2 are 3.9 $\times 10^{21}$ and 6.7 $\times 10^{21}$ cm^-2. The peak column densities of SMM1 and SMM2 are 1.3 $\times 10^{21}$ and 3.0 $\times 10^{22}$ cm^-2, which are comparable to within a factor of a few with those estimated by Ward-Thompson et al. (2006).

We estimated magnetic field strengths in the Horsehead Nebula using the angular dispersion function (Hildebrand et al., 2009) implementation of the DCF method (Davis, 1951; Chandrasekhar & Fermi, 1953). The DCF method is widely used to estimate magnetic field strengths in star-forming regions from dust polarization observations (e.g. Pattle et al., 2022, and refs. therein). DCF estimates magnetic field strength from the gas density, velocity dispersion and polarization angle dispersion in a given region. The fundamental assumption of DCF is that the distortion of magnetic field lines by small-scale turbulent gas motions is indicative of the Alfvén Mach number of the gas, and can be estimated as a polarization angle dispersion in the region. In traditional implementations of DCF, polarization angle dispersion is measured under the assumption of a uniform underlying magnetic field direction. However, magnetic field lines can show ordered variation and curvature due to gravitational collapse, rotation, shocks or outflows. In order to estimate magnetic field strengths using traditional DCF, it is thus essential to assume a model for the structure of the mean magnetic field.

Hildebrand et al. (2009) suggested a method to separate the ordered and turbulent components of the magnetic field, a structure-function based approach to DCF known as the angular dispersion function. To calculate the angular dispersion function, we calculate the differences in polarization angles between individual pairs of polarization segments, $\Delta\theta(l)\equiv\theta(x)-\theta(x+l)$, where $\theta(x)$ is the polarization angle of a segment at position $x$, and $\theta(x+l)$ is the polarization angle of a segment separated from $x$ by a distance $l$. If the number of pairs is given by $N(l)$, the angular dispersion function is then given by

The angular dispersion function can, at small $l$, be fitted with the quadratic function

where $ml$ is the mean magnetic field (referred to as the large-scale field by Hildebrand et al. (2009)), $b$ is the turbulent dispersion about the mean magnetic field and $\sigma_{M}(l)$ is the measurement uncertainty. Hildebrand et al. (2009) expressed polarization angle dispersion as the ratio of the turbulent to large-scale magnetic field strength, which is given $b/\sqrt{2-b^{2}}$. When the turbulent component of the magnetic field is much smaller than the ordered component, i.e. $b\ll 1$ rad, $b/\sqrt{2-b^{2}}\to b/\sqrt{2}$. In this case, the DCF equation for magnetic strength can be expressed as

where $B_{pos}$ is the magnetic field strength in the plane of the sky, $\rho$ is the mass density of gas coupled to the magnetic field and $\sigma_{v}$ is the velocity dispersion of the gas. The mass density is given by $\rho=\mu m_{\text{H}}n({\text{H}}_{2})$, where $\mu$ is the mean molecular weight, taken to be 2.8 (Crutcher, 2004a; Kauffmann et al., 2008), $m_{\text{H}}$ is the mass of a hydrogen atom, and $n({\text{H}}_{2})$ is the volume number density of the gas.

We estimated each of the parameters in equation (8) in order to measure magnetic field strengths in SMM1 and SMM2. We fitted the quadratic model of the angular dispersion function, as given in equation (7), to the dispersion of polarization angle differences in SMM1 and SMM2 as a function of the distance between pairs of polarization segments, as shown in Figure 5. We fitted the data in the range $l<36^{\prime\prime}$. The fitted $b$ values are $21.2^{\circ}\pm 3.4^{\circ}$ and $15.8^{\circ}\pm 2.6^{\circ}$ in SMM1 and SMM2, respectively. Note that the units of $b$ in equation (8) are radians.

We estimated volume density values from the column density values obtained in Section 3. We assumed effective radii for SMM1 and SMM2 of 0.15 and 0.10 pc, which we took to be the radii of circles with areas equal to the areas within the orange contours shown in Figure 2. The volume density is estimated using

where $n(\text{H}_{2})$ is the volume density and $r$ is the effective radius. We thus estimate mean volume densities in SMM1 and SMM2 of 6.4$\times$10³ and 1.7$\times$10⁴ cm^-3.

We estimated non-thermal velocity dispersion values from our C¹⁸O spectral line data. We fitted each pixel of the observed line data with a single Gaussian profile. The non-thermal velocity dispersion is given by

where $\sigma_{obs}$ is the estimated dispersion of the Gaussian profile (Figure 3), $k$ is the Boltzmann constant, $T_{k}$ is the kinetic temperature and $m_{\text{C}^{18}\text{O}}$ is the mass of C¹⁸O molecule. We took $T_{k}$ to be equal to the dust temperatures of 22 and 15 K in SMM1 and SMM2, respectively. The mean observed velocity dispersions in SMM1 and SMM2 are 0.11 and 0.10 km s^-1. These dispersions are dominated by non-thermal motions; the non-thermal velocity dispersions given by equation (10) agree with these values to two decimal places.

By substituting the above values into equation (8), we estimated magnetic field strengths of 56$\pm$9 and 129$\pm$21 $\mu$G in SMM1 and SMM2, respectively. The uncertainties on these magnetic field strengths were obtained by propagating the fitting uncertainties on $b$ for the two sources.

The relative energetic importance of magnetic fields in the Horsehead Nebula can be estimated using the mass-to-flux ratio and Alfvén Mach number to compare magnetic fields to self-gravity and internal turbulent motions. We derive mass-to-flux ratios and Alfvén Mach numbers for both SMM1 and SMM2. Additionally, we estimate the magnetic, gravitational and turbulent energies in SMM2, and consider their relative importance. From these results, we discuss the role of magnetic fields in the evolution of the Horsehead Nebula.

Broadly, both simulations (e.g. Arthur et al., 2011) and previous observations suggest that magnetic fields in molecular gas typically run plane-parallel to the interface with an HII region, whether in a PDR (Ward-Thompson et al., 2017), or in columns sheltered behind PDRs (Pattle et al., 2018). In the heads of pillars, this is predicted to result in ‘hairpin’ magnetic field structures (e.g. Arthur et al., 2011), in which the pre-shock magnetic field is bent back on itself as the pillar forms from its parent molecular cloud.

The SMM1 PDR appears to be a relatively extended, bar-like PDR sitting at the head of a broad and short pillar. The plane-of-sky magnetic field morphology which we observe is broadly perpendicular to the projected major axis of the PDR. Such a projected geometry could result from a plane-parallel magnetic field if either (1) we observe hairpin magnetic field lines, with the pre-shock magnetic field lines having been folded back on themselves along the line of sight, and/or (2) the Horsehead pillar is slightly inclined with respect to the plane of the sky, and we observe in projection magnetic fields which are running parallel to the plane of the SMM1 PDR, which is itself extended along the line of sight. Although we appear to be observing the $\sigma$ Ori-Horsehead interaction nearly edge-on (Abergel et al., 2003), $\sigma$ Ori may be located slightly behind the Horsehead (Pound et al., 2003), which may lend weight to the latter of these scenarios. In both of these cases, we would expect the average orientation of the magnetic field to be inherited from that in the pre-shock molecular cloud, and so it is not surprising that the average magnetic field in the densest parts of SMM1 is similar to that in the apparently less disturbed SMM2.

Figure 6 shows a schematic view of the Horsehead Nebula which illustrates our proposed magnetic field geometry. The gray lines show magnetic field lines in the HII region, which likely trace the initial (pre-shock) magnetic field direction, while the magnetic field lines in SMM1 and SMM2 are shown as red lines. The O star, $\sigma$ Ori, is located to the southwest of the pillar. The observer views the region either edge-on, or slightly inclined along the line of sight. The average field direction in the HII region is 161 degrees east of north, the mean of the magnetic field orientation angles obtained from Palomar observations. In SMM1, our proposed hairpin magnetic field geometry is shown: the field lines are aligned along the pillar, and folded or curved back on themselves in the PDR at the head of the pillar (e.g., Arthur et al. 2011; Pattle et al. 2018). The magnetic field in SMM2 is linear, with a mean field direction of 115 degrees east of north, as estimated from JCMT observations. The field is not aligned with the major axis of the pillar structure, and differs by about 45 degrees from that in the HII region. As discussed above, we hypothesize that SMM2 is embedded within the pillar and has inherited its magnetic field from the pre-shock molecular cloud, although we cannot rule out the scenario in which the core has been forced into collapse by photon pressure or gas flows due to the formation of the pillar.