<figure><img src="image_1.png" /></figure>

HAL
open science

Magnetic Fields of the Starless Core L 1512

Sheng-Jun Lin, Shih-Ping Lai, Kate Pattle, David Berry, Dan P. Clemens, Laurent Pagani, Derek Ward-Thompson, Travis J. Thieme, Tao-Chung Ching

► To cite this version:

Sheng-Jun Lin, Shih-Ping Lai, Kate Pattle, David Berry, Dan P. Clemens, et al., Magnetic Fields of the Starless Core L 1512. The Astrophysical Journal, 2024, 961, 10.3847/1538-4357/ad0cc1. insu-04822495

HAL Id: insu-04822495

https://insu.hal.science/insu-04822495v1

Submitted on 6 Dec 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Magnetic Fields of the Starless Core L 1512

Sheng-Jun Lin<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> Ⓓ, Shih-Ping Lai<fnref n="1" />,<fnref n="3" /> Ⓓ, Kate Pattle⁴ Ⓓ, David Berry<fnref n="5" /> Ⓓ, Dan P. Clemens<fnref n="6" /> Ⓓ, Laurent Pagani<fnref n="7" /> Ⓓ, Derek Ward-Thompson⁸ Ⓓ, Travis J. Thieme<fnref n="2" />,<fnref n="3" /> Ⓓ, and Tao-Chung Ching<fnref n="9" /> Ⓓ

Received 2022 December 24; revised 2023 November 13; accepted 2023 November 13; published 2024 January 18

## Abstract

We present JCMT POL-2 850 µm dust polarization observations and Mimir H-band stellar polarization observations toward the starless core L 1512. We detect the highly ordered core-scale magnetic field traced by the POL-2 data, of which the field orientation is consistent with the parsec-scale magnetic fields traced by Planck data, suggesting the large-scale fields thread from the low-density region to the dense core region in this cloud. The surrounding magnetic field traced by the Mimir data shows a wider variation in the field orientation, suggesting there could be a transition of magnetic field morphology at the envelope-scale. L 1512 was suggested to be presumably older than 1.4 Myr in a previous study via time-dependent chemical analysis, hinting that the magnetic field could be strong enough to slow the collapse of L 1512. In this study, we use the Davis-Chandrasekhar-Fermi method to derive a plane-of-sky magnetic field strength ($B_{pos}$) of $18 \pm 7$ µG and an observed mass-to-flux ratio ($\lambda_{obs}$) of $3.5 \pm 2.4$, suggesting that L 1512 is magnetically supercritical. However, the absence of significant infall motion and the presence of an oscillating envelope are inconsistent with the magnetically supercritical condition. Using a virial analysis, we suggest the presence of a hitherto hidden line-of-sight magnetic field strength of $\sim 27$ µG with a mass-to-flux ratio ($\lambda_{tot}$) of $\sim 1.6$, in which case both magnetic and kinetic pressures are important in supporting the L 1512 core. On the other hand, L 1512 may have just reached supercriticality and will collapse at any time.

*Unified Astronomy Thesaurus concepts: Interstellar magnetic fields (845); Interstellar medium (847); Molecular clouds (1072); Polarimetry (1278); Submillimeter astronomy (1647); Star forming regions (1565); Star formation (1569); Starlight polarization (1571)*

## 1. Introduction

Magnetic fields (B-fields) may play a key role in the formation of starless cores and the evolution of protostellar cores. In a magnetic-field-dominated scenario, the magnetic field tends to be uniform, and a starless core will collapse into a flattened structure (a so-called pseudo-disk) and pull the magnetic field into an hourglass shape with its symmetry axis aligned with the minor axis of the pseudo-disk and with the rotation/outflow axes (e.g., Shu et al. 1987; Galli & Shu 1993a, 1993b; Li & Shu 1996; McKee & Ostriker 2007). On the other hand, magnetohydrodynamic (MHD) simulations suggest that turbulence can cause the magnetic field to become disordered and chaotic, and either compress the gas to form stars or dissipate gas to suppress star-forming activities (e.g., Padoan 1999; Mac Low & Klessen 2004; Federrath & Klessen 2012).

There are numerous observations of the magnetic field morphology in the protostellar stage. For scales larger than cores ($\sim 0.1$ pc), both hourglass-like patterns and chaotic B-

field patterns are found in the protostellar envelopes of young Class 0 and hot massive cores (e.g., Girart et al. 2006, 2009; Ching et al. 2017). At pseudo-disk scales ($\sim 0.01$ pc), Chapman et al. (2013) found a positive correlation between the mean B-field direction and the modeled pseudo-disk symmetry axis in six Class 0 sources that have outflows near the plane of the sky. On the other hand, Hull et al. (2013) showed that the angle between the B-field direction and outflow is likely statistically random at smaller scales of $\sim 1000$ au in a sample of 16 cores hosting Class 0 or I protostars. Hence, in order to understand the origin of the variation in the B-field morphology, it is crucial to measure dust polarization at the early, starless stage.

Unlike protostellar cores, the magnetic fields of starless cores are rarely measured due to their low polarized thermal emission intensity. Observations of background starlight polarimetry in the near-infrared (NIR) have been performed toward a few starless cores. These results showed significant nonuniformity in the plane-of-sky B-field structure around the cores (e.g., Clemens et al. 2016; Kandori et al. 2017, 2020a, 2020b). While background starlight polarimetry provides insights into these peripheral regions, it has limitations when it comes to probing the magnetic fields in the large-$A_V$ central regions. To address this, dust continuum emission polarimetry in the submillimeter wavelengths has emerged as a key tool for investigating the

<figure><img src="image_4.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Institute of Astronomy, National Tsing Hua University, No. 101, Section 2, Kuang-Fu Road, Hsinchu 30013, Taiwan; shengjunlin@asiaa.sinica.edu.tw
slai@phys.nthu.edu.tw</footnote>

<footnote n="2">² Academia Sinica Institute of Astronomy and Astrophysics, No. 1, Section 4, Roosevelt Road, Taipei 10617, Taiwan;</footnote>

<footnote n="3">³Center for Informatics and Computation in Astronomy, National Tsing Hua University, No. 101, Section 2, Kuang-Fu Road, Hsinchu 30013, Taiwan
⁴Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="5">⁵East Asian Observatory, 660 N. A'ohoktu Place, University Park, Hilo, HI 96720, USA</footnote>

<footnote n="6">⁶Institute for Astrophysical Research, Boston University, 725 Commonwealth Avenue, Boston, MA 02215, USA</footnote>

<footnote n="7">⁷LERMA & UMR8112 du CNRS, Observatoire de Paris, PSL University, Sorbonne Universités, CNRS, F-75014 Paris, France
⁸Jeremiah Horrocks Institute, University of Central Lancashire, Preston PR1 2HE, UK</footnote>

<footnote n="9">⁹Jansky Fellow, National Radio Astronomy Observatory, 1003 Lopezville Road, Socorro, NM 87801, USA</footnote>
magnetic fields in starless cores. However, starless cores, typically extended objects without compact substructures, tend to be entirely resolved out by interferometers (Dunham et al. 2016; Kirk et al. 2017; Caselli et al. 2019; Tokuda et al. 2020). As a result, single-dish telescopes, such as the James Clerk Maxwell Telescope (JCMT) and the Atacama Pathfinder Experiment (APEX) telescope, remain the preferred tools for observing starless cores.

Mapping the magnetic field morphology across entire cores is largely restricted to nearby star-forming regions due to low surface brightness and small core sizes. So far, only ten nearby individual starless cores have had 850/870 µm dust polarization detection: L 183, L 1544, L 43 (Ward-Thompson et al. 2000; Crutcher et al. 2004), L 1498, and L 1517B (Kirk et al. 2006) were observed with JCMT-SCUPOL; FeSt 1-457 (Alves et al. 2014) was observed with APEX-PolKa; and Oph C (Liu et al. 2019), IRAS 16293E, L 1689 SMM-16, and L 1689B (Pattle et al. 2021) were observed with JCMT-POL2. Their core sizes (the FWHM of the submillimeter intensity) range from approximately 0.04 to 0.1 pc, for distances spanning from approximately 110–150 pc (Kirk et al. 2005; Alves et al. 2014; Liu et al. 2019; Karoly et al. 2020; Pattle et al. 2021). With resolutions of 14" and 20" for JCMT and APEX, respectively, a spatial resolution of ~0.01 pc can be achieved at these distances to resolve core structure. In addition, a number of smaller starless cores embedded in clumps and filaments have also had polarization detection with JCMT-POL2 (Pattle & Fissel 2019). However, these studies mostly focused on either the role of B-fields of the parent clumps themselves (e.g., Oph A & B clumps; Kwon et al. 2018; Soam et al. 2018) or the mean B-field orientations of the cores compared to large-scale B-fields (e.g., L 1495 filament; Eswaraiah et al. 2021; Ward-Thompson et al. 2023), rather than the resolved B-field morphologies in the individual cores, due to sensitivity limitations. The protostellar cores have brighter surface brightness and internally complex substructures (e.g., disks, outflows), making them not just easier to be observed with single-dishes (e.g., Chapman et al. 2013; Pattle et al. 2021) but also more commonly observed with interferometers with much higher resolutions (see Hull & Zhang 2019, and references therein).

The aforementioned nearby isolated starless cores often exhibit magnetic fields that are relatively smooth and well-ordered (Pattle et al. 2023). Of the first three starless core detected with submillimeter polarizations (L 183, L 1544, and L 43; Ward-Thompson et al. 2000), Crutcher et al. (2004) smoothed these data from 14" to 21" and applied the Davis-Chandrasekhar-Fermi (DCF; Davis 1951; Chandrasekhar & Fermi 1953) method, but could not conclude whether the cores were magnetically subcritical (i.e., magnetically supported) due to the unknown orientations/inclinations of the cores. One unexpected result from the submillimeter single-dish observations is that the mean B-field direction and the core minor axes are not perfectly aligned. Ward-Thompson et al. (2009) found an angular offset of ~20°–50° among the five polarization-detected cores (Ward-Thompson et al. 2000; Kirk et al. 2006) available at that time. This is opposite to what magnetically regulated star formation models predict. Although this offset may be explained by the projection effect of tri-axial objects (Basu 2000), observations with a higher signal-to-noise ratio (S/N) are needed to critically compare to theoretical models. Recently, L 183 and L 43 were reobserved with JCMT-POL2

(Karoly et al. 2020, 2023), and the results suggest that L 183 is magnetically subcritical throughout the entire core, and L 43 is also subcritical. The Karoly et al. studies also demonstrate the sensitivity of JCMT-POL2 is greatly improved compared with the previous polarimeter JCMT-SCUPOL.

