The JCMT BISTRO Survey: Studying the Complex Magnetic Field of L43

We observed L43 at 450 and 850 μm using the Submillimetre Common-User Bolometer Array 2 (SCUBA-2)/POL-2 on JCMT; however, this work focuses only on the 850 μm observations, and we leave the 450 μm data for a potential future multiwavelength study. The observations were taken between 2020 February and 2021 March as part of the BISTRO large survey program (Project ID: M20AL018) in its third generation of observations, "BISTRO-3." The observations consisted of 27 repeats of ∼31 minutes with one observation listed as questionable which we omitted from the data reduction. The data were observed in Band 1 weather conditions with the atmospheric opacity at 225 GHz (τ₂₂₅) less than 0.05. The JCMT has a primary dish diameter of 15 m and a beam size of 146 at 850 μm when approximated with a two-component Gaussian (Dempsey et al. 2013). The observations were performed using a modified SCUBA-2 DAISY mode (Holland et al. 2013) optimized for POL-2 (Friberg et al. 2016), which produces a central 3' region with uniform coverage with noise and exposure time increasing and decreasing, respectively, to the edge of the map. The 3' region covers most of the L43 molecular cloud including the starless core and RNO 91. This mode has a scan speed of 8'' s⁻¹ with a half-wave plate rotation frequency of 2 Hz (Friberg et al. 2016).

To reduce the data, we used the Submillimetre User Reduction Facility (SMURF) package (Chapin et al. 2013) from the Starlink software (Currie et al. 2014). The SMURF package contains the data reduction routine for SCUBA-2/POL-2 observations named pol2map. ¹⁰¹

In the first step, the raw bolometer time streams are separated into Stokes I , Q, and U time streams. The command makemap (Chapin et al. 2013) is then called to create an initial Stokes I map from the Stokes I time streams. The Stokes I map from the first step is used to create an AST and PCA mask at a fixed signal-to-noise ratio (S/N), which is then used to mask out nonastronomical signals (we refer to these as the auto-generated masks). The AST mask is used to define background regions that are forced to zero after each iteration in order to prevent growth of spurious structures in the map. The PCA mask is used to define regions that are excluded from the principle component analysis (PCA), which removes correlated large-scale noise components from the bolometer time streams. These excluded regions are generally source regions since these contain uncorrelated data and would disrupt the PCA process. The second step of the reduction creates the final Stokes I , Q, and U maps and a polarization half-vector catalog. The polarization vectors are often referred to as "half-vectors" due to the 180° ambiguity since direction is not known (e.g., 45° is treated the same as 225°). We included the parameter skyloop and followed the same data reduction technique as Pattle et al. (2021). We also set the parameter pixsize to a value of 8 to reduce the data with 8'' pixels as explained in Section 3.3.

We corrected for instrumental polarization (IP) in the Stokes Q and U maps based on the final Stokes I map and the "August 2019" IP polarization model. ¹⁰² The 850 μm Stokes I, Q, and U maps were also multiplied by a flux conversion factor (FCF) of 748 Jy beam⁻¹ pW⁻¹ to convert from pW to Jy beam⁻¹ and account for loss of flux from POL-2 inserted into the telescope. This value was calculated using the post-2018 June 30 SCUBA-2 FCF of 495 Jy beam⁻¹ pW⁻¹ (Mairs et al. 2021) multiplied by a factor of 1.35 to account for the additional losses in POL-2 (Friberg et al. 2016) and then multiplied by a factor of 1.12 to account for the 8'' pixels. This extra factor was determined from SCUBA-2 calibration plots. ¹⁰³ The final Stokes I map has an rms noise of ≈2.2 mJy beam⁻¹. To further increase the S/N of our polarization half-vectors and attempt to account for the JCMT beam size, we binned them to a resolution of 12''. The polarization half-vectors are also debiased as described in Wardle & Kronberg (1974) to remove statistical bias in regions of low S/N (see Equation (1)).

The values for the debiased polarization fraction P were calculated from

$$\begin{eqnarray}&&P=\displaystyle \frac{1}{I}\sqrt{{Q}^{2}+{U}^{2}-\displaystyle \frac{{Q}^{2}\delta {Q}^{2}+{U}^{2}\delta {U}^{2}}{{Q}^{2}+{U}^{2}}},\end{eqnarray} \tag{ 1 }$$

where I , Q, and U are the Stokes parameters, and δ Q, and δ U are the uncertainties for Stokes Q and U . The uncertainty δ P of the polarization degree was obtained using

$$\begin{eqnarray}&&\delta P=\sqrt{\displaystyle \frac{({Q}^{2}\delta {Q}^{2}+{U}^{2}\delta {U}^{2})}{{I}^{2}({Q}^{2}+{U}^{2})}+\displaystyle \frac{\delta {I}^{2}({Q}^{2}+{U}^{2})}{{I}^{4}}},\end{eqnarray} \tag{ 2 }$$

with δ I being the uncertainty for the Stokes I total intensity.

The polarization position angles θ, measured from north to east in the sky projection (north is 0°), were measured using the relation

$$\begin{eqnarray}&&\theta =\displaystyle \frac{1}{2}{\tan }^{-1}\displaystyle \frac{U}{Q}.\end{eqnarray} \tag{ 3 }$$

The corresponding uncertainties in θ were calculated using

$$\begin{eqnarray}&&\delta \theta =\displaystyle \frac{1}{2}\displaystyle \frac{\sqrt{{Q}^{2}\delta {U}^{2}+{U}^{2}\delta {Q}^{2}}}{({Q}^{2}+{U}^{2})}\times \displaystyle \frac{180^\circ }{\pi }.\end{eqnarray} \tag{ 4 }$$

The plane-of-sky orientation of the magnetic field is inferred by rotating the polarization angles by 90° (assuming that the polarization is caused by elongated dust grains aligned perpendicular to the magnetic field).

Standard reductions of SCUBA-2/POL-2 observations are done with a 4'' pixel size (e.g., Pattle et al. 2021). Nearly all of these reductions have been done on bright, high S/N sources. As we approach the limit of the POL-2 polarimeter, we can explore the use of larger pixel sizes to attempt to boost the S/N in these dense, starless sources. The use of different pixel sizes in reductions of JCMT SCUBA-2 data has been explored before, such as in the Gould Belt Survey (Ward-Thompson et al. 2007) where originally 6'' (see Sadavoy et al. 2013) and 3'' (see Mairs et al. 2015) pixels were used with the latter being chosen in order to recover small-scale structure. The current pixel size of 4'' was picked in order to properly sample the Gaussian beam and allow the mapmaking algorithm to converge in a reasonable time (Chapin et al. 2013), as well as avoid smoothing of larger pixels. However, with faint sources such as starless cores, we need to investigate the potential of using larger pixel sizes.

