# An early transition to magnetic supercriticality in star formation

### Data reduction

The FAST Zeeman observations towards the HINSA column density peak in L1544 were carried out on five days between August and November 2019 with a total integration time of 7.6 h. The HINSA spectra were obtained with the central beam of the L-band 19-beam receiver^{31}. The central beam has an average system temperature of 24 K, a main beam efficiency of 0.63 and a main beam diameter at the half-power point of 2.9′ with a pointing accuracy of 7.9″. The 19-beam receiver had orthogonal linear polarization feeds followed by a temperature-stabilized noise injection system and low noise amplifiers to produce the X and Y signals of the two polarization paths. The XX, YY, XY and YX correlations of the signals then were simultaneously recorded using the ROACH backend with 65,536 spectral channels in each polarization. The spectral bandwidth was 32.75 MHz centred at the frequency of the H I 21-cm line for a channel spacing of 500 Hz, and the *V*(*v*) spectrum presented in this work was Hanning smoothed, which produced a spectral resolution of 0.21 km s^{−1}.

The data reduction, including gain and phase calibrations of the two polarization paths, bandpass calibrations of the four correlated spectra and polarization calibrations to generate the Stokes *I*, *Q*, *U* and *V* spectra, was carried out using the IDL RHSTK package written by C. Heiles and T. Robishaw, which is widely used for Arecibo and GBT polarization data. The 19-beam receiver is rotatable from −80° to +80° with respect to the line of equatorial latitude. The polarization calibrations used drifting scans of the continuum source 3C286 at rotation angles of −60°, −30°, 0°, 30° and 60° over 1.5 h surrounding its transit. The details of the polarization calibration procedure are provided in ref. ^{32}. We performed polarization calibrations once per month during the observations. The calibrated polarization of 3C286 of the three epochs were 8.9% ± 0.1%, 8.7% ± 0.2% and 9.0% ± 0.1% for polarization degrees and 30.4° ± 0.3°, 33.8° ± 0.5° and 29.4° ± 0.3° for polarization angles. Considering that the ionosphere can generate a Faraday rotation of 1°–3° in polarization angle at the L band^{33}, our results were consistent with the intrinsic polarization degree of 9.5% and polarization angle of 33° of 3C286 at 1,450 MHz (ref. ^{34}). In addition to the polarization observations of L1544 and 3C286, we observed the circularly polarized OH maser source IRAS02524+2046^{35} to verify that our procedures produced consistent *B*_{los}, including the sign or direction of the magnetic field, as had been obtained previously.

The convolutions of the sidelobes of the Stokes *V* beam with the spatial gradient of the Stokes *I* emission may generate a false ‘S curve’ in the *V* spectrum^{27}. To check the credibility of our Zeeman detections, we measured the Stokes *V* beam of FAST and convolved the beam with the Galactic Arecibo L-band Feed Array (GALFA) Stokes *I* cube^{36} of L1544. The convolved *V* spectrum showed a profile with a shape similar to the *I* spectrum and a strength less than 0.03% of the *I* spectrum, different from the ‘S curve’ patterns in the observed *V* spectrum. Meanwhile, the 19-beam receiver was rotated to −45°, 0° and 45° in the three epochs of the L1544 observations, and all of the three epochs showed ‘S curve’ patterns in the *V* spectra, indicating that our Zeeman results were true detections.

Although the data of the 19 beams of the FAST L-band receiver were simultaneously taken in our observations, only the polarization of the central beam was commissioned at the time of writing. The results represented in this work were made with only the central beam pointing towards the HINSA column density peak in Fig. 1. The Zeeman results of the 18 off-central beams will be published in the future.

### Multiple Gaussians and radiative transfer fitting to *I*(*v*) and *V*(*v*)

We adopt the least-squares fits of multiple Gaussians with radiative transfer^{26} to decompose the *I*(*v*) into the HINSA, CNM and WNM components. The expected profile of *I*(*v*) consists of multiple CNM components providing opacity and also brightness temperature and a WNM component providing only brightness temperature:

$$I(v)={I}_{{rm{CNM}}}(v)+{I}_{{rm{WNM}}}(v).$$

(1)

The *I*_{CNM}(*v*) is an assembly of *N* CNM components

$${I}_{{rm{CNM}}}(v)=mathop{sum }limits_{n=1}^{N}{I}_{{rm{peak}},n}(1-{{rm{e}}}^{-{tau }_{n}(v)}){{rm{e}}}^{-({sum }_{m=0}^{M}{tau }_{m}(v)+{tau }_{0})},$$

(2)

where the subscript *m* with its associated optical depth profile *τ*_{m}(*v*) represents each of the *M* CNM clouds that lie in front of cloud *n*. The optical depth of the *i*th component is