Another challenge for starless core polarimetry is that the power-law index $\alpha$ of the polarization-fraction to total intensity relation ($p \propto I^{-\alpha}$) has been found to be close to 1 (e.g., Alves et al. 2015; Liu et al. 2019), indicating that grain alignment with the magnetic field inside the core could be lost (e.g., Goodman et al. 1995; Andersson et al. 2015). Typically, the polarization fraction measurements are debiased with Gaussian noise and filtered with S/N criteria to fit the above single power-law relation in the conventional approach. However, Pattle et al. (2019) found that, in cases of low polarized intensity (e.g., starless cores) or partial grain alignment loss ($0 < \alpha < 1$), a higher-S/N threshold in Stokes *I*, which discards more data, is necessary to reliably recover the $\alpha$ index (see their Figure 1). Below this S/N threshold, the fitted $\alpha$ index will tend to be 1 because the actual polarization signal is weaker than the non-Gaussian-distributed noise, which has an approximate $I^{-1}$ dependence. Therefore, without a suitable S/N threshold, the conventional approach (fitting with a single power-law model of $\alpha$) would overestimate the $\alpha$ index. Pattle et al. (2019) found that $\alpha$ can be better estimated with a Ricean noise model without debiasing the polarization; i.e., favoring a nondebiased-polarization to total intensity relation ($p'-I$). The authors thus revised the $\alpha$ index of the starless core Oph C to $0.6 \sim 0.7$ from the value of 1.03 previously found by Liu et al. (2019). In FeSt 1-457, another starless core, $\alpha$ was also revised from 0.92 (Alves et al. 2015) to a lower value of 0.41 (Kandori et al. 2020c) using this method. Wang et al. (2019) also adopted the same noise model but used a Bayesian approach to revise the $\alpha$ index in IC 5146 from 1.03 ~ 1.08 to 0.56. In addition, an equivalent relation in NIR is $p/A_v \propto A_v^{-\alpha}$ (Andersson et al. 2015; Pattle et al. 2019). It has been found that $\alpha < 1$ in the dense clouds, even with debiased data (see Wang et al. 2017, 2019, and references therein). Therefore, dust grains remain aligned at higher densities than previously expected, allowing investigation of the magnetic fields within starless cores.

L 1512 (Lynds 1962) is a nearby isolated dark cloud harboring a starless core located near the edge of the Taurus–Auriga molecular cloud complex (Myers et al. 1983; Lombardi et al. 2010; Launhardt et al. 2013) at a nominal distance of 140 pc (Kenyon et al. 1994; Torres et al. 2009; Roccatagliata et al. 2020). Lin et al. (2020) modeled the physical and chemical structure of the L 1512 core using infrared dust extinction measurements and multiline observations of $N_2H^+$, $N_2D^+$, DCO$^+$, $C^{18}O$, and $H_2D^+$ (1$_{10}$–1$_{11}$) using non-LTE radiative transfer. With the high-density tracer $N_2H^+$, they found a low central temperature of $8 \pm 1$ K, a small 3D isotropic nonthermal velocity dispersion of $0.080 \text{ km s}^{-1}$ (corresponding to a 1D nonthermal velocity dispersion, $\sigma_{v,NT}$, of $0.046 \text{ km s}^{-1}$, a nonthermal FWHM linewidth, $\Delta v_{NT}$, of $0.109 \text{ km s}^{-1}$, and a Mach number, $M = \sigma_{v,NT}/c_s$, of 0.27, where $c_s$ is the isothermal sound speed of $0.17 \text{ km s}^{-1}$ at 8 K), and the absence of inward motions. Lin et al. (2020) concluded that L 1512 is chemically evolved and older than 1.4 Myr, suggesting that the dominant core formation mechanism could be a slow process, such as ambipolar diffusion (e.g., Mouschovias 1991; Tassis & Mouschovias 2004; Mouschovias et al. 2006). In addition, Falgarone et al. (2001) found that the
transverse velocity gradients (i.e., gradients perpendicular to filaments) change their sign periodically within the CO filaments surrounding L 1512, and they suggested an MHD instability may be developing in a helical B-field within the filaments with ~10 µG at a density of 500 cm⁻³. They also found inward motions toward the L 1512 core along a 1 pc long north-south filament with a low accretion rate of ~4 × 10⁻⁶ M☉ yr⁻¹, likely mediated by B-field forces. These results suggest that magnetic fields in the L 1512 core could be strong enough to slow its evolution.

In this paper, we aim to examine the role of magnetic fields in L 1512 by studying the multi-scale B-field morphology and evaluating the B-field strength using linear polarization data. We describe these polarization observations in Section 2, and our reduction recipes in Section 3, and we present our result and analysis in Section 4. In Section 5, we discuss the energy budgets of L 1512. Lastly, we summarize our results in Section 6.

## 2. Observations

### 2.1. James Clerk Maxwell Telescope (JCMT) Observations

The Submillimeter Common-User Bolometer Array 2 (SCUBA-2; Holland et al. 2013) photometric 850 µm observations of L 1512 were retrieved from the JCMT archive.<fnref n="10" /> They are part of the M13BC01 project (PI: James di Francesco) observed in 2013 November. The raw data from six observations and one reduced image were downloaded. Each observation was integrated for 40 minutes with the PONG-900 scan pattern, which fully samples a 15' diameter circular region. The telescope beam size at 850 µm is 14". The downloaded image is presented in Figure 1(a), showing some possibly spurious extended structures. This image was reduced by the pipeline using the iterative map-making routine *makemap* provided by the SMURF package<fnref n="11" /> in the *Starlink* software suite (Chapin et al. 2013). We applied another reduction routine on the raw data to eliminate the spurious structures, and the details are described in Section 3.1.

For 850 µm polarization, we observed L 1512 with the POL-2 polarimeter (Friberg et al. 2016) mounted on the SCUBA-2 camera on JCMT. L 1512 was observed 27 times, for a total of 14 hr, between 2020 August to 2020 November in Band 2 weather (0.05 < τ<sub>225GHz</sub> < 0.08) under project code M20BP046 (PI: Sheng-Jun Lin). Each observation was integrated for 31 minutes with the POL-2-DAISY scan pattern, which fully samples a 3' diameter circular region. The reduction details are described in Section 3.2.

### 2.2. Mimir Polarization Observations

We conducted H-band polarization observations toward L 1512 on 2019 December 16, 17, and 20, and on 2020 January 6 and 9, and February 9 with the Mimir instrument (Clemens et al. 2007) mounted on the 1.8 m Perkins telescope outside Flagstaff, Arizona, U.S., operated by Boston University. The seeing conditions were between 1"/6 and 3"/4 for all observations. The data were calibrated with the Basic Data Processor (MSP-BDP) and Photo-Polarimetry Analysis Tool (MSP-PPOL) from the Mimir Software Package (MSP).<fnref n="12" /> The

data selection criteria were that the H-band magnitude was less than or equal to 13.5 mag, the polarization fraction uncertainty (σ<sub>p</sub>) was less than or equal to 0.9% (or σ<sub>p</sub> ≤ 0.009 in a decimal equivalent), and an S/N of the polarization fraction (p/σ<sub>p</sub>) ≥ 2 or σ<sub>PA</sub> ≤ 15°, yielding a sample of 31 stars out of the 194 total observed.

### 2.3. Other Observations

Herschel 500 µm data with the beam size of 36" (observation ID: 1342191182, quality: level 2 processed) were retrieved from the NASA/IPAC Infrared Science Archive.<fnref n="13" /> Planck 353 GHz (850 µm) polarization data were retrieved from the Planck Legacy Archive.<fnref n="14" /> We employed the Planck 2018 GNLIC Stokes I, Q, and U foreground thermal dust emission maps at variable resolution (5'–80') over the sky (Planck Public Data Release 3 (PR3); Planck Collaboration IV 2020; Planck Collaboration XII 2020), where the Galactic foreground emission is estimated by the GNLIC foreground dust model (Remazeilles et al. 2011), with reduced contamination from the cosmic infrared background, cosmic microwave background, and instrumental noise. The data set we used for L 1512 has an effective beam size of 15'.

## 3. Data Reduction

We describe our JCMT data reduction recipes and compare them to the standard processes in the following.

### 3.1. JCMT SCUBA-2 Photometric Map

Figure 1(a) shows the archival SCUBA-2 image with the possibly spurious extended structures (see Section 2.1). In order to improve the map quality, we reduced the raw SCUBA-2 data using the *skyloop* routine with a configuration file optimized for extended emission,<fnref n="15" /> provided by the SMURF package. With this configuration file, the *skyloop* routine performed principal component analysis (PCA) to remove correlated noise (see Chapin et al. 2013). Figure 1(b) shows this reprocessed map, where the spurious extended structures are now absent. The map was gridded to 4" pixels and calibrated using a flux conversion factor (FCF) of 537 Jy pW⁻¹ (Dempsey et al. 2013).

### 3.2. Missing Flux of JCMT POL-2 for Extended Structures

The POL-2 polarization data were first reduced with the standard two-stage reduction process<fnref n="16" /> using the *pol2map* routine in the SMURF package, which includes a PCA model to remove correlated noise. Pattle et al. (2021) contains a detailed description of the standard POL-2 reduction process. Figure 1(c) presents the POL-2 total intensity map reduced with the above standard procedure, while Figure 1(d) presents the total POL-2 intensity map reduced from the same POL-2 data but with a different reduction procedure that excluded the default PCA model (described later in this section). In contrast to the obvious detection of L 1512 in the SCUBA-2 map (Figure 1(b), peak at 70 mJy beam⁻¹), we found that the dust emission of L 1512 is undetected in Stokes I (Figure 1(c)), and

<footnote n="10">¹⁰ http://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/en/jcmt/</footnote>

<footnote n="11">¹¹ http://starlink.eao.hawaii.edu/docs/sun258.htm/sun258.html</footnote>

<footnote n="12">¹² https://people.bu.edu/clemens/mimir/software.html</footnote>

<footnote n="13">¹³ https://irsa.ipac.caltech.edu/data/Herschel/HHLI/index.html</footnote>

<footnote n="14">¹⁴ http://pla.esac.esa.int/pla/#maps</footnote>

<footnote n="15">¹⁵ https://www.eaobservatory.org/jcmt/2019/04/a-new-dimmconfig-that-uses-pca/</footnote>

<footnote n="16">¹⁶ http://starlink.eao.hawaii.edu/docs/sc22.htm/sc22.html</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 1.** L 1512 JCMT 850 µm total intensity maps sampled on a 4" grid. (a) Archival SCUBA-2 Pong-15' map from the *makemap* routine (see Section 2.1). (b) Reprocessed SCUBA-2 map from the *skyloop* routine with a PCA model (see Section 3.1). (c) POL-2 Stokes *I* map with the standard *pol2map* reduction process, including a PCA model, and (d) the reprocessed POL-2 Stokes *I* map from the *makemap* routine with the PCA model disabled (see Section 3.2). The uniform noise fields with diameters of 15' for SCUBA-2 and 3' for POL-2 are indicated by white circles. Panels other than (c) have the same color scales and same relative contour levels (20% and 60%) with respect to their peak intensities (74, 70, and 69 mJy beam⁻¹ for panels (a), (b), and (d)) for straightforward comparison. The color scale of panel (c) is [-3σ, 3σ], where σ = 3.1 mJy beam⁻¹ is measured in the central 3' field.</figcaption></figure>

neither *Q* nor *U* are detected with the standard POL-2 data reduction. Given the rms noise of 3.1 mJy beam⁻¹ measured in the central 3' field in the POL-2 map using the standard reduction, the dust emission peak of L 1512 should be detected with an S/N of ~20 based on the SCUBA-2 map. Therefore, we suspected that the POL-2 scanning pattern or other instrument/reduction issues might have filtered out extended emission when observing low-surface-brightness objects.

One difference between observations with and without POL-2 is the scan speed. POL-2 contains a half-wave plate (HWP) and a grid analyzer to measure polarized light. The standard POL-2 observing mode (POL-2-DAISY) is a “scan and spin” mode, in which the telescope is continuously moving while the HWP spins (Friberg et al. 2016). The POL-2-DAISY scan speed is slow (8" s⁻¹), which is limited by the integration time for a full polarization cycle (i.e., the HWP rotation speed of

2 Hz). Thus, the data can be sufficiently sampled, and the Stokes parameters *Q* and *U* can be registered to the maps with at least 4" pixels, which is the default pixel size (about one-third of the 850 µm beam size of 14"; i.e., oversampled by a factor of ~1.5 relative to a standard Nyquist sampling rate) adopted by the *pol2map* routine. Such a slow speed could result in less extended structure recovery because the sensitivity of bolometers is limited by the low-frequency noise dominated by slow atmospheric variation between bolometer readouts (Chapin et al. 2013; Friberg et al. 2016).

One approach to evaluate this possibility would be to work with the observatory to perform POL-2 observations with the other scan patterns. An alternate approach is to disable the PCA model included in the standard two-stage reduction process. PCA can remove the correlated components of the bolometer readouts. However, it may be not a good solution for producing
<figure><img src="image_6.png" /><figcaption>**Figure 2.** L 1512 850 µm total intensity profiles. The curves are averaged radial profiles while the filled areas show the ranges of the maximum and minimum intensities. The red curve and filled area show the profile of the SCUBA-2 map in Figure 1(b). The blue dashed curve shows the profile of the POL-2 Stokes *I* map in Figure 1(c). The blue solid curve and filled area show the profile of the POL-2 Stokes *I* map in Figure 1(d).</figcaption></figure>

the extended structures because the real signal could be also
correlated among the bolometers and hence is removed, as
discussed by Chapin et al. (2013). The POL-2 map size is about
half the size of the SCUBA-2 PONG-900 map, so the relative
coverage of extended structures on the POL-2 map is larger
than that of the SCUBA-2 map. As a result, these extended
structures are sampled more with POL-2, producing more
correlated signals among the bolometers for POL-2. This may
explain why the PCA model can help to improve the SCUBA-2
map (Figure 1(b)) but failed for the POL-2 map (Figure 1(c)).
We rereduced the POL-2 raw data using the *makemap* routine
without PCA and set *flagslow*=0.01 to take into account
the slow scan speed used by POL-2 so that the data taken with
the POL-2 scan speed will not be flagged and ignored when
producing the map. This allows comparing the results of the
standard two-stage process with and without PCA, to see what
the impact of the PCA model is. The new reduction result is
shown in Figure 1(d). Without the default PCA model, L 1512
appears in the new, but noisier, map with a similar peak
intensity as seen in our reprocessed SCUBA-2 map
(Figure 1(b)).

Figure 2 shows the azimuthally averaged intensity profiles of our improved SCUBA-2 map (Figure 1(b)) compared to the (no PCA) POL-2 total intensity map (Figure 1(d)). The intensity profiles are similar, showing evidence that we did detect the signal from L 1512 with POL-2, although we cannot rule out the possibility that parts of the emission were still filtered out. The standard POL-2 data reduction procedure incorporates the PCA model for the purpose of removing background signals. This step is particularly essential due to the nonuniform and

radially dependent background variation within the POL-2
observing field. Given that all bolometers in POL-2 observe the
same background but with time delays in the time-stream data,
these signals are correlated and identified as such by PCA. The
PCA approach does effectively address the background signal
issue. However, the presence of extended structures also leads
to correlated total intensity signals between bolometer readouts,
and these are subsequently removed by the PCA model. The
above experiment of not using the PCA model indicates that
the Stokes *I* intensity from L 1512 could be recovered in the
polarization data despite the limitations in map quality.

Therefore, the POL-2 missing flux issue in the total intensity
of L 1512 is likely due to several factors: (1) the scan speed
of 8" s⁻¹ is slow, (2) the background variation is radially
dependent in the POL-2 observations, and (3) L 1512 is faint
and extended, resulting in a low surface brightness contrast.

On the other hand, the HWP rotation speed of 2 Hz provides
a fast 8 Hz modulation of linear polarized intensity, making the
signals of Stokes Q and U extractable (Friberg et al. 2016).
However, POL-2 data reduction requires a total intensity
(Stokes I) map to define the source region with a fixed S/N to
enable reducing the Stokes Q and U signals. We thus used the
reprocessed SCUBA-2 map (Figure 1(b) and Figure 3(a) for a
zoomed-in view) as the input Stokes I model for the POL-2
data reduction instead of the POL-2 Stokes I map. The reduced
Q and U maps were gridded to 4" pixels and calibrated using a
POL-2 FCF of 725 Jy pW⁻¹ (Dempsey et al. 2013; Friberg
et al. 2016). This POL-2 data reduction, involving the use of an
external SCUBA-2 map, was recommended upon the initial
release of POL-2 for open use (Friberg et al. 2016). This
procedure was employed in the first POL-2 850 µm polariza-
tion studies by Pattle et al. (2017), Ward-Thompson et al.
(2017), which focused on the Orion A filament, marking the
initial publication from the B-fields In STar-forming Region
Observations (BISTRO) survey. The newer standard proce-
dure, utilizing the POL-2 Stokes I map, was adopted in
subsequent BISTRO publications shortly after Kwon et al.
(2018), as the targeted objects typically exhibit notable
brightness, thereby facilitating effective reduction based on
the POL-2 Stokes I data.

