Magnetic Fields of the Starless Core L 1512

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 ± 7 μG and an observed mass-to-flux ratio (λ obs) of 3.5 ± 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 ∼27 μG with a mass-to-flux ratio (λ tot) of ∼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.


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 soshengjunlin@asiaa.sinica.edu.tw,slai@phys.nthu.edu.twcalled 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,b;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 arXiv:2311.08026v1[astro-ph.GA] 14 Nov 2023 Lin et al. star-forming activities (e.g., Padoan & Nordlund 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 (∼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. 2006Girart et al. , 2009;;Ching et al. 2017).At pseudodisk scales (∼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 ∼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 non-uniformity in the planeof-sky B -field structure around the cores (e.g., Clemens et al. 2016;Kandori et al. 2017Kandori et al. , 2020a,b),b).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 (submm) wavelengths has emerged as a key tool for investigating the 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 µm/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) with APEX-PolKa; and Oph C (Liu et al. 2019), IRAS 16293E, L 1689 SMM-16, and L 1689B (Pattle et al. 2021) with JCMT-POL2.Their core sizes (the full width at half maximum, FWHM, of the submm intensity) range from approximately 0.04 pc to 0.1 pc, for distances spanning from approximately 110 pc to 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 at the parent clumps themselves (e.g., Oph A & B clumps; Kwon et al. 2018;Soam et al. 2018) or the mean Bfield 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 submm 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 submm 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 (SNR) are needed to critically compare to theoretical models.Recently, L 183 and L 43 were re-observed with JCMT-POL2 (Karoly et al. 2020(Karoly et al. , 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 α of the polarization-fraction to total-intensity relation (p ∝ I −α ) 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 SNR 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 < α < 1), a higher SNR threshold in Stokes I , which discards more data, is necessary to reliably recover the α index (see their Fig. 1).Below this SNR threshold, the fitted α 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 SNR threshold, the conventional approach (fitting with a single power-law model of α) would overestimate the α index.Pattle et al. (2019) found that α can be better estimated with a Ricean noise model without debiasing the polarization; i.e., favoring a non-debiased-polarization to total-intensity relation (p ′ -I ).The authors thus revised the α index of the starless core Oph C to 0.6∼0.7 from the value of 1.03 previously found by Liu et al. (2019).In FeSt 1-457, another starless core, α 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 α index in IC 5146 from 1.03∼1.08 to 0.56.In addition, an equivalent relation in NIR is p/A V ∝ A −α V (Andersson et al. 2015;Pattle et al. 2019).It has been found that α < 1 in the dense clouds, even with debiased data (see Wang et al. 2017, 2019, andreferences 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 multi-line observations of N 2 H + , N 2 D + , DCO + , C 18 O, and H 2 D + (1 10 -1 11 ) using non-LTE radiative transfer.With the high-density tracer N 2 H + , they found a low central temperature of 8±1 K, a small threedimensional isotropic non-thermal velocity dispersion of 0.080 km s −1 (corresponding to a one-dimensional nonthermal velocity dispersion, σ v,NT , of 0.046 km s −1 , a non-thermal FWHM linewidth, ∆v NT , of 0.109 km s −1 , and a Mach number, M = σ v,NT /c s , of 0.27, where c s is the isothermal sound speed of 0.17 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 −3 .They also found inward motions toward the L 1512 core along a one-parseclong north-south filament with a low accretion rate of ∼ 4 × 10 −6 M ⊙ yr −1 , 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 Sec. 2, our reduction recipes in Sec. 3, and present our result and analysis in Sec. 4. In Sec. 5, we discuss the energy budgets of L 1512.Lastly, we summarize our results in Sec. 6.

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 1 .They are part of the M13BC01 project (PI: James di Francesco) observed in November 2013.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 Fig. 1a, 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 2 in the Starlink software suite (Chapin et al. 2013).We applied another reduction routine on the raw data to eliminate 2 http://starlink.eao.hawaii.edu/docs/sun258.htx/sun258.htmlthe spurious structures and the details are described in Sec.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 hours, between August 2020 to November 2020 in Band 2 weather (0.05 < τ 225GHz < 0.08) under project code M20BP046 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 Fig. 1b.The blue dashed curve shows the profile of the POL-2 Stokes I map in Fig. 1c.The blue solid curve and filled area show the profile of the POL-2 Stokes I map in Fig. 1d.
(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 Sec.3.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.Seeing conditions were between 1.6-3.4arcsec for all observations.Data were calibrated with the Basic Data Processor (MSP-BDP) and Photo-Polarimetry Analysis Tool (MSP-PPOL) from the Mimir Software Package 3 .The data selection criteria were that the H band magnitude was less than or equal to 13.5 mag, the polarization fraction uncertainty (σ p ) was less than or equal to 0.9% (or σ p ≤ 0.009 in a decimal equivalent), and an SNR of the polarization fraction (p/σ p ) ≥ 2 or σ PA ≤ 15 • , yielding a sample of 31 stars out of the 194 total observed.

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 Archive4 .Planck 353 GHz (850 µm) polarization data were retrieved from the Planck Legacy Archive5 .We employed the Planck 2018 GNLIC Stokes I , Q, and U foreground thermal dust emission maps at variable resolution (5 ′ -80 ′ ) over the sky (PR3; Planck Collaboration IV 2020; Planck Collaboration XII 2020), where the Galatic 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 dataset we used for L 1512 has an effective beam size of 15 ′ .

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

JCMT SCUBA-2 Photometric Map
Figure 1a shows the archival SCUBA-2 image with the possibly spurious extended structures (see Sec. 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 emission6 , 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 1b 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 −1 (Dempsey et al. 2013).

Missing Flux of JCMT POL-2 for Extended Structures
The POL-2 polarization data were first reduced with the standard two-stage reduction process7 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 1c presents the POL-2 total intensity map reduced with the above standard procedure, while Figure 1d 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 (Fig. 1b, peak at 70 mJy beam −1 ), we found that the dust emission of L 1512 is undetected in Stokes I (Fig. 1c) and neither Q nor U are detected with the standard POL-2 data reduction.Given the rms noise of 3.1 mJy beam −1 measured in the central 3 arcmin field in the POL-2 map using the standard reduction, the dust emission peak of L 1512 should be detected with an SNR 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 halfwave 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 −1 ), 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., over-sampled 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 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 (Fig. 1b) but failed for the POL-2 map (Fig. 1c).We re-reduced 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 Fig. 1d.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 (Fig. 1b).
Figure 2 shows the azimuthally-averaged intensity profiles of our improved SCUBA-2 map (Fig. 1b) compared to the (no PCA) POL-2 total intensity map (Fig. 1d).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 non-uniform 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 −1 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 SNR to enable reducing the Stokes Q and U signals.We thus used the reprocessed SCUBA-2 map (Fig. 1b and Fig. 3a 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 −1 (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 polarization studies by Ward-Thompson et al. (2017); Pattle 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 procedure, 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  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 submm wavelengths is perpendicular to the plane-of-sky magnetic field, while the dust-extincted starlight polarization in the NIR is parallel to the planeof-sky magnetic field (e.g., Andersson et al. 2015;Pattle et al. 2023).Thus the vectors overlaid on Fig. 3a and 3b were rotated 90 • so they may represent the planeof-sky B -field direction.The vectors have been then filtered, based on the SNR of the polarization fraction being larger than three (p/σ p ≥ 3) and the SNR of the total intensity larger than ten (I/σ I ≥ 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 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 Sec. 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.

Polarization Properties
The formulae for computing the polarization properties from the reduced data are as follows: the nondebiased polarization fraction (p ′ ) is given by where I , Q, U are the Stokes parameters measured.We let σ I , σ Q , and σ U be their uncertainties.Uncertainty of the polarization fraction (σ p ) is given by The polarization fraction is debiased using the asymptotic estimator (Vaillancourt 2006;Montier et al. 2015), and is rearranged by adopting the aforementioned σ Q = σ U = 0.79 mJy beam −1 and ignoring the I −4 term of the radicand in Equation 2 (Kwon et al. 2018) as where p is the debiased polarization fraction.The polarization angle (PA) is given by and the uncertainty of the polarization angle (σ PA ) by 4. RESULTS AND ANALYSIS

Foreground Star Census
The dust-extincted starlight polarization in the Mimir H band from the background stars can trace the planeof-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 (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 σ d = σ π /π 2 (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 5a shows the Mimir H band polarization vectors superimposed on the CFHT H band image.A grey 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. = 5 h 04 m 07.s 5 and decl.= 32 • 43 ′ 25. ′′ 0, J2000; Lin et al. 2020) and continue out to 410 ′′ .Therefore, their polarization angles 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.The Mimir vectors mostly trace the magnetic field morphology in the diffuse envelope, with a spatial scale of ∼0.56 pc down to ∼0.18 pc (see Sec. 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 over-sampled for visualization, as the entire field of view (11.5 ′ ×12 ′ ) of Fig. 5a fits within one effective Planck beam size of 15 ′ .The Planck data indicate the large-scale mean field angle of θ Planck = −30 • over a spatial scale of ∼0.6 pc at the distance of 140 pc.

Magnetic Field Morphology
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 polarization angles around the Planck mean field angle (i.e., |θ H − θ Planck |) 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-to envelope-scales.These deviations of the NIR vectors seem to show the effects of a relative motion be-  tween 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 Sec.3.2, Figure 3b 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 (θ POL2 ≈ 0 • ).We note that the pattern exhibits a mostly smooth change in position angle from θ POL2 ≈ 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 5c shows a spatial comparison between POL-2 vectors and Mimir vectors, where the POL-2 vectors are identical to those in Fig. 3b but undersampled by a factor of two for clarity.Their distinct different spatial coverages indeed show that the Mimir data 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 (p.d.f.) are derived to account for the uncertainties (σ PA ) 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 p.d.f. 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 Fig. 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.
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 polarimetric 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 (Figs.5b and 5c) 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 Fig. 5c.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 pc 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 Fig. 6b) 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.

Davis-Chandrasekhar-Fermi Analysis
The Davis-Chandrasekhar-Fermi (DCF;Davis 1951;Chandrasekhar & Fermi 1953) method is widely used to derive the plane-of-sky B -field strength (B 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 non-thermal gas motions result in a distortion of the magnetic field.By measuring the dispersions in the non-thermal gas velocities and in the polarization position angles, the field strength is (Crutcher et al. 2004): where ⟨ρ⟩ is the mean volume mass density, ⟨n H2 ⟩ is the mean H 2 number density (⟨ρ⟩ = µ H2 ⟨n H2 ⟩m H , and µ H2 = 2.8), σ v,NT is the 1D non-thermal velocity dispersion, ∆v NT is the non-thermal FWHM linewidth (∆v NT = σ v,NT √ 8 log 2), δθ 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 pos , if δθ < 25 • .

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 θ POL2 ≈ 0 • to θ POL2 ≈ −30 • in the southwestern side (see Fig. 3b).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 Figs.7  and 8. Figure 7a shows the B -field position angle (θ) 12 ′′ -pixel map, where the vectors are shown with uniform lengths for clarity.Figure 7d shows the corresponding B -field-angle uncertainty (σ PA ) map.We make a map of the underlying B -field by smoothing Fig. 7a  with a 3-by-3-pixel boxcar filter (i.e., averaging over the nearby nine pixels).This smoothed B -field angle (⟨θ⟩ local ) map is shown in Fig. 7b. Figure 7c shows the residual (∆θ = θ − ⟨θ⟩ local ) map, in which three isolated pixels from the θ and ⟨θ⟩ local maps are rejected and 104 pixels remained.
On the residual angle map (∆θ), the underlying ordered B -field geometry has been removed.In order to estimate a representative intrinsic B -field angular dispersion (δθ) 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 δθ is via the inverse-variance-weighted quadratic mean of the residual angles, where w i = 1/σ 2 PA,i (e.g., Clemens et al. 2016;Wang et al. 2019).Hence the well-characterized residual angles with smaller σ PA dominate the estimation.We find δθ = 8.8 • ± 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 δθ can be better recovered, closer to the input δθ value in their simulations (their socalled "true value"), via averaging a set of standard deviation estimates of the residual angles selected with different maximum allowed PA uncertainties (σ PA,max ); i.e., ⟨σ(∆θ|uncertainties ≤ σ PA,max )⟩. Figure 8 shows these standard deviation estimates as the function of σ PA,max , with black dots, and the cumulative number of the residual angles limited by σ PA,max , with red dots.Figure 7e shows the angular uncertainty of ∆θ, denoted  2017) ran a set of Monte Carlo simulations by setting different "true values" of the intrinsic B -field angular dispersion (δθ) to generate their simulated datasets.They found the above procedure can approach the input δθ when σ PA,max is small and the cumulative number of ∆θ is still sufficient.They took the mean of such well-characterized standard deviation estimates as the representative estimation of δθ for the whole region.In our analysis, we collect the standard deviation estimates computed with 6 • ≤ σ PA,max ≤ 9.3 • (two vertical dashed black lines in Fig. 8) because these σ(∆θ) estimates remain relatively similar values in contrast to the sharp σ(∆θ) decline at σ PA,max < 6 • , 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 σ(∆θ) estimates, we obtain an intrinsic B -field angular dispersion (δθ) of 8.2 • ± 1.4 • (horizontal blue lines in Fig. 8) for the whole region.Compared to the conventional approach, both δθ estimations are consistent, and we adopt δθ = 8.2 • ±1.4 • for the following analysis.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 ) η ), where n 0 is the central density, R 0 is the characteristic radius, and η is the power-law index of the density profile at r ≫ R 0 .The best-fit Plummer parameters fitted to N 2 H + spectral line observations and dust extinction measurements were found to be n 0 = 1.1 × 10 5 cm −3 , R 0 = 0.022 pc, and η = 2.0 for the east side, and n 0 = 1.1 × 10 5 cm −3 , R 0 = 0.027 pc, and η = 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 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 edge to obtain a column density (N H2 ) map, and then (2) divide the column densities by the corresponding sightline depth inside the sphere to obtain the line-of-sightaveraged H 2 number density (n H2,los ) map (i.e., for a pixel indexed j, n H2,los,j = N H2,j /(2 R 2 edge,cm − d 2 j,cm ), Note-Quantities shown are the dispersion in magnetic field position angles (δθ), the non-thermal FWHM linewidth (∆vNT), the mean H2 number density (⟨nH 2 ⟩), the mean H2 column density (⟨NH 2 ⟩), the plane-of-sky B -field strength (Bpos), the observed mass-to-flux ratio (λ obs ), the statistically corrected total B -field strength (|Bcor|), the statistically corrected mass-to-flux ratio (λcor), and the total/lineof-sight B -field strength (Btot, B los ) and the corresponding mass-to-flux ratio (λtot) are derived by assuming the L 1512 core is virially stable (see Sec. 5.1).The uncertainties are their measured dispersion or computed by the standard error propagation.For ∆vNT, we adopt the N2H + (1-0) spectral resolution from Lin et al. (2020) as the uncertainty.For Btot, B los , and λtot, we do not estimate uncertainties because their derivation involves energy budgets (Table 2) that may be accurate only to order of magnitude.

Number Density and Velocity Dispersion
where R edge,cm and d j,cm are in the units of centimeter, and d j,cm is the projected distance to the core center).
Figure 9 shows the N H2 map and n H2,los map.For selfconsistency, we compute the mean, and standard deviation, of n H2,los and N H2 values on the same set of pixels, of which the residual angles (Fig. 7c) are used to estimate the representative B -field angular dispersion δθ.We find ⟨n H2 ⟩ = ⟨n H2,los ⟩ = (2.1 ± 0.9) × 10 4 cm −3 , and ⟨N H2 ⟩ = (8.3± 4.7) × 10 21 cm −2 .Lin et al. (2020) found that a non-thermal FWHM linewidth of ∆v NT = 0.109 km s −1 can reproduce their N 2 H + (1-0) spectral observations.We take their N 2 H + (1-0) spectral resolution of 0.031 km s −1 as the uncertainty.We adopt these values for the DCF analysis and list them in Table 1.

Magnetic Field Strength and Mass-to-flux Ratio
Using the DCF method (Equation 8), and the aboveestimated values (see Table 1), we calculate the mean plane-of-sky B -field strength (B pos ) across the core to be 18 ± 7 µG, where the uncertainty is computed with the standard error propagation from the dispersions of δθ, ∆v NT , and n H2 and hence the uncertainty of B 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/Φ) crit .We use the formula from Crutcher et al. (2004), The core is magnetically supercritical if λ > 1 (i.e., unstable to collapse) and magnetically subcritical if λ < 1 (i.e., magnetically supported).With our derived ⟨N H2 ⟩ and B pos , we calculate the observed mass-to-flux ratio (λ obs ) as 3.5 ± 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 λ 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 λ 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 non-thermal energies are not considered.Crutcher et al. (2004) showed that B pos and λ obs can be statistically corrected by averaging the cases of all Bfield 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π subtended from the core center.The correction is In this case, the total magnetic field would have a strength of 23 ± 9 µG.For our observed mass-to-flux ratio, λ obs , Crutcher et al. (2004) derived the correction as This correction also takes into account the projection effect of N H2 , assuming that the magnetic flux tube is aligned with the core minor axis.This yields a mass-toflux ratio of λ cor = 1.2 ± 0.8, still suggesting the core is approximately magnetically critical or slightly super-  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.
critical.We will further discuss the B -field strength in the context of the Virial analysis in Sec. 5.

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;Jones et al. 2015;Andersson et al. 2015;Wang et al. 2019;Pattle 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 relation8 , is sought, where 0 ≤ α ≤ 1 (e.g., Jones et al. 2015;Alves et al. 2015).In the α = 0 case, the polarization fraction is constant at every depth along each sightline, indicating a perfectly aligned case.For the α = 1 case, the polarized intensity (= p true •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 α, 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 SNR criterion for fitting α.In fact, these polarization estimators are biased estimators but they work better in the high-SNR regime because the bias is minimized (Montier et al. 2015); therefore, the low-SNR data should be removed in order to better determine α.For faint starless cores, such SNR criteria may remove too much data, and thus α might not be well-constrained.
Instead of debiasing data and removing the low-SNR ones, Pattle et al. (2019) assumed that the non-biased polarization fraction measurements follow a Rice distribution, including both low-and high-SNR 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, where L 1 2 is a Laguerre polynomial of order 1/2, σ QU is the average of σ Q and σ U , α is the power-law index, and p σ QU is the expected polarization fraction when I = σ QU .We note that technically, p σ QU and σ QU are the free scaling factors of the assumed underlying p − I relation of p true = p σ QU (I/σ QU ) −α , where the I scaling factor is chosen to be the practical noise (σ QU ) while the polarization fraction scaling factor (p σ QU ) is left for fitting.This Ricean estimator allows for a better estimation of the index α by being able to properly include the low-SNR data.
We use our data to determine the index α in L 1512 with the above Ricean-mean model.We collect data within a 3 ′ uniform noise field (σ QU = 0.79 mJy beam −1 ), and supply σ p as a data weighting in order to perform a fitting with the Python scipy.optimize.curvefit function.
We obtain a best-fit α of 0.48 ± 0.05, and p σ QU of 0.60 ± 0.10.Figure 10 shows the best-fit Ricean-mean model in the solid red curve.The true value of α = 0.48 ± 0.05 is recoverable above the critical intensity (I crit ) of 3.22 mJy beam −1 (corresponding to a SNR of I crit /σ QU = 4.05).In contrast, owing to the low SNR, the true α is not recoverable below I crit but apparently approaches the value of unity expected for disaligned grains, shown as the dashed black line.
The index of α = 0.48 ± 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.

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 1 2 where I is the moment of inertia, E thermal is the thermal energy in the core, E surface is the surface kinetic term due to external pressure, E turb is the turbulent energy, E rot is the rotational energy, E grav is the gravitational potential energy, and E mag is the magnetic energy calculated by the total B -field strength (B tot ).The magnetic energy can be divided into the plane-of-sky and line-ofsight components such that E mag = E mag,pos + E mag,los since B 2 tot = B 2 pos + B 2 los .The total kinetic energy can be defined as T = E thermal − E surface + E turb + E rot .By defining the virial energy of E vir = 2T + E grav , we rewrite Equation 16as 1 2 The sign of d 2 I/dt 2 determines whether the core is virially unbound with a positive net energy (d 2 I/dt 2 > 0) or virially bound with a negative net energy (d 2 I/dt 2 < 0), and d 2 I/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 2 H + multi-transition 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 2 H + 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 Figs. 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 Fig. 5c from Lin et al. (2020) for the density, temperature, and rotational velocity profiles.In addition, the constant one-dimensional non-thermal velocity dispersion, σ v,NT , used in the onion model is 0.046 km s −1 .We calculate the plane-of-sky component magnetic energy (E mag,pos ) by using the DCF-derived plane-of-sky magnetic field strength (B pos ) of 18 µG (Table 1).Table 2 summarizes the results for each energy of the L 1512 core and Appendix A 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 vir ) will change from −0.51 to −0.39 and the plane-of-sky component magnetic energy (E mag,pos ) will change from 0.16 to 0.11, in units of |E grav |.Accordingly, if E 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 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 grav is dependent on the square of mass.
The rotational energy (E rot ) in Table 2 is one of the interpretations from the N 2 H + (1-0) observation toward L 1512, which reveals a uniform velocity gradient (2.26±0.04km s −1 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 solidbody rotation (Caselli et al. 2002;Fig. 5c in Lin et al. 2020), of which the corresponding rotation period of 2.72±0.05Myr (Ω = (7.3± 0.1) × 10 −14 rad s −1 and v rot ≤ 0.07 km s −1 ) is about twice as long as the core lifetime of ≳1.4 Myr (Lin et al. 2020).The corresponding E rot is only 1% of |E grav |.This contribution does not affect the sign of d 2 I/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 blue-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 Fig. 3b.
The plane-of-sky magnetic energy (E mag,pos ) in Table 2 is calculated with the DCF-derived planeof-sky B -field strength (B pos ) of 18 µG.Because E vir + E mag,pos = −0.35|Egrav | < 0, the L 1512 core is virially bound if E 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 2 H + (1-0), Lin et al. (2020) performed non-LTE radiative transfer modeling and found no significant in-fall 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 −1 by analyzing their N 2 H + data, indicating that the core region is relatively quiescent.Lin et al. also found that fitting their multi-line observations of N 2 H + , N 2 D + , DCO + , and o-H 2 D + 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 Bfield alone can support the L 1512 core (i.e., 2T < |E grav | and E mag,pos < |E 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-to-sight magnetic field provides additional support, making L 1512 nearly stable.While Zeeman measurements of L 1512 would help to estimate the line-to-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 vir value, leading the L 1512 core to be virially stable (i.e., 0 = 1 2 d 2 I dt 2 = E vir + E mag implies E mag = 0.51|E grav |; see Appendix A), the total B -field strength (B tot ) would need to be ∼32 µG.This strength is within the range measured in other starless cores derived via the DCF technique (e.g., Kirk et al. 2006;Pattle et al. 2021;Myers & Basu 2021) and via the Zeeman effect technique (e.g., Crutcher & Troland 2000;Troland & Crutcher Lin et al. 2008).In this case, the line-of-sight B -field component (B los = B 2 tot − B 2 pos ) is estimated to be ∼27 µG and the inclination angle (i = sin −1 (B pos /B tot )) of total B -field direction is ∼34 • with respect to the line of sight.We note that the derived values of B tot and B los , as well as the DCF-derived B pos strength and energy budgets, carry certain uncertainties.It is better to consider that if B los is present with a similar order of magnitude to B pos = 18 ± 7 µG, these values suggest an approximate virial stability of L 1512.By adopting N H2 = ⟨N H2 ⟩ cos(i) and B = B tot in Equation 10, the corresponding mass-to-flux ratio (λ tot ) is ∼1.6, suggesting an approximately magnetically critical or sightly supercritical condition.Although B tot does not make the L 1512 core magnetically subcritical, the magnetic pressure and the kinetic pressure are of comparable importance (2T ∼ E mag and 2T + E mag ∼ |E 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 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.

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 µ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 Figs. 3 and 5, the magnetic field orientation of the L 1512 envelope (with an average field angle of θ H = −15 •+40 • −39 • ) appears to be inherited from that of the large-scale B -field (θ Planck = −30 • ) 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 θ POL2 bends from ≈0 • in the core center to ≈−30 • , aligning with the nearby envelope-scale field at θ H ≈ −30 • (see Sec. 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 lowestdensity 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 (λ obs ∼ 0.8 at the scale of ∼0.8 pc; Sharma et al. 2022), the L 1512 core (λ obs = 3.5 ± 2.4 and λ tot ∼ 1.6 at the scale of ∼0.1 pc) may have undergone a sub-to-supercritical transition, through such as ambipolar diffusion.
As mentioned in Sec. 1, only a few starless cores have resolved submm 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 non-uniform B -field structures (Clemens et al. 2016;Alves et al. 2014), 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 tran-sition 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) suggests 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 corescale B -fields are aligned with the parsec-scale B -field.However, the L 183 cloud has a total mass of ∼80 M ⊙ (Pagani et al. 2004), which is considerably larger than the total mass of the L 1544 cloud, estimated to be ∼10 M ⊙ (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×10 6 cm −3 ; Pagani et al. 2007) comparable to that of L 1544 (8.6×10 6 cm −3 ; 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.

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 principal component analysis (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 Bfield 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 non-debiased polarization fraction data yielded a power-law index α of 0.48±0.05 in the p ′ ∝ I −α relation, indicating the dust grains retain substantial alignment with the magnetic field at the higher densities within the core.
4. A Davis-Chandrasekhar-Fermi analysis revealed a plane-of-sky B -field strength of 18±7 µG, and a mass-to-flux ratio of λ obs = 3.5 ± 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 λ 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.
where V is the core volume.Here we assume that P ext = n i=N k B T kin,i=N is equal to the thermal pressure in the outermost layer.
In our onion model, the one-dimensional non-thermal velocity dispersion, σ v,NT , is a constant value of 0.046 km s −1 .The turbulent energy is calculated by where M is the total core mass.The rotational energy is calculated by where I i and ω i are the moment of inertia and angular velocity of the ith layer, respectively.The magnetic energy is calculated in the cgs units by )

Figure 1 .
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 Sec. 2.1).(b) Reprocessed SCUBA-2 map from the skyloop routine with a PCA model (see Sec. 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 Sec. 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 −1 for panels (a), (b), and (d )) for straightforward comparison.The color scale of panel (c) is [−3σ, 3σ], where σ = 3.1 mJy beam −1 is measured in the central 3 ′ field.
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 Fig. 1b.The blue dashed curve shows the profile of the POL-2 Stokes I map in Fig. 1c.The blue solid curve and filled area show the profile of the POL-2 Stokes I map in Fig. 1d.

Figure 3 .
Figure 3. L 1512 JCMT 850 µm Stokes I , Q, and U maps.(a) SCUBA-2 4 ′′ -sampled map (a zoomed-in view of Fig. 1b) used as the input total intensity model (see Sec. 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 SNR (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 −1 in panels (a)-(d ), respectively.

Figure 4 .
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/σp ≥ 3, and in blue for 3 > p/σp ≥ 2. The black dashed line indicates the nominal distance to L 1512 of 140 pc.

Figure 5 .
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 grey vector as it is a foreground star, identified using Gaia data in Fig. 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 Fig. 1b 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.

Figure 5
Figure5shows 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 sen-

Figure 6 .
Figure 6.Probability density function (p.d.f.) 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 Fig. 5c) resulting from the 2D Gaussian fitting to the 850 µm intensity is plotted as a black dashed line in panel (b).

Figure 7 .
Figure 7. Position angles of the B -field orientations and their uncertainties.(a) B -field position angle map (θ) inferred from the polarization angles (PA) rotated by 90 • .(b) Mean B -field angle map (⟨θ⟩ 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 (∆θ = θ − ⟨θ⟩ local ) and blue vectors of these residual angles are overlaid.(d) Uncertainty map (σPA) of B -field angles.(e) Map of the local largest PA uncertainty (σ PA,local ), on which each pixel value represents the largest PA uncertainty (σ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 σ PA,local (see Sec. 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 σPA, and 4.7 • , 7.6 • , and 9.3 • for σ 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.
Figure 8. Residual angle (∆θ) cumulative function of the maximum allowed B -field angle uncertainty (σPA,max), in that the standard deviation of ∆θ with the uncertainty (σ PA,local ) equal to or less than σPA,max is evaluated.The standard deviation estimate, σ(∆θ), is shown as black dots, while the cumulative number of ∆θ is shown as red dots.The intrinsic B -field angular dispersion (δθ) of 8.2 • ±1.4 • (horizontal solid and dashed blue lines) is determined by computing the mean and standard deviation of the σ(∆θ) evaluated from 6 • ≤ σPA,max ≤ 9.3 • (vertical dashed black lines).

Figure 9 .
Figure 9. Density maps for the L 1512 core with R edge = 108 ′′ .(a) NH 2 map.(b) Line-of-sight-averaged nH 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.

Figure 10 .
Figure10.JCMT 850 µm non-debiased polarization fraction (p ′ ) is plotted against the total intensity (Stokes I ).The best-fit Ricean-mean model (Equation14) and the grain disalignment model (Equation15) are shown as solid red and dashed black lines, respectively.The red-and green-colored regions are the 68% and 95% confidence intervals, respectively.

Table 1 .
Estimated properties in the DCF analysis from the 850 µm polarimetry.

Table 2 .
Energy budgets in L 1512.Note-Quantities shown are the gravitational energy (Egrav), the thermal energy (E thermal ), the surface kinetic term due to external pressure (E surface ), the turbulent energy (E turb ), the rotational energy (Erot), the total kinetic energy (T = E thermal − E surface + E turb + Erot), the virial energy (Evir = 2T + Egrav), and the plane-of-sky component magnetic energy (Emag,pos).