$${tau }_{i}(v)={tau }_{i}{{rm{e}}}^{-{[(v-{v}_{0,i})/{sigma }_{v,i}]}^{2}}$$

(3)

in which *τ*_{0} represents the HINSA providing only opacity and no brightness temperature. For the WNM in the background

$${I}_{{rm{WNM}}}(v)={I}_{{rm{peak}},{rm{WNM}}}{{rm{e}}}^{-{[(v-{v}_{0,{rm{WNM}}})/{sigma }_{v,{rm{WNM}}}]}^{2}}{{rm{e}}}^{-{sum }_{i=0}^{N}{tau }_{i}(v)}.$$

(4)

The fitting of *I*(*v*) thus yields values for the intrinsic peak Stokes *I* emission (*I*_{peak}), *τ*, *v*_{0} and the Gaussian dispersion (*σ*_{v}) of the components.

We consider the radiative transfer of *V*(*v*) in terms of right circular polarization (RCP) and left circular polarization (LCP). The Zeeman effect states that with the existence of *B*_{los}, the frequency of RCP shifts from its original frequency *ν*_{0} to *ν*_{0} + *ν*_{z} and the frequency of LCP shifts to *ν*_{0} − *ν*_{z} with *ν*_{z} = (*Z*/2) × *B*_{los}, where *Z* is the Zeeman splitting factor (2.8 Hz μG^{−1} for the H I 21-cm line). As the RCP and LCP are orthogonal components of radiation, the radiative transfer processes of RCP and LCP are independent of each other. For RCP, equation (1) becomes

$${T}_{{rm{RCP}}}={T}_{{rm{RCP}},{rm{CNM}}}(v,{tau }_{{rm{RCP}},i})+{T}_{{rm{RCP}},{rm{WNM}}}(v,{tau }_{{rm{RCP}},i}),$$

(5)

where for the *i*th component, ({T}_{{rm{RCP}},i}={I}_{{rm{peak}},i}/2), *τ*_{RCP,i} is the optical depth in the RCP radiation to substitute the *τ*_{i} in equation (3) with ({tau }_{{rm{RCP}},i}={tau }_{i}({nu }_{0}+{nu }_{{rm{z}},i})) for *B*_{los,i} of the component, and the parameters of *v*_{0} and *σ*_{v} remain the same. Similarly, for LCP

$${T}_{{rm{LCP}}}={T}_{{rm{LCP}},{rm{CNM}}}(v,{tau }_{{rm{LCP}},i})+{T}_{{rm{LCP}},{rm{WNM}}}(v,{tau }_{{rm{LCP}},i})$$

(6)

with ({T}_{{rm{LCP}},i}={I}_{{rm{peak}},i}/2) and ({tau }_{{rm{LCP}},i}={tau }_{i}({nu }_{0}-{nu }_{{rm{z}},i})). The fitting of (V(v)={T}_{{rm{RCP}}}-)({T}_{{rm{LCP}}}+cI(v)), which includes a *c* term accounting for leakage of *I*(*v*) into *V*(*v*), thus yields values for *B*_{los} of the components. In Extended Data Table 1, we list the parameters of the components obtained from least-squares fits to *I*(*v*) and *V*(*v*). The leakage of our HINSA Zeeman observations is *c* = 0.034%.

### HINSA and H_{2} column density maps

L1544 is a low-mass prestellar core in the Taurus molecular cloud complex at a distance of about 140 pc. The core has a size of about 0.1 pc (ref. ^{37}), presumably formed out of a parsec-long elongated molecular ridge^{7}, which, for simplicity, we refer as the molecular envelope. We show the HINSA and H_{2} column density maps of L1544 in Fig. 1a, and we use the H_{2} column density map to calculate the *N*_{H2} of the envelope and core at the beams of FAST and Arecibo observations in Table 1. The HINSA column density map is a revision of Fig. 8 in ref. ^{4}. To derive the H_{2} column density map, we retrieved the level-2.5 processed, archival Herschel images that were taken at 250 μm/350 μm/500 μm using the SPIRE instrument^{38} (observation ID 1342204842). We smoothed the Herschel images to a common angular resolution of the 36″ beam at 500 μm and regridded the images to the same pixel size of 6″. We performed least-squares fits of the 250 μm/350 μm/500 μm spectral energy distributions weighted by the squares of the measured noise levels to derive the pixel-to-pixel distributions of dust temperature *T*_{d} and dust optical depth *τ*_{ν} using ({S}_{nu }={{Omega }}_{m}{B}_{mu }({T}_{{rm{d}}})(1-{{rm{e}}}^{-{tau }_{nu }})), where *S*_{ν} is the flux density at frequency *ν*, *Ω*_{m} is the is the solid angle of the pixel, *B*_{μ}(*T*_{d}) is the Planck function at *T*_{d} and ({tau }_{nu }={tau }_{230}{(nu ({rm{GHz}})/230)}^{beta }) with a dust opacity index *β* of 1.8. Next, we obtained the H_{2} column density with ({N}_{{rm{H2}}}=g{tau }_{230}/({kappa }_{230},{mu }_{{rm{m}}}{m}_{{rm{H}}})), where *g* = 100 is the gas-to-dust mass ratio, *κ*_{230} = 0.09 cm^{2} g^{−2} (ref. ^{39}) is the dust opacity at 230 GHz, *μ*_{m} = 2.8 is the mean molecular weight and *m*_{H} is the atomic mass of hydrogen. To estimate the uncertainties in the H_{2} column density, we used a Monte Carlo technique. For each pixel, we created artificial 250 μm/350 μm/500 μm flux densities by adding the original flux densities with normal-distributed errors taking account the uncertainty in the measured flux and a 10% correlation for the calibration uncertainty in SPIRE^{40}. We then estimated the uncertainty in each pixel with 1,000 fittings of the H_{2} column density. The *N*_{H2} and its uncertainty in Table 1 were obtained from the convolutions of the H_{2} column density map and uncertainty map with the FAST and Arecibo beams.