Figure 3 shows the Stokes *I*, *U*, and *Q* maps, where the Stokes *I* maps are overlaid with *B*-field vectors. Panel (a) shows a zoomed-in view of 4" sampled Figure 1(b), while panels (b), (c), and (d) are binned and sampled on 12" grid. The polarization vectors were computed from the Stokes *I*, *Q*, and *U* maps on the 12" grid to make each vector a nearly independent measurement (close to the resolution of 14"). Assuming dust grains are aligned with the magnetic field, the dust emission polarization at submillimeter wavelengths is perpendicular to the plane-of-sky magnetic field, while the dust-extincted starlight polarization in the NIR is parallel to the plane-of-sky magnetic field (e.g., Andersson et al. 2015; Pattle et al. 2023). Thus, the vectors overlaid on Figures 3(a) and (b) were rotated 90°, so they may represent the plane-of-sky *B*-field direction. The vectors have been then filtered, based on the S/N of the polarization fraction being larger than 3 (*p*/*σ*<sub>*p*</sub>≥3) and the S/N of the total intensity larger than 10 (*I*/*σ*<sub>*I*</sub>≥10). We note that a few vectors appear in the periphery with very high polarization fractions of about 50%. This is probably because the total intensity is underestimated, as the whole point of using a SCUBA-2 PONG observation is to
<figure><img src="image_7.png" /><figcaption>Figure 3. L 1512 JCMT 850 µm Stokes I, Q, and U maps. (a) SCUBA-2 4" sampled map (a zoomed-in view of Figure 1(b)) used as the input total intensity model (see Section 3.2) and (b) the same map but smoothed and sampled on a 12" grid. Both maps are overlaid with contours at 0%, 20%, 40%, 60%, and 80% of their peak intensities. The B-field vectors (polarization vectors rotated by 90°) are overlaid on the 12" grid, where vectors with polarization fraction S/N (p/σp) larger than 6 are shown in red, and vectors with 6 > p/σp ≥ 3 in blue. Vector lengths are scaled according to polarization fraction in panel (a), where a scale bar of 10% (or p = 0.1 in a decimal equivalent) is denoted in the bottom right corner, and shown with uniform lengths in panel (b). (c) Stokes Q and (d) Stokes U 12" sampled maps overlaid with the contours from panel (b). The beam sizes of 14" are denoted in the bottom left corners. The 3' POL-2 uniform noise fields are indicated by white and black circles, in which the rms noises are measured as 4.35, 1.03, 0.79, and 0.79 mJy beam⁻¹ in panels (a)–(d), respectively.</figcaption></figure>

maximize large-scale structure recovery, but the flux could still be missed, as L 1512 is fainter and more extended (i.e., low-contrast surface brightness) compared to the other starless cores, which have had 850 µm polarization detection (see Section 1). The average rms noise of the Stokes Q and U maps over the central 3' field is ~0.79 mJy for 12" pixels. For the total intensity map, the average rms noise is ~1.03 mJy for 12" pixels. We note that the total intensity rms noise is not related to the Stokes Q and U rms noise because we used the total intensity map observed with SCUBA-2 instead of POL-2.

### 3.3. Polarization Properties

The formulae for computing the polarization properties from the reduced data are as follows: the nondebiased polarization fraction (p') is given by

$$p' = \frac{1}{I} \sqrt{Q^2 + U^2}, \qquad (1)$$

where I, Q, U are the Stokes parameters measured. We let σI, σQ, and σU be their uncertainties. The uncertainty of the polarization fraction (σp) is given by

$$\sigma_p = \sqrt{\frac{Q^2\sigma_Q^2 + U^2\sigma_U^2}{I^2(Q^2 + U^2)} + \frac{(Q^2 + U^2)\sigma_I^2}{I^4}}. \qquad (2)$$

The polarization fraction is debiased using the asymptotic estimator (Vaillancourt 2006; Montier et al. 2015),

$$p = \sqrt{p'^2 - \sigma_p^2}, \qquad (3)$$

and is rearranged by adopting the aforementioned σQ = σU = 0.79 mJy beam⁻¹ and ignoring the I⁻⁴ term of the radicand in Equation (2) (Kwon et al. 2018) as

$$p \approx \frac{1}{I} \sqrt{Q^2 + U^2 - \frac{1}{2}(\sigma_Q^2 + \sigma_U^2)}, \qquad (4)$$
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Mimir *H*-band polarization fraction (*p*), expressed in percentage, is plotted against the stellar distance for 30 stars with positive parallax measurements (π) from Gaia DR3 (Gaia Collaboration 2023), out of the total 31 stars. Dots with their uncertainties are shown in red for *p*/σ<sub>p</sub> ≥ 3, and in blue for 3 > *p*/σ<sub>p</sub> ≥ 2. The black dashed line indicates the nominal distance to L 1512 of 140 pc.</figcaption></figure>

where *p* is the debiased polarization fraction. The polarization
angle (PA) is given by

$$
\mathrm{PA} = \frac{1}{2} \arctan \left( \frac{U}{Q} \right), \qquad (5)
$$

and the uncertainty of the PA ($\sigma_{\text{PA}}$) by

$$
\sigma_{\text{PA}} = \frac{\sqrt{Q^2\sigma_U^2 + U^2\sigma_Q^2}}{2(Q^2 + U^2)}. \quad (6)
$$

4. Results and Analysis

4.1. Foreground Star Census

The dust-extincted starlight polarization in the Mimir *H* band from the background stars can trace the plane-of-sky magnetic field toward L 1512. We aim to identify and exclude foreground stars from our *H*-band polarimetric detection. We match our 31 Mimir NIR stellar polarization measurements with the Gaia parallax (π) measurements (Gaia Data Release 3 (DR3); Gaia Collaboration 2023). Figure 4 shows the *H*-band polarization fraction (*p*) plotted against their stellar distance, except for one star with a negative parallax (−0.0480 ± 0.2224 mas). The stellar distance (*d*) is calculated by *d* = 1/π, and the corresponding uncertainty is approximated by *σ*<sub>d</sub> = *σ*<sub>π</sub>/π<sup>2</sup> (Luri et al. 2018). We adopt the nominal distance of 140 pc to L 1512 (Kenyon et al. 1994; Torres et al. 2009; Roccatagliata et al. 2020) and identify one foreground star with a distance of 107 ± 4 pc. This foreground star indeed has the lowest *H*-band polarization fraction among our data (0.330% ± 0.077%, expressed in percentage). The other 30 stars are categorized as background stars due to their distances being greater than 183 ± 1 pc. Accordingly, they trace the magnetic field toward L 1512 because CO isotopologue line observations exhibit only one gas component along the sightlines (Falgarone et al. 1998). Figure 5(a) shows the Mimir *H*-band polarization vectors superimposed on the CFHT *H*-band image. A gray vector corresponds to the foreground star, while the remaining 30 vectors are associated with background stars. The sightlines toward these background stars sample from about 130" away

from the core center (R. A. = 5h04m07s5, and decl. =
32°43′25″ 0, J2000; Lin et al. 2020) and continue out to 410″.
Therefore, their PAs mostly trace the morphology of the plane-
of-sky magnetic field in the diffuse envelope, spanning from a
spatial scale of ~0.56 pc down to ~0.18 pc, which is the outer
layers of the L 1512 cloud surrounding the dense core.

4.2. Magnetic Field Morphology

Figure 5 shows Mimir *H*-band polarization (B-field) vectors overlaid on the continuum maps. These images trace different spatial structures of L 1512. Both the JCMT 850 µm and Herschel 500 µm dust emission trace the central cold dust. Herschel had much greater sensitivity to large-scale emission, while SCUBA-2 mostly traces the cold dust in the L 1512 core. The *H*-band dust scattered light (cloudshine; Foster & Goodman 2006) traces the L 1512 envelope (the outer diffuse region surrounding the L 1512 core; Saajasto et al. 2021).

The Mimir vectors mostly trace the magnetic field morph-
ology in the diffuse envelope, with a spatial scale of ~0.56 pc
down to ~0.18 pc (see Section 4.1). The envelope-scale B-field
pattern seems to generally follow the large-scale magnetic field
traced by the Planck 353 GHz (850 µm) data. These Planck
vectors were rotated 90° to represent the plane-of-sky B-field
directions (θ = PA + 90°) and oversampled for visualization, as
the entire field of view (11′5×12′) of Figure 5(a) fits within
one effective Planck beam size of 15′. The Planck data indicate
the large-scale mean field angle of θ<sub>Planck</sub> = −30° over a
spatial scale of ~0.6 pc at the distance of 140 pc.

In the southern region of L 1512, there are significant
deviations of the NIR vector orientations compared to the
Planck vector orientations. For these NIR vectors, the absolute
deviation of their PAs around the Planck mean field angle (i.e.,
|θ<sub>H</sub> − θ<sub>Planck</sub>|) ranges from a minimum of 6° ± 13° to a
maximum of 80° ± 8°, with the 68th percentile absolute
deviation at 51°. In contrast, the northern NIR vectors show
comparatively less deviation from the Planck mean field angle,
with the 68th percentile absolute deviation at 18°. This
indicates a *B*-field orientation transition from large scale to
envelope-scales. These deviations of the NIR vectors seem to
show the effects of a relative motion between the surrounding
medium and L 1512, and most of the interaction with the
surrounding medium occurs in the south-southeast.

In highly extincted regions, NIR background starlight is
faint, and the corresponding NIR polarization is difficult to
detect. To investigate the magnetic fields in the core,
polarimetry at longer, submillimeter wavelengths is necessary.
In Section 3.2, Figure 3(b) has shown the JCMT POL-2 *B*-field
vectors overlaid on the 850 µm continuum map, zoomed to the
core region. Unlike the envelope-scale *B*-field, the core-scale *B*-
field morphology in L 1512 shows a much more ordered
pattern in a nearly vertical orientation (θ<sub>POL2</sub> ≈ 0°). We note
that the pattern exhibits a mostly smooth change in position
angle from θ<sub>POL2</sub> ≈ 0° in the northwest to ≈−30° in the
southeast, similar to the large-scale mean field angle of −30°
seen in the Planck data. The POL-2 data may reveal a twist or
kink in the plane-of-sky orientations in the core region, blended
with the large-scale magnetic field. Figure 5(c) shows a spatial
comparison between POL-2 vectors and Mimir vectors, where
the POL-2 vectors are identical to those in Figure 3(b) but
undersampled by a factor of 2 for clarity. Their distinct
different spatial coverages indeed show that the Mimir data
<figure><img src="image_9.png" /><figcaption>**Figure 5.** L 1512 Mimir H-band background starlight polarimetry compared to dust maps. Mimir polarizations are shown as red, magenta, or white vectors on each panel with scale bars of 1% polarization fraction (or *p* = 0.01 in a decimal equivalent). One star shows a gray vector as it is a foreground star, identified using Gaia data in Figure 4. (a) CFHT H-band image from Lin et al. (2020) and yellow Planck 353 GHz B-field vectors. The Mimir vectors are displayed in magenta for $3>p/σ_p≥2$, and in red for $p/σ_p≥3$. (b) Herschel 500 µm map with contours at 10, 20, 30, 40, and 50 MJy sr$^{-1}$. (c) JCMT SCUBA-2 850 µm map from Figure 1(b) with contours at 10%, and 50% of the peak intensity at 69.5 mJy beam$^{-1}$. The 2D Gaussian intensity fitting result is denoted by a green ellipse. The JCMT 850 µm B-field vectors with $p/σ_p > 3$ are displayed in blue with uniform lengths. The central cross in each panel indicates the center of L 1512 (Lin et al. 2020). The scale bars of 0.1 pc and submillimeter-wavelength beam sizes are denoted in the top right and bottom left corners, respectively.</figcaption></figure>

trace the envelope-scale *B*-field while the POL-2 data trace the core-scale *B*-field.

To further compare the magnetic field orientations at different scales, Figure 6 shows the distributions of the plane-of-sky *B*-field position angle (θ) estimated from our Mimir and POL-2 data plus the mean field angle inferred from the Planck data. These *B*-field angle probability density functions (pdf) are derived to account for the uncertainties (σ<sub>PA</sub>) in the individual PA measurements via accumulation of the representative Gaussian distributions (Clemens et al. 2020). For POL-2 and Mimir data, the mean *B*-field angle and the 68% highest-density interval (68% HDI; equivalent to one standard deviation of a Gaussian distribution) calculated from the pdf are −7°+18°/−12° and −15°+40°/−39°, respectively. In addition, the mean *B*-field angle with standard deviation of −23°±17° inferred from the AIMPOL *R*-band polarization of 94 stars (Sharma et al. 2022) is also overplotted. The optical *R*-band polarization traces the magnetic field beyond the region probed by our NIR *H*-band polarization because the dust extincts more in the *R* band.

Both the AIMPOL data and Planck data trace the large-scale plane-of-sky magnetic field of L 1512, but with different resolutions. The Planck measurements, with their low spatial resolution of 15' and small optical depth at 850 µm (Planck Collaboration XII 2020), not only average the polarization on a large scale toward the L 1512 cloud but also integrate it along the line of sight, weighting it by dust emission intensity. In contrast, the *R*-band measurements could not be made toward the central ~10' extincted area of the L 1512 cloud (see Figure 4 of Sharma et al. 2022). However, the AIMPOL data coverage has a diameter of ~20' (~0.8 pc at the distance of 140 pc). As a result, the magnetic field traced in *R* band spans a wide field, sampling primarily the outer low-density region, but with relatively higher resolution, depending on the *R*-band background star distribution. Given that foreground and background stars can be identified, one can confidently associate the measured polarization with the cloud. Consequently, both the Planck and AIMPOL measurements effectively trace the large-scale magnetic field pattern and are consistent with each other toward the low-density periphery of the L 1512 cloud.

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Probability density function (pdf) of *B*-field position angles (θ). The upper panel shows the Mimir data (red), while the lower panel shows the JCMT POL-2 data (blue). The mean and the 68% highest-density interval (HDI) of the *B*-field angle distributions are displayed with vertical lines and filled areas. The Planck mean field angle of −30° is plotted as a purple line in panel (a). The mean field angle, with standard deviation, of AIMPOL *R*-band polarizations Sharma et al. (2022) is represented as a Gaussian distribution on panel (a), vertically offset by a probability of 1% for clarity. The position angle of the core major axis (PA = −14°;2; see Figure 5(c)) resulting from the 2D Gaussian fitting to the 850 µm intensity is plotted as a black dashed line in panel (b).</figcaption></figure>
The magnetic field toward L 1512 shows a consistent average orientation from the large scale down to the core scale. The largest angular dispersion is seen in the H-band data, related to the previously mentioned deviation of the NIR vectors in the southern region. In addition, Sharma et al. (2022) found that the R-band polarization fraction is lower toward the southern region compared to the other surrounding regions. They suggested that it could be due to the lack of aligned grains owing to the different grain size distribution, collisional disalignment with gas, or the depolarization caused by a tangled B-field. Therefore, with these multiwavelength polari- metric data, magnetic fields are suggested to thread from the large scale to the dense core scale in L 1512, while the magnetic field could interact with the diffuse medium in the envelope-scale.

The dust emission of L 1512 (Figures 5(b) and (c)) shows an elongated, cometary morphology. We performed a 2D Gaussian intensity fitting on the 850 µm map. The fitting result is denoted by the green ellipse in Figure 5(c). The PA of the major axis is −14°2 ± 0°4, and the major and minor FWHM axes are 141′′2 ± 2′′1 and 86′′2 ± 1′′3 (~0.10 and ~0.06 pc at the distance of 140 pc), respectively. The major axis of the projected core shape is parallel to the B-field orientation (see Figure 6(b)) instead of being perpendicular to the magnetic field as suggested by the B-field-dominated core formation scenario at the core scale of ~0.1 pc (e.g., Galli & Shu 1993a; Li & Shu 1996; Ciolek & Basu 2000; Myers et al. 2018). This discrepancy can be attributed to the projection effect of a tri- axial core, as suggested by Basu (2000). Chen & Ostriker (2018) conducted a comprehensive analysis of turbulent MHD simulations, and they discovered that dense cores do tend to be tri-axial, unlike the idealized oblate cores assumed in classical theories. Moreover, they observed that environmental factors, such as ram pressure or magnetic pressure, also play a crucial role in shaping dense cores.

### 4.3. Davis-Chandrasekhar-Fermi Analysis

The DCF (Davis 1951; Chandrasekhar & Fermi 1953) method is widely used to derive the plane-of-sky B-field strength ($B_{\text{pos}}$) using linear polarization data. This method assumes that the turbulence is sub-Alfvénic, and the magnetic field is frozen into the gas, so the nonthermal gas motions result in a distortion of the magnetic field. By measuring the dispersions in the nonthermal gas velocities and in the polarization position angles, the field strength is as follows (Crutcher et al. 2004):

$$B_{\text{pos}} = \text{Q}\sqrt{4\pi\langle\rho\rangle} \frac{\sigma_{v,\text{NT}}}{\delta\theta} \quad (\text{cgs units}) \quad (7)$$

$$\approx 9.3 \sqrt{\langle n_{\text{H}_2} \rangle / \text{cm}^{-3}} \frac{\Delta v_{\text{NT}} / \text{km} \cdot \text{s}^{-1}}{\delta\theta / \text{degree}} \mu\text{G}, \quad (8)$$

where $\langle\rho\rangle$ is the mean volume mass density, $\langle n_{\text{H}_2} \rangle$ is the mean H$_2$ number density ($\langle\rho\rangle = \mu_{\text{H}_2}\langle n_{\text{H}_2} \rangle m_{\text{H}}$, and $\mu_{\text{H}_2}=2.8$), $\sigma_{v,\text{NT}}$ is the 1D nonthermal velocity dispersion, $\Delta v_{\text{NT}}$ is the nonthermal FWHM linewidth ($\Delta v_{\text{NT}} = \sigma_{v,\text{NT}}\sqrt{8\log 2}$), $\delta\theta$ is the intrinsic dispersion in B-field angles, and Q is a correction factor. Based on calibrations from numerical simulations (Heitsch et al. 2001; Ostriker et al. 2001), Q = 0.5 provides a good estimate of $B_{\text{pos}}$, if $\delta\theta < 25°$.

### 4.3.1. Magnetic Field Angular Dispersion

We use the DCF method to estimate the plane-of-sky magnetic field strength at the core scale with the POL-2 data. The DCF method requires the intrinsic B-field angular dispersion measured from the random component of the magnetic field perturbed by local turbulence. However, the PA measurements include not just the random B-field component but also the underlying ordered B-field component and the PA measurement uncertainties. Particularly, the core-scale B-field shows a slightly curved pattern by changing the B-field orientation from $\theta_{\text{POL2}} \approx 0°$ to $\theta_{\text{POL2}} \approx -30°$ in the southwestern side (see Figure 3(b)). Thus, we adopt the “unsharp-masking method” (Pattle et al. 2017) to estimate the intrinsic angular dispersion of the random B-field component in L 1512. In this method, a smoothed B-field position angle map is created by convolving the original B-field angle map with a 3-by-3-pixel boxcar filter (i.e., a width of ~2.6 independent beams). This smoothed map represents the underlying ordered B-field component. The 3-by-3-pixel boxcar filter minimizes the influence of the curvature of the field pattern. After subtracting this smoothed map from the original map, the residual angles represent the random B-field component and the PA measurement uncertainties. The intrinsic B-field angular dispersion can be estimated by the standard deviation of the residual angles.

Our analysis results for 850 µm are shown in Figures 7 and 8. Figure 7(a) shows the B-field position angle ($\theta$) 12"-pixel map, where the vectors are shown with uniform lengths for clarity. Figure 7(d) shows the corresponding B-field-angle uncertainty ($\sigma_{\text{PA}}$) map. We make a map of the underlying B-field by smoothing Figure 7(a) with a 3-by-3-pixel boxcar filter (i.e., averaging over the nearby nine pixels). This smoothed B-field angle ($\langle\theta\rangle_{\text{local}}$) map is shown in Figure 7(b). Figure 7(c) shows the residual ($\Delta\theta = \theta - \langle\theta\rangle_{\text{local}}$) map, in which three isolated pixels from the $\theta$ and $\langle\theta\rangle_{\text{local}}$ maps are rejected, and 104 pixels remained.

On the residual angle map ($\Delta\theta$), the underlying ordered B-field geometry has been removed. In order to estimate a representative intrinsic B-field angular dispersion ($\delta\theta$) across L 1512, we present two approaches. One is the conventional approach, and the other one is a subsequent approach introduced in the unsharp-masking method. The conventional approach for estimating $\delta\theta$ is via the inverse-variance-weighted quadratic mean of the residual angles,

$$\delta\theta = \sqrt{\frac{\sum_i w_i \Delta\theta_i^2}{\sum_i w_i}}, \quad (9)$$

where $w_i = 1/\sigma_{\text{PA},i}^2$ (e.g., Clemens et al. 2016; Wang et al. 2019). Hence, the well-characterized residual angles with smaller $\sigma_{\text{PA}}$ dominate the estimation. We find $\delta\theta = 8.8 \pm 0.5$ using this approach, where the uncertainty is calculated with the standard error propagation. On the other hand, Pattle et al. (2017) ran Monte Carlo simulations and demonstrated that a representative value of $\delta\theta$ can be better recovered, closer to the input $\delta\theta$ value in their simulations (their so-called “true value”), via averaging a set of standard deviation estimates of the residual angles selected with different maximum allowed PA uncertainties ($\sigma_{\text{PA,max}}$); i.e., $\langle\sigma(\Delta\theta|\text{uncertainties} \le \sigma_{\text{PA,max}})\rangle$. Figure 8 shows these standard deviation estimates as the
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Position angles of the B-field orientations and their uncertainties. (a) B-field position angle map ($\theta$) inferred from the polarization angles (PA) rotated by 90°. (b) Mean B-field angle map ($\langle\theta\rangle_{\text{local}}$) from convolving panel (a) by a boxcar filter with a size of 3-by-3 pixels. The black vectors overlaid on panels (a) and (b) show the B-field angles, while the magenta vectors on panel (b) show the mean B-field angles. (c) Residual map ($\Delta\theta = \theta - \langle\theta\rangle_{\text{local}}$) and blue vectors of these residual angles are overlaid. (d) Uncertainty map ($\sigma_{\text{PA}}$) of B-field angles. (e) Map of the local largest PA uncertainty ($\sigma_{\text{PA,local}}$), on which each pixel value represents the largest PA uncertainty ($\sigma_{\text{PA}}$) among the neighboring nine pixels within a 3-by-3 pixel box centered on that pixel. The “unsharp-masking method” sets the selection criteria based on $\sigma_{\text{PA,local}}$ (see Section 4.3.1 and also Pattle et al. 2017). For panels (d) and (e), the minimum, mean, and maximum values are 3°5, 5°9, and 9°3 for $\sigma_{\text{PA}}$, and 4°7, 7°6, and 9°3 for $\sigma_{\text{PA,local}}$, respectively. The beam sizes are denoted in the bottom left corners. The 3' POL-2 uniform noise field is indicated by cyan circles.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Residual angle ($\Delta\theta$) cumulative function of the maximum allowed B-field angle uncertainty ($\sigma_{\text{PA,max}}$), in that the standard deviation of $\Delta\theta$ with the uncertainty ($\sigma_{\text{PA,local}}$) equal to or less than $\sigma_{\text{PA,max}}$ is evaluated. The standard deviation estimate, $\sigma(\Delta\theta)$, is shown as black dots, while the cumulative number of $\Delta\theta$ is shown as red dots. The intrinsic B-field angular dispersion ($\delta\theta$) of $8.2^\circ \pm 1.4^\circ$ (horizontal solid and dashed blue lines) is determined by computing the mean and standard deviation of the $\sigma(\Delta\theta)$ evaluated from $6^\circ \le \sigma_{\text{PA,max}} \le 9.3^\circ$ (vertical dashed black lines).</figcaption></figure>

function of $\sigma_{\text{PA,max}}$, with black dots, and the cumulative number of the residual angles limited by $\sigma_{\text{PA,max}}$, with red dots. Figure 7(e) shows the angular uncertainty of $\Delta\theta$, denoted by $\sigma_{\text{PA,local}}$, adopted by Pattle et al. (2017) for this step. At each pixel, the $\sigma_{\text{PA,local}}$ is the largest $\sigma_{\text{PA}}$ of the neighboring nine pixels previously selected to compute $\langle\theta\rangle_{\text{local}}$. Pattle et al. (2017) ran a set of Monte Carlo simulations by setting different “true values” of the intrinsic B-field angular dispersion ($\delta\theta$) to generate their simulated data sets. They found the above procedure can approach the input $\delta\theta$ when $\sigma_{\text{PA,max}}$ is small, and the cumulative number of $\Delta\theta$ is still sufficient. They took the mean of such well-characterized standard deviation estimates as the representative estimation of $\delta\theta$ for the whole region. In our analysis, we collect the standard deviation estimates computed with $6^\circ \le \sigma_{\text{PA,max}} \le 9.3^\circ$ (two vertical dashed black lines in Figure 8) because these $\sigma(\Delta\theta)$ estimates remain relatively similar values in contrast to the sharp $\sigma(\Delta\theta)$ decline at $\sigma_{\text{PA,max}} < 6^\circ$, which is due to the sample of the residual angles being too small (fewer than 6 vectors for this case). By taking the mean and standard deviation of these $\sigma(\Delta\theta)$ estimates, we obtain an intrinsic B-field angular dispersion ($\delta\theta$) of $8.2^\circ \pm 1.4^\circ$ (horizontal blue lines in Figure 8) for the
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Density maps for the L 1512 core with $R_{\text{edge}} = 108''$. (a) $N_{\text{H}_2}$ map. (b) Line-of-sight-averaged $n_{\text{H}_2}$ map, on which each pixel is calculated as the column density divided by sightline depth. The 3' POL-2 uniform noise field is indicated by white circles for reference.</figcaption></figure>

whole region. Compared to the conventional approach, both $\delta\theta$ estimations are consistent, and we adopt $\delta\theta = 8.2 \pm 1.4$ for the following analysis.

### 4.3.2. Number Density and Velocity Dispersion

Lin et al. (2020) built an onion-like volume density model of the L 1512 core, which is comprised of different eastern and western hemispheres. They modeled the density structure with Plummer-like profiles of $n_0/(1 + (r/R_0)^\eta)$, where $n_0$ is the central density, $R_0$ is the characteristic radius, and $\eta$ is the power-law index of the density profile at $r \gg R_0$. The best-fit Plummer parameters fitted to $N_2H^+$ spectral line observations and dust extinction measurements were found to be $n_0 = 1.1 \times 10^5 \text{ cm}^{-3}$, $R_0 = 0.022 \text{ pc}$, and $\eta = 2.0$ for the east side, and $n_0 = 1.1 \times 10^5 \text{ cm}^{-3}$, $R_0 = 0.027 \text{ pc}$, and $\eta = 3.1$ for the west side. Since the density required by the DCF method should be the average density along the line of sight, we assume the polarized light is mainly contributed by the L 1512 core, which has an $R_{\text{edge}}$ of 108'' (= 0.073 pc; Lin et al. 2020). With the above density profiles, we (1) integrate the volume density profiles within the sphere defined by $R_{\text{edge}}$ to obtain a column density ($N_{\text{H}_2}$) map, and then (2) divide the column densities by the corresponding sightline depth inside the sphere to obtain the line-of-sight-averaged $H_2$ number density ($n_{\text{H}_2,\text{los}}$) map (i.e., for a pixel indexed $j$, $n_{\text{H}_2,\text{los},j} = N_{\text{H}_2,j}/(2\sqrt{R_{\text{edge},\text{cm}}^2 - d_{j,\text{cm}}^2})$, where $R_{\text{edge},\text{cm}}$ and $d_{j,\text{cm}}$ are in the units of centimeter, and $d_{j,\text{cm}}$ is the projected distance to the core center). Figure 9 shows the $N_{\text{H}_2}$ map and $n_{\text{H}_2,\text{los}}$ map. For self-consistency, we compute the mean, and standard deviation,

<table><caption>**Table 1**

Estimated Properties in the DCF Analysis from the 850 µm Polarimetry</caption><thead><tr><th>Property</th><th>Value</th></tr></thead><tbody><tr><td>δθ (degree)</td><td>8°2° ± 1°4</td></tr><tr><td>Δv<sub>NT</sub> (km s<sup>-1</sup>)</td><td>0.109 ± 0.031</td></tr><tr><td>⟨n<sub>H2</sub>⟩ (cm<sup>-3</sup>)</td><td>(2.1 ± 0.9) × 10<sup>4</sup></td></tr><tr><td>⟨N<sub>H2</sub>⟩ (cm<sup>-2</sup>)</td><td>(8.3 ± 4.7) × 10<sup>21</sup></td></tr><tr><td>B<sub>pos</sub> (μG)</td><td>18 ± 7</td></tr><tr><td>λ<sub>obs</sub></td><td>3.5 ± 2.4</td></tr><tr><td>|B<sub>cor</sub>| (μG)</td><td>23 ± 9</td></tr><tr><td>λ<sub>cor</sub></td><td>1.2 ± 0.8</td></tr><tr><td>B<sub>tot</sub> (μG)</td><td>~32</td></tr><tr><td>B<sub>los</sub> (μG)</td><td>~27</td></tr><tr><td>λ<sub>tot</sub></td><td>~1.6</td></tr></tbody></table>

**Note.** Quantities shown are the dispersion in magnetic field position angles ($\delta\theta$), the nonthermal FWHM linewidth ($\Delta v_{\text{NT}}$), the mean $H_2$ number density ($\langle n_{\text{H}_2} \rangle$), the mean $H_2$ column density ($\langle n_{\text{H}_2} \rangle$), the plane-of-sky B-field strength ($B_{\text{pos}}$), the observed mass-to-flux ratio ($\lambda_{\text{obs}}$), the statistically corrected total B-field strength ($|B_{\text{cor}}|$), and the statistically corrected mass-to-flux ratio ($\lambda_{\text{cor}}$); and the total/line-of-sight B-field strength ($B_{\text{tot}}$, $B_{\text{los}}$) and the corresponding mass-to-flux ratio ($\lambda_{\text{tot}}$) are derived by assuming the L 1512 core is virially stable (see Section 5.1). The uncertainties are their measured dispersion or computed by the standard error propagation. For $\Delta v_{\text{NT}}$, we adopt the $N_2H^+$ (1–0) spectral resolution from Lin et al. (2020) as the uncertainty. For $B_{\text{tot}}$, $B_{\text{los}}$, and $\lambda_{\text{tot}}$, we do not estimate uncertainties because their derivation involves energy budgets (Table 2) that may be accurate only to order of magnitude.

of $n_{\text{H}_2,\text{los}}$ and $N_{\text{H}_2}$ values on the same set of pixels, of which the residual angles (Figure 7(c)) are used to estimate the representative B-field angular dispersion $\delta\theta$. We find $\langle n_{\text{H}_2} \rangle = \langle n_{\text{H}_2,\text{los}} \rangle = (2.1 \pm 0.9) \times 10^4 \text{ cm}^{-3}$, and $\langle n_{\text{H}_2} \rangle = (8.3 \pm 4.7) \times 10^{21} \text{ cm}^{-2}$. Lin et al. (2020) found that a nonthermal FWHM linewidth of $\Delta v_{\text{NT}} = 0.109 \text{ km s}^{-1}$ can reproduce their $N_2H^+$ (1–0) spectral observations. We take their $N_2H^+$ (1–0) spectral resolution of $0.031 \text{ km s}^{-1}$ as the uncertainty. We adopt these values for the DCF analysis and list them in Table 1.

### 4.3.3. Magnetic Field Strength and Mass-to-flux Ratio

Using the DCF method (Equation (8)), and the above-estimated values (see Table 1), we calculate the mean plane-of-sky B-field strength ($B_{\text{pos}}$) across the core to be $18 \pm 7 \mu\text{G}$, where the uncertainty is computed with the standard error propagation from the dispersions of $\delta\theta$, $\Delta v_{\text{NT}}$, and $n_{\text{H}_2}$, and hence, the uncertainty of $B_{\text{pos}}$ represents the dispersion of its distribution.

It is also important to understand if the magnetic field could support the core against gravity, which can be determined by the mass-to-flux ratio in units of the critical ratio, $(M/\Phi)_{\text{crit}}$. We use the formula from Crutcher et al. (2004),

$$ \lambda = \frac{(M/\Phi)_{\text{obs}}}{(M/\Phi)_{\text{crit}}} = 7.6 \times 10^{-21} \frac{N_{\text{H}_2}/\text{cm}^{-2}}{B/\mu\text{G}}. \quad (10) $$

The core is magnetically supercritical if $\lambda > 1$ (i.e., unstable to collapse) and magnetically subcritical if $\lambda < 1$ (i.e., magnetically supported). With our derived $\langle n_{\text{H}_2} \rangle$ and $B_{\text{pos}}$, we calculate the observed mass-to-flux ratio ($\lambda_{\text{obs}}$) as $3.5 \pm 2.4$ or a range of 1.1–5.9, suggesting that the magnetic field alone may not fully support the core against gravity. However, the mass-to-flux ratio could be overestimated by $\lambda_{\text{obs}}$, due to the unknown
inclination of the B-field and the 3D geometry of the core (Crutcher et al. 2004). On the other hand, the $\lambda_{obs}$ lower limit of 1.1 does not entirely rule out the possibility of the magnetically critical condition, owing to various sources of uncertainty in the DCF method (Crutcher 2012; Pattle & Fissel 2019). Therefore, our estimation of the mass-to-flux ratio suggests that the L 1512 core is approximately magnetically critical or supercritical. We note that the mass-to-flux ratio calculated here is considered as a qualitative indicator rather than as a precise measure of the L 1512 core's stability against gravity. We also note that, when computing the mass-to-flux ratio, contributions from thermal and nonthermal energies are not considered.

Crutcher et al. (2004) showed that $B_{pos}$ and $\lambda_{obs}$ can be statistically corrected by averaging the cases of all B-field inclinations with respect to the line of sight. Under this methodology, $B_{pos}$ is a statistical average of the total magnetic field, $B_{cor}$, which can point over a total solid angle of $2\pi$ subtended from the core center. The correction is

$$B_{pos} = \frac{\pi}{4} |B_{cor}|. \qquad (11)$$

In this case, the total magnetic field would have a strength of $23 \pm 9$ µG. For our observed mass-to-flux ratio, $\lambda_{obs}$, Crutcher et al. (2004) derived the correction as

$$\lambda_{cor} = \frac{\lambda_{obs}}{3}. \qquad (12)$$

This correction also takes into account the projection effect of $N_{H_2}$, assuming that the magnetic flux tube is aligned with the core minor axis. This yields a mass-to-flux ratio of $\lambda_{cor} = 1.2 \pm 0.8$, still suggesting the core is approximately magnetically critical or slightly supercritical. We will further discuss the B-field strength in the context of the virial analysis in Section 5.

### 4.4. Grain Alignment

A long-standing debate about submillimeter polarimetry is whether dust grains remain aligned with magnetic fields inside starless cores (e.g., Alves et al. 2014; Andersson et al. 2015; Jones et al. 2015; Pattle et al. 2019; Wang et al. 2019). From a theoretical standpoint, the main alignment mechanism is thought to be radiative alignment torques (B-RATs), in which an anisotropic radiation field is a key to aligning the dust grains with the magnetic field (e.g., Dolginov & Mitrofanov 1976; Lazarian & Hoang 2007; Andersson et al. 2015). In the center of starless cores, the absence of an internal radiation source producing an anisotropic radiation field might lead to an outcome where dust grains are not aligned with the magnetic field. Therefore, it is crucial to investigate whether the dust grains are aligned at the centers of starless cores.

To evaluate grain alignment in the submillimeter regime, a power-law index in the $p-I$ relation,<fnref n="17" />

$$p_{\text{true}}(I) \propto I^{-\alpha}, \qquad (13)$$

is sought, where $0 \le \alpha \le 1$ (e.g., Alves et al. 2015; Jones et al. 2015). In the $\alpha = 0$ case, the polarization fraction is constant at every depth along each sightline, indicating a perfectly aligned case. For the $\alpha = 1$ case, the polarized intensity ($=p_{\text{true}} \cdot I$) is constant. In such a case, one interpretation is that only the outer shell of the starless core contributes to the polarized intensity, and no contribution comes from the inner region, implying a completely disaligned case (Pattle et al. 2019). We can only probe the B-field morphology at the core surface in this case.

In order to determine the index $\alpha$, we should take into account that the polarization fraction measurements do not follow a Gaussian distribution. Typically, the polarization fraction measurements are statistically debiased with a polarization estimator (e.g., the asymptotic estimator defined by Equation (3); please refer to Montier et al. 2015, for the detailed discussion of other estimators) and filtered with an S/N criterion for fitting $\alpha$. In fact, these polarization estimators are *biased* estimators, but they work better in the high-S/N regime because the *bias* is minimized (Montier et al. 2015); therefore, the low-S/N data should be removed in order to better determine $\alpha$. For faint starless cores, such S/N criteria may remove too much data, and thus, $\alpha$ might not be well-constrained.

Instead of debiasing data and removing the low-S/N ones, Pattle et al. (2019) assumed that the nondebiased polarization fraction measurements follow a Rice distribution, including both low- and high-S/N data, and took the mean of the Rice distribution (Rice 1945; Serkowski 1958) as the estimator of the true polarization fraction at a given $I$. They referred to this estimator as the Ricean-mean model,

$$p'(I) = \sqrt{\frac{\pi}{2}} \left(\frac{I}{\sigma_{QU}}\right)^{-1} \mathcal{L}_{\frac{1}{2}} \left(-\frac{p_{\sigma_{QU}}^2}{2} \left(\frac{I}{\sigma_{QU}}\right)^{2(1-\alpha)}\right), \qquad (14)$$

where $\mathcal{L}_{\frac{1}{2}}$ is a Laguerre polynomial of order $1/2$, $\sigma_{QU}$ is the average of $\sigma_Q$ and $\sigma_U$, $\alpha$ is the power-law index, and $p_{\sigma_{QU}}$ is the expected polarization fraction when $I = \sigma_{QU}$. We note that, technically, $p_{\sigma_{QU}}$ and $\sigma_{QU}$ are the free scaling factors of the assumed underlying $p-I$ relation of $p_{\text{true}} = p_{\sigma_{QU}}(I/\sigma_{QU})^{-\alpha}$, where the $I$ scaling factor is chosen to be the practical noise ($\sigma_{QU}$) while the polarization fraction scaling factor ($p_{\sigma_{QU}}$) is left for fitting. This Ricean estimator allows for a better estimation of the index $\alpha$ by being able to properly include the low-S/N data.

We use our data to determine the index $\alpha$ in L 1512 with the above Ricean-mean model. We collect data within a 3' uniform noise field ($\sigma_{QU} = 0.79$ mJy beam$^{-1}$), and supply $\sigma_p$ as a data weighting in order to perform a fitting with the Python `scipy.optimize.curve_fit` function. We obtain a best-fit $\alpha$ of $0.48 \pm 0.05$, and $p_{\sigma_{QU}}$ of $0.60 \pm 0.10$. Figure 10 shows the best-fit Ricean-mean model in the solid red curve. The true value of $\alpha = 0.48 \pm 0.05$ is recoverable above the critical intensity ($I_{\text{crit}}$) of 3.22 mJy beam$^{-1}$ (corresponding to a S/N of $I_{\text{crit}}/\sigma_{QU} = 4.05$). In contrast, owing to the low S/N, the true $\alpha$ is not recoverable below $I_{\text{crit}}$ but apparently approaches the value of unity expected for disaligned grains,

$$p'(I) = \sqrt{\frac{\pi}{2}} \left(\frac{I}{\sigma_{QU}}\right)^{-1}, \qquad (15)$$

<footnote n="17">¹⁷ We note that the polarization fraction discussed in this power law is the true/intrinsic value ($p_{\text{true}}$), while both of the nondebiased and debiased polarization fractions ($p'$ and $p$ defined by Equations (1) and (3)) are measurements. The measurement $p$ is an attempt made by the asymptotic estimator to find $p_{\text{true}}$.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 10.** JCMT 850 µm nondebiased polarization fraction ($p'$) is plotted against the total intensity (Stokes I). The best-fit Ricean-mean model (Equation (14)) and the grain disalignment model (Equation (15)) are shown as solid red and dashed black lines, respectively. The red- and green-colored regions are the 68% and 95% confidence intervals, respectively.</figcaption></figure>

shown as the dashed black line. The index of $\alpha = 0.48 \pm 0.05$ suggests that some of the dust grains remain aligned with the magnetic field at higher densities (Andersson et al. 2015; Pattle & Fissel 2019). Therefore, our data are able to probe the magnetic field in the L 1512 core.

# 5. Discussion

## 5.1. Gravitational Stability

The gravitational stability of the L 1512 core can be assessed using the virial theorem (McKee & Zweibel 1992; McKee & Ostriker 2007), which can be written as

$$ \frac{1}{2} \frac{d^2 I}{dt^2} = 2(E_{\text{thermal}} - E_{\text{surface}} + E_{\text{turb}} + E_{\text{rot}}) + E_{\text{grav}} + E_{\text{mag}}, \quad (16) $$

where I is the moment of inertia, $E_{\text{thermal}}$ is the thermal energy in the core, $E_{\text{surface}}$ is the surface kinetic term due to external pressure, $E_{\text{turb}}$ is the turbulent energy, $E_{\text{rot}}$ is the rotational energy, $E_{\text{grav}}$ is the gravitational potential energy, and $E_{\text{mag}}$ is the magnetic energy calculated by the total B-field strength ($B_{\text{tot}}$). The magnetic energy can be divided into the plane-of-sky and line-of-sight components such that $E_{\text{mag}} = E_{\text{mag.pos}} + E_{\text{mag.los}}$, since $B_{\text{tot}}^2 = B_{\text{pos}}^2 + B_{\text{los}}^2$. The total kinetic energy can be defined as $T = E_{\text{thermal}} - E_{\text{surface}} + E_{\text{turb}} + E_{\text{rot}}$. By defining the virial energy of $E_{\text{vir}} = 2T + E_{\text{grav}}$, we rewrite

<table><caption>**Table 2**
Energy Budgets in L 1512</caption><thead><tr><th>Energy Compared with E<sub>grav</sub></th><th>Value</th></tr></thead><tbody><tr><td>E<sub>grav</sub> (erg)</td><td>-3.9 × 10<sup>42</sup></td></tr><tr><td>E<sub>thermal</sub>/|E<sub>grav</sub>|</td><td>0.61</td></tr><tr><td>E<sub>surface</sub>/|E<sub>grav</sub>|</td><td>0.41</td></tr><tr><td>E<sub>turb</sub>/|E<sub>grav</sub>|</td><td>0.04</td></tr><tr><td>E<sub>rot</sub>/|E<sub>grav</sub>|</td><td>0.01</td></tr><tr><td>T/|E<sub>grav</sub>|</td><td>0.25</td></tr><tr><td>E<sub>vir</sub>/|E<sub>grav</sub>|</td><td>-0.51</td></tr><tr><td>E<sub>mag.pos</sub>/|E<sub>grav</sub>|</td><td>0.16</td></tr><tr><td>(E<sub>vir</sub> + E<sub>mag.pos</sub>)/|E<sub>grav</sub>|</td><td>-0.35</td></tr></tbody></table>

**Note.** Quantities shown are the gravitational energy ($E_{\text{grav}}$), the thermal energy ($E_{\text{thermal}}$), the surface kinetic term due to external pressure ($E_{\text{surface}}$), the turbulent energy ($E_{\text{turb}}$), the rotational energy ($E_{\text{rot}}$), the total kinetic energy ($T = E_{\text{thermal}} - E_{\text{surface}} + E_{\text{turb}} + E_{\text{rot}}$), the virial energy ($E_{\text{vir}} = 2T + E_{\text{grav}}$), and the plane-of-sky component magnetic energy ($E_{\text{mag.pos}}$).

Equation (16) as

$$ \frac{1}{2} \frac{d^2 I}{dt^2} = E_{\text{vir}} + E_{\text{mag.pos}} + E_{\text{mag.los}}. \quad (17) $$

The sign of $d^2I/dt^2$ determines whether the core is virially unbound with a positive net energy ($d^2I/dt^2 > 0$) or virially bound with a negative net energy ($d^2I/dt^2 < 0$), and $d^2I/dt^2 = 0$ means that the core is virially stable.

To assess the core stability, the physical structure of L 1512 is needed. Lin et al. (2020) focused on $N_2H^+$ multitransition spectra modeling and visual extinction measurement to build an onion model to describe the volume density, kinetic temperature, and rotational velocity profiles for the case of a constant turbulent velocity in the L 1512 core. Their onion model can well reproduce the $N_2H^+$ and the other line data observed along a horizontal (RA) cut and a vertical (Dec) cut across the L 1512 core. The onion model is comprised of eastern (nine layers) and western (six layers) hemispheres, because the west extent of L 1512 is a factor of ~2/3 shorter than along the eastern side. Here, we adopt the eastern hemisphere model to represent the entire core because the eastern side of L 1512 is more spherically symmetric with respect to the core center than the western side. Thus, the eastern hemisphere onion model could provide a better description for the core (also see Figures 1 and 2 from Lin et al. 2020).

We use the physical parameters in the eastern onion model from Lin et al. (2020) to calculate each energy term in the virial equation (Equation (16)). Please refer to Table C.1 and Figure 5(c) from Lin et al. (2020) for the density, temperature, and rotational velocity profiles. In addition, the constant 1D nonthermal velocity dispersion, $\sigma_{v,NT}$, used in the onion model is 0.046 km s$^{-1}$. We calculate the plane-of-sky component magnetic energy ($E_{\text{mag.pos}}$) by using the DCF-derived plane-of-sky magnetic field strength ($B_{\text{pos}}$) of 18 µG (Table 1). Table 2 summarizes the results for each energy of the L 1512 core, and Appendix shows the formulae we used for deriving each energy term. If we use the physical parameters in the western onion model from Lin et al. (2020), the virial energy ($E_{\text{vir}}$) will change from -0.51 to -0.39, and the plane-of-sky component magnetic energy ($E_{\text{mag.pos}}$) will change from 0.16 to 0.11, in units of $|E_{\text{grav}}|$. Accordingly, if $E_{\text{mag.los}} = 0$, the value of the
inertia term on the left-hand side of Equation (17) will change from −0.35 to −0.28 in units of $|E_{\text{grav}}|$. We note that these energies could also have significant uncertainties; thus, these derived values are accurate only to order of magnitude. However, as the dominant uncertainty is mass, all the derived energies are linearly dependent on the mass, except that $E_{\text{grav}}$ is dependent on the square of mass.

The rotational energy ($E_{\text{rot}}$) in Table 2 is one of the interpretations from the $N_2H^+$ (1–0) observation toward L 1512, which reveals a uniform velocity gradient ($2.26 \pm 0.04 \text{ km s}^{-1} \text{pc}^{-1}$) along roughly the north–south direction across an extent of ~0.1 pc (Caselli et al. 2002; Lin et al. 2020). An interpretation is a slow solid-body rotation (Caselli et al. 2002; Figure 5(c) in Lin et al. 2020), of which the corresponding rotation period of $2.72 \pm 0.05 \text{ Myr}$ ($\Omega = (7.3 \pm 0.1) \times 10^{-14} \text{ rad s}^{-1}$ and $v_{\text{rot}} \le 0.07 \text{ km s}^{-1}$) is about twice as long as the core lifetime of $\gtrsim 1.4 \text{ Myr}$ (Lin et al. 2020). The corresponding $E_{\text{rot}}$ is only 1% of $|E_{\text{grav}}|$. This contribution does not affect the sign of $d^2I/dt^2$. Another interpretation involves the projection of tilted inward motions along filaments threaded by magnetic flux tubes (Balsara et al. 2001), which can be related to the slow, subparsec-scale accretion flows toward the core along the north–south CO filament (Falgarone et al. 2001). The blueshifted and redshifted inward motions may lead to a kink of magnetic fields that contributes to the deviations of POL-2 polarization vectors from being uniform, as shown in Figure 3(b).

The plane-of-sky magnetic energy ($E_{\text{mag.pos}}$) in Table 2 is calculated with the DCF-derived plane-of-sky B-field strength ($B_{\text{pos}}$) of $18 \mu\text{G}$. Because $E_{\text{vir}} + E_{\text{mag.pos}} = -0.35|E_{\text{grav}}| < 0$, the L 1512 core is virially bound if $E_{\text{mag.los}} = 0$, and further contraction could happen. However, the molecular spectral line observations toward L 1512 do not suggest that L 1512 is a contracting core but instead suggest it may have an oscillating envelope. Using the high-density tracer $N_2H^+$ (1–0), Lin et al. (2020) performed non-LTE radiative transfer modeling and found no significant infall in the L 1512 core region. Employing the same radiative transfer model, we estimate an upper limit for the radial infall velocity of ~0.04 km s⁻¹ by analyzing their $N_2H^+$ data, indicating that the core region is relatively quiescent. Lin et al. also found that fitting their multiline observations of $N_2H^+$, $N_2D^+$, DCO+, and $\text{O}_2\text{H}^+$ does not require an infall velocity field in the radiative transfer model, even though the hyperfine structures were carefully considered. Additionally, the envelope of L 1512 was suggested to be oscillating because a mixture of blue and red asymmetric spectral line profiles was observed across the entire cloud in the CS (2–1) line (Lee & Myers 1999; Lee et al. 2001; Lee & Myers 2011), which is a low-density envelope tracer significantly depleted in the central core region. Another envelope tracer, HCN (1–0), which suffers fewer depletion effects compared with CS, also shows a mixed spectral feature across the observing area (Sohn et al. 2007; Kim et al. 2016). Although Schnee et al. (2013) found a red asymmetric feature from their HCO⁺ (3–2) observation, indicating an outward motion, their single-pointing observation does not conflict with the oscillation notion. Based on these spectral observations, L 1512 is likely a long-lived starless core, which is consistent with the core lifetime estimated to be longer than 1.4 Myr according to deuteration chemical modeling (Lin et al. 2020). Therefore, we expect the L 1512 core to be approximately virially stable. However, neither the kinetic pressure nor

magnetic pressure of the plane-of-sky B-field alone can support the L 1512 core (i.e., $2T < |E_{\text{grav}}|$ and $E_{\text{mag.pos}} < |E_{\text{grav}}|$). Given that the magnetic field seems to thread from the large scale to the core scale in L 1512, and the plane-of-sky magnetic field maintains a well-ordered field pattern, we speculate that the magnetic field is not entirely dynamically unimportant in the L 1512 core. It can be the case that the magnetic field does not lie close to the sky plane, and a hitherto hidden line-of-sight magnetic field provides additional support, making L 1512 nearly stable. While Zeeman measurements of L 1512 would help to estimate the line-of-sight B-field strength, such a measurement is currently unavailable.

If the magnetic field is strong enough for the total magnetic energy to compensate for the negative $E_{\text{vir}}$ value, leading the L 1512 core to be virially stable (i.e., $0 = \frac{1}{2}\frac{d^2I}{dt^2} = E_{\text{vir}} + E_{\text{mag}}$ implies $E_{\text{mag}} = 0.51|E_{\text{grav}}|$; see Appendix), the total B-field strength ($B_{\text{tot}}$) would need to be ~32 $\mu\text{G}$. This strength is within the range measured in other starless cores derived via the DCF technique (e.g., Kirk et al. 2006; Myers & Basu 2021; Pattle et al. 2021) and via the Zeeman effect technique (e.g., Crutcher & Troland 2000; Troland & Crutcher 2008). In this case, the line-of-sight B-field component ($B_{\text{los}} = \sqrt{B_{\text{tot}}^2 - B_{\text{pos}}^2}$) is estimated to be ~27 $\mu\text{G}$, and the inclination angle ($i = \sin^{-1}(B_{\text{pos}}/B_{\text{tot}})$) of total B-field direction is ~34° with respect to the line of sight. We note that the derived values of $B_{\text{tot}}$ and $B_{\text{los}}$, as well as the DCF-derived $B_{\text{pos}}$ strength and energy budgets, carry certain uncertainties. It is better to consider that, if $B_{\text{los}}$ is present with a similar order of magnitude to $B_{\text{pos}} = 18 \pm 7 \mu\text{G}$, these values suggest an approximate virial stability of L 1512. By adopting $N_{\text{H}_2} = \langle N_{\text{H}_2} \rangle \cos(i)$, and $B = B_{\text{tot}}$ in Equation (10), the corresponding mass-to-flux ratio ($\lambda_{\text{tot}}$) is ~1.6, suggesting an approximately magnetically critical or sightly supercritical condition. Although $B_{\text{tot}}$ does not make the L 1512 core magnetically subcritical, the magnetic pressure and the kinetic pressure are of comparable importance ($2T \sim E_{\text{mag}}$ and $2T + E_{\text{mag}} \sim |E_{\text{grav}}|$) in supporting the core. On the other hand, if the magnetic fields in the L 1512 core do not have a line-of-sight component or just have a small $B_{\text{los}}$, L 1512 would be magnetically supercritical and should be collapsing. However, the aforementioned spectral observations show this is not the case. Therefore, either L 1512 has just recently reached supercriticality and will collapse at any time and we happened to observe it in this special state, or L 1512 is nearly stable and there is an as-yet unseen line-of-sight B-field.

## 5.2. Relationship between Large- to Core-scale Magnetic Fields

The H-band polarimetry is an important tool to trace the magnetic field at scales between those observed with POL-2 and with Planck. Our POL-2 850 $\mu$m, Mimir H band, and the Planck polarization observations enable efficient plane-of-sky B-field characterization across the small, intermediate, and large scales of the L 1512 cloud. As shown in Figures 3 and 5, the magnetic field orientation of the L 1512 envelope (with an average field angle of $\theta_H = -15^{\circ}+40^{\circ}$) appears to be inherited from that of the large-scale B-field ($\theta_{\text{Planck}} = -30^{\circ}$) and reveals a twist in the B-field morphology between the core and the envelope. The twist occurs in the southwestern core region, where the field angle $\theta_{\text{POL2}}$ bends from ~0° in the core center
to ≈−30°, aligning with the nearby envelope-scale field at $\theta_H$ ≈ −30° (see Section 4.2).

While the 3D field geometry of the field lines remains unknown, the observed twisted field may hint that the matter altered the field at the core scale of ~0.1 pc, with the field orientation still close to the initial large-scale field. This suggests that L 1512 is likely in an intermediate phase, transitioning from the magnetically dominated phase (i.e., magnetically subcritical phase) to the matter-dominated phase (i.e., magnetically supercritical phase). This intermediate phase was proposed by Ward-Thompson et al. (2023) based on their polarimetric observations of nine starless cores embedded within the L 1495A-B10 filaments. The authors found that the plane-of-sky core-scale B-field orientations of these cores are roughly perpendicular to the filaments. However, they are not correlated with the large-scale B-field orientations measured by Planck, except for the lowest-density and possibly youngest core, where the core-scale field is still close to the large-scale field. In this case, a twisted field may be due to the early mass accumulation in the core, and the local field would become perpendicular to the major axis of the core if the gravitational instability is further enhanced.

On the other hand, twisted magnetic fields are also observed in protostellar sources such as filamentary gas flows on the subparsec scale in Serpens South (Pillai et al. 2020) and a twist of ~45° within the inner ~0.005 pc region of L 483, a Class 0 protobinary core (Cox et al. 2022). These authors suggested that these twists result from the gas flow feeding onto the nearby cluster-forming regions and interactions involving binary systems. Thus, the aforementioned protostellar twists are formed in the magnetically supercritical phase, where gravity can efficiently influence the local magnetic field orientations. In contrast, the twist observed in L 1512 is likely developed during the intermediate phase.

In terms of kinematics, the L 1512 core exhibits quiescent gas motions. With the lack of significant gravitational instability, the magnetic field may be dynamically important in the core evolution of L 1512. Moreover, since the L 1512 cloud was reported to be magnetically subcritical based on R-band data ($\lambda_{obs}$ ~0.8 at the scale of ~0.8 pc; Sharma et al. 2022), the L 1512 core ($\lambda_{obs}$ = 3.5 ± 2.4 and $\lambda_{tot}$ ~ 1.6 at the scale of ~0.1 pc) may have undergone a sub-to-supercritical transition, through a process such as ambipolar diffusion.

As mentioned in Section 1, only a few starless cores have resolved submillimeter polarization detections. Among these cores, L 183 (Clemens 2012; Karoly et al. 2020), FeSt 1-457 (Alves et al. 2014; Kandori et al. 2017), and L 1544 (Clemens et al. 2016) also had polarization detection in the H band. In the cases of L 1544 and FeSt 1-457, despite the presence of some nonuniform B-field structures (Alves et al. 2014; Clemens et al. 2016), the orientations of their core-, envelope-, and large-scale magnetic fields remain roughly consistent, similar to L 1512. In contrast, L 183 is the only core in this sample where the orientation of the core-scale B-field tends to be perpendicular to that of the large-scale B-field; the transition of the B-field in the L 183 from its envelope to the core was captured by the H-band polarization measurements (Clemens 2012; Karoly et al. 2020). The apparently similar orientations in the large- and core-scale magnetic fields in L 1512, L 1544, and FeSt 1-457 (unless due to a projection effect) suggest that material can accrete onto these cores along magnetic field lines more efficiently than in the case in L 183. Such core formation scenario has been

demonstrated by Chen et al. (2020) with their turbulent MHD simulations, suggesting that dense cores could accumulate more mass when the core-scale B-fields are aligned with the parsec-scale B-field. However, the L 183 cloud has a total mass of ~80 $M_\odot$ (Pagani et al. 2004), which is considerably larger than the total mass of the L 1544 cloud, estimated to be ~10 $M_\odot$ (Kim et al. 2022). It could be possible that the L 183 cloud is older, enabling it to accumulate sufficient mass, or that the L 183 cloud was supplied with a substantial ancient mass reservoir.

The potential transition from subcritical to supercritical magnetic conditions through ambipolar diffusion is a shared characteristic among L 1512, L 1544, and FeSt 1-457. This phenomenon has been proposed specifically for L 1544 (Ciolek & Basu 2000; Li et al. 2002) and FeSt 1-457 (Kandori et al. 2020c; Bino & Basu 2021). In terms of kinematics, L 1544 exhibits extended inward motions (Tafalla et al. 1998), whereas L 1512 and FeSt 1-457 have quiescent cores and oscillating envelopes (Aguti et al. 2007; Lee & Myers 2011; Juárez et al. 2017; Lin et al. 2020). Moreover, FeSt 1-457 is suggested to be supported by both kinetic pressure and magnetic pressure (Kandori et al. 2018), similar to L 1512. In contrast, the entire L 183 core is found to be subcritical according to POL-2 observations (Karoly et al. 2020), suggesting the dominance of ambipolar diffusion in its core evolution. This could be consistent with the absence of inward motions for the L 183 core (Pagani et al. 2007) and that its surrounding envelope is suggested to be oscillating (Schnee et al. 2013). Despite that, the L 183 core has developed a central density ($2.3 \times 10^6$ cm⁻³; Pagani et al. 2007) comparable to that of L 1544 ($8.6 \times 10^6$ cm⁻³; Keto et al. 2015; Sipilä et al. 2022). Further observations of more starless cores remain vital to shed light on their core formation process and to better understand their environmental influences.

## 6. Conclusions

We present JCMT POL-2 850 µm dust continuum polarization observations and Mimir H-band NIR polarization observations toward L 1512. Our observations reveal an ordered core-scale B-field morphology in L 1512. From our analysis, we find the following:

1. The L 1512 850 µm data, as obtained, likely suffer from missing large-scale flux for the POL-2 data collection. We found that PCA in the standard reduction process removed extended emission, resulting in apparent non-detection of the total intensity. By including a SCUBA-2 Stokes *I* map in the reduction procedure, POL-2 Stokes *Q* and *U* maps could be correctly recovered.

2. The magnetic field traced by POL-2 850 µm, Mimir H band, AIMPOL R band, and Planck polarization data are in agreement as to the average field orientation, suggesting that the large-scale B-field threads the L 1512 cloud down into the dense core region. The largest angular dispersion, found in Mimir H-band data, indicates that a transition of B-field morphology could be happening at the envelope-scale.

3. Ricean-mean modeling of the nondebiased polarization fraction data yielded a power-law index $\alpha$ of 0.48 ± 0.05 in the $p' \propto I^{-\alpha}$ relation, indicating the dust grains retain substantial alignment with the magnetic field at the higher densities within the core.
4. A DCF analysis revealed a plane-of-sky B-field strength of 18 ± 7 µG, and a mass-to-flux ratio of $\lambda_{obs} = 3.5 \pm 2.4$ or a range of 1.1–5.9, suggesting that L 1512 is magnetically supercritical; however, the true mass-to-flux ratio may be being overestimated by $\lambda_{obs}$, and the magnetically critical condition is not entirely ruled out. Given the absence of significant inward motions, and the presence of a well-ordered core-scale B-field and an oscillating envelope, it is likely that L 1512 is supported by both magnetic and kinetic pressures. By assuming L 1512 is virially stable and including the kinetic energy, we estimated that a total B-field strength of ~32 µG could support the L 1512 core against gravity, suggesting a corresponding mass-to-flux ratio of ~1.6. This requires a hitherto hidden line-of-sight B-field component of ~27 µG, which could be sought using Zeeman effect techniques.

5. Alternatively, if there is little to no line-of-sight B-field, then L 1512 should be collapsing. In this case, L 1512 may have just recently reached supercriticality and will collapse at any time.

## Acknowledgments

We thank the referee for providing thoughtful and constructive feedback that helped to improve this article. We also thank Tyler Bourke for his comments on this manuscript that helped to improve it. S.J.L. and S.P.L. acknowledge grants from the National Science and Technology Council (NSTC) of Taiwan 109-2112-M-007-010-MY3, 112-2112-M-007-011, 111-2124-M-001-005 and 112-2124-M-001-014. K.P. is a Royal Society University Research Fellow, supported by grant no. URF\R1\211322. D.W.-T. acknowledges UK Science and Technology Facilities Council (STFC) support under grant No. ST/R000786/1. This work used high-performance computing facilities operated by the Center for Informatics and Computation in Astronomy (CICA) at National Tsing Hua University. This equipment was funded by the Ministry of Education of Taiwan, the Ministry of Science and Technology of Taiwan, and National Tsing Hua University. Nawfel Bouflous and Patrick Hudelot (TERAPIX data center, IAP, Paris, France) are warmly thanked for their help in preparing the CFHT/WIRCAM observation scenario and for performing the data reduction. This work was supported by the Programme National, Physique et Chimie du Milieu Interstellaire (PCMI), of CNRS/INSU with l'Institut de Chimie (INC)/l'Institut de Physique (INP) cofunded by le Commissariat à l'énergie atomique et aux énergies alternatives (CEA) and le Centre national d'études spatiales (CNES) and by Action Fédératrice Astrochimie de l'Observatoire de Paris. The James Clerk Maxwell Telescope is operated by the East Asian Observatory on behalf of The National Astronomical Observatory of Japan; Academia Sinica Institute of Astronomy and Astrophysics; the Korea Astronomy and Space Science Institute; the National Astronomical Research Institute of Thailand; Center for Astronomical Mega-Science (as well as the National Key R&D Program of China with No. 2017YFA0402700). Additional funding support is provided by the Science and Technology Facilities Council of the United Kingdom and participating universities and organizations in the United Kingdom and Canada. Additional funds for the construction of SCUBA-2 were provided by the Canada Foundation for

Innovation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. This research was conducted in part using the Mimir instrument, jointly developed at Boston University and Lowell Observatory and supported by NASA, NSF, and the W.M. Keck Foundation. Analysis software for Mimir data was developed under NSF grants AST 06-07500 and 09-07790 to Boston University. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency. This research has made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

**Facilities:** James Clerk Maxwell Telescope (JCMT), Perkins Telescope Observatory.

**Software:** Starlink (Currie et al. 2014), Mimir Software Package (Clemens et al. 2012), APLpy (Robitaille & Bressert 2012; Robitaille 2019), Astropy (Astropy Collaboration et al. 2013, 2018), SciPy (Virtanen et al. 2020).

## Appendix
Virial Analysis

Here, we present the equations for calculating each energy term in the virial equation. For the ith layer in our onion model, we denote the H₂ number density as $n_{H_2,i}$, the temperature as $T_{kin,i}$, the inner and outer radii as $r_{i-1}$ and $r_i$, and the midpoint radius of the layer as $r_{i-1/2} = \frac{1}{2}(r_{i-1} + r_i)$, where $i = 1,...,N$, $r_0 = 0$, and $r_N = R$. The layer width is denoted as $\Delta r$.

We assume the gas is composed of molecular hydrogen, helium, and metals with the mass fractions of $X = 0.7110$, $Y = 0.2741$, and $Z = 0.0149$, respectively (Lodders 2003). The total number density, $n$, can be expressed as

$$ 
\begin{aligned}
n &= n_{H_2} + n_{He} + n_{metal} = \frac{\rho X}{2m_H} + \frac{\rho Y}{4m_H} + \frac{\rho Z}{15.5m_H} \\
&= \frac{\rho}{m_H} \left( \frac{X}{2} + \frac{Y}{4} + \frac{Z}{15.5} \right) = \frac{\rho}{m_H} \frac{1}{\mu},
\end{aligned} 
\quad (A1) $$

where $\rho$ is volume mass density, $m_H$ is the mass of the hydrogen atom, and $\mu$ is the mean molecular weight per particle ($\mu = 2.35$). The mean molecular weight per molecular hydrogen ($\mu_{H_2} = 2.81$) can be found by

$$ n_{H_2} = \frac{\rho X}{2m_H} = \frac{\rho}{m_H} \frac{1}{\mu_{H_2}}. \quad (A2) $$

The volume mass density, $\rho$, can be expressed by either the total number density, $n$, or the H₂ number density, $n_{H_2}$, with

$$ \rho = \mu n m_H \quad (A3) $$

$$ = \mu_{H_2} n_{H_2} m_H. \quad (A4) $$

The total number density, $n$, can be expressed by the H₂ number density, $n_{H_2}$, with

$$ n = \frac{\mu_{H_2}}{\mu} n_{H_2}. \quad (A5) $$
Another widely assumed gas composition is $n_{H_2} = 5n_{He}$, where metals are negligible (e.g., Myers 1983; and also the DCF formula of Equation (8)). The corresponding mean molecular weights under this assumption are $\mu = 2.33$, and $\mu_{H_2} = 2.8$. Although we adopt a different but more realistic gas composition in the virial analysis from the one Crutcher et al. (2004) adopted in the DCF formula, the $\mu_{H_2}$ from the two assumptions are similar, and thus, the derived volume mass density (Equation (A4)) would not have any significant difference.

We calculate the gravitational energy by

$$E_{\text{grav}} = -4\pi \int_{0}^{R} GM(r)\rho(r)rdr \quad (A6)$$

$$\approx -4\pi G \mu_{H_2} m_{H_2} \sum_{i=1}^{N} M(r_{i-\frac{1}{2}}) n_{H_2,i} r_{i-\frac{1}{2}} \Delta r, \quad (A7)$$

where G is the gravitational constant, and $M(r)$ is the fractional core mass inside a radius of r.

We calculate the thermal energy by

$$E_{\text{thermal}} = 4\pi \int_{0}^{R} \frac{3}{2} n(r) k_{B} T_{\text{kin}}(r) r^{2} dr \quad (A8)$$

$$\approx 4\pi \frac{3\mu_{H_2} k_B}{2\mu} \sum_{i=1}^{N} n_{H_2,i} T_{\text{kin},i} r_{i-\frac{1}{2}}^2 \Delta r, \quad (A9)$$

where $k_B$ is the Boltzmann constant.

The surface kinetic term due to external energy, $P_{\text{ext}}$, is

$$E_{\text{surface}} = \frac{3}{2} P_{\text{ext}} V = 2\pi P_{\text{ext}} R^3, \quad (A10)$$

where V is the core volume. Here, we assume that $P_{\text{ext}} = n_{i=N}k_B T_{\text{kin},i=N}$ is equal to the thermal pressure in the outermost layer.

In our onion model, the 1D nonthermal velocity dispersion, $\sigma_{v,NT}$, is a constant value of 0.046 km s⁻¹. The turbulent energy is calculated by

$$E_{\text{turb}} = \frac{3}{2} M \sigma_{v,NT}^2, \quad (A11)$$

where M is the total core mass.

The rotational energy is calculated by

$$E_{\text{rot}} = \frac{1}{2} \sum_{i=1}^{N} I_i \omega_i^2 = \frac{1}{2} \sum_{i=1}^{N} \frac{2}{5} (M(r_i) - M(r_{i-1})) \frac{r_i^5 - r_{i-1}^5}{r_i^3 - r_{i-1}^3} \omega_i^2, \quad (A12)$$

where $I_i$ and $\omega_i$ are the moment of inertia and angular velocity of the i-th layer, respectively.

The magnetic energy is calculated in the cgs units by

$$E_{\text{mag}} = \frac{B^2}{8\pi} V. \quad (A13)$$

## ORCID iDs

Sheng-Jun Lin https://orcid.org/0000-0002-6868-4483

Shih-Ping Lai https://orcid.org/0000-0001-5522-486X

Kate Pattle https://orcid.org/0000-0002-8557-3582

David Berry https://orcid.org/0000-0001-6524-2447

Dan P. Clemens https://orcid.org/0000-0002-9947-4956

Laurent Pagani https://orcid.org/0000-0002-3319-1021

Derek Ward-Thompson https://orcid.org/0000-0003-1140-2761

Travis J. Thieme https://orcid.org/0000-0003-0334-1583

Tao-Chung Ching https://orcid.org/0000-0001-8516-2532

## References

Aguti, E. D., Lada, C. J., Bergin, E. A., Alves, J. F., & Birkinshaw, M. 2007, ApJ, 665, 457

Alves, F. O., Frau, P., Girart, J. M., et al. 2014, A&A, 569, L1

Alves, F. O., Frau, P., Girart, J. M., et al. 2015, A&A, 574, C4

Andersson, B. G., Lazarian, A., & Vaillancourt, J. E. 2015, ARA&A, 53, 501

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Balsara, D., Ward-Thompson, D., & Crutcher, R. M. 2001, MNRAS, 327, 715

Basu, S. 2000, ApJL, 540, L103

Bino, G., & Basu, S. 2021, ApJ, 911, 15

Caselli, P., Benson, P. J., Myers, P. C., & Tafalla, M. 2002, ApJ, 572, 238

Caselli, P., Pineda, J. E., Zhao, B., et al. 2019, ApJ, 874, 89

Chandrasekhar, S., & Fermi, E. 1953, ApJ, 118, 116

Chapin, E. L., Berry, D. S., Gibb, A. G., et al. 2013, MNRAS, 430, 2545

Chapman, N. L., Davidson, J. A., Goldsmith, P. F., et al. 2013, ApJ, 770, 151

Chen, C. Y., Behrens, E. A., Washington, J. E., et al. 2020, MNRAS, 494, 1971

Chen, C. Y., & Ostriker, E. C. 2018, ApJ, 865, 34

Ching, T. C., Lai, S. P., Zhang, Q., et al. 2017, ApJ, 838, 121

Ciolek, G. E., & Basu, S. 2000, ApJ, 529, 925

Clemens, D. P. 2012, ApJ, 748, 18

Clemens, D. P., Cashman, L. R., Cerny, C., et al. 2020, ApJS, 249, 23

Clemens, D. P., Pavel, M. D., & Cashman, L. R. 2012, ApJS, 200, 21

Clemens, D. P., Sarcia, D., Grabau, A., et al. 2007, PASP, 119, 1385

Clemens, D. P., Tassis, K., & Goldsmith, P. F. 2016, ApJ, 833, 176

Cox, E. G., Novay, G., Sadavoy, S. I., et al. 2022, ApJ, 932, 34

Crutcher, R. M. 2012, ARA&A, 50, 29

Crutcher, R. M., Nutter, D. J., Ward-Thompson, D., & Kirk, J. M. 2004, ApJ, 600, 279

Crutcher, R. M., & Troland, T. H. 2000, ApJL, 537, L139

Currie, M. J., Berry, D. S., Jenness, T., et al. 2014, in ASP Conf. Ser. 485, Astronomical Data Analysis Software and Systems XXIII, ed. N. Manset & P. Forshay (San Francisco, CA: ASP), 391

Davis, L. 1951, PhRv, 81, 890

Dempsey, J. T., Friberg, P., Jenness, T., et al. 2013, MNRAS, 430, 2534

Dolginov, A. Z., & Mitrofanov, I. G. 1976, Ap&SS, 43, 291

Dunham, M. M., Offner, S. S. R., Pineda, J. E., et al. 2016, ApJ, 823, 160

Eswaraiah, C., Li, D., Furuya, R. S., et al. 2021, ApJL, 912, L27

Falgarone, E., Panis, J. F., Heithausen, A., et al. 1998, A&A, 331, 669

Falgarone, E., Pety, J., & Phillips, T. G. 2001, ApJ, 555, 178

Federrath, C., & Klessen, R. S. 2012, ApJ, 761, 156

Foster, J. B., & Goodman, A. A. 2006, ApJL, 636, L105

Friberg, P., Bastien, P., Berry, D., et al. 2016, Proc. SPIE, 9914, 991403

Gaia Collaboration 2023, A&A, 674, A1

Galli, D., & Shu, F. H. 1993a, ApJ, 417, 220

Galli, D., & Shu, F. H. 1993b, ApJ, 417, 243

Girart, J. M., Beltrán, M. T., Zhang, Q., Rao, R., & Estalella, R. 2009, Sci, 324, 1408

Girart, J. M., Rao, R., & Marrone, D. P. 2006, Sci, 313, 812

Goodman, A. A., Jones, T. J., Lada, E. A., & Myers, P. C. 1995, ApJ, 448, 748

Heitsch, F., Zweibel, E. G., Mac Low, M. M., Li, P., & Norman, M. L. 2001, ApJ, 561, 800

Holland, W. S., Bintley, D., Chapin, E. L., et al. 2013, MNRAS, 430, 2513

Hull, C. L. H., Plambeck, R. L., Bolatto, A. D., et al. 2013, ApJ, 768, 159

Hull, C. L. H., & Zhang, Q. 2019, FrASS, 6, 3

Jones, T. J., Bagley, M., Krejny, M., Andersson, B. G., & Bastien, P. 2015, AJ, 149, 31

Juárez, C., Girart, J. M., Frau, P., et al. 2017, A&A, 597, A74

Kandori, R., Tamura, M., Kusakabe, N., et al. 2017, ApJ, 845, 32

Kandori, R., Tamura, M., Saito, M., et al. 2020a, ApJ, 890, 14

Kandori, R., Tamura, M., Saito, M., et al. 2020b, PASJ, 72, 8

Kandori, R., Tomisaka, K., Saito, M., et al. 2020c, ApJ, 888, 120

Kandori, R., Tomisaka, K., Tamura, M., et al. 2018, ApJ, 865, 121

Karoly, J., Soam, A., Andersson, B. G., et al. 2020, ApJ, 900, 181
Karoly, J., Ward-Thompson, D., Pattle, K., et al. 2023, ApJ, 952, 29

Kenyon, S. J., Dobrzycka, D., & Hartmann, L. 1994, AJ, 108, 1872

Keto, E., Caselli, P., & Rawlings, J. 2015, MNRAS, 446, 3731

Kim, G., Lee, C. W., Gopinathan, M., Jeong, W. S., & Kim, M. R. 2016, ApJ, 824, 85

Kim, S., Lee, C. W., Tafalla, M., et al. 2022, ApJ, 940, 112

Kirk, H., Dunham, M. M., Di Francesco, J., et al. 2017, ApJ, 838, 114

Kirk, J. M., Ward-Thompson, D., & André, P. 2005, MNRAS, 360, 1506

Kirk, J. M., Ward-Thompson, D., & Crutcher, R. M. 2006, MNRAS, 369, 1445

Kwon, J., Doi, Y., Tamura, M., et al. 2018, ApJ, 859, 4

Launhardt, R., Stutz, A. M., Schmiedeke, A., et al. 2013, A&A, 551, A98

Lazarian, A., & Hoang, T. 2007, MNRAS, 378, 910

Lee, C. W., & Myers, P. C. 1999, ApJS, 123, 233

Lee, C. W., & Myers, P. C. 2011, ApJ, 734, 60

Lee, C. W., Myers, P. C., & Tafalla, M. 2001, ApJS, 136, 703

Li, Z. Y., Shematovich, V. I., Wiebe, D. S., & Shustov, B. M. 2002, ApJ, 569, 792

Li, Z. Y., & Shu, F. H. 1996, ApJ, 472, 211

Lin, S. J., Pagani, L., Lai, S. P., Lefèvre, C., & Lique, F. 2020, A&A, 635, A188

Liu, J., Qiu, K., Berry, D., et al. 2019, ApJ, 877, 43

Lodders, K. 2003, ApJ, 591, 1220

Lombardi, M., Lada, C. J., & Alves, J. 2010, A&A, 512, A67

Luri, X., Brown, A. G. A., Sarro, L. M., et al. 2018, A&A, 616, A9

Lynds, B. T. 1962, ApJS, 7, 1

Mac Low, M. M., & Klessen, R. S. 2004, RvMP, 76, 125

McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565

Mckee, C. F., & Zweibel, E. G. 1992, ApJ, 399, 551

Montier, L., Plaszczynski, S., Levrier, F., et al. 2015, A&A, 574, A136

Mouschovias, T. C. 1991, ApJ, 373, 169

Mouschovias, T. C., Tassis, K., & Kunz, M. W. 2006, ApJ, 646, 1043

Myers, P. C. 1983, ApJ, 270, 105

Myers, P. C., & Basu, S. 2021, ApJ, 917, 35

Myers, P. C., Basu, S., & Auddy, S. 2018, ApJ, 868, 51

Myers, P. C., Linke, R. A., & Benson, P. J. 1983, ApJ, 264, 517

Ostriker, E. C., Stone, J. M., & Gammie, C. F. 2001, ApJ, 546, 980

Padoan, P., & Nordlund, AA. 1999, ApJ, 526, 279

Pagani, L., Bacmann, A., Cabrit, S., & Vastel, C. 2007, A&A, 467, 179

Pagani, L., Bacmann, A., Motte, F., et al. 2004, A&A, 417, 605

Pattle, K., & Fissel, L. 2019, FrASS, 6, 15

Pattle, K., Fissel, L., Tahani, M., Liu, T., & Ntormousi, E. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al., 193

Pattle, K., Lai, S. P., Di Francesco, J., et al. 2021, ApJ, 907, 88

Pattle, K., Lai, S. P., Hasegawa, T., et al. 2019, ApJ, 880, 27

Pattle, K., Ward-Thompson, D., Berry, D., et al. 2017, ApJ, 846, 122

Pillai, T. G. S., Clemens, D. P., Reissl, S., et al. 2020, NatAs, 4, 1195

Planck Collaboration IV 2020, A&A, 641, A4

Planck Collaboration XII 2020, A&A, 641, A12

Remazeilles, M., Delabrouille, J., & Cardoso, J. F. 2011, MNRAS, 418, 467

Rice, S. O. 1945, BSTJ, 24, 46

Robitaille, T., 2019 APLpy v2.0: The Astronomical Plotting Library in Python, Zenodo, doi:10.5281/zenodo.2567476

Robitaille, T., & Bressert, E. 2012, APLpy: Astronomical Plotting Library in Python, Astrophysics Source Code Library, ascl:1208.017

Roccatagliata, V., Franciosini, E., Sacco, G. G., Randich, S., & Sicilia-Aguilar, A. 2020, A&A, 638, A85

Saajasto, M., Juvela, M., Lefèvre, C., Pagani, L., & Ysard, N. 2021, A&A, 647, A109

Schnee, S., Brunetti, N., Di Francesco, J., et al. 2013, ApJ, 777, 121

Serkowski, K. 1958, AcA, 8, 135

Sharma, E., Gopinathan, M., Soam, A., Lee, C. W., & Seshadri, T. R. 2022, MNRAS, 517, 1138

Shu, F. H., Adams, F. C., & Lizano, S. 1987, ARA&A, 25, 23

Sipilä, O., Caselli, P., Redaelli, E., & Spezzano, S. 2022, A&A, 668, A131

Soam, A., Pattle, K., Ward-Thompson, D., et al. 2018, ApJ, 861, 65

Sohn, J., Lee, C. W., Park, Y. S., et al. 2007, ApJ, 664, 928

Tafalla, M., Mardones, D., Myers, P. C., et al. 1998, ApJ, 504, 900

Tassis, K., & Mouschovias, T. C. 2004, ApJ, 616, 283

Tokuda, K., Fujishiro, K., Tachihara, K., et al. 2020, ApJ, 899, 10

Torres, R. M., Loinard, L., Mioduszewski, A. J., & Rodríguez, L. F. 2009, ApJ, 698, 242

Troland, T. H., & Crutcher, R. M. 2008, ApJ, 680, 457

Vaillancourt, J. E. 2006, PASP, 118, 1340

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Wang, J. W., Lai, S. P., Eswaraiah, C., et al. 2017, ApJ, 849, 157

Wang, J. W., Lai, S. P., Eswaraiah, C., et al. 2019, ApJ, 876, 42

Ward-Thompson, D., Karoly, J., Pattle, K., et al. 2023, ApJ, 946, 62

Ward-Thompson, D., Kirk, J. M., Crutcher, R. M., et al. 2000, ApJL, 537, L135

Ward-Thompson, D., Pattle, K., Bastien, P., et al. 2017, ApJ, 842, 66

Ward-Thompson, D., Sen, A. K., Kirk, J. M., & Nutter, J. M. 2009, MNRAS, 398, 394