Mapping the magnetic field in the Taurus/B211 filamentary cloud with SOFIA HAWC + and comparing with simulation

ABSTRACT

1 INTRODUCTION

Filamentary structures have been found at almost all size scales in the Galaxy. Massive, long filamentary dark clouds are commonly found inside giant molecular clouds (GMCs; e.g. Bergin & Tafalla 2007; André et al. 2014, and references therein), such as the dark clouds L1495 in the Taurus cloud complex (e.g. Chapman et al. 2011) and the Serpens South cloud in the Serpens region (e.g. Dhabal et al. 2018). Filamentary clouds of 4–6 pc length are common, and possibly longer than 10 pc. Some of these clouds are dark at infrared wavelengths. The line–width size relation observed for molecular gas indicates that the thermal Mach number would exceed 10 at such size scales. The long-term survival of these filamentary structures requires a reinforcing mechanism. As shown in the ideal magnetohydrodynamical (MHD) simulations of Li & Klein (2019), a moderately strong, large-scale magnetic field (Alfv|$\acute{\rm e}$|n Mach number, |${{\cal M}_{\rm A}}\sim 1$|⁠) can provide such a mechanism. In the weak-field model with |${{\cal M}_{\rm A}}=10$|⁠, the appearance of molecular clouds is clumpy, rather than the long and slender filamentary clouds seen in moderately strong field models. High-resolution images of massive molecular clouds from the Herschel space telescope reveal complex filamentary substructures (e.g. André et al. 2014). The characteristic inner width of molecular filaments found with Herschel is about ∼0.1 pc (Arzoumanian 2011; Arzoumanian et al. 2019). Dense cores, where stars form, are located along or at the intersections of some of these fine substructures (e.g. Könyves et al. 2015; Tafalla & Hacar 2015). From these observations of molecular cloud structures at different size scales, one can visualize an evolutionary sequence of star formation starting from highly supersonic, magnetized GMCs, continuing on to filamentary dark clouds that form within them, and then on to finer filamentary substructures. Fragmentation of these filamentary structures and substructures leads to the clumps and dense cores that form protostellar clusters and protostars. Knowing the physical conditions inside filamentary clouds would provide crucial information on the formation of filamentary substructures and dense cores, and on the origin of the initial mass function (IMF) and the star formation rate. Particularly important is the characterization of the physical properties of transcritical filamentary structures whose mass per unit length is within a factor of ∼2 of the critical line mass |${M_{\rm crit,\, th,\, \ell }}=2\, c_{\rm s}^2/G$| of nearly isothermal cylindrical filaments (e.g. Ostriker 1964; Inutsuka & Miyama 1997), where c_s is the isothermal sound speed. Indeed, Herschel observations suggest that transcritical filamentary structures dominate the mass function of star-forming filaments and that their fragmentation may set the peak of the prestellar core mass function and perhaps ultimately the peak of the IMF (André et al. 2019). In this paper, we report the results of polarimetric observations of the pristine section B211 of one such transcritical filament, the Taurus B211/B213 filament, using the High-resolution Airborne Wideband Campera plus (HAWC+) onboard Stratospheric Observatory For Infrared Astronomy (SOFIA). We determine the magnetic field structure inside a filamentary cloud with filamentary substructures.

The filamentary cloud L1495 is located in the Taurus molecular cloud at a distance of about 140 pc (Elias 1978). Using their H-band polarization observation and the Davis–Chandrasekhar–Fermi (DCF) method (Davis 1951; Chandrasekhar & Fermi 1953), Chapman et al. (2011) estimated the plane-of-sky (POS) magnetic field strength to be 10–17 μG in the low density regions near the L1495 cloud, including the B211 region, and 25–28 μG inside the cloud. From observations of ¹²CO and ¹³CO, they find that the velocity dispersion is 0.85–1.16 km s⁻¹. Their observations have a resolution of 0.135 pc. The mean surface density is about N(H₂) ∼ 1.45 × 10²² cm⁻² (Palmeirim et al. 2013). Using the density estimated by Hacar et al. (2013) and the measured velocity dispersion ∼1 km s⁻¹ cited above, the Alfv|$\acute{\rm e}$|n Mach number of the long filamentary cloud is about 2.7. Combining the polarization observations of Heiles (2000), Heyer et al. (2008), and Chapman et al. (2011), Palmeirim et al. (2013) found that the large-scale mean field direction is almost orthogonal to the cloud axis in B211/B213 and roughly parallel to faint striations seen in both CO and Herschel data. There is also some kinematic evidence that the B211/B213 filament is embedded in a sheet- or shell-like ambient cloud and in the process of accreting mass from this ambient cloud (Shimajiri et al. 2019), perhaps through the magnetically aligned striations. Is it possible that the magnetic field pierces straight through the cloud? If so, then the picture of the formation of filamentary clouds is simple: gas is simply gathered into the cloud along approximately straight field lines.

A portion of the Taurus/B213 filament was recently mapped with JCMT-POL2 as part of the BISTRO project (Eswaraiah et al. 2021), but the corresponding |$850\, \mu$|m polarization data only constrained the magnetic field towards the dense cores within the filament. Other recent high-resolution polarization observations of magnetic field structures within more massive molecular clouds, such as Vela C (Soler et al. 2013, 2017; Dall’Olio et al. 2019) and M17 SWex (Sugitani et al. 2019), show that magnetic fields inside dense filamentary clouds with complex substructures are not simple. Magnetic fields inside clouds can have large deviations from the large-scale field orientation. Within the ∼0.8 pc outer diameter measured on Herschel data, the B211/B213 filament system has a characteristic half-power width of ∼0.1 pc and exhibits complex filamentary substructures (Hacar et al. 2013; Palmeirim et al. 2013). Hacar et al. (2013) found that the B211 region has a mass of 138 M_⊙ and is roughly 2 pc long. From their C¹⁸O intensity map, the width of the B211 region encompassed by the lowest C¹⁸O contour is about 0.3 pc. Hacar et al. (2013) identified multiple velocity components in C¹⁸O at different locations in B211–B213 with separations as large as about 2 km s⁻¹. Tafalla & Hacar (2015) find that the relative velocities between filamentary substructures in the filamentary cloud range over 2.2 km s⁻¹, possibly implying that the substructures are converging at high velocity. B211 is very bright in both C¹⁸O and SO and has intense dust millimeter emission. The gas in B211 has an unusually young chemical composition and lack of young stellar objects, indicating that this region is at a very early state of evolution (Hacar et al. 2013). This region is therefore particularly suitable for the study of magnetic field structures in filamentary clouds without any confusion from protostellar activity.

In the high-resolution infrared dark cloud (IRDC) simulation by Li & Klein (2019) using the adaptive mesh refinement code ORION2, a long filamentary cloud is formed in a moderately strong magnetic field environment, even though the thermal Mach number was 10. The long filamentary cloud created in the simulation has filamentary substructures similar to those in L1495. The simulation may therefore provide unique information on the physical environment inside filamentary clouds and on how they form. In the simulation, the large-scale magnetic field is approximately perpendicular to the cloud axis, similar to the case in L1495. However, the small-scale magnetic field inside the simulated cloud, which has a width similar to that of B211, has a chaotic structure. Until now, there has never been a polarimetric observation with a resolution and a sensitivity high enough to peer into a filamentary cloud with no star formation. This motivates us to map a portion of L1495 to determine the field morphology inside the cloud and thereby gain an understanding of the physical environment inside such a cloud.

We report in this paper our observations of the filamentary cloud L1495/B211 using the recently optimized HAWC + polarimeter on SOFIA to probe the complex magnetic field inside a slender filamentary cloud with complex filamentary substructures. The observation using HAWC + polarimeter, the data reduction method, and the results are presented in Section 2. In Section 3, we investigate the physical state of the B211 region from observation. Using the DCF method, we estimate the magnetic field strength in Section 3.1 with the aid of recent C¹⁸O (1-0) line emission data from the Institut de Radioastronomie Millimétrique (IRAM) 30-m telescope. In Section 3.2, we study the relation between the inferred magnetic field from HAWC + observation and the surface density contours. In Section 4, results of our numerical simulation of filamentary clouds are used to provide insights into the physical state of B211. Our conclusions are presented in Section 5.

2 OBSERVATIONS

2.1 SOFIA HAWC + mapping observations and data reduction methods

L1495 was observed (ID: 07_0017, PI: Li, P.S.) at 214 |$\mu$|m (Δλ = 44 |$\mu$|m, full width at half maximum, FWHM) using HAWC + (Vaillancourt et al. 2007; Dowell et al. 2010; Harper et al. 2018) on the 2.7-m SOFIA telescope. HAWC + polarimetric observations simultaneously measure two orthogonal components of linear polarization arranged in two arrays of 32 × 40 pixels each, with a detector pixel scale of 9|${_{.}^{\prime\prime}}$|37 pixel⁻¹ and beam size (FWHM) of 18|${_{.}^{\prime\prime}}$|2 at 214 |$\mu$|m. At 214 |$\mu$|m, HAWC + suffers of vignetting where five columns cannot be used for scientific analysis (Harper et al. 2018), therefore the field of view (FOV) of the polarimetric mode is 4.2 × 6.2 sqarcmin. We performed observations using the on-the-fly-map (OTFMAP) polarimetric mode. This technique is an experimental observing mode performed during SOFIA Cycle 7 observations as part of the shared-risk time to optimize the polarimetric observations of HAWC+. Although we will focus on the scientific results of L1495, we here describe the high-level observational steps used in these these observations, where Sections 2.2 and 2.3 describe the details of the OTFMAP polarimetric mode.

We performed OTFMAP polarimetric observations in a sequence of four Lissajous scans, where each scan has a different halfwave plate (HWP) position angle in the following sequence: 5°, 50°, 27.5°, and 72.5°. This sequence is called ‘set’ hereafter (Table 1 – column 9). In this new HAWC + observing mode, the telescope is driven to follow a parametric curve at a non-repeating period whose shape is characterized by the relative phases and frequency of the motion. Each scan is characterized by the scan angle, scan amplitude, scan rate, scan phase, and scan duration. The scan angle is the relative angle of the cross-elevation direction of the FOV of the scan with respect to north, where 0° is North and positive increase is in the east of north direction (Table 1 – column 7). An example of the OTFMAP for total intensity observations of NGC 1068 using HAWC + /SOFIA is shown by Lopez-Rodriguez et al. (2018, fig. 1). The OTFMAP polarimetic mode using HAWC + /SOFIA at 89 |$\mu$|m has recently been successfully applied to the Galaxy Centaurus A (Lopez-Rodriguez 2021). A summary of the observations at 214 |$\mu$|m are shown in Table 1. We performed square scans (Table 1 – column 8) at three different positions as shown in Table 1 (columns 5 and 6). After combining all scans, the full FOV is 20 × 20 sqarcmin. Although Table 1 lists all performed observations for this program with a total executed time of 6.37 h, only a final executed time of 4.73 h (with a total on-source time of 4.40 h) was used for scientific analysis. The removed sets listed in Table 1 were not used due to loss of tracking during observations.