Note that the equivalent H_{2} column density *N*_{H2} of the CNM1 are derived from H I data towards 3C132 and 3C133, a method different to the *N*_{H2} of the L1544 envelope and core that are derived from dust emission. Therefore, in addition to the statistical errors listed in Table 1, there is a systematic difference between the *N*_{H2} derived from the two methods. Considering that the regime traced by dust emission can be different from those traced by HINSA or OH, which is particularly noticeable from the different spatial extents of dust and HINSA in Fig. 1a, we expect that the systematic difference could be as large as a factor of a few. As the values of *λ* between CNM1 and L1544 are different at least by a factor of 13, the systematic difference between the two methods should not change the qualitative conclusion of this work.

In Fig. 1a, the peak of HINSA column density appears to be shifted from the centre of L1544 by 0.15 pc and the 70% and 90% contours of the peak HINSA column density do not enclose L1544. We note that such offset has also been seen for other dense gas tracers in prestellar cores^{41}. The core geometry may not be as simple as envisioned in idealized theories, where the dense core sits near the centre of a lower-density molecular envelope. In particular, the L1544 core appears to sit near one end of an elongated molecular (and dust) ridge, which roughly coincides with the region traced by the HINSA. Such an offset can result from complexities in chemistry and formation history, but does not affect the main science result of this work, namely, the HINSA Zeeman probes the magnetic fields of the current molecular ridge that is the progenitor of the dense core.

### Maximum likelihood

We adopt the analysis of maximum likelihood^{18} to study the uniformity of magnetic fields in the envelope of L1544. Assuming that the true *B*_{los} follows a Gaussian distribution with mean *B*_{0} and intrinsic spread *σ*_{0}, the likelihood *l*_{j} for a single observation in a set of *N* measurements (*j* = 1, …, *N*) to measure *B*_{j} with Gaussian error *σ*_{j} is proportional to the convolution of the probability (exp [-{(B-{B}_{0})}^{2}/2{sigma }_{0}^{2}]/sqrt{2{rm{pi }}}{sigma }_{0}) for the magnetic field to have a true value of *B* with the probability (exp [-{(B-{B}_{j})}^{2}/2{sigma }_{j}^{2}]/sqrt{2{rm{pi }}}{sigma }_{j}) of observing a value *B*_{j} of the field. Therefore, *l*_{j} is the integral over all possible true values of the magnetic field

$${l}_{j}={int }_{-infty }^{infty }{rm{d}}Bfrac{exp [-{(B-{B}_{j})}^{2}/2{sigma }_{j}^{2}]}{sqrt{2{rm{pi }}}{sigma }_{j}}frac{exp [-{(B-{B}_{0})}^{2}/2{sigma }_{0}^{2}]}{sqrt{2{rm{pi }}}{sigma }_{0}}.$$

(7)

Although the overall likelihood ( {mathcal L} ) for a set of observations is the product of individual likelihoods of the observations (( {mathcal L} ={prod }_{J=1}^{N}{l}_{j})), the *B*_{0} and *σ*_{0} can be estimated by maximizing the likelihood ( {mathcal L} ). After performing the integration in equation (7) and some algebraic manipulations

$$ {mathcal L} ({B}_{0},{sigma }_{0})=left(mathop{prod }limits_{J=1}^{N}frac{1}{sqrt{{sigma }_{0}^{2}+{sigma }_{j}^{2}}}right)exp left[-frac{1}{2}mathop{sum }limits_{j=1}^{N}frac{{({B}_{j}-{B}_{0})}^{2}}{{sigma }_{0}^{2}+{sigma }_{j}^{2}}right]$$

(8)

Extended Data Fig. 1 shows the distribution of ( {mathcal L} ) as functions of *B*_{0} and *σ*_{0} and the probability distributions of *B*_{0} and *σ*_{0} by integrating ( {mathcal L} ) along the *B*_{0} axis and the *σ*_{0} axis, respectively. The probability distribution of *B*_{0} is similar to a normal distribution with a mean value of +4.1 μG and a standard deviation of 1.6 μG. The probability distribution of *σ*_{0} is highly asymmetric as the values of *σ*_{0} cannot be negative. The first, second and third quartiles of the *σ*_{0} distribution are 0.6 μG, 1.2 μG and 2.4 μG. We therefore suspect that the Zeeman measurements in the L1544 envelope can be explained by a magnetic field with *B*_{0} = +4.1 ± 1.6 μG and ({s}_{0}={1.2}_{-0.6}^{+1.2},{rm{mu }}{rm{G}}).

### Inclination angle of magnetic field

Given the uniformity of magnetic fields in the envelope of L1544 and CNM1 is well constrained by the maximum likelihood analysis, the coherent *B*_{los} suggests that the inclination angles of magnetic fields in the CNM1 and L1544 envelope are likely to be similar, or a special geometry of magnetic field structure across multi-scales and multi-phases of the interstellar medium is needed. In contrast, the *B*_{los} of the L1544 envelope and core differ by a factor of 2.6. There are two physical explanations for the 2.6-times difference between the HINSA and OH Zeeman measurements. First, the OH measurement probably samples a denser gas than the HINSA measurement, as the column density along the OH sightline is twice that along the HINSA sightline (Table 1). As the magnetic field strength in molecular clouds tends to increase with number density^{15}, the stronger field is naturally expected in the denser core. Alternatively, the inclination angle of the L1544 core magnetic field could differ substantially from that of the coherent field. As we cannot rule out the second possibility, an assumption of similar inclination angles of the magnetic fields in the L1544 envelope and core thus is required to calculate the relative values of *λ*.

We note that dust polarization observation may give some clues as dust polarization traces the position angle of the plane-of-sky component of the magnetic field. The near-infrared polarization observations of L1544^{42} indicate that the mean position angle of the magnetic field towards the core location of the Arecibo beam is 29.0°–36.9°, and the mean position angles of magnetic fields towards the four envelope locations of the GBT beams are 30.5°–55.8°. The difference in the position angles between the core and envelope thus may be about 10°–20°.

We perform Monte Carlo simulations^{43} to study whether the 2.6-times difference between the *B*_{los} of the L1544 envelope and core can be explained by different inclination angles. The simulations randomly generate two unit vectors in three dimensions, and then measure the difference between the inclination angles and the difference between the position angles of the two vectors. The probability of the cases that the line-of-sight length of one vector is 2.6-times larger than that of the other is roughly 0.19. For those cases, the mean difference between the inclination angles and the mean difference between the position angles of the two vectors are about 38° and 45°, respectively. As the probability of 0.19 is small and the difference of about 45° between the simulated position angles is about a factor of two to four times larger than the difference of about 10°–20° between the observed position angles, it is less likely that the 2.6-times difference between the *B*_{los} of the L1544 envelope and core can be solely explained by different inclination angles.

### CCS Zeeman measurements

Ref. ^{20} reported a CCS Zeeman detection of 117 ± 21 μG in the dense core of TMC-1 that has an estimated H_{2} column density of 3 × 10^{22} cm^{−2}, which is four-times higher than that probed by OH Zeeman measurements in L1544 and nearly one order of magnitude higher than that probed by our HINSA measurements. It appears to provide further support to the evolutionary scenario suggested by our HINSA measurements: namely, once the gas loses its magnetic support during the transition from the CNM to the molecular envelope (or ridge) and becomes magnetically supercritical, there is no longer any need to lose magnetic flux further (relative to the mass) for a piece of the envelope/ridge to condense into a (magnetically supercritical) core (for example, the L1544 core probed by OH) and for the core to evolve further by increasing its column density (for example, the TMC-1 core probed by CCS).

Technically, we note that one potential source of significant uncertainty in frequency shift, namely the uncertainty of beam squint, was not included in the CCS result, which may affect the level of significance. In comparison, the HINSA measurement is robust with a greater than 10*σ* significance with the beam squint and velocity gradient being taken into account by convolving the FAST Stokes *V* beam with the Stokes *I* cube of L1544 (see the third paragraph in the data reduction section in Methods).

https://www.nature.com/articles/s41586-021-04159-x