One issue with using larger pixels is that the larger pixel size tends to produce masks that cover a larger area of the sky. Doubling the pixel size from 4'' to 8'' typically causes the number of bolometer samples falling in each pixel to increase by a factor of 4, thus increasing the S/N of each pixel value by a factor near to 2. Since each mask is defined by a fixed S/N cutoff, this causes a larger fraction of the map to be covered by the mask. An increase in the size of the AST mask is potentially problematic, as it can encourage the growth of artificial large-scale structures within the masked areas (see Chapin et al. 2013). Distinguishing such artificial structures from real astronomical signals requires care.

Our solution to this problem is to re-use the 4'' pixel auto-generated masks when creating externally masked maps with 8'' pixels, rather than using new masks based on the auto-masked 8'' maps. The smaller 4'' masks will then restrict the growth of artificial extended structures giving us more confidence in the remaining extended structure. To do this, we ran the entire reduction using the standard 4'' pixel size. We then regridded the AST and PCA masks from that reduction to 8'' using the command compave from the KAPPA package (Currie & Berry 2014). We then ran the second step of the reduction using the regridded AST and PCA masks to create the externally masked Stokes I, Q, and U maps as well as the polarization vector catalogs, using a pixel size of 8''. This resulted in a molecular cloud that looked similar to the original 4'' reduction but with better S/N and therefore more polarization half-vectors (vector catalog increased from 98 vectors to 133 vectors at the same S/N cut). This is the reduction that we present in this work.

As a further check, we performed a jackknife test by dividing our observations into two populations and comparing the Stokes I, Q, and U maps from both populations. The results from this test can be found in the Appendix. We saw a more significant difference between the populations when using the auto-masked 8'' maps. This difference occurred mainly in the areas where emission was present in the 8'' maps but not present in the 4'' maps, raising further doubt as to the validity of the new extended emission in the auto-masked 8'' maps. Any differences seen in 8'' reduction done using the regridded masks were the same as differences seen in 4'' reduction, just smoothed due to larger pixel sizes.

We used archival observations of the CO J = 1–0 line carried out with the Berkely Illinois Maryland Array (BIMA) 10 antenna interferometry array. The CO J = 1–0 data were obtained from Lee et al. (2002), and details of the observations and data reduction can be found therein. BIMA has a similar beam size to that of JCMT at 128.

We also used archival observations of the CO J = 3–2 line carried out with the Heterodyne Array Receiver Program (HARP) to remove CO contribution from the 850 μm Stokes I map. This is discussed further in Section 4.3. The data were accessed from the Canadian Astronomy Data Centre database ¹⁰⁴ (Project ID: M07AU11) and were downloaded as reduced spectral cubes that were then mosaicked using the PICARD recipe MOSAIC_JCMT_IMAGES. ¹⁰⁵

Figure 2 shows the 850 μm dust contours in green overlaid on the Herschel SPIRE 250 μm where the 850 μm dust traces the densest part of the filament. The filament does continue to the east and west, but this may be more extended structure and is therefore lost by SCUBA-2/POL-2. The 850 μm dust traces the northern edge of the CO outflow cavity, which is discussed in the next section. The 850 μm emission is peaked in the main starless core (L43) and then the dust region surrounding RNO 91.

Figure 3 shows the column density map, which is discussed later in Section 4.3 but it has the 850 μm dust contours overlaid with more levels to better show the emission structure. In the main starless core, the densest emission peaks toward the center, but then there are two lobes that extend to the northwest and southeast. A small peak can be seen in southeast lobe in Figure 3. We do not have resolved kinematic or significant magnetic field data between these three areas (the center part and the two lobes) so it is not possible to tell if they are fragmenting. However, the 850 μm emission shows structure suggesting these could be on the way to fragmentation.

Zoom In Zoom Out Reset image size

We can model these three regions as Bonnor–Ebert (BE) spheres (Ebert 1955; Bonnor 1956) and estimate their critical BE masses. We take the sound speed c_s to be ∼0.19 km s⁻¹, which was calculated assuming a dust temperature of 12.1 K (Planck Collaboration et al. 2016c). The critical BE mass can be calculated using the relation (Equation (3.9), Bonnor 1956),

$$\begin{eqnarray}&&{M}_{\mathrm{BE},\mathrm{crit}}=3.3\displaystyle \frac{{c}_{{\rm{s}}}^{2}}{G}{R}_{\mathrm{crit}},\end{eqnarray} \tag{ 5 }$$

where G is the gravitational constant and R_crit is the critical radius of the core. We estimate R_crit from the observed flux structure (where the flux drops before peaking again in the center) and use it to also calculate total flux for total mass estimates. Values for R_crit are given in Table 1 along with locations of the three potentially fragmenting regions in the submillimeter core (northwest, main, and southeast). The estimated critical BE masses of the two lobes are ${M}_{\mathrm{BE},\mathrm{crit}}^{\mathrm{NW}}\sim 0.21\,{M}_{\odot }$ and ${M}_{\mathrm{BE},\mathrm{crit}}^{\mathrm{SE}}\sim 0.20\,{M}_{\odot }$ for the northwest and southeast lobes, respectively. For the central (main) peak, the estimated critical BE mass is ${M}_{\mathrm{BE},\mathrm{crit}}^{\mathrm{main}}\sim 0.32\,{M}_{\odot }.$

We can estimate the total mass from the 850 μm dust emission using the relation from Hildebrand (1983),

$$\begin{eqnarray}&&{M}_{\mathrm{TOT}}=\displaystyle \frac{{F}_{\nu }{D}^{2}}{{\kappa }_{\nu }{B}_{\nu }({T}_{{\rm{d}}})},\end{eqnarray} \tag{ 6 }$$

and see also Ward-Thompson & Whitworth (2011), where

$$\begin{eqnarray}&&{\kappa }_{\nu }={\kappa }_{o}{\left(\displaystyle \frac{\nu }{{\nu }_{o}}\right)}^{\beta },\end{eqnarray} \tag{ 7 }$$