We reduced the data using the Comprehensive Reduction Utility for SHARP II v.2.42-1 (crush; Kovács et al. 2006, 2008) and the HAWC_DRP_v2.3.2 pipeline developed by the data-reduction pipeline group at the SOFIA Science Center. Each scan was reduced by crush, which estimates and removes the correlated atmospheric and instrumental signals, solves for the relative detector gains, and determines the noise weighting of the time streams in an iterated pipeline scheme. Each reduced scan produces two images associated with each array. Both images are orthogonal components of linear polarization at a given HWP position angle. We estimated the Stokes |$I,\, Q,\, U$| parameters using the double difference method in the same manner as the standard chop-nod observations carried by HAWC + described in section 3.2 by Harper et al. (2018). The degree (P) and position angle of polarization were corrected by instrumental polarization (IP) estimated using OTFMAP polarization observations of planets. We estimated an IP of Q/I = −1.0 per cent and U/I = −1.4 per cent at 214 |$\mu$|m, respectively, with an estimated uncertainty of |$\sim 0.8{{\ \rm per\ cent}}$|⁠. The IP using OTFMAP observations are in agreement with the estimated IP using chop-nod observations. To ensure the correction of the position angle of polarization of the instrument with respect to the sky, we took each set with a fixed line of sight (LOS) of the telescope. For each set, we rotated the Stokes Q and U from the instrument to the sky coordinates. The polarization fraction was debiased (Wardle & Kronberg 1974) and corrected by a polarization efficiency of 97.8 per cent at 214 |$\mu$|m. The final Stokes I and its associated errors were calculated and downsampled to the beam size (18|${_{.}^{\prime\prime}}$|20). The final Stokes |$Q,\, U,\, P$|⁠, position angle, polarized intensity (PI), and their associated errors were calculated and re-sampled to a super-pixel of 3 × 3 detector pixel size, which corresponds to a re-sampled pixel size of 28|${_{.}^{\prime\prime}}$|1 (or 0.019 pc at the distance of the cloud). This super-pixel was chosen to optimize the signal-to-noise ratio (SNR), obtain statistically independent measurements and significant polarization measurements without compromising spatial resolution for the data analysis.

2.2 HAWC + OTFMAP polarization: advantages and limitations

Several advantages and limitations are found with the OTFMAP polarization mode. The advantages are the reduction of overheads and radiative offsets when compared with the chop-nod technique. The overheads of the OTFMAP are estimated to be 1.1 in comparison with the typical overheads of 2.7 by the chop-nod technique, which shows an improvement by a factor ≥2. This improvement is due to OTFMAP constantly integrating with the source on the FOV while covering off-source regions to estimate the background levels, and observing overheads. For the OTFMAP method, the telescope is always on-axis, without chopping the secondary mirror as it is in the chop-nod technique. Thus, the radiative offset is not present and the sensitivity of the observations was estimated to improve by a factor of 1.6. The OTFMAP technique provides a larger map area when compared to the chop-nod technique. Our observations were taken at three different positions covering an FOV of 10 × 10 and 11 × 11 sqarcmin, yield a final FOV of 20 × 20 sqarcmin. Note the advantage of the large FOV by the OTFMAP when compared with the single 4.2 × 6.2 sqarcmin by the chop-nod technique.

The limitation of the OTFMAP technique lies in the recovering of large-scale diffuse and faint emission from the astrophysical objects. This is a result of the finite size of the array, variable atmosphere conditions, variable detector temperatures, and the applied filters in the reduction steps to recover extended emission. We applied several filters using crush to recover large-scale emission structures of L1495 while paying close attention to any change that may compromise the intrinsic polarization pattern of the astrophysical object. We conclude that the faint filter with a number of 30 iterations from crush can recover large-scale emission structures larger than the Band E FOV from our observations of L1495. The faint option applies filtering to the timestreams and extended structures to recover fluxes with SNR <10 in a single scan. In addition, the number of rounds are such as that the iterative pipeline is able to recover large-scale structures without introducing additional artificial structures not identified in the Herschel images. In general, the noise increases as a function of the length, L, of the extended emission as ∼L². We force each individual scan produced by crush to have a pixel scale of 3 × 3 detector pixels (28|${_{.}^{\prime\prime}}$|1), which increase the SNR of each scan by a factor of 3 helping to recover larger and fainter structures.

2.3 HAWC + OTFMAP polarization: zero-level background

An important step is the estimation of the zero-level background of the observations. We remind the reader that HAWC + measures the power of the emissive and variable atmosphere and the astrophysical object. The data reduction scheme described above produces regions of negative fluxes in areas of extended and low surface brightness due to the similar levels of noise and astrophysical signal. Thus, it is of great importance to characterize and estimate the zero-level background across the observations of L1495, because there is a potential loss of flux that requires to be estimated and added to the full image.

We have determined and corrected the zero-level background of our observations as follows. Using Herschel images at 160 and 250 |$\mu$|m from the Herschel Archive,¹ we identified a region in the sky where the fluxes of an individual pixel of size 28|${_{.}^{\prime\prime}}$|1 are below the sensitivity of HAWC + at 214 |$\mu$|m. This area is shown in Fig. 1, which is located in a common region for all scans across the multiple flights. The size of this region is chosen to be equal to the HAWC + FOV at Band E, i.e. 4.2 × 6.2 sqarcmin. The size of the background region was chosen to be the same as if the observations were performed using the chop-nod observing mode. The size of this region does not influence the estimation of the zero-background level, rather the location and the surface brightness do. Then, for both arrays and each HWP position angle produced from the first step by crush, we estimate the mean and standard deviation within the zero-level region. To remove negative values across the image, the mean is added to all pixels in each scan and HWP position angle. After this step, the same reduction procedures as described in Section 2.1 are followed. Using archival chop-nod and OTFMAP observations of well-known objects, e.g. 30 Doradus and OMC-1, and applied the same methodology, we reached similar conclusions and methodologies, while the polarization pattern was shown to be consistent between reduction schemes. Finally, we computed the SED of the source using 70–500 |$\mu$|m Herschel images and estimated the expected flux at 214 |$\mu$|m. We estimated that the fluxes from the zero-level background corrected image are within 8 per cent from the expected flux from the Herschel SED, which is within the flux calibration uncertainty of HAWC + of ≤15 per cent provided by the SOFIA Science Center.

Figure 1.

Total surface brightness at 214 |$\mu$|m of L1495 within a 1200 × 1200 sqarcsec region using the OTFMAP observations. Contours start at 45σ and increase in steps of 5σ, with σ = 0.032 mJy sqarcsec⁻¹. Although quality cuts have been performed to cut edge effects, there are still structural artefacts at the corners of the map (specially in the West region) due to the low coverage at the edges of the map. Polarization measurements (black lines) have been rotated by 90° to show the inferred magnetic field morphology. The length of the polarization vector is proportional to the degree of polarization. Only vectors with P/σ_P ≥ 2 are shown. A legend with a |$5{{\ \rm per\ cent}}$| polarization and the beam size (18|${_{.}^{\prime\prime}}$|2) are shown. The zero-level background region (white dashed line) described in Section 2.3 is shown.

Open in new tabDownload slide

Although the zero-level background region has low surface brightness, the polarization may be high and contaminate the astrophysical signal after the zero-level background correction. Here, we estimate the contribution of the zero-level background to the polarization measurements. As mentioned above, the mean and standard deviation within the zero-level background region was estimated for each array and HWP position angle. Using the double difference method (section 3.2 by Harper et al. 2018), we estimate the Stokes Q and U and their uncertainties by spatial averaging within full FOV of the zero-level background region. Then, the Stokes Q and U were corrected by instrumental polarization, and P and position angle were estimated and corrected by bias and polarization efficiency. Finally, the P and position angle of the zero-level background region were estimated to be 8.5 ± 3.5 per cent and 42 ± 8°, respectively. The minimum detectable flux from Stokes I is estimated to be 3σ_I = 0.096 mJy sqarcsec⁻¹, which corresponds to a polarized flux of 3σ_PI = 0.008 mJy sqarcsec⁻¹ using P = 8.5 per cent. From our polarization measurements with P/σ_P ≥ 2, we estimate a median polarized flux of 0.055 mJy sqarcsec⁻¹. Thus, the zero-level background correction contributes a median of ∼14 per cent to the polarized flux in our measurements.

2.4 HAWC + polarization map and orientation of magnetic fields

Fig. 1 shows polarization measurements projected on to the total surface brightness at 214 |$\mu$|m of the 1200 × 1200 sqarcsec region of L1495/B211 that we observed. The polarization measurements have been rotated by 90° to show the inferred magnetic field morphology. All polarization angles (PAs) cited in this paper have been rotated in this manner. Only polarization measurements with P/σ_P ≥ 2 are shown (Wardle & Kronberg 1974). The length of the polarization lines are proportional to the degree of polarization, where a 5 per cent polarization measurement is shown as reference. Images had edges artefacts due to the sharp changes in fluxes and limited number of pixels. The final 214 |$\mu$|m HAWC + polarization measurements contain pixels with the following quality cuts: (1) pixels with a scan coverage ≥30 per cent of all observations, (2) pixels which Stokes I measurements have an uncertainty ≥2.5σ_I, where σ_I is the minimum uncertainty in Stokes I, (3) pixels with a surface brightness of ≥1 mJy sqarcsec⁻¹, (4) pixels with P ≤ 30 per cent given by the maximum polarized emission found by Planck observations (Planck Collaboration XII 2013), and (5) polarization measurements with P/σ_P ≥ 2. We find that 14 per cent (40 out of 282) of the measurements are within 2 ≤ P/σ_P ≤ 3.

In Fig. 2, we overplotted the magnetic field orientations from near-infrared H-band observations obtained by Chapman et al. (2011). Note that the inferred magnetic field from the H band arises from dichroic absorption, while our 214-|$\mu$|m measurements arise from dichroic emission. We detect many multiple structures of magnetic field over just 0.82-pc region inside the 2-pc-long B211 region from the HAWC + polarization mapping at smaller scales 28|${_{.}^{\prime\prime}}$|1 compared to the lower resolution observation from the H band and Planck observations (see Fig. 5). We note that the magnetic field of the lower half of the observed area is more uniform and close to the perpendicular direction of the filamentary cloud. From the Herschel intensity map, this part of B211 appears to have two filamentary substructures crossing each other in an x-shape appearance near RA of 4^h18^m20^s and Dec. of 27|${_{.}^{\circ}}$|29. The two structures may be spatially nearby and appeared to be overlapped along the LOS. The magnetic field could be a combined result of the overlapping projection. In the upper half, the magnetic field structure appears as a highly non-uniform chaotic state. It is consistent with the turbulent appearance of the underlining intensity map that there could be three tangling filamentary substructures at this location as identified by Hacar et al. (2013). The deviation of polarization angle is large from the uniform large-scale field direction indicated by the near-infrared H band and Planck observations. Fig. 3 is a line integral convolution (Cabral & Leedom 1993) plot of the inferred magnetic field from the HAWC + polarization observations.

Figure 2.

The 250 |$\mu$|m total surface brightness from Herschel/SPIRE (colour scale) with the magnetic field morphology as inferred from the 214 |$\mu$|m HAWC + (black lines) and H-band (grey lines; Chapman et al. 2011) observations. Contours start at 0.5 mJy sqarcsec⁻¹ and increase in steps of 0.25 mJy sqarcsec⁻¹.

Open in new tabDownload slide

Figure 3.

The inferred magnetic field orientation from HAWC + polarization observations at 214 |$\mu$|m of the B211 filament is shown using the linear integral convolution algorithm (LIC; Cabral & Leedom 1993). Same polarization measurements as Fig. 1, a resample scale of 20, and a contrast of 4 were used. The colourscale shows the |$250 \,\mu$|m total intensity image from the Herschel Gould Belt survey as shown in Fig. 2.

Open in new tabDownload slide

The histogram of the B-field PA distribution of all the polarization measurements detected with P/σ_P > 2 is shown in Fig. 4(a). The peak is near 30°, which is very close to the large-scale mean magnetic field orientation of 26° ± 18° and is nearly orthogonal to the L1495 filamentary cloud axis at 118° ± 20° (Palmeirim et al. 2013). If the PAs have a range approaching 180°, then the dispersion can depend on the choice of θ = 0° since PAs near +90° can be flipped to −90° by a change in the orientation of 0°. In this paper, we evaluate the dispersion in PAs by choosing 0° to be consistent with the minimum dispersion in PAs. We find that the PAs in Fig. 4 have a dispersion of 36.1°, indicating that the small-scale magnetic field is strongly perturbed inside the cloud compared with the large-scale field. Basically, the inferred B field are pointing at all directions in the northwestern side region of 0.82 pc in size. The large-scale field PA distribution from Planck is also plotted in Fig. 4(a) for direct comparison. Note that most of the magnetic field orientations from Planck are located far from the B211 region (see Fig. 5) and have a dispersion of 37°. The several Planck polarization orientations inside the observed B211 region (indicated by a black dash-line box in Fig. 5) are all inside the two bins between 24° and 57° at the peak of the Planck distribution and close to the peak of the mean magnetic field inside B211 from HAWC+. The resolution of the large-scale field from Planck is ∼0.4 pc, which is about 21 times the size of the super-pixel that we adopted for the HAWC + results. In Section 4.3, we use a numerical simulation to discuss how the resolution of a polarization map can affect the interpretation of magnetic field structures.

Figure 4.

(a) The histograms of all the inferred magnetic field orientations from HAWC + polarization observations of the B211 filament at resolution of 28.1 arcsec shown in Fig. 2 (empty) and from Planck polarization observations (green) of the entire Taurus/B211 region at resolution of 10 arcmin shown in Fig. 5. The six Planck polarization measurements within the FOV of our observations are inside the two bins 24°–57° marked by the two blue arrows. The bin width of 16.4° is chosen to be larger than the measurement uncertainty of the PA. The vertical dashed line indicates the orientation of the B211 filament at PA = −62° equivalent to +118° (Palmeirim et al. 2013). The angle differences of the peaks of the two sets of histograms from the PA of the B211 filament are shown. (b) The histograms of inferred magnetic field orientations from HAWC + polarization observations of sub-regions SR1 (yellow) and SR2 (empty).

Open in new tabDownload slide

Figure 5.

Left: Column density map from Herschel Gould Belt survey data at 18.2 arcsec resolution (André et al. 2010; Palmeirim et al. 2013), with magnetic field vectors derived from Planck polarization data at 10 arcmin resolution (Planck Collaboration XII 2013) displayed in orange/yellow (and spaced every 10 arcmin). Contours are N(H₂) = 3 × 10²¹, 6.7 × 10²¹, and 10²² cm⁻². Right: Blow-up of the left image in the area mapped with SOFIA/HAWC+. Yellow vectors are from Planck and are here spaced by half a beam (5 arcmin). Smaller black segments show the magnetic field vectors derived from HAWC+ at 28 arcsec resolution. The solid red circles mark positions where both significant HAWC+ polarization measurements and C¹⁸O line data from the IRAM 30-m telescope are available. The two sub-regions (SR1 and SR2) for which a DCF analysis has been carried out are marked by white dotted rectangles; the two components of SR1 (SR1a and SR1b) are outlined by green dashed contours.

Open in new tabDownload slide

2.5 IRAM 30-m observations

C¹⁸O(1–0) mapping observations of a portion of the B211 field imaged with HAWC+ were carried out with the Eight MIxer Receiver (EMIR) receiver on the IRAM-30-m telescope at Pico Veleta (Spain) in 2016 April, as part of another project (Palmeirim et al. in preparation). At 109.8 GHz, the 30-m telescope has a beam size of ∼23 arcsec (HPBW), a forward efficiency of 94 per cent, and a main beam efficiency of 78 per cent.² As backend, we used the VESPA autocorrelator providing a frequency resolution of 20 kHz, which corresponds to a velocity resolution of ∼0.05 km s⁻¹ at 110 GHz. The standard chopper wheel method³ was used to convert the observed signal to the antenna temperature |$T_{\rm A}^{*}$| in units of K, corrected for atmospheric attenuation. During the observations, the system noise temperatures ranged from ∼85 to ∼670 K. The telescope pointing was checked every hour and found to be better than ∼3 arcsec throughout the run.

3 PHYSICAL STATE OF B211 REGION INFERRED FROM OBSERVATIONS

3.1 Magnetic field strength in B211

We derive magnetic field strengths using the DCF method based on the observed velocity dispersion, surface density (which provides an estimate of the gas density), and the dispersion of polarization angles. We used IRAM 30 m C¹⁸O data to derive the velocity dispersion and a HAWC + polarization map to derive the dispersion in polarization angles. The density is based on Herschel column density data published in the literature. We describe the detailed methods below. The observed and derived parameters are summarized in Table 2.

3.1.1 Davis–Chandrasekhar–Fermi method

A key step in the DCF method is to infer the dispersion in the field, |$\sigma _{\delta B_\perp }$|⁠, from the dispersion in PAs, σ_θ. A complication is that the field angles (FAs) can range over −180° to +180°, whereas the PAs are restricted to the range −90° to +90°. As a result, the direction of the implied field depends on the choice of zero angle for the PAs. A field angle of 60° if 0° is vertical becomes −30°if the coordinate system is rotated 90° counterclockwise. Since the magnitude of the PA depends on the choice of coordinate system, it follows that the value of σ_θ does also. In this paper, we adopt the convention that we choose the coordinate system that minimizes the dispersion, σ_θ, as recommended by Padoan et al. (2001). This becomes relevant only if some of the field angles differ by more than 180°, which generally occurs only if σ_θ is not small.

3.1.2 IRAM 30 m C¹⁸O data and velocity dispersion

Thanks to their high sensitivity, the IRAM 30 m molecular line data highlight the kinematic complexity of the region mapped by HAWC+, with the presence of multiple velocity components. These multiple velocity components are consistent with the presence of filamentary substructures in this region as discussed by Hacar et al. (2013). The variety of observed C¹⁸O(1–0) spectra is illustrated in Fig. 6, which shows clear changes in the number of velocity components and in overall centroid velocity as a function of position within and around the B211 filament.

Figure 6.

C¹⁸O(1–0) integrated intensity map over all channel velocities from 4.5 to 7 km s⁻¹. The contours correspond to 30, 50, and 70 per cent of the maximum integrated intensity (3 K km s⁻¹). Markers (x) in green indicate positions where statistically significant polarization measurements were obtained with HAWC + . Representative C¹⁸O(1–0) spectra observed with the IRAM 30-m telescope at selected positions in the field are shown to the left, bottom, and right of the map.

Open in new tabDownload slide

C¹⁸O(1–0) molecular line data trace the kinematics of the gas and can be used to estimate the level of non-thermal motions due to turbulence in the region. As the C¹⁸O(1–0) transition is usually optically thin, multiple peaks in the spectra, when present, likely trace the presence of independent velocity components as opposed to self-absorption effects. For better characterization of the different velocity peaks, we performed multiple Gaussian fits which allowed us to identify the centroid position of each velocity component where multiple components are observed. Comparing all of the C¹⁸O(1–0) spectra observed in a given sub-region, it was possible to identify a dominant velocity component in each case. Table 2 provides the centroid velocity and velocity dispersion of the dominant and total velocity component in each sub-region for the significant HAWC+ detection where P/σ_P ≥ 2, where P represents the polarization degree (first row for each sub-region in Table 2).

The centroid velocities of the relevant velocity components range from 5.4 to 5.9 km s⁻¹, and the associated line-of-sight velocity dispersion range from ∼0.2 to ∼0.3 km s⁻¹ for the dominant components and from ∼0.4 to ∼0.5 km s⁻¹ if all velocity components are considered.

3.1.3 Dispersion in polarization angles from HAWC+

While SR2 is dominated by one C¹⁸O velocity component [#12 with |$V_{\rm LSR} = 5.6\,$| km s⁻¹ in Hacar et al. (2013)], SR1 consists of two parts, SR1a to the north-east and SR1b to the south-west, where two distinct velocity components dominate (#11 with |$V_{\rm LSR} = 6.7\,$| km s⁻¹ and #12 with |$V_{\rm LSR} = 5.6\,$| km s⁻¹, respectively). These two components may be interacting with one another, increasing the dispersion in polarization angles. It may therefore seem justified to subdivide SR1 into these two parts (cf. SR1a and SR1b in Fig. 5) when estimating the field strength with the DCF method. Doing this results in a polarization angle dispersion of 65° ± 9° in SR1a for 45 independent HAWC+ measurements and a dispersion of 42° ± 5° in SR1b for 73 independent measurements. In both parts of SR1, the dispersion of polarization angles remains significantly higher than in SR2. For both SR1 and SR2, we also analysed the data using the structure function (SF) variant of the DCF method [Hildebrand et al. (2009); Section 3.1.1, Appendix  A].⁴ In Fig. 7, we fit 〈ΔΦ(ℓ)²〉^1/2 for the two sub-regions using the SF method after correcting for measurement error by computing the error-weighted ΔΦ²(ℓ) as in equation (6). It is common practice to restrict |ΔΦ| to be less than 90°; that is, whenever |ΔΦ| is found to be larger than 90°, it is replaced by |180° − ΔΦ|. As discussed in the Appendix, this often results in an underestimate of the dispersion in field angles and a corresponding overestimate of the field; in some cases, however, it can improve the accuracy of the field determination. We therefore provide both values. The intercepts of the fits at ℓ = 0 in Fig. 7 are ΔΦ₀ = 50.3 ± 4.4° and ΔΦ₀ = 60.5 ± 3.0° for the restricted and unrestricted approaches, respectively. The difference is only 10°. The corresponding angular dispersions contributed from the turbulence (⁠|$\sigma _{\theta } \sim \Delta \Phi _0/\sqrt{2}$|⁠) (equation A27) are about 36° and 43°, respectively. From the fitting, the turbulent correlation length-scale of SR1 and SR2 is about 4 to 5 superpixels, corresponding to 0.075–0.095 pc, about the size scale of the filamentary substructures. We divide the SR1 into four smaller portions, each with 30 polarization measurements. The root-mean-square of the angular dispersion in these four smaller portions is 30.3°, close to the angular dispersion of the turbulence from the SF analysis. The intercepts of fitting at ℓ = 0 in Fig. 7 are 23.2 ± 3.5° and 24.0 ± 3.6° from the restricted and un-restricted approaches, almost the same. The summary of all the measured parameters, and the results of the DCF and SF analysis are listed in Table 2. The estimation for the ambient cloud around the entire L1495 is also provided in the table for comparison.

Figure 7.

Structure functions and fits for sub-regions SR1 (red) and SR2 (blue) using ΔΦ restricted to less than 90° (circles and solid curve) and no restriction on ΔΦ (squares and dash curve) as functions of scale in units of HAWC + superpixels. For SR1, the fitting is from ℓ = 5 to 10. The fitted intercepts with restriction and no restriction are 50.9 ± 4.3° and 60.9 ± 3.0°, respectively. Error bars are the standard deviations of angle differences at a given distance. For SR2, the fitting is from ℓ = 4 to 18. The fitted intercepts with restriction and no restriction are 23.2 ± 3.5° and 24.0 ± 3.6°, respectively.

Open in new tabDownload slide

3.1.4 Volume density from Herschel column density data

The average volume density in each of the two portions of the B211 filament marked in Fig. 5 was estimated using the surface density map at 18.2 arcsec resolution published by Palmeirim et al. (2013) and Marsh et al. (2016) from Herschel Gould Belt survey (HGBS) data.⁵ To do so, we assumed that the depth of each sub-region along the LOS is the same as the mean projected outer width. This is a very reasonable assumption, especially for SR2 which corresponds to a segment of the filament, since there is good observational evidence that B211 is a true cylinder-like filament as opposed to a sheet seen edge-on (Li & Goldsmith 2012). The mean outer width was obtained using the projected area of pixels above the minimum surface density with a detected polarization signal [log N(H₂) = 21.59], divided by the length of each sub-region. The average surface density 〈N(H₂)〉 above log N(H₂) = 21.59 in each sub-region was derived from the Herschel column density map, and the resulting value was divided by the mean outer width, namely L ∼ 0.15 pc for SR2 and L ∼ 0.3 pc for SR1. This provided the average density, |$\langle n({\rm H_2})\rangle \, =\, \langle N({\rm H_2})\rangle /L$|⁠, given in Table 2.

3.1.5 Magnetic field strength

Using equation (3) with the volume densities, velocity dispersions, and dispersions in polarization angles estimated in Sections 3.1.3–3.1.5, we can determine the field strengths for the two sub-regions marked by white dashed rectangles in Fig. 5. The results are summarized in Table 2. We begin with SR2, which has a relatively smooth field with a small dispersion, σ_θ = 20°. The field strength ranges between 43 and 66 μG from the standard DCF method and 79–82 μG with the SF variant. Knowing the magnetic field, it is possible to determine the POS mass-to-flux ratio relative to the critical value, μ_{Φ, POS}, and the Alfv|$\acute{\rm e}$|n Mach number (Appendix  B). SR2 is trans-Alfv|$\acute{\rm e}$|nic, with |${{\cal M}_{\rm A}}\simeq 1.0\!-\!1.3$| and magnetically critical to mildly supercritical, μ_{Φ, POS} ≃ 0.9–1.7, depending on the method of analysis that is adopted. The critical mass per unit length, M_{crit, ℓ}, is that value of M_ℓ such that the pressure and magnetic forces are in balance with gravity (Appendix  B). The SR2 filament segment is slightly subcritical, with M_ℓ ≃ 0.7M_{crit, ℓ} – i.e. it is gravitationally stable against radial collapse. In the absence of perpendicular magnetic fields, filaments that are moderately subcritical (⁠|$0.9\gtrsim M_\ell /M_{{\rm crit},\ell }\gtrsim 0.2$|⁠, with an optimum value of M_ℓ/M_{crit, ℓ} ∼ 0.5) are subject to fragmentation into prestellar cores (i.e. starless cores with M ≥ M_BE) since gas can flow along the filament (Nagasawa 1987; Fischera & Martin 2012; see Appendix  B). Perpendicular fields suppress fragmentation for μ_Φ < 1. SR2 contains at least 5 candidate prestellar cores (Marsh et al. 2016), which suggests that the lower estimates of the field in Table 2 are more accurate.

By contrast, SR1 has a chaotic field with a large dispersion in polarization angles, σ_θ = 54° ± 5°. This dispersion substantially exceeds the upper limit of applicability of the DCF method recommended by Ostriker et al. (2001), as well as the less stringent criterion in Appendix A1. We note that the same remains true even if we subdivide SR1 into the two parts SR1a and SR1b considered in Section 3.1.3. None the less, the large dispersion implies a small, albeit uncertain, field: The standard method yields B₀ ∼ 7–13 μG, depending on whether the velocity dispersion is estimated from the dominant velocity component (⁠|$B_0^d$|⁠) or the total line width (⁠|$B_0^m$|⁠). For the two sub-components of SR1, the DCF method gives B₀ ∼ 10–11 μG for SR1a and B₀ ∼ 22–28 μG for SR1b. The large dispersion in angles is due in part to large-scale variations in the field structure that are allowed for in the DCF/SF analysis. Using that method with the total line width, the estimated magnetic field strength |$B_{0}^{m}$| is 23 or 16 μG, depending on whether |ΔΦ| is restricted to be less than 90° or not (see Appendix A2.1). The turbulent magnetic field strength |$\delta B_0^{m} = 18$| μG, comparable to |$B_{0}^{m}$|⁠. We note that Marsh et al. (2016) found only one candidate prestellar core in the sub-region SR1. Comparing the inset of fig. 12 of Hacar et al. (2013) with the Herschel column density map suggests that SR1 may be the location where material from the ambient cloud is presently being accreted on to B211. In particular, the fiber #11 in Hacar et al. (2013) is not straight and part of it is parallel to the striations seen in CO and Herschel data; it matches a ‘spur’ or ‘strand’ [in the terminology of Cox et al. (2016)] and may correspond to the tip of a striation where it meets and interacts with the main B211 filament (Shimajiri et al. 2019). This suggests that the flow velocities in the plane of the sky could be substantial, so that the observed LOS velocity is smaller than the POS velocities that determine σ_θ. In fact, Shimajiri et al. (2019) estimated that the inclination angle of the northeastern accretion flow to the line of sight is 70°, corresponding to a POS velocity 2.75 times larger than the LOS velocity. If so, the DCF value of the field there is an underestimate.

Myers & Goodman (1991) and Houde et al. (2009) have pointed out that if the turbulent correlation length, δ, is less than the thickness of the region being observed along the LOS, w, then the dispersion in PAs will be reduced. Houde et al. (2009) found that the reduction factor is [w/(⁠|$2\pi$|⁠)^1/2δ]^−1/2 when δ is much larger than the beamwidth. From Fig. 7, we find that δ is about 3 super pixels in size for SR2 and 5 super pixels for SR1, significantly greater than the beamwidth, which is less than one super pixel. In both cases, the turbulent correlation length is about w/3, so the reduction factor is of order unity. This is to be expected in a filament that forms in a turbulent medium. Since this effect is small compared to the uncertainties in the observations and in the method, we ignore it.

We also estimated the field strength of a larger area of Taurus/B211 using Planck polarization data from Planck Collaboration XII (2013) (at an effective HPBW resolution of 10 arcmin). The independent polarization measurements from Planck in this area (displayed as orange vectors in Fig. 5) indicate a dispersion in polarization angles of about 24° at 10 arcmin resolution. The average velocity dispersion in this extended environment around almost the entire L1495/B213 filament is ∼0.85 km s⁻¹ as estimated by Chapman et al. (2011) from ¹³CO(1–0) observations. We estimated the average volume density, |$n_{\rm H_2}\simeq 10^3$| cm⁻³, following the same approach as described in Section 3.1.5 but adopting a characteristic depth of ∼0.5 pc for the ambient cloud around Taurus/B211 (see Shimajiri et al. 2019). Applying the DCF formula of equation (3) with these values lead to a field strength of |$\sim 41\, \mu$|G.

3.2 Polarization vectors and surface density contours

As discussed in Soler et al. (2017) and references therein, the gas that feeds a cloud appears to be gathered along the magnetic field direction. Physically, it is easier for gas to flow along the field than perpendicular to the field when the field is dynamically important. Furthermore, a long, slender filament can accrete gas much more easily on its sides than at its ends. This accounts for the observation that the dense regions in many molecular clouds show magnetic fields that tend to be perpendicular to contours of the surface density (Planck Collaboration XXXV 2016).

Gas flows near the cloud determine how gas is accreted on to the cloud and thus how the cloud forms (Shimajiri et al. 2019). However, observations provide only the LOS velocity information, which can be very different from the POS velocity and thereby give a misleading idea of the true spatial gas movement (Li & Klein 2019). As noted above, fields with a substantial component normal to a filament can facilitate accretion of gas on to the filament. To assess the importance of magnetic fields in B211, we present two complementary plots of the data. In Fig. 8(a), we plot the orientations of the gradient of the surface density from Herschel data against the PAs from the HAWC + observations. Note that the gradient of the surface density is normal to the contours of surface density, so that fields perpendicular to the filament are parallel to the gradient. In Fig. 8(b), we plot a histogram of the ϕ distribution (i.e. the distribution of angles between the PAs and the tangents of the surface density contours). Because of the relatively small number of detected pixels and the limited dynamic range of the SOFIA polarization data in terms of column density, we cannot meaningfully apply a tool such as the histogram of relative orientations (HROs) as a function of surface density to the HAWC + observations in B211. Therefore, we show only one HRO in Fig. 8(b) from all the detected pixels of the observed B211 region. The histogram shape parameter for B211 is ξ = −0.28. The negative value is primarily due to SR2, which has ξ = −0.48; the chaotic field in SR1 has ξ ∼ 0. It is clear from Fig. 8(b) that there are more pixels at 90° than at 0°. The distribution of angles in this figure is similar to the high surface density Centre-Ridge region in the Vela C molecular complex. As noted above, a negative value of ξ is consistent with gas accretion along field lines that thread the cloud.

Figure 8.

(a) Orientations of the inferred magnetic field (black) from HAWC + and of the surface density gradient vectors (magenta) are plotted over the Herschel surface density. (b) Histogram of the distribution of angles between the inferred magnetic field directions and the tangents of the contours of surface density.

Open in new tabDownload slide

4 COMPARISON WITH SIMULATION

Above we used observational data from HAWC + and the IRAM 30-m telescope to obtain the LOS velocity, the magnetic field orientation, and an estimate of the field strength. In this section, we shall compare these observations with a numerical simulation that was designed not to simulate L1495 in particular, but rather to simulate the formation of filamentary structures in a typical supersonically turbulent, magnetized interstellar molecular cloud (Li & Klein 2019). Although there are some differences between the simulated filamentary cloud and L1495, such as the mass per length and probably the overall magnetic field strength in the regions, the filamentary substructures in the simulated cloud are similar to those in L1495 (Hacar et al. 2013). In fact, the results of our simulation inspired this high-resolution polarization observation of the L1495/B211 region with the aim of understanding the 3D structure of the magnetic field inside filamentary clouds.

To compare the HAWC + observational results in Section 2.1 with simulation, we use our high-resolution simulation results of the formation of filamentary molecular clouds described in detail in Li & Klein (2019). This simulation used our multiphysics, adaptive mesh refinement (AMR) code orion2 (Li et al. 2012). Since the purpose of the simulation was to study the formation of filamentary structures prior to the onset of star formation, radiation transport, and feedback physics were ignored. The ideal MHD simulation begins with turbulent driving but without gravity for two crossing times in order to reach a turbulent equilibrium state. The entire simulation region is 4.55 pc in size with a base grid of 512³. Two levels of refinement were imposed to refine pressure jumps, density jumps, and shear flows to reach a maximum resolution of 2.2 × 10⁻³ pc, which was chosen to be sufficient to study filamentary substructures with a width of order 0.1 pc. Turbulence was driven throughout the simulation at a 3D thermal Mach number |${\cal M}= 10$| on the largest scales, with wavenumber k = 1–2. Gravity was turned on after two crossing times. After gravity was turned on, we included an additional refinement requirement, the Jeans condition (Truelove et al. 1997). We adopted a Jeans number of 1/8, which means that the Jeans length is resolved by at least 8 cells. We adopted periodic boundary conditions and assumed an isothermal equation of state for the entire simulation at a temperature of 10 K. Using the turbulent line-width–size relation (McKee & Ostriker 2007), setting the Alfv|$\acute{\rm e}$|n Mach number to be 1, and setting the virial parameter to be 1, implies that the total mass of the entire cloud is |$M = 3110 \, \mathrm{ M}_{\odot }$| and the initial magnetic field is 31.6 μG. A long, massive filamentary cloud formed after gravity was turned on, and at a time of 700 000 yr, it had a length of 4.42 pc and a mass of about |$471 \, \mathrm{ M}_{\odot }$|⁠. The moderately strong large-scale field was found to be crucial in maintaining the integrity of the long and slender filamentary cloud. Details of the physical properties of the filamentary cloud can be found in Li & Klein (2019).

4.1 Simulation parameters and methods

In our simulation, even the base grid has resolution of ∼0.009 pc per cell, higher than HAWC+ superpixels. To produce the same resolution map for direct comparison with HAWC + or Planck data, we first integrate LOS quantities, such as volume density to obtain the surface density, over the base grid to create a 2D map at 512² resolution. We compute the Stokes parameters following Zweibel (1996). Density weighting is used when computing the Stokes parameters and the LOS velocity dispersion. We then coarsen the 2D map to the resolution of a HAWC + superpixel or of the Planck data by computing the mean of the corresponding number of pixels.

The surface density map of the entire simulated region is shown in Fig. 9(a). The polarization field indicating the density-weighted large-scale magnetic field at a resolution of 0.4 pc, which is the best resolution that Planck can achieve at the distance of L1495, is superimposed on the map. In the other two panels of Fig. 9, the large-scale polarization field are shown at 0.2-pc resolution. The main filamentary cloud in between the two white lines is enlarged in Fig. 9(b). The simulated filamentary cloud is composed of rich filamentary substructures along the entire length, similar to L1495 and other filamentary clouds. To study the magnetic field structures of filamentary clouds at the early stage of the formation, it is helpful to observe a cloud before the formation of protostars because powerful protostellar outflows can disrupt the magnetic field structures within filamentary substructures. The region B211 in L1495 has no protostars but contains filamentary substructures (Hacar et al. 2013) and prestellar cores (Marsh et al. 2016). Therefore, the simulated cloud is suitable for comparison with B211. Due to the collision of two filamentary clouds in our simulation at x ∼ 3.3 pc, our comparison with observations will be in the range of x = 0–2.7 pc, i.e. up to the left of the vertical yellow line in Fig. 9(b). The length of B211 with signal detected by HAWC+ is about 0.82 pc. We can create a projection of the cloud of the same length within this range for comparison. An example of a projected window, the white box in Fig. 9(b), is shown in Fig. 9(c). The small-scale magnetic field structures at the resolution of 0.019 pc, corresponding to the super-pixel resolution in the HAWC+ observation, are shown together with the low-resolution magnetic field. All the following comparisons between the simulation and HAWC + observations will be at this resolution. For clarity, we show only vectors at pixels with surface density log₁₀N(H₂) ≥ 21.59, corresponding to the minimum surface density with detected polarization signal in the observed B211 region by HAWC+. We can see the small-scale magnetic fields inside the cloud have large deviations from the low-resolution large-scale fields surrounding the dense substructures, as shown in Li & Klein (2019).

Figure 9.

(a) Surface density map of the entire simulated region as viewed along the y-axis; the mean field is in the z direction. Magenta lines indicate the large-scale magnetic field at a resolution of 0.4 pc. The filamentary cloud lies between the two long white lines. (b) Enlargement of the surface density map of the filamentary cloud with contours of log₁₀N(H₂) ranging from 21.875 to 23, separated by Δlog₁₀N(H₂) = 0.125. Due to a collision of two filamentary clouds near x = 3.3 pc, only the range x = 0–2.7 pc (up to the yellow vertical line) of the cloud is used for comparison with observation. Magenta lines indicate the large-scale magnetic field at 0.2-pc resolution. (c) Zoom in around the 0.82 pc × 0.69 pc FOV window in panel (b) showing the highly perturbed magnetic field at 28|${_{.}^{\prime\prime}}$|1 resolution (the same as the HAWC + observation). The local orientation of the field is indicated by the short black lines (shown only at pixels with log₁₀N(H₂) ≥ 21.59). As in (b), magenta lines indicate the large-scale magnetic field at 0.2-pc resolution. The surface density contours start from log₁₀N(H₂) of 21.59 and separated by Δlog₁₀N(H₂) = 0.125.

Open in new tabDownload slide

4.2 PA distribution

The HAWC + observations of B211 show a larger dispersion of PAs than the lower resolution Planck observations of the large-scale field as discussed in Section 2.4. The results indicate that small-scale perturbations of the magnetic field are present in B211. In Fig. 10, we show the PA distributions of three FOVs in the simulation. They have a length of 0.82 pc, which corresponds the HAWC + map, and height of 0.69 pc, which is large enough to include the width of the filament. The distribution in Fig. 10(a) is a single group peaking at about −15°. In Fig. 10(b), the distribution becomes double-humped, with peaks at −25° and 15°. These two FOVs along the filamentary cloud have quite different PA distributions even though they are offset by only 0.6 pc. In Fig. 10(c), which is the white coloured FOV shown in Fig. 9(b), the distribution returns to a single group again and peaks near 35°. We see that the PA distribution and the mean PA vary along the simulated cloud. In Palmeirim et al. (2013), the mean PA of the extended optical and infrared polarization vectors also changes along the filamentary cloud L1495. More polarization mapping in different parts of the filamentary cloud Taurus/B211 will be needed to find out if the PA distribution would change as in Fig. 10.

Figure 10.

PA distributions of magnetic field of three 0.82 pc long FOVs of the simulated cloud. The angle θ is measured relative to the mean direction of the field in the simulation box. (a) Single group PA distribution from an FOV starting at x = 0.16 pc; the distribution peaks at −25°. (b) Double-hump PA distribution from an FOV starting at x = 0.75 pc. (c) A single group PA distribution from an FOV starting at x = 1.69 pc, with a peak at 35°, at the location of the FOV window in Fig. 9(b).

Open in new tabDownload slide

In Fig. 11, we compare two PA distributions in the simulation by viewing the simulated cloud at the same distance of L1495, one in a small region at the HAWC + superpixel resolution of 28.1 arcsec and one in the whole simulated box at the Planck resolution of 10 arcmin. For the small region, we choose the FOV outlined in Fig. 9 since the PA distribution of this segment of the simulated cloud is similar to that of the observed B211 region. The other two FOV windows are quite different from B211, so we shall not discuss them further. The PA distribution at the Planck 10-arcmin resolution (red histogram in Fig. 11) is obtained from all the vectors in Fig. 9(a). At this resolution the dispersion is only 6.6°, much smaller than the dispersion of 29.2° of the polarization at the HAWC + superpixel scale (see Table 3). The resolution effect on the dispersion of PAs is clear both in simulation and observation (Fig. 4a). Some of this reduction in dispersion is likely due to a much lower dispersion in the low column-density gas that fills much of the Planck field: We found a dispersion of only 6.8° in a low-column region above the FOV window in the simulation. Since the polarization in the high-column LOSs is dominated by emission in the filament whereas that in the low-column LOSs is spread more uniformly over the entire LOS, the dispersion in the low-column directions is reduced by averaging along the LOS. In other words, the longer effective path-length in the low-column directions leads to an LOS resolution effect.

Figure 11.

Comparison of the PA distributions of magnetic field in the simulation FOV window of Fig. 9(b) at HAWC + 28.1 arcsec resolution (yellow histogram) and the entire simulation box (Fig. 9a) at Planck 10-arcmin resolution (red histogram; overlapping points are in orange). The 10 arcmin low resolution field vectors inside the simulation FOV window are all within the two bins between 10°–30°, marked by the blue arrows.

Open in new tabDownload slide

In addition to this observational effect, the dispersion of PAs inside a molecular cloud is increased by the combined results of differential motions of dense substructures during cloud formation (Li & Klein 2019) and small-scale local gravity-driven motion as seen in numerical simulations of molecular cloud formation (e.g. Chen, King & Li 2016; Li & Klein 2019; Seifried et al. 2020). These motions stretch the magnetic field locally, causing large changes in the direction of the magnetic field, as shown in fig. 10 of Li & Klein (2019) and Fig. 9(c) in this paper.

4.3 DCF field estimates in the simulation

Here, we apply the DCF and DCF/SF methods at HAWC + super-pixel resolution to the FOV outlined in white in Fig. 9(b). The velocity dispersion is density-weighted along the LOS. The mean width of the simulated filament, w = 0.33 pc, was computed by dividing the projected area of pixels above log N(H₂) = 21.59 (the minimum surface density of the observed B211 region with a detected polarization signal) by the 0.82-pc length of the segment. Following the procedure used in analysing the observations, we then estimated the density from the column density by assuming that the mean depth of the cloud is the same as the mean width, w.

In Table 3, we compare the results for the simulated cloud at HAWC + super-pixel resolution using the DCF and the DCF/SF methods. The turbulent correlation length, δ, in the simulated cloud segment is about 5–6 super-pixels, similar to that in SR1 and SR2 (see Fig. A1). This length is resolved by more than 40 cells at the highest resolution, so the DCF/SF results should be reliable. The magnetic field strength estimated using the DCF/SF method is in the range 41.0–43.6 μG, a little larger than the estimated value using DCF method and slightly closer to the true value.

The mean volume density, n(H₂), LOS velocity dispersion, σ_V, and dispersion of polarization PA, σ_θ, are listed in Table 3, and the Alfv|$\acute{\rm e}$|n Mach number based on the POS field in this window, |${{\cal M}_{\rm A}}/\cos \gamma$|⁠, is given in Table 4; all these values are intermediate between the values for SR1 and SR2 given in Table 2. The primary difference between the simulated and observed regions is that the simulation has a higher mass per unit length, M_ℓ, and correspondingly a smaller value of the filament virial parameter, |$\alpha _{\rm vir,f}=2\sigma _V^2/(GM_\ell)$| (Fiege & Pudritz 2000) (Table 4); the magnetic properties are similar.

To determine the true mean POS magnetic field strength above the surface density threshold in the simulation, the volume-means of the two projected magnetic field vector components for each pixel above the threshold were computed first. These vector fields were then averaged along the line of sight to a depth of 0.8 pc to obtain the mean vector field, |${\boldsymbol {\rm B}}_{0,\rm true}$|⁠. Since the value of the true field is volume weighted whereas the DCF field is based on the density-weighted polarization (which is quite different from a density-weighted field), we do not expect the true field to exactly agree with the DCF field. The true 3D field is 57.1 μG, which implies that the mean field is at an angle of 12° with respect to the plane of the sky. Table 3 shows that the DCF and DCF/SF estimates of the mean POS field are about 70 per cent of the true value. The SF variant is slightly closer to the true value than the standard DCF method, and the restricted value for the SF variant is slightly closer than the unrestricted value, but given the uncertainty in f_DCF and the fact that the true field and the DCF field have different weightings, it is not clear that these differences are significant.

The turbulent field, |$\delta B_{0,\rm true}$|⁠, is the root-mean-square of the vector difference of the field vectors from the volume-mean field vector in all the cells in the volume corresponding to the FOV. The true values of δB are significantly greater than the true values of δB_⊥ (the rms value of component δB perpendicular to the mean POS field) due to the substantial parallel component of δB in δB. The value of δB_⊥ is the same for the DCF and DCF/SF methods since both are based on equation (A2). The approximations made in the DCF/SF method imply that the value of δB calculated in this method is actually δB_⊥. The fact that the true value of δB_⊥ is about 1.7 times larger than the DCF value, B₀tan σ_θ, accounts for most of the difference between the DCF value of B₀ and the true value.

The derived physical parameters of the observed sub-regions in B211 and the segment of the simulated cloud are summarized in Table 4. The values of α_vir,f are based on the total velocity dispersions (including thermal motions) of SR1, SR2, and the segment of the simulated cloud. The normalized mass-to-flux ratio in the simulation based on the DCF field is μ_{Φ, POS} = 2.7, comparable to the observed values in SR1. Since the true value of the total POS field is 79 μG, the true value of μ_{Φ, POS} is 1.5, intermediate between the values for SR1 and SR2 and comparable to the initial value of 1.62 of the entire box viewed normal to the initial field. It is this value that we have used in determining μ_{Φ, POS} in the simulation. (Note that the use of the POS field to estimate μ_Φ leads to a slight overestimate of M_ℓ/M_{crit, ℓ} – see equation B6.) Including the effect of a perpendicular magnetic field, the ratio of the line mass to the critical line mass, M_ℓ/M_{crit, ℓ} (equation B8), is less than unity for both observed subregions and for the simulated cloud segment, as shown in the table. These structures are therefore gravitationally stable against radial collapse on the scale at which this ratio is determined. In the absence of perpendicular magnetic fields, the critical mass and the virial mass are the same, and the filaments would be subject to fragmentation (Nagasawa 1987; Fischera & Martin 2012; see Appendix B2). Perpendicular fields would stabilize the cloud against fragmentation if the normalized mass-to-flux ratio is μ_Φ < 1. The results in Table 4 show that SR1 and the simulated filament should be subject to fragmentation; since SR2 has μ_Φ ∼ 1 it is marginally susceptible. There is some evidence that pre-stellar cores are forming in both SR1 and SR2 (e.g. Marsh et al. 2016); as noted above, the fact that pre-stellar cores are observed in SR2 favours a lower estimate for the field there than the value given in Table 4. Equation (B16) shows that fragments that form in SR1 are near the critical mass and could collapse, whereas those in SR2 appear to be stable against collapse; however, given the uncertainty in the parameters in Table 4, these conclusions are tentative. The FOV window in the simulation is magnetically supercritical, and dense cores are forming along some filamentary substructures. On small scales, we expect the velocity dispersion to be primarily thermal. At a temperature of 10 K, the thermal values of the virial parameter are |$\alpha _{\rm vir,f}= 0.30,\, 0.46,\, \rm{and}\,\, 0.15$|⁠, respectively.

5 CONCLUSIONS

In this work, we have used HAWC + on-board SOFIA to observe the L1495/B211 region in Taurus to investigate the magnetic field morphology in thin filamentary clouds. This observation is challenging because of the low surface brightness of the filamentary cloud. We needed to re-sample 3 × 3 detector pixel data to a super-pixel of 28|${_{.}^{\prime\prime}}$|1 to optimize the SNR. We have total 282 independent measurements that have P/σ_P ≥ 2. The morphology of the observed polarization map clearly reveals two sub-regions, designated SR1 and SR2, in the observed B211 region. With IRAM 30-m C¹⁸O (1-0) data, we estimate the magnetic field strengths using the standard DCF method and the alternative DCF method using a structure function. We then compared the physical states of the two sub-regions with a simulated filamentary cloud.

1. Polarization morphology of the two sub-regions in B211. The chaotic appearance of the polarization vectors in SR1 indicates a strongly perturbed region, in contrast to SR2, which has a well organized magnetic field structure mostly, but not entirely, perpendicular to the filamentary cloud axis. The organized field in SR2 matches the large-scale field from Planck observation very well. The dispersion of the PAs in SR1 is |$54\deg$|⁠, almost 3 times of that in SR2.

2. Filamentary substructures in B211. The IRAM 30 m C¹⁸O (1-0) data reveals multiple velocity components in the observed B211 region, similar to what has been reported in a previous study of L1495 (e.g. Hacar et al. 2013). There are at least 3 velocity components in SR1 and 2 velocity components in SR2. Multiple filamentary substructures are also clearly seen in the high resolution Herschel map. The chaotic appearance of substructures, polarization vectors, and the multiple-component line profiles in SR1 may indicate strong interaction among substructures.

3. Magnetic fields of the two sub-regions and of the simulated filament. Using the DCF and DCF/SF methods, the estimated field strength based on the total LOS velocity dispersion in SR1 is in the range 13–23 μG. Because of the very large dispersion of the polarization angles (σ_θ = 54°), the field estimate in this region is very uncertain, but it is clear that the field is small. By contrast, the estimated field strength of SR2 is from 66 to 82 μG, significantly larger than that in SR1. These estimates are based on the assumption that the numerical coefficient introduced to correct for the approximations in the DCF method, f_DCF, is 0.5 (Ostriker et al. 2001). In the part of the simulated filament that we analyzed in detail, the field strength is intermediate between that of SR1 and of SR2. The measured value of f_DCF was slightly larger than 0.5, but consistent with that value within the expected statistical uncertainties.

4. Comparison with Zeeman field estimates. Based on the Zeeman data summarized by Crutcher et al. (2010), Li et al. (2015) concluded that the average 3D magnetic field in molecular clumps in the ISM is |$66n_{{\rm H_2},4}^{0.65}$| μG. For a typical inclination with respect to the plane of the sky of 30°, this corresponds to a POS field |$B_0=57n_{{\rm H_2},4}^{0.65}$| μG. This is several times larger than the DCF estimate of the field in SR1, and we suggested that this could be due to the measured LOS velocity dispersion being less than the POS velocity dispersion. The POS field (for γ = 30°) corresponding to the average interstellar Zeeman field agrees reasonably well with the DCF field in SR2 and with the true POS field in the simulation.

5. Resolution effect on the magnetic field dispersion. The dispersion in polarization angles from the low resolution Planck data is significantly smaller than that of the high resolution HAWC + data. Heitsch et al. (2001) and Falceta-Gonçalves et al. (2008) found this resolution effect in their simulations, and we do also. The simulation shows that the angle dispersion in low-column regions is less than in high-column regions, which contributes to the observed resolution effect.

6. Polarization vectors and surface density gradients. The relative distribution of the inferred magnetic field vectors and the tangent of surface density contours in the observed B211 region shows that the magnetic field has a tendency to be normal to the contours of surface density. This can be quantified by the histogram shape parameter, ξ, defined in equation (8). In B211, we find ξ = −0.28, meaning that the number of pixels with a magnetic field nearly normal to the contours of the surface density is about 1.8 times that with the magnetic field nearly parallel to the contours. The tendency for the field to be normal to the contours is primarily due to SR2, which has ξ = −0.48 and an average angle between the contours and the field of 〈ϕ〉 = 55°. The fact that there is some correlation between the orientation of the field and the column density contours of the gas indicates that the magnetic field is at least marginally dynamically important there. In the chaotic region SR1, the fact that ξ = −0.03 indicates that the magnetic field is dynamically sub-dominant, in agreement with the large value of the projected Alfv|$\acute{\rm e}$|n Mach number there (Table 4).

7. Physical states of the two subregions and of the simulated filament. From the mass-to-flux ratios and Alfv|$\acute{\rm e}$|n Mach numbers, SR1 is magnetically supercritical and slightly super-Alfv|$\acute{\rm e}$|nic, although we have suggested that the DCF method underestimates the field in SR1. The magnetic field in SR2 is significantly greater than that in SR1. Both the standard DCF analysis and the DCF/SF method suggest that SR2 is approximately magnetically critical and that it is trans-Alfv|$\acute{\rm e}$|nic. The segment of the simulated filament we have analysed is magnetically supercritical like SR1, although it has a significantly smaller dispersion of PA angles; it has an Alfv|$\acute{\rm e}$|n Mach number of about unity. The ratio of the line mass to the critical line mass is slightly less than unity for SR1, SR2, and the simulated filament if the full velocity widths of the filaments are used to estimate the virial parameters. Pre-stellar cores are suggested in both SR1 and SR2 (Marsh et al. 2016). There are two cores forming in the segment of the simulated filament that we have analysed.

8. The DCF method. In Appendix  A, we present derivations of both the standard DCF method and the structure function (SF) variant that are not restricted to small values of the polarization angles. We show that the standard DCF result often applies for the case of equipartition even if the perturbed field is not due to Alfv|$\acute{\rm e}$|n waves. Our simulation confirms that δ, the correlation length of the turbulent magnetic field, is small, as assumed in the derivation of the DCF/SF method (Hildebrand et al. 2009). For SR1, SR2, and our simulation, we find that δ ∼ FWHM of the filament ∼0.1 pc, consistent with the formation of a filament in a turbulent medium. We discuss the restriction procedure often used in the DCF/SF method in which differences in angles that exceed 90° are converted to |180–90°| and suggest that restriction provides a lower limit on the structure function and is significant only when the dispersion in PAs is large, so that the DCF method is of questionable accuracy.

9. Different versions of the DCF method. In Appendix  A, we also compare the standard DCF method with the DCF/SF version (Hildebrand et al. 2009) and the parallel-δB version (Skalidis & Tassis 2021; Skalidis et al. 2021). In most cases, the three methods agree within the uncertainties for both the observed regions, SR1 and SR2, and for the simulation. The exception is the standard DCF method, which gives a low value for the mean field in the highly tangled region SR1, probably because this method does not allow for spatial variation of B₀.

10. Equilibrium filaments and their fragmentation. In Appendix  B, we give analytical estimates of the fragment mass and the condition for the formation of a pre-stellar core in an unmagnetized filament.

ACKNOWLEDGEMENTS

Support for this research is based on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy (SOFIA) under the 07_0017 Program. SOFIA is jointly operated by the Universities Space Research Association, Inc. (USRA), under NASA contract NNA17BF53C, and the Deutsches SOFIA Institut (DSI) under DLR contract 50 OK 0901 to the University of Stuttgart. We thank Che-Yu Chen, Amitava Bhattacharjee, Martin Houde, Alex Lazarian, and Junhao Liu for a number of very helpful conversations. We also thank the two anonymous referees for their many helpful suggestions that greatly improve the paper. Support for this research was provided by NASA through a NASA ATP grant NNX17AK39G (RIK, CFM, and PSL) and the US Department of Energy at the Lawrence Livermore National Laboratory under contract DE-AC52-07NA 27344 (RIK). CFM acknowledges the hospitality of the Center for Computational Astronomy of the Flatiron Institute in New York, where he was a visiting scholar. JR thanks partial support from SOFIA program 07_0017 and support from the SOFIA program 07_0047 and NASA Astrophysics Data Analysis grant (80NSSC20K0449). PhA acknowledges support from ‘Ile de France’ regional funding (DIM-ACAV + Pro- gram) and from the French national programs of CNRS/INSU on stellar and ISM physics (PNPS and PCMI). This work is partly based on observations carried out under project number 129-15 with the IRAM 30-m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain). The present study also made use of data from the Herschel Gould Belt survey (HGBS) project (http://gouldbelt-herschel.cea.fr). This work used computing resources from an award from the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation grant number ACI-1548562, through the grant TG-MCA00N020, computing resources provided by an award from the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center, and an award of computing resources from the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231.

DATA AVAILABILITY

The processed HAWC + data in FITS format is available at CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via https://cdsarc.unistra.fr/viz-bin/cat/J/MNRAS.

Footnotes

REFERENCES

APPENDIX A: THE DAVIS–CHANDRASEKHAR–FERMI METHOD

Davis (1951) and Chandrasekhar & Fermi (1953) proposed a method for estimating magnetic field strengths in the ISM based on the assumptions that the medium is isotropic and that variations in the orientation of the field are due to Alfv|$\acute{\rm e}$|n waves. Hereafter, we refer to this as the DCF method. However, there are different approximations used and assumptions made in the literature, particularly in dealing with a large dispersion in the polarization angles (PAs). Therefore in this appendix we give a more rigorous derivation of the DCF result based on the method of Hildebrand et al. (2009). We then discuss two variants of the DCF method, the structure function method (DCF/SF) of Hildebrand et al. (2009) and the parallel-δB version of Skalidis & Tassis (2021). In applications of the DCF/SF method, differences between the PAs at different points are often restricted to be less than 90°, and we show how that can be problematic.

Isotropy is an important assumption in the DCF method. Observations by Heyer et al. (2008) of the Taurus molecular cloud show that the turbulence there is anisotropic; it is not known if this is typical for molecular clouds. Their simulations for |$\beta =2c_{\rm s}^2/v_{\rm A}^2=0.02$| were strongly anisotropic, with 1D velocities perpendicular to the mean field 2–4 times greater than those along the field. As noted in the discussion of Alfv|$\acute{\rm e}$|nic turbulence above, the DCF method can still be applied in the presence of such anisotropy if the mean 3D field is close to the POS. (The median value of the inclination γ is 30°, for which cos γ ≃ 0.87–1.) For larger values of β, the simulations of Heyer et al. (2008) for β ≥ 0.2 and those of Heitsch et al. (2001) for β ≥ 0.05 showed approximately isotropic turbulence.

Another important assumption that went into the derivation of equation (A2) for Alfv|$\acute{\rm e}$|n waves and equation (A3) for the case of equipartition is that a single turbulent region dominates the signal along the LOS; if there is one dominant object along the LOS, its depth must be smaller than, or at most comparable to, the turbulent correlation length. If there are multiple turbulent regions, then σ_V includes the differences in mean velocities of the regions and |$\sigma _{\delta B_\perp }$| includes the differences in the mean field orientation along the LOS. Such effects have been analysed by Zweibel (1990), Myers & Goodman (1991), and Houde et al. (2009). As discussed in Section 3.1.5, possible effects of this type of inhomogeneity in the region we have observed are small.

The ratio |$\sigma _{\delta B_\perp }/B_0$| is estimated from fluctuations in the orientation of the field as revealed by polarization observations. We now discuss two methods of doing this, the standard method and the structure function method developed by Hildebrand et al. (2009). Bear in mind that a basic assumption of the DCF method is that the polarization traces an appropriately weighted (including by the density) integral of the direction of the magnetic field along the LOS. It must be borne in mind that the polarization angles (PAs), θ_i, are limited to the range −90° ≤ θ_i ≤ 90°, whereas the field angles (FAs), θ_{FA, i}, extend over the range −180° ≤ θ_{FA, i} ≤ 180°, so that that there is a 180° ambiguity in the relation between the FAs and the PAs. Martin (1974) showed that the PA traces the mean FA in the simple case in which the FA is a linear function of position and the density is constant; for variable density, the conclusion holds if the FA is a linear function of surface density.

A1 The standard DCF method

Different assumptions lead to different approximations for |$\sigma _{\delta B_\perp }$|⁠. For example, Zweibel (1996) assumed that the FAs are identical to the PAs at the outset. She therefore excluded the possibility that individual FAs could exceed 90°, in contrast to our approach. With δB_i = B₀tan θ_i, averaging δB² over different lines of sight gives |$\sigma _{\delta B_\perp }^2/B_0^2={\langle \tan ^2\theta \rangle }$|⁠; she obtained the same result through an analysis using the Stokes parameters. Heitsch et al. (2001) recognized that this is problematic for flows with Alfv|$\acute{\rm e}$|n Mach numbers |$\gtrsim 1$| since the average of tan ²θ is dominated by angles near 90°, and they suggested several approximations to overcome this. As noted above, Falceta-Gonçalves et al. (2008) suggested replacing the average of the tangent by the tangent of the average, which we derived above; this also overcomes this problem.

A2 The structure function version of the DCF method (DCF/SF)

Hildebrand et al. (2009) improved on the standard DCF approach by allowing the direction of the mean magnetic field to be a slowly varying function of position, B₀(x); the magnitude of the field was assumed to be constant, however. A strength of their method is that the unknown direction of the mean field is not needed in order to determine its magnitude. Furthermore, it is relatively independent of the dispersion of PAs on large scales and can therefore handle cases in which large dispersions on large scales cause the standard DCF method to break down. Houde et al. (2009) extended this method to allow for variations along the line of sight and across the telescope beam, but at the expense of adding an additional parameter that must be fit from the data. Here, we follow the simpler approach of Hildebrand et al. (2009). We include possible effects of integration along the line of sight, in addition to the effects of other approximations made in the method, in a numerical factor f_DCF, as in equation (A9).

Hildebrand et al. (2009) assumed that δ, the maximum scale of the turbulent velocity correlations, was of order 1 mpc, well below the resolution of the observations they were fitting. Subsequent work using the method of Houde et al. (2009) obtained larger values: for example, Guerra et al. (2021) found δ ∼ 10−100 mpc for OMC-1. Our analysis of SR1 and SR2 also gives δ ∼ 100 mpc, as does our simulation (see below). The observations we have analysed and our simulation are consistent with the turbulent correlation length δ being of order the FWHM of the filament, which is plausible for a filament formed in a turbulent medium and consistent with the results of Palmeirim et al. (2013). Note that the value of δ does not enter; all that is required is that it be small enough that there is a range of ℓ over which ΔΦ(ℓ)² is accurately fit by the first two terms in equation (A21).

How well does our simulation agree with the SF variant of the DCF method? The actual values of the field strength are compared with the DCF/SF values in Table 3. Here, we focus on the validity of the SF relation between the dispersion in field angles and the dispersion in field strength, equation (A21). For simplicity, we adopt the small angle approximation; had we not done that, the results would have changed by only 7 per cent. Fig. A1 compares the dispersion in angle, ΔΨ(ℓ), measured in the FOV window of the simulation with values from equation (A21). The blue curve plots the first two terms in the equation, using the measured value of |$\sigma _{\delta B_\perp }/B_{\rm rms}$| and treating m as a free parameter. It provides an excellent fit to the data for |$\ell \gtrsim 5$| superpixels, or about 0.1 pc for the parameters we have adopted. This is the turbulent correlation length, δ, and is about 1/3 the width of the filament, as noted above. It is well resolved in the simulation, with more than 40 grid cells at the highest resolution. The red curve shows the last term in equation (A21). As assumed by Hildebrand et al. (2009), it is negligible except at small scales, ℓ < δ ∼ 0.1 pc. The figure shows that the approximations made in deriving equation (A21) are reasonably good in this case, and that the turbulent velocity correlations extend to scales large enough that they must be taken into account.

Figure A1.

The structure function of magnetic field vectors in the simulation window (blue circles) as a function of scale in units of HAWC + superpixels (28.1 arcsec), with the best fit (solid curve) and the negative of the last term of equation (A21) (red squares).

Open in new tabDownload slide

A2.1 Restriction of ΔΦ

A fundamental problem with determining the field strength from polarization observations is that the field angles (FAs), Ψ, can range over 360° whereas the PAs, Φ, are limited to a range of 180°. For a given choice of the direction corresponding to 0°, let the orientation of the PAs lie in the range ±90°. FAs lying outside that range will have PAs in the opposite direction – i.e. such an FA will differ from the corresponding PA by 180°. As a result, the measured value of ΔΦ, based on the PAs, will differ from the actual value, which is based on the FAs. The error depends on the number of FAs that are flipped in direction, which in turn depends on the choice of the 0° direction; we choose that to be the direction that gives the minimum dispersion in the PA angles, Φ_i. In an attempt to reduce this error in the DCF/SF method, it is common to restrict the difference between angles, ΔΦ, to be less than 90° by replacing |ΔΦ| with |180° − ΔΦ| when |ΔΦ| > 90° (e.g. Davidson et al. 2011). Under what conditions is this valid? If the angles are restricted when they should not be, then the dispersion will be underestimated and the field overestimated.

First assume that the PAs are an accurate reflection of the FAs, up to an ambiguity of 180°. If the FAs are confined to a range less than 180°, then the FAs and PAs can be in alignment and restricting ΔΦ would lead to an error. If the FAs extend beyond that range, but the mean field has a constant direction, then the dispersion in the PAs will be less than that in the FAs and the the field strength will be overestimated. The error will only grow larger if ΔΦ is restricted. It follows that restriction should never be used if the mean field has a constant direction.

If B₀ changes direction as a function of position, the situation becomes more complicated. If the average angle differs significantly from the FAs in a local region, then it is possible that some of those FAs will be flipped by 180° when converted to PAs, thereby increasing the dispersion relative to neighbouring PAs. Restriction corrects this by significantly reducing ΔΦ for the PAs in that region. On the other hand, it also reduces ΔΦ for the PAs that were initially quite different from the initial FA direction. This discussion suggests that restriction provides a lower limit on ΔΦ₀ and that it is significant only when the dispersion in angles is large, when the DCF method is of questionable accuracy. For SR2, with σ_θ = 20°, restriction reduced ΔΦ₀ by 4 per cent and therefore increased B₀ by the same factor; for the simulation, with σ_θ = 29°, restriction increased the inferred field by 5 per cent; and for SR1, with σ_θ = 54°, restriction increased B₀ by 20 per cent. An approach that reduces the uncertainties due to restriction is to map the field locally (Guerra et al. 2021), so that there is less variation of the mean field in each region.

A3 The parallel-δB version of the DCF method

In our simulation, we find that |${\langle |\delta {\hat{{\boldsymbol B}}}\cdot {\hat{{\boldsymbol B}}}_0| \rangle }\simeq 0.8$|⁠, so the parallel component of δB is indeed significant. On the other hand, the positive and negative values nearly cancel so that |${\langle \delta {\hat{{\boldsymbol B}}}\cdot {\hat{{\boldsymbol B}}}_0 \rangle }\simeq 0.06$| – i.e. rarefactions, which have |$\delta {\hat{{\boldsymbol B}}}\cdot {\hat{{\boldsymbol B}}}_0\lt 0$|⁠, nearly cancel the effect of compressions, which have |$\delta {\hat{{\boldsymbol B}}}\cdot {\hat{{\boldsymbol B}}}_0\gt 0$|⁠. This effect is not included in the model of Skalidis & Tassis (2021) since they assumed that 〈δB · B₀〉 can vanish only for incompressible turbulence and proceeded to make the incorrect assumption that 〈δB · B₀〉 = 〈|δB|〉B₀. They attempted to justify this step by appealing to Bhattacharjee, Ng & Spangler (1998), although that work applies only to very subsonic turbulence and has 〈δB〉 = 0. Skalidis et al. (2021) argued that the maximum kinetic energy in fluctuations, a second-order quantity, is in equipartition with the maximum magnetic energy in the fluctuations, a first order quantity; as shown by Zweibel & McKee (1995), however, it is the second-order energies that are in equipartition. In agreement with equation (A11), they note that a non-zero polarization angle is possible only in the presence of a perpendicular component of the field. As a result their method requires δB_∥ ≃ δB_⊥, which they find to be satisfied to within a factor 2 in the simulations they analyse.

For tan σ_θ ≃ σ_θ, the ratio of their result to the standard DCF one is (σ_θ/2)^1/2/f_DCF. For f_DCF = 0.5, the two values of the field agree for σ_θ = 27°; since that is close to the values we find in our simulation, we are not able to determine whether their result is more accurate than the standard DCF method. It should be noted, however, that their result has no free parameters, whereas f_DCF is a free parameter for the standard method. In view of the questionable assumptions underlying their method, more work is needed to understand the physical basis for the method and the circumstances under which it works.

Applying their method to the observed fields in B211 with the full line width |$\sigma _V^m$| gives B₀ = 27 μG for SR1, about twice the value with the standard DCF method but only slightly larger than the 23 μG with the DCF/SF method with restriction. If σ_θ is replaced by tan σ_θ in their formula, their result would be 22 μG for this region, not that much larger than the DCF/SF value without restriction, 16 μG. For SR2 they find B₀ = 57 μG, somewhat less than the 66 μG with the standard DCF method and ∼80 μG with the DCF/SF method. In most cases the three methods agree within the uncertainties for both observation and simulation. The exception is the standard DCF method, which gives a low value for the tangled-field SR1, most likely because it includes variation in B₀ in its determination of σ_θ. The parallel-δB method includes such variation as well, but the result is less sensitive since it enters only as the square root.

APPENDIX B: EQUILIBRIUM AND FRAGMENTATION OF FILAMENTS

B1 Equilibria of cylindrical filaments

B2 Fragmentation of filaments stable against radial collapse

The growth rate of the instability is within 25 per cent of 0.3(⁠|$4\pi$|Gρ_c)^−1/2 for all R/H based on Fischera & Martin (2012)’s fit to the results of Nagasawa (1987). Thus, unmagnetized filaments that are stable against radial collapse always fragment, although the amplitude of the perturbation could be small, as we shall now see.

Motiei et al. (2021) studied the stability of polytropic filaments with a field parallel to the filament. For B = 0 and |$\gamma =\frac{3}{4}$| (the case closest to γ = 1), they found that the most unstable wavelength and the critical wavelength are within about 30 per cent of the values for an isothermal filament (equation B11). This suggests that there is no significant difference between the unstable wavelengths of isothermal filaments and polytropic filaments with γ → 1 despite the fact that isothermal filaments have a much steeper density gradient at large radii, ρ ∝ r⁻⁴ (e.g. Ostriker 1964), than polytropic filaments, ρ ∝ r^{−2/(2 − γ)} (Viala & Horedt 1974). The results discussed above should therefore apply to both isothermal filaments and polytropic ones with γ → 1. We note that observationally, the polytropic solution appears to be favoured – for example, Palmeirim et al. (2013) found ρ ∝ r⁻² with γ just below 1.