and F_ν is the total measured flux density at the observed frequency ν, B_ν (T_d) is the Planck function for a dust temperature T_d (T_d = 12.1 K Planck Collaboration et al. 2016c), and κ_ν is the monochromatic opacity per unit mass of dust and gas. κ_ν ∼ 0.0125 cm² g⁻¹ assuming κ_o = 0.1 cm² g⁻¹, ν_o = 10¹² Hz (Beckwith et al. 1990) and β = 2. We should note that κ_ν can have a systematic uncertainty of up to 50% (Roy et al. 2014).

Assuming a distance of 125 pc, we estimate total masses from the 850 μm dust emission of ${M}_{\mathrm{TOT}}^{\mathrm{NW}}\sim 0.10\,{M}_{\odot }$ and ${M}_{\mathrm{TOT}}^{\mathrm{SE}}\sim 0.12\,{M}_{\odot }$ for the northwest and southeast lobes, respectively, and ${M}_{\mathrm{TOT}}^{\mathrm{main}}\sim 0.37\,{M}_{\odot }$ for the main core. This suggests that if they are indeed fragmented, the central part of the starless core may be undergoing gravitational collapse (i.e., ${M}_{\mathrm{TOT}}^{\mathrm{main}}$/${M}_{\mathrm{BE},\mathrm{crit}}^{\mathrm{main}}\gt 1$) while the two smaller lobes are not, rather than a coherent collapse of the whole core. All of the masses are summarized in Table 1.

We also estimated the envelope mass from the 850 μm dust emission of the two T-Tauri sources RNO 90 and RNO 91 using Equation (6). We find a total estimated envelope mass for RNO 90 of ${M}_{\mathrm{TOT}}^{\mathrm{RNO}\ 90}$ = 0.033 ± 0.016 M_⊙ assuming a distance of 115.7 ± 1 pc and a radius of 144. This is orders of magnitude greater than the dust mass of the disk as seen by the Atacama Large Millimeter/submillimeter Array (ALMA), which is closer to 2 × 10⁻⁵ M_⊙ (Garufi et al. 2022), but we do not resolve this structure and are more likely seeing the remaining dusty envelope. For RNO 91, we estimate a total mass from the 850 μm dust emission of ${M}_{\mathrm{TOT}}^{\mathrm{RNO}\ 91}$ = ∼0.335 ± 0.169 M_⊙ assuming a distance of 125 pc and that the dusty envelope is an ellipse with dimensions 320 × 180 rotated 60° east of north. This estimated total mass value is in good agreement with Young et al. (2006) who found a mass of 0.3 ± 0.1 M_⊙. Assuming just a uniform sphere of radius 156, we get an estimated total mass of 0.215 ± 0.109 M_⊙.

As mentioned in Section 2, there is a weak CO outflow driven by the embedded Class I protostar in RNO 91. The southern outflow traces the southern edge of the L43 starless core and forms a limb-brightened U shape (Lee et al. 2002), which is seen in all transitions. The southern outflow is heavily blueshifted, indicating the outflow is tilted toward us, potentially by up to 60° (Lee & Ho 2005). Weintraub et al. (1994) found that RNO 91 is not very deeply embedded in the L43 molecular cloud and rather sits nearer to us than the main submillimeter starless core.

However, a very clear dust cavity can be seen in Figure 2, which the outflow traces nearly perfectly. This dust cavity sits along the filament, and it appears that the filament has been disrupted by the outflow, as material is cleared to the south and potentially pushed north to form the kink in the filament, though there is not much redshifted CO emission to the north. This morphology of the dust in the filament suggests some sort of interaction with or influence by the outflow. This does contradict the above claim that the source is not deeply embedded. The 850 μm dust emission also shows this U-shaped bend to the south, suggesting even the densest part of the filament is affected by, or was initially affected by, the outflow. Bence et al. (1998) did suggest that the outflow has been weakened over time by a UV radiation field, so the current outflow we observe may not be the original morphology or strength.

One possibility then is that the source was previously embedded and cleared out the dust cavity we see including along the LOS so that it presently sits in the foreground of the dense filament. This was also suggested by Mathieu et al. (1988) who determined that RNO 91 was once associated with the dense molecular core, but has since blown through the dense gas with the outflow. They also suggested that the outflow energy is only coupled with a small fraction of the core mass. So the majority of the dense starless core L43 is undisturbed, though as discussed later, our observations of the magnetic field suggest that we are either only tracing affected foreground dust or that the outflow has influenced some of the dense material.

Regardless, the fact that the dust appears to be heavily influenced by the outflow suggests we must be careful in our analysis of the magnetic field, which is traced by the dust. A more in-depth discussion of the interaction of the outflow with the magnetic field and potential CO emission or polarization contribution is presented in Section 5.2.

We used archival Herschel Photodetector Array and Camera Spectrometer (PACS) 160 μm, SPIRE 250, 350, and 500 μm dust emission maps, ¹⁰⁶ along with the JCMT 850 μm dust emission map from this work to create column density maps. We filtered the Herschel maps in order to remove the large-scale structure that SCUBA-2/POL-2 is not sensitive to. We followed the method from Sadavoy et al. (2013) and Pattle et al. (2015) of introducing the Herschel maps into the Stokes I time stream and repeating the reduction process from Section 3.2, using 4'' pixels. We then subtracted the original 850 μm only SCUBA-2/POL-2 Stokes I emission from the map, which included the Herschel maps in the reduction and the resulting map was the filtered Herschel map.

We also attempted to "correct" the 850 μm Stokes I maps by removing potential contamination from the ¹²CO (J = 3–2) line. The ¹²CO (J = 3–2) line, which sits at 345.796 GHz, is within the SCUBA-2 850 μm bandpass filter, so it can contribute to total intensity continuum. We followed the method of Parsons et al. (2018) using the HARP data mentioned in Section 3.4. We used the regular 4'' Stokes I map because this correction method is best characterized for 4'' maps before and for the purposes of fitting the blackbody spectrum, we do not require the increase in S/N that is helpful for our polarization vectors. The contribution to total intensity from CO was ∼5%–10%, getting up to ∼20% directly around RNO 91. We should note that the reduction produced slight negative bowling to the north of the L43 emission, though not in a region of any emission.

We then fit the five maps with a modified blackbody (Hildebrand 1983)

$$\begin{eqnarray}&&{F}_{\nu }={\mu }_{{{\rm{H}}}_{2}}{m}_{{\rm{H}}}{N}_{{{\rm{H}}}_{2}}{B}_{\nu }({T}_{{\rm{d}}}){\kappa }_{\nu },\end{eqnarray} \tag{ 8 }$$

where again F_ν is the measured flux density at the observed frequency ν, B_ν (T_d) is the Planck function for a dust temperature T_d, ${\mu }_{{{\rm{H}}}_{2}}$ is the mean molecular weight of the hydrogen gas in the cloud, m_H is the mass of an hydrogen atom, ${N}_{{{\rm{H}}}_{2}}$ is the column density, and κ_ν is the dust opacity (see Equation (7)). We use a value of 2.8 for ${\mu }_{{{\rm{H}}}_{2}}$, and κ_ν was calculated for each frequency observed using Equation (7), where β is the emissivity spectral index of the dust and is taken to be 1.8 (an approximate value in starless cores; Shirley et al. 2005; Schnee et al. 2010; Sadavoy et al. 2013), and we again assume κ_o = 0.1 cm² g⁻¹ and ν_o = 10¹² Hz (Beckwith et al. 1990). We used temperature values from previously derived dust temperature maps using just the nonfiltered SPIRE maps and 850 μm maps, where we had convolved the data to the 500 μm resolution of ∼35'' and then regridded to the 850 μm maps.

We see column densities in the main starless core on the order of 10^22.8 cm⁻², which is ∼6 × 10²² cm⁻², with a maximum column density of ∼3 × 10²³ cm⁻².

In Figure 4 we plot polarization fraction versus intensity of the nondebiased polarization half-vectors in L43. We focus only on the central 3' region of L43, as this is where the exposure and noise are roughly uniform. A very clear decrease in polarization fraction can be seen toward the regions of high intensity. Within starless cores, this depolarization occurs in the highest density regions, due to some combination of field tangling and the loss of grain alignment at high enough A_V 's (≈20 mag) as predicted by RAT theory (Andersson et al. 2015). A common method to study the grain alignment efficiency in molecular clouds is to determine the relationship between polarization efficiency and visual extinction, where polarization fraction and total intensity can be substituted for those two quantities, respectively, at submillimeter wavelengths (see Pattle et al. 2019 and references therein).

Zoom In Zoom Out Reset image size

The relationship between polarization and intensity should follow a power law, p ∝ ${A}_{V}^{-\alpha }$, where an α of 1 indicates a loss of alignment and an α of 0 would indicate perfect alignment. We follow the methods of Pattle et al. (2019) and use the Ricean fitting technique to fit the data. We get α = 0.83 ± 0.06 for the 12'' vectors and see an obvious offset from the null, which would indicate we retain some alignment. Additionally, the ordered polarization geometry suggests that we are continuing to trace the magnetic field to high A_V 's. We performed the Ricean fitting for polarization vectors binned from 8'' up to 32'' and see α values from 0.89–0.70, respectively.

The vectors chosen for analysis have an S/N cut of I/δ_I > 10 and p/δ_p > 2. An S/N cut of p/δ_p > 2 can be quite poor in polarization so we must proceed with caution with those vectors. In Figure 5 we plot the lower-S/N vector distribution (dashed histogram) and the higher-S/N vectors with p/δ_p > 3 (solid histogram). Within most of the molecular cloud, the two S/N cuts agree well with the lower-S/N vectors following the same orientation as the higher-S/N vectors. The polarization angle distributions follow the same shape between the two S/N cuts, and they agree well with that found by Matthews et al. (2009) in the SCUPOL legacy survey (blue histogram). Using the lower-S/N cut, we get more data points to then use when calculating magnetic field strength (see Section 4.6), which could increase the spread of the position angles but can also increase the statistical confidence in the calculated dispersion.

Zoom In Zoom Out Reset image size

As can be seen in Figure 5, there is no clear single morphology of the magnetic field, and it instead must be considered as either randomized or a multiple-component field. There is a rather distinct peak around 150° with then more scatter toward the lower magnetic field polarization angles. Some of this scatter is structured fields in other parts of the molecular cloud. As will be discussed later in Section 5.2, we suspect that the magnetic field is partially influenced by the CO outflow from RNO 91. This was discussed as well in Ward-Thompson et al. (2000), where they suggested the western edge of the field they observed was being influenced by RNO 91. With the more sensitive POL-2 observations and a larger field of view, we can actually see the overlap of the CO emission with some of the magnetic field vectors.

We split the magnetic field inferred from the 850 μm polarized emission of L43 into three parts, the two labeled regions seen in Figure 6 and then the magnetic field vectors, which spatially (in the plane-of-sky) overlap with the CO outflow or are nearby and follow the same orientation. We list the mean field orientation, 〈θ_B〉, and standard deviation from a Gaussian fit of the magnetic field position angle distributions in Table 2. Regions 1 and 2 show different magnetic field orientations, although both are rather scattered. Region 2, which corresponds to the northern half of the starless core, has a magnetic field that has an average orientation of 63° east of north, which is roughly parallel to the filament (≈67°) and Planck magnetic field orientations (≈60°). It also lies roughly perpendicular to the local core elongation axis, which is something seen across starless and prestellar cores (see Pattle et al. 2022a, for a recent review). This is further discussed in the context of the region's evolution in Section 5.3. There is more scatter toward the center of the region, which is what causes the spread we see in the position angles, but the structured component can be seen on either side.

Zoom In Zoom Out Reset image size

Table 2. DCF+ADF Values

Region	Outflow	Reg 1	Reg 2
〈θB〉 a (°)	146(22)	140(25)	63(24)
R.A. (J2000)	16:34:34.9	16:34:41.7	16:34:35.6
Decl. (J2000)	−15:47:30.0	−15:47:49.8	−15:46:37.1
a ('')	62.9	32.3	63.0
b ('')	17.4	24.8	28.8
PA b (°)	130	0	120
c ('')	33.1	25.7	42.6
N(H2) (×1021 cm−2)	33.0(7.0)	4.0(1.0)	52.0(9.0)
n(H2) (×105 cm−3)	4.0(0.9)	0.57(0.15)	4.9(0.8)
ΔvNT c (km s−1)	0.35(0.02)	0.35(0.02)	0.35(0.02)
b (°)	14.0(1.5)	15.3(6.2)	23.8(4.0)
σθ (°)	10.1(1.1)	11.0(4.5)	17.6(3.0)
Bpos (μG)	116(20)–257(42)	40(17)–88(38)	73(14)–162(32)
${{ \mathcal M }}_{{\rm{A}}}$	0.4(0.04)–0.9(0.1)	0.4(0.2)–1.0(0.4)	0.7(0.1)–1.6(0.3)
vA (km s−1)	0.3(0.04)–0.7(0.1)	0.3(0.1)–0.7(0.3)	0.2(0.03)–0.4(0.1)
λ	⋯	0.3(0.2)–0.8(0.4)	2.4(0.6)–5.4(1.4)
EB (×1035 J)	0.5(0.2)–2.6(0.9)	⋯	0.5(0.2)–2.2(0.9)

Notes.

^a Standard deviation of the Gaussian fit is in parentheses. ^b PA of ellipses is counterclockwise from north. ^c Δv_NT values from Caselli et al. (2002).

Download table as:  ASCIITypeset image

Region 1 has a slightly more coherent magnetic field structure that is oriented ≈140° east of north, nearly perpendicular to the filament direction and parallel to the CO outflow. Considering it is still near to the CO outflow and is a less-dense region, the magnetic field could still be influenced by the CO outflow, or we are simply seeing another component of the complex magnetic field. The mean field orientation of the vectors spatially overlapping with the CO outflow is 146°, which is well aligned with the outflow direction, which we have taken to be ∼150° ± 10° due to it curving slightly.

We also detect a few B-field vectors in the dust envelope of RNO 90 and in a very diffuse "blob," isolated to the west. RNO 90 is shown in the inset of Figure 6, and the magnetic field is oriented roughly north–south. The magnetic field in the diffuse blob to the west appears to still follow the large-scale Planck field, something that has been seen in diffuse cores (Ward-Thompson et al. 2023) and other isolated starless cores (L1689B; Pattle et al. 2021). The fact that this more diffuse region still follows the Planck field while Region 1 does not suggests that Region 1 may indeed be, or have been, affected by the outflow.

We estimated the magnetic field strength in L43 using the Davis–Chandrasekhar–Fermi (DCF) method (Davis 1951; Chandrasekhar & Fermi 1953). The DCF method (see Equation (11)) assumes that the geometry of the mean magnetic field is uniform in each region. It then assumes that deviations from this uniformity are Alfvénic such that the deviations are due to nonthermal gas motions. The Alfvénic Mach number of the gas is given by

$$\begin{eqnarray}&&{{ \mathcal M }}_{{\rm{A}}}=\displaystyle \frac{{\sigma }_{\mathrm{NT}}}{{v}_{{\rm{A}}}}=\displaystyle \frac{{\sigma }_{\theta }}{Q},\end{eqnarray} \tag{ 9 }$$

where the nonthermal deviations are quantified by a dispersion in magnetic field position angles, σ_θ . σ_NT is the 1D nonthermal velocity dispersion of the gas and Q is a correction factor that accounts for variations of the magnetic field on scales smaller than the beam and along the LOS where 0 < Q < 1 (Ostriker et al. 2001). v_A is the Alfvén velocity of the magnetic field, and so ${{ \mathcal M }}_{{\rm{A}}}\lt 1$ suggests the magnetic field is more important than turbulent motions (sub-Alfvénic) while ${{ \mathcal M }}_{{\rm{A}}}\gt 1$ means the turbulent motions are more important (super-Alfvénic). The Alfvén velocity is given by

$$\begin{eqnarray}&&{v}_{{\rm{A}}}=\displaystyle \frac{B}{\sqrt{4\pi \rho }}=Q\displaystyle \frac{{\sigma }_{\mathrm{NT}}}{{\sigma }_{\theta }},\end{eqnarray} \tag{ 10 }$$

where B is the magnetic field strength and ρ is the gas density. Since the dispersion in position angles, σ_θ , is for POS observations, we can only calculate the plane-of-sky magnetic field strength, B_pos, which is given by

$$\begin{eqnarray}&&{B}_{\mathrm{pos}}\approx Q\sqrt{4\pi \rho }\displaystyle \frac{{\sigma }_{\mathrm{NT}}}{{\sigma }_{\theta }}.\end{eqnarray} \tag{ 11 }$$

This can then be simplified to

$$\begin{eqnarray}&&{B}_{\mathrm{pos}}(\mu {\rm{G}})\approx 18.6\,Q\sqrt{n({{\rm{H}}}_{2})({\mathrm{cm}}^{-3})}\displaystyle \frac{{\rm{\Delta }}{v}_{\mathrm{NT}}(\mathrm{km}\,{{\rm{s}}}^{-1})}{{\sigma }_{\theta }(\mathrm{degree})}.\end{eqnarray} \tag{ 12 }$$

Typically Q is taken to be 0.5 (see Ostriker et al. 2001; Crutcher et al. 2004) but we will consider a range of Q values, 0.28 < Q < 0.62 from Liu et al. (2022; see their Table 3) to obtain upper and lower limits of the B-field strength. Then n(H₂) is the volume density of molecular hydrogen where n(H₂) = ρ/${\mu }_{{{\rm{H}}}_{2}}{m}_{{\rm{H}}}$ and ${\mu }_{{{\rm{H}}}_{2}}$ = 2.8 and m_H is the mass of hydrogen. Δ v_NT is the FWHM of the nonthermal gas velocity calculated by Δ v_NT = ${\sigma }_{\mathrm{NT}}\sqrt{8\mathrm{ln}2}$. As mentioned above, σ_θ is the dispersion of the position angles of the magnetic field vectors, which we calculated using an angular dispersion function (ADF) as discussed later in this section. It should be noted that Crutcher et al. (2004) found on average B_pos/B ≈ π/4, but since this is a general statistical correction, we do not use this when calculating the magnetic field strength.

We can rewrite Equations (9) and (10) as

$$\begin{eqnarray}&&{{ \mathcal M }}_{{\rm{A}}}=1.74\times {10}^{-2}\sqrt{2}\ \displaystyle \frac{{\sigma }_{\theta }(\mathrm{degree})}{Q},\end{eqnarray} \tag{ 13 }$$

and

$$\begin{eqnarray}&&{v}_{{\rm{A}}}(\mathrm{km}\,{{\rm{s}}}^{-1})=24.2Q\ \sqrt{2}\ \displaystyle \frac{{\rm{\Delta }}{v}_{\mathrm{NT}}(\mathrm{km}\,{{\rm{s}}}^{-1})}{{\sigma }_{\theta }(\mathrm{degree})},\end{eqnarray} \tag{ 14 }$$

respectively. We have also included a $\sqrt{2}$ factor in Equations (13) and (14), which is suggested by Heiles & Troland (2005) to account for the velocity line width assumptions, specifically converting the 1D LOS velocity measurements to an approximate value suitable for estimating the POS magnetic field strength.

We calculated the magnetic field strength in Regions 1 and 2 (shown in Figure 6) using Equation (12). We treated the regions as ellipses with semimajor and semiminor axes a and b (see Table 2), and assumed the depth of those regions to be the geometric mean, $c=\sqrt{{ab}}$. We used column density values from Figure 3 to calculate the volume density, n(H₂), in each region. We used N₂H⁺ (1–0) velocity line profiles from Caselli et al. (2002), which have a resolution of ∼0.063 km s⁻¹. We corrected them to account for the thermal component (since Equation (12) uses nonthermal velocity line widths), which was calculated using the excitation temperature, T_ex = 7 ± 1 K (also from Caselli et al. 2002), giving 0.35 ± 0.02 km s⁻¹. It should be noted that the velocity line profile observations are from the main starless core (Region 2). These observations do not necessarily extend to Region 1, but we do not have observations of Region 1 specifically so elect to use the same line width value as Region 2. Line widths of other tracers in the main starless core vary with some larger than and some smaller than the 0.35 km s⁻¹ value we use (see Chen et al. 2009), but N₂H⁺ (1–0) traces dense regions of molecular clouds, which should coincide with the depths we are observing at 850 μm as well. There are also NH₃ observations of the starless core from Jijina et al. (1999) and Fehér et al. (2022), with a spread of line width values from 0.273 km s⁻¹ (HFS fitting from Fehér et al. 2022) to 0.718 km s⁻¹ (Gaussian fitting from Fehér et al. 2022) and then 0.32 km s⁻¹ (Jijina et al. 1999). Considering the largest line width, the B-field strength could be up to two times larger.

We determined the dispersion of position angles in each region using the ADF (Hildebrand et al. 2009; see their Equations (1) and (3)). This assumes that there is a large-scale structured field and a smaller turbulent or random component. We assume that the length scales we are observing with JCMT/POL-2 (ℓ) are greater than the turbulent correlation length and much smaller than the large-scale field (such as Planck). While the former statement may be more difficult to accurately determine, the latter is true in our situation as we see a very structured large-scale B-field from Planck (see Figure 2) with no variations in the region we are looking at (admittedly a single Planck beam nearly covers the whole region). The contribution of both the turbulent and large-scale components to the angular dispersion of B-field vectors 〈ΔΦ²(ℓ)〉^1/2 is given by the terms b and m ℓ, respectively, and the relation (Hildebrand et al. 2009),

$$\begin{eqnarray}&&\langle {\rm{\Delta }}{{\rm{\Phi }}}^{2}({\ell }){\rangle }_{\mathrm{tot}}\simeq {b}^{2}+{m}^{2}{{\ell }}^{2}+{\delta }_{\theta }^{2}({\ell }),\end{eqnarray} \tag{ 15 }$$

where ${\delta }_{\theta }^{2}({\ell })$ is the additional contribution to the dispersion from measurement uncertainties (see Equation (4)). In each region, we calculated 〈ΔΦ²(ℓ)〉_tot and then fit Equation (15) to determine values for m and b. The plots are seen in Figure 7 where the best-fit parameters are shown as well. In all three regions, limiting the fitting to the first three bins (36'') provides the best fit (as determined by χ² < 1), though for the outflow region, extending the fit to 48'' still provided χ² < 1. For this region we then take the mean of the b values from both fits (14.2 for 36'' and 13.7 for 48'') to get b ≈ 14.0 ± 1.5. Generally we see the dispersion slowly increase with distance (Hildebrand et al. 2009; Hwang et al. 2023; see also the top panel of Figure 7), but in the bottom two panels of Figure 7 we see the dispersion growing and then at a large distance, suddenly drop again. In the case of Region 2 (bottom panel of Figure 7), the distance this happens at, ∼120'', is the distance between the two structured components mentioned in Section 4.5, further justifying a structured field approximately oriented at 56° and parallel with the Planck field and the filament direction. The dispersion in magnetic field position angles, σ_θ , can then be calculated by

$$\begin{eqnarray}&&{\sigma }_{\theta }=\displaystyle \frac{b}{\sqrt{2-{b}^{2}}}\displaystyle \frac{180^\circ }{\pi },\end{eqnarray} \tag{ 16 }$$

where b (in radians) is found from the ADF fit (see Hildebrand et al. 2009; see also Figure 7).

Zoom In Zoom Out Reset image size

The magnetic field strength in Regions 1 and 2 is ∼40–90 μG and ∼70–160 μG, respectively. All of the values calculated when considering the magnetic field strength are listed in Table 2. In our situation, the difference in magnetic field strength between regions is due to variations in density and angular dispersion of the vectors since we use a constant velocity line width value across the region. For example, despite the larger σ_θ in Region 2, it is also denser, which is what increases the magnetic field strength and makes it comparable, if not greater than, the magnetic field in Region 1. We calculate an upper and lower bound for each of our magnetic field strengths based on the variation in Q of 0.28 < Q < 0.62 (Liu et al. 2022). The value in Region 2 is approximately that found by Crutcher et al. (2004), although they treated the region as a sphere and had a smaller dispersion angle of 12°. We find a slightly larger column and volume density, by a factor of ∼1.3.

The magnetic field strength calculated for the outflow is ∼120–260 μG. However, we note that this value may be severely overestimated and its use for interpretation limited. This is because we do not generally apply the DCF method to regions interacting with outflows since we assume the deviations in the magnetic field to be small nonthermal gas motions in the region, and an outflow is a much stronger disruptive force. In our case, we are arguing that some of the dense gas and dust has been affected by the outflow and hence that the outflow may have dragged the field (which is flux-frozen into the gas), aligning it with the outflow walls and giving us the low position angle dispersion. In that case, the low position angle dispersion may be due to a strong outflow rather than a strong field; although this would require more outflow modeling to determine if the outflow would preferentially align and order the field, or disorder the field. Additionally, we consider all of the vectors spatially aligned with the outflow, but take just the dust density in the southern half of the L43 starless core. If we were to assume a much lower volume density, one that is more likely associated with outflow material, then our field strength would be lower.

We acknowledge the limitations of using the DCF method, especially when determining the position angle dispersion and acknowledge in such a low-S/N environment that these uncertainties are increased. However, we do find magnetic field strengths that are on the order of strengths seen in other starless cores (Karoly et al. 2020; Pattle et al. 2021), including values that agree within error with those found in Crutcher et al. (2004).

The mass-to-flux ratio λ is a parameter that can be used to quantify the importance of magnetic fields relative to gravity (Crutcher 2004). It compares the critical value for the mass, which can be supported by the magnetic flux Φ, ${M}_{{B}_{\mathrm{crit}}}$ = Φ/2 $\pi \sqrt{G}$ (Nakano & Nakamura 1978) to the observed mass and flux values. If the column density N and magnetic field strength B can be measured, the observed value for the ratio between mass and flux is (M/Φ)_obs = mNA/BA and then the mass-to-flux ratio λ is,

$$\begin{eqnarray}&&\lambda =\displaystyle \frac{{(M/{\rm{\Phi }})}_{\mathrm{obs}}}{{(M/{\rm{\Phi }})}_{\mathrm{crit}}}=7.6\times {10}^{-21}\displaystyle \frac{{N}_{{{\rm{H}}}_{2}}({\mathrm{cm}}^{-2})}{{B}_{\mathrm{pos}}(\mu {\rm{G}})}\end{eqnarray} \tag{ 17 }$$

where m = 2.8 m_H, ${N}_{{{\rm{H}}}_{2}}$ is the molecular hydrogen column density and B_pos is the plane-of-sky magnetic field strength. When λ < 1, then the magnetic field is strong enough to support against gravity; this is referred as a "magnetically subcritical" regime. Alternatively, if λ > 1, then the magnetic field is insufficient by itself to oppose gravity, and the cloud is instead "magnetically supercritical."

Using Equation (17), we calculate the mass-to-flux ratio to be 0.3–0.8 in the low-density periphery of the core, Region 1. We calculate a mass-to-flux ratio of 2.4–5.4 in the denser part of the core, Region 2. According to Crutcher (2004) these ratios can be overestimated due to geometric biases, and they suggest it can be overestimated by a factor up to 3, although it is a statistical correction and its application to individual measurements is unclear. If we consider this statistical correction, the mass-to-flux ratios are ∼0.1–0.3 and ∼0.8–1.8 in Regions 1 and 2, respectively. These are then the lower limits for the mass-to-flux ratio in each region.

We also get Alfvén Mach numbers of 0.4–1.0 and 0.7–1.6 for Regions 1 and 2, respectively, suggesting that both regions are roughly transcritical, and magnetic field and turbulence may play equal parts in support. As noted previously, the velocity information we are using is for the main starless core, which is near Region 2, but farther from Region 1, which is the low-density area on the periphery of the main dense core; it does not resolve individual parts of the molecular cloud. We also find Alfvén velocities in the range of 0.2–0.7 km s⁻¹ throughout the cloud.

For the lower-density Region 1, the region is entirely magnetically subcritical, indicating that the region may still be sufficiently supported against gravitational collapse by the magnetic field. It is also slightly sub-Alfvénic, meaning the magnetic field may play the dominant role. In the denser Region 2, we obtain a more definitively supercritical mass-to-flux ratio (note the large uncertainties with this value, as well as the results of the statistical correction), with a lower limit approaching trans- to subcritical. This gradient of a subcritical envelope (Region 1) transitioning into a trans- to supercritical core (Region 2) at sufficient densities is described in Crutcher (2004). It does suggest that the magnetic field is not sufficiently strong to support against gravitational collapse in the main starless core. However, it is worth mentioning that there are still many other processes in the molecular cloud such as turbulence and the influence of the CO outflow that could prevent gravitational collapse into a stellar object. While we do not yet see a protostar-forming, as mentioned in Section 4.1, we do see some possible fragmentation within the starless core where the densest part may be undergoing gravitational collapse.

Myers (2017) suggested from modeling that L43 has formed all of the stars it will form in its lifetime and does not contain sufficient amounts of dense gas for further star formation. However, we find higher column density values than they use for their modeling and do see potential fragmentation in the core. Chen et al. (2009) found that the main L43 starless core has observed DCO⁺ and HCO⁺ abundances that are higher and lower than modeled abundances, respectively, for an assumed amount of CO depletion. They suggest this indicates more CO depletion in the core and that the L43 starless core is spending a longer time at the higher-density pre-protostellar core phase. If this is the case, additional supports such as turbulence may be needed to continue support against gravitational collapse in Region 2 since it seems to be moving beyond the stage where magnetic fields are critical. But findings of Region 2 in near equilibrium (within errors) also validates the long-lived age of the core that Chen et al. (2009) found. The local magnetic field appears to still be significant in the more diffuse Region 1, but this region is beyond the area considered by Chen et al. (2009). Additionally, we must consider that the CO outflow has potentially altered the structure of the magnetic field in the cloud, even within the starless core, potentially weakening or strengthening the field.

The spatial alignment (in the plane-of-sky) of the magnetic field in L43 with the outflow cavity walls can be seen in Figure 6, where the cyan outflow spatially overlaps with many of the magnetic field vectors. The magnetic field vectors that coincide with the CO outflow show a uniform distribution and strong peak at 146°. This coincides well with the outflow direction, which we have taken to be ∼150° ± 10° due to it curving slightly.

Weintraub et al. (1994) suggested that RNO 91 sits in the foreground of the general L43 molecular cloud and so the possibility exists that there is in fact no physical association between the magnetic field and outflow. However, as was mentioned in Section 4.1 and as can be clearly seen in Figure 2, the outflow appears to have carved a cavity out of the molecular cloud, indicating it is to some degree embedded. This was also suggested by Mathieu et al. (1988).

Alternatively, we may be tracing magnetic fields in the outflow cavity walls. In this case, as mentioned in Section 4.3, some of the Stokes I emission we see, especially in the regions coincident with the outflow cavity walls, may be CO features (Drabek et al. 2012), especially the isolated emission to the south of the main cloud. We expect some contribution to the measured Stokes I emission from CO, but such contributions are typically less than 20% of the total emission observed (Drabek et al. 2012; Pattle et al. 2015; Coudé et al. 2016) and in our case, we see contributions of ∼5%–15%. However, considering we do have some CO contribution, we cannot rule out the possibility that some fraction of the polarized emission in this region arises from CO polarization, polarized through the Goldreich-Kylafis effect (Goldreich & Kylafis 1981, 1982). This would add a further ±90° ambiguity on the magnetic field orientation.

On the other hand, we can assume the polarization and emission are not purely CO based, as the emission features are also seen in all of the Herschel bands and persist after CO subtraction in 850 μm, indicating that there is a real dust feature present. A similar "hollow shell" morphology of dust emission in the presence of outflows is discussed in Moriarty-Schieven et al. (2006). Additionally, Bence et al. (1998) suggested that the mere presence of CO suggests some amount of dust shielding (from the UV field) in the outflow region. Thus, we still consider it probable that we are tracing dust polarization in the outflow cavity walls.

Additionally, the relationship between magnetic fields with outflows has been observed on numerous occasions in other sources, on both JCMT and ALMA scales (see Hull et al. 2017, 2020; Lyo et al. 2021; Pattle et al. 2022b). Hull et al. (2017, 2020) and Pattle et al. (2022b) also saw a similar alignment between the magnetic field and the cavity wall of the outflow to that which we see in L43/RNO 91. The benefit of our larger field of view here, when compared to ALMA, is that we can compare the magnetic field of the outflow to that in the surrounding regions and see that there is not just a preferential direction northwest to southeast, but rather that the magnetic field in the outflow region is actually different to that in the rest of the cloud. It should be noted that in regions observed by the JCMT, a preferred misalignment of 50° between magnetic fields and outflows has previously been identified in a larger statistical sample (Yen et al. 2021). The outflow observations in that study are largely on envelope- or small scales (i.e., tens of arcseconds) rather than large-scale outflows like we see here. So while there is a statistically preferred misalignment between magnetic fields as observed by JCMT and outflows, we do not see a clear indication of this occurring in L43 and we instead see good alignment between outflow and magnetic fields. This is perhaps because we have such distinct large-scale cavity outflow walls, which is what JCMT may be preferentially tracing. RNO 91 would also be interesting to follow up with ALMA polarization observations since there are smaller-scale outflows in the envelope as well (Lee & Ho 2005; Arce & Sargent 2006). This could be more directly compared to observations by Hull et al. (2017, 2020) and to the statistical sample in Yen et al. (2021).

In RNO 91, the CO outflow was found to have a lower limit of its energy at ≈10²⁹ J (Bence et al. 1998). However, this assumes a Class II source lifetime when considering how long the CO has been exposed to UV radiation, though a Class I lifetime would still be longer than the un-shielded CO lifetime of a few hundred years (Bence et al. 1998). So the CO outflow likely had a larger energy once and may have been able to influence the magnetic field orientation. Other studies have found the CO outflow energy to be ∼5 × 10³⁵ J (Myers et al. 1988) and ∼1.4 × 10³⁵ J (Arce & Sargent 2006). These values are comparable to the magnetic energy values we see in L43, which is what we would expect. We can calculate the magnetic energy in Region 2 and the outflow region using

$$\begin{eqnarray}&&{E}_{{\rm{B}}}({\rm{J}})={10}^{-20}\displaystyle \frac{{B}^{2}(\mu {{\rm{G}}}^{2})V({{\rm{m}}}^{3})}{2{\mu }_{o}({\rm{N}}\,{{\rm{A}}}^{-2})}\end{eqnarray} \tag{ 18 }$$

where B is the magnetic field strength in microGauss (as measured in the plane-of-sky), V is the volume of the region in cubic meters, and μ_o is the permeability of free space, giving us the magnetic energy in Joules. Since we consider the outflow to have affected the dense gas and dust and dragged the magnetic field, we would expect the outflow energy to be at least equal to the magnetic energy. Assuming an ellipsoid shape when calculating the volume of both regions (see Table 2 for ellipse parameters), we find magnetic energies of ≈0.5–2.5 × 10³⁵ J. This suggests that we can help further place a lower limit on the outflow energy of ≈0.5–2.5 × 10³⁵ J, which is comparable to the values stated above, though we do remain cautious of the magnetic field strength derived in the outflow region as mentioned before.

One point of interest in this region is that there is a very clear evolutionary gradient from southwest to northeast. Initially there is the evolved T-Tauri star RNO 90, which is the oldest source and currently has no known large-scale outflows. It has also formed a protostellar disk (Pontoppidan et al. 2010). RNO 91, which is farther along the filament, is a protostellar source that drives the now familiar CO outflow. Then finally the starless core sits ≈10,000 au farther along the filament. The orientation of this filament is such that it extends roughly radially away from Sco OB2, with RNO 90 the closest to Sco OB2. The filament sits ≈42 pc away from Sco OB2 (in the plane-of-sky) assuming a general distance of 125 pc. While this may be merely a coincidence, it is interesting that the evolutionary track in such an isolated filament starts in the part of the filament pointing directly toward Sco OB2 (in the plane-of-sky). The Ophiuchus region and star formation within has previously been suggested to be shaped and driven by Sco OB2 (Loren 1989).

Additionally, we can picture two evolutionary scenarios for the starless core, scenarios that with present observations we cannot distinguish between and that may very well be happening at the same time. The starless core L43 has formed with its long axis parallel to the outflow cavity wall. This could suggest that material has been funnelled down the filament that is also parallel with the large-scale Planck field and is building up on the outflow cavity walls. Build-up has not occurred so readily along the western wall of the outflow cavity because there is less material available for accretion to the west since RNO 90 has already been formed. However, it could also be the case that the dense core already existed and fragmented to form both the starless core and RNO 91 and the starless core has not been compressed along the filament orientation by the outflow. It is difficult to differentiate between these two scenarios and the possibility of course remains that they could both be true, with an initial fragmentation that has become denser over time. Kim et al. (2020) suggested that the starless core is a "late" or chemically evolved starless core as determined by a high N(DNC)/N(HN¹³C) ratio and line detection in N₂D⁺. So the core may have formed at a similar time to RNO 91 but has since had its evolution slightly delayed due to injected turbulence by RNO 91 as well as less readily available material to form a star with.