The magnetic properties of the protostellar core IRAS 15398-3359

Context. Magnetic fields can affect significantly the star formation process. The theory of the magnetically-driven collapse in a uniform field predicts that initially the contraction happens along the field lines. When the gravitational pull grows strong enough, the magnetic field lines pinch inwards, giving rise to a characteristic hourglass shape. Aims. We investigate the magnetic field structure of a young Class 0 object, IRAS 15398-3359, embedded in the Lupus I cloud. Previous observations at large scales suggest that this source evolved in an highly magnetised environment. This object thus appears an ideal candidate to study the magnetically driven core collapse in the low-mass regime. Methods. We have performed polarisation observations of IRAS 15398-3359 at 214$\mu$m using the SOFIA/HAWC+ instrument, thus tracing the linearly polarised thermal emission of cold dust. Results. Our data unveil a significant bend of the magnetic field lines due to the gravitational pull. The magnetic field appears ordered and aligned with the large-scale B-field of the cloud and with the outflow direction. We estimate a magnetic field strength of $B= 78 \mu$G, expected to be accurate within a factor of two. The measured mass-to-flux parameter is $\lambda= 0.95$, indicating that the core is in a transcritical regime.


Introduction
Magnetic fields (B) are expected to play an important role in the star formation process, for instance providing a source of nonthermal pressure against the gravitational pull (see e.g. McKee & Ostriker 2007). High-density molecular clouds are threaded by magnetic fields aligned perpendicularly to the clouds' main axes (Planck Collaboration et al. 2016), which is the expected configuration when B-fields regulate the star formation process more effectively than turbulence (Mouschovias & Spitzer 1976;Nakamura & Li 2008;Li et al. 2013). At smaller scales (≈ 0.1 − 1 pc), the collapse of spherical cores in uniform fields have been widely investigated (Mouschovias 1991;Galli & Shu 1993;Allen et al. 2003). Theory predicts the formation of a flattened structure (pseudodisc) because the collapse happens preferentially along the field lines. Given that the interstellar gas is often slightly ionised, matter is expected to be coupled with the magnetic field lines at envelope scales. Therefore, magnetic lines pinch inwards due to the gravitational pull, exhibiting a characteristic hourglass shape. In low-mass star forming regions this feature has been detected only in 30% of Young Stellar Objects in polarisation (9 sources out of 32, Hull & Zhang 2019), suggesting that this is not a universal picture. Furthermore, out of these nine detections only two show a clear hourglass shape, namely IRAS 4A (Girart et al. 2006) and L1448 (Kwon et al. 2018). Polarisation observations in the far-infrared (FIR)/submillimeter wavelength are an effective way to investigate the magnetic properties of clouds and cores at intermediate/high visual extinctions (A V 10 mag). In fact, asymmetric dust grains illuminated by an external radiation field are expected to develop a magnetic moment and thus to align with their minor axes par-allel to the magnetic field direction (radiative torque alignment, or RAT, Lazarian 2007). As a result, the background starlight at near-infrared/optical frequency is absorbed and linearly polarised in the B-field direction. Similarly, since the dust grains emit preferentially in the direction of their major axes, dust thermal emission is linearly polarised perpendicularly to the local B-field. Cold dust (T ≈ 10 − 30 K) emission happens typically at FIR/sub-mm wavelengths.
IRAS 15398-3359 (hereafter IRAS15398) is a low-mass Class 0 protostar (Andre et al. 2000) located in the Lupus I molecular cloud, at a distance of 156 pc (Dzib et al. 2018). The protostellar mass, estimated from modelling of the envelope rotation, is subsolar (M * < 0.1 M , Oya et al. 2014;Yen et al. 2017). Lupus I is the least evolved cloud of the Lupus complex (Rygl et al. 2013), and optical polarisation observations have shown that it is threaded by a very ordered magnetic field, perpendicular to its whole filamentary extension (Franco & Alves 2015). IRAS15398 has a luminosity of 1.8 L and an envelope mass of 1.2 M . It powers a bipolar outflow detected in CO lines both with single dish and interferometric observations (Tachihara et al. 1996;Bjerkeli et al. 2016).
Using the SOFIA telescope, we have observed the polarised dust thermal emission arising from the core and the filamentary structure in which IRAS15398 is embedded. For the first time, we can compare the large scale field with the small scale field in a Class 0 source In this paper, we report our findings, which unveil an ordered magnetic field, with hints of pinching of the field lines due to the gravitational pull. Our results suggest that IRAS15398 evolved in an highly magnetised environment, and experienced a magnetically-driven collapse.

Polarimetric data
The polarimetry observations were carried out with the SOFIA telescope, using the HAWC+ instrument (Harper et al. 2018) in summer 2018. We used the band E, with a nominal wavelength of 214 µm (1.4 THz). At this frequency, the SOFIA beam full-width-half-maximum is FWHM ≈ 19 (corresponding to ≈ 0.01 pc at the source's distance), and the field of view is 5.0 × 6.3 . The data were processed with the standard pipeline (HAWC DRP, version 1.3.0), and we used the flux-calibrated Level 4 data products (already flux calibrated and mosaiced). The total integration time was ≈ 1h15min, and the sensitivity achieved in the Stokes I is 65 mJy/beam. In order to increase the signal-to-noise ratio (SNR) of the data, and to have significant detections over a large fraction of the source, we smoothed the Stokes parameter I, Q, and U to 42 resolution (≈ 0.03 pc), before producing maps of polarization intensity and position angle (see Sec. 3). The final mean sensitivity is rms = 48 mJy/beam (Stokes I), 77 mJy/beam (Stokes Q), and 77 mJy/beam (Stokes U). Since we found no indication of spatially correlated noise, the smoothing technique as a tool to increase SNR is justified.

H 2 column density maps
To obtain the gas column density N(H 2 ) and the visual extinction A V maps of the source, we use the archive data of the Gould Belt Survey, performed with the Herschel space telescope (André et al. 2010). They provide directly the H 2 column density map of Lupus I, obtained with a spectral fitting of the dust emission at 5 wavelengths (70, 160, 250, 350, and 500µm). This map has a resolution of ≈ 38 , which allows a fair comparison with the SOFIA smoothed data. We derived the visual extinction map using the standard relation A V /mag = 1.06×10 −21 ·N(H 2 )/cm −2 , (Bohlin et al. 1978).

Results
Appendix A reports the maps of the three Stokes parameters I, Q, and U. From these, we derive the polarised flux (I P ), the polarised fraction (P) and the polarisation angle (PA), following the standard equations: In order to debias the polarised intensity, in Eq. (1) we have removed the contribution of the flux uncertainty (σ I P ). This procedure is necessary since I P is biased to positive values, while the Stokes U and Q parameters can be positive or negative (Vaillancourt 2006). To avoid oversampling, after smoothing the data we re-gridded the maps to a pixel size of 9.1 (≈ 4 pixels per smoothed beam). Table 1 summarises the peak and mean values for the derived parameters. Figure 1 shows the Stokes I emission (in colorscale), with polarization vectors already rotated by 90 • (in black) to indicate the plane-of-sky component of the magnetic field. We mask data points with a signal-to-noise ratio in polarisation fraction lower than 3. Notes. (a) The mean values are computed over positions satisfying σ P /P > 3.0 and P < 30%

Analysis and discussion
4.1. The magnetic field direction Figure 1 shows the magnetic field direction traced by the SOFIA data. The field at sub-parsec scale is aligned with the large (parsec) scale one observed in starlight optical polarisation. The histograms of the two position angle distributions are presented in Figure 2. The mean values of the two distributions differ by ≈ 5 • , less than the mean uncertainty on the SOFIA position angles ( PA = 7 • ). Since the two datasets are sensitive to two different regimes, the optical data tracing large cloud scale, whilst the SOFIA ones being sensitive to the core scales, this means that the uniform magnetic field of the cloud has been inherited by the core, and that the gravitational collapse was magnetically-driven.
In Fig. 1 we also show that the outflow direction (PA = 35 • , Bjerkeli et al. 2016) lies almost parallel to the magnetic field, which has a mean direction of PA = (45 ± 7) • . This behaviour is consistent with the prediction of the theory of the magnetically driven collapse, according to which a strong magnetic field can efficiently remove the excess of angular momentum from rotating cores, thus aligning the rotation axis with the magnetic one (see e.g. Li et al. 2014, and references therein).
Another prediction of the theory is the hourglass shape of the magnetic field lines. In IRAS15398 this is not clearly visible, even though we do detect a bending of the field lines which may hint to a partial hourglass shape. In fact, on the side of the core facing South-East, the field lines pinch inwards toward the centre of the object. The mean position angle North-East of the source is PA = (55 ± 7) • , whilst PA = (36 ± 7) • is derived South-West of it. The change in angle is thus 19 • , significantly larger than the uncertainty. We do not detect the North-West side of the hourglass, which can be due to two main reasons. First of all, on that side we lack the sensitivity to detect a bend in the field lines, since the polarised flux is less prominent. This can be due to the presence of the filamentary structure that extends towards West. This filament provides more shielding from the external radiation field with respect to the other side of the core, which is more exposed. Since this radiation is responsible of the grain alignment, according to the RAT theory, grains are less efficiently aligned on the North-West side of IRAS15398, which results in a lower polarised flux. Moreover, the presence of the filament can affect the detected morphology of the field lines in a second way. In fact, the presence of the extended dust emission in the West direction can have played a role in disturbing the spherical collapse of the core. This scenario is supported by the few detections along the filament, which reveal a significant twist in the B direction (PA ≈ 69 • ). This suggests that in this portion of the source the magnetic field morphology has been perturbed, likely by accretion motions towards the central object.  The SOFIA data present a three peaked distribution, with peaks being at ≈ 40 • , ≈ 55 • , and ≈ 77 • . These correspond to vectors coming from the two halves of the presumed hourglass shape and the filament, respectively. Figure 1 shows that the emission is less polarised as the column density gets higher, towards the centre of the core. This effect has been widely observed in the literature (Matthews et al. 2009;Alves et al. 2014;Jones et al. 2015Jones et al. , 2016. In Figure 3 we show the scatterplot of the polarisation fraction as a function of the visual extinction. We also report the polarisation efficiency for the optical data, which is defined as the polarised fraction normalised by the visual extinction (see the Introduction in Santos et al. 2017). In this way, one correctly takes into account the attenuation of the background stellar light. We use the data from Franco & Alves (2015), limiting to the Galactic latitude range 16.2 • < b < 17.0 • (the "Middle" region in Franco & Alves 2015).

Depolarisation at high column densities
There are a number of possible explanations for the depolarisation at high column densities. In general, if the alignment is dominated by radiative torque due to the interstellar radiation field, at higher visual extinction less and less radiation penetrates the core and dust grains are less aligned. Moreover, the alignment efficiency is highly sensitive to the dust grain population (Bethell et al. 2007;Pelkonen et al. 2009;Brauer et al. 2016). In particular, Brauer et al. (2016) showed that polarization can decrease by up to 10% when grains have sub-mm/mm sizes. Grain growth is expected in Class 0 objects, as shown by Jørgensen et al. (2007); Kwon et al. (2009);Chiang et al. (2012). Finally, geometrical smearing is possible. Close to the central object, a highly perturbed magnetic field is likely to be present due to the dominant effect of gravity at the small scales which are unresolved by the SOFIA data. Higher angular resolution observations are needed to disentangle these scenarios.
We fit a linear relation in the log-log space for the polarization efficiency vs visual extinction in each data-set, deriving the slope of the relation P eff ∝ A −α V (Figure 3). For the optical Article number, page 3 of 7  data, tracing the cloud scales, we find α = 0.57 ± 0.07, which is consistent with modelling of RAT alignment in molecular clouds (see e.g. Whittet et al. 2008, who also report the measurements in Taurus and Ophiucus clouds). For the FIR data, we find a steeper slope (α = 1.21 ± 0.12), implying a strong depolarisation for 10 A V (mag) 50. Values of α ≈ 1 are often found in literature for prestellar and young Class 0 objects (Alves et al. 2014;Jones et al. 2015Jones et al. , 2016. The different slopes exhibited by the optical and the FIR data suggest that the two datasets are probing two different regimes in the grain alignment. This in turn can be related to two different grain populations, with the larger grains in the core less efficiently aligned to the B-field than the ones in the cloud (see e.g. Sect. 3 in Andersson et al. 2015), or it can be due to a less efficient RAT, with the grains traced by SOFIA more shielded from the interstellar radiation field with respect to those traced by the optical data. Another possibility is a more entangled magnetic field lines, as mentioned before. According to recent results, a slope of α ≈ 1 can be due to biases introduced by the SNR cut (Wang et al. 2019). To exclude this possibility, we fit the linear relation lowering the threshold to SNR > 0.5, obtaining α = 1.26 ± 0.09, in agreement the previous result, thus suggesting that no significant bias is introduced by the data selection.

Angular dispersion function and field strength
The angular dispersion function (ADF) method, introduced by It allows to derive quantitative information on the B-field without depending on assumption on the field morphology. The method assumes that the magnetic field is composed by an ordered component (B 0 ), slowly changing spatially, and a turbulent or random component (B t ), which is characterised by a coherent length scale δ. We now consider the autocorrelation function of the position angles ∆Φ, i.e. the difference in angle of every pair of vectors separated by the distance l: Hildebrand et al. (2009) showed that on scales d such that δ < d << ∆, where ∆ is the physical scale of the analysed source, the structure function has the form: where σ M is the uncertainty on the measurements. The parameter a is linked to the large-scale variations of the B-field, while b is related to the ratio of the turbulent and uniform magnetic field Houde et al. (2009) further expanded the analysis, in particular including the effects introduced by the telescope beam size. In case of a Gaussian beam with standard deviation W, they found the following expression: where ∆ is the cloud effective thickness, and m is a parameter related to the large-scale structure of the magnetic field, which does not involve turbulence. In order to fit Eq. (6) to derive the B t /B 0 parameter, we compute ∆Φ 2 for all the available pairs of points with 0 < l < 180 , divided in 9 bins 20 wide. The choice of binning is determined by the angular resolution of the observations (we do not want to oversample the beam size) and by the map size (each bin should contain a fairly constant number of points). We compute the uncertainty σ M in each bin propagating the uncertainties on the position angles, assumed to be uncorrelated. They are usually within 3-5%. The full-with-half-maximum of SOFIA smoothed beam (42 ) corresponds to a standard deviation of W = 18 = 13 mpc.
Due to the limited number of data-points available for the fit, we fix the turbulence correlation scale on δ = 20 mpc. This value is in good agreement with previous estimation of this quantity in other star forming regions: for instance Houde et al. (2009) found δ = 16 mpc in OMC-1, Frau et al. (2014) derived δ = 13−33 mpc in the high-mass star forming region NGC 7538, and Coudé et al. (2019) reported δ = 7 mpc in Barnard 1. Furthermore, we assume the cloud thickness ∆ to be equal to the source effective radius (r eff ), defined as the radius of a circular region of identical surface area. To compute the latter, we considered only positions with Stokes I > 0.05 Jy/pix. This threshold roughly corresponds to the lowest contour in N(H 2 ) in Figure  1 (1 × 10 22 cm −2 ), and it is the first closed contour that comprises both the central core and the filamentary structure towards West. We obtain r eff = 0.1 pc, a result which is also in agreement with the typical size of cores and filamentary structures in molecular clouds (see e.g. Arzoumanian et al. 2011;André et al. 2014). We explore later in the section the effects of these assumptions on the results.
In Figure 4 we show the resulting data points and their uncertainties. We also show the best-fit solution found for Eq. (6). The obtained best-fit parameters are reported in Table 2. We emphasise that while figures and tables are expressed in deg/arcsec units, the actual fit is performed in rad/pc units.
In order to estimate the strength of the magnetic field, we apply a modified version of the method proposed by Chandrasekhar & Fermi (1953) (CF). The CF approach assumes the equipartition of kinetic and perturbed magnetic field energy to link the magnetic field strength on the plane of sky (B pos ) to the velocity dispersion of the gas (σ V ) and to the dispersion of the polarization angle (δφ), which in turn corresponds to the turbulent ratio (δφ ≡ B t /B 0 ). Further development of the theory led to the equation: which is valid in the assumption of small polarization angle (δφ < 25 • , Ostriker et al. 2001;Padoan et al. 2001;Crutcher et al. 2004). In Eq. (7), µ is the gas mean molecular weight per hydrogen molecule, assumed to be µ = 2.8 (Kauffmann et al. 2008), m H the hydrogen mass (m H = 1.67 10 −24 g), and n H 2 is the gas volume density.
To compute B pos , we need to determine the quantities n H 2 and σ V . For the gas volume density, we consider the H 2 column density map, integrated over the region used to derive the effective radius, obtaining the total number of H 2 molecules H tot in the considered region. Assuming a uniform distribution of the gas, we can thus derive: For the gas velocity dispersion, we use the results of Benedettini et al. (2012), who observed several molecular tracers in Lupus I with the MOPRA telescope, which at 3 mm has a beam size similar to the one of our data. At the position of IRAS15398, they report full width at half maximum (FWHM) values in the range 0.31 − 0.72 km s −1 . We adopt the mean of their results (excluding the CS (2-1) line, which is most likely optically thick), FWHM = 0.41 km s −1 , corresponding to a velocity dispersion σ V = 0.17 km s −1 . Inserting these values in Eq. (7), we obtain B pos = 78 µG. Due to the strong assumptions of the method, such as the core uniform density, we expect this value to be accurate within a factor of two. Furthermore, we want to highlight that it is intrinsically a lower limit, since it neglects the line-of-sight component of the field. A fundamental parameter is the mass-to-flux ratio (M/Φ), which is determined by the ratio between the gravitational and the magnetic energy. It provides the dynamical state of the core (equilibrium or collapse). In particular, it is interesting to compute its observed value with respect to the critical value, defined by the parameter λ (Crutcher et al. 2004): where N(H 2 ) is the mean column density of the source, computed in the region over which B is measured. Focusing on the same region used to derive n H 2 , we obtain N(H 2 ) = 9.8 × 10 21 cm −2 , and therefore λ = 0.95. Our derived value of λ suggests that the core is in transitional state between being subcritical and supercritical. Our analysis of the mass-to-flux ratio considers only the gas mass of the source, and neglects the protostellar mass. However, the total mass within the contour Stokes I > 0.05 Jy/pix is ≈ 6 M , significantly larger than the the protostellar mass (the upper limit is M * < 0.1M , Yen et al. 2017). This approximation therefore causes an error of less than 2%, well below our uncertainties.

Influence of the assumed parameters
The analysis described above depends on the assumption of several key parameters. We explore in this section how the results are affected by changing the assumed values. The first assumption regards the turbulence correlation scale, δ. The literature results seem to suggest that δ is lower for low-mass starforming regions, with respect to high-mass ones: both Liu et al. (2019) and Coudé et al. (2019) found δ < 10 mpc in low-mass cores. We therefore tested other 3 values (5, 10, 15 mpc), obtaining B pos = 16, 40, 62 µG respectively. The magnetic field strength hence increases with increasing δ. The turbulent-touniform component ratio is lower than one in all cases but the extreme one (δ = 5 mpc), for which is marginally higher than one (B t /B 0 = 1.3). On the other hand, the mass-to-flux ratio for these values is λ > 1.0, pushing the core's state more towards the supercritical regime. A second assumption that is made regards the cloud effective thickness (∆ ), which may not necessary correspond to the real thickness ∆. In particular, according to Houde et al. (2009), ∆ ≤ ∆, where the equality holds only if the large-scale magnetic field is completely uniform. We therefore repeated the analysis using ∆ = r eff /2. The results vary of less than 30%, below our uncertainties.
The magnetic field strength can also be calculated directly with the modified Chandrasekhar-Fermi method (Eq. (7)), computing directly the standard deviation of the position angles (δφ), hence avoiding the assumptions made in the ADF method. In this case, Eq. (7) must be corrected for a numerical factor (Q), in order to take into account that its original version tends to overestimate the magnetic field (Ostriker et al. 2001;Crutcher et al. 2004). A usually adopted value is Q = 0.5. From the SOFIA data, we measure δφ obs = 14 • . This value is however biased up-ward by the observational uncertainties on the position angles ( PA ), and it must be corrected using δφ 2 = δφ 2 obs − 2 PA = (12 • ) 2 (Crutcher et al. 2004). Hence, introducing δφ = 12 • in Eq. (7), we obtain B pos = 47 µG, which is in well agreement with the result of the full ADF analysis, given the uncertainties. We conclude that the results are robust with respect to our assumptions.

Conclusions
We have performed polarimetric observations of the dust thermal emission at 1.4 THz in the class 0 protostar IRAS15398. Our data unveil a very ordered magnetic field. The direction of the magnetic field is consistent within uncertainties with the one found by Franco & Alves (2015) using optical data, which trace large cloud scales. This suggests that during its evolution the core preserved the B-field morphology inherited from the parental cloud, and experienced a magnetically driven collapse.
The magnetic field lines on the South-East side of the core pinch inwards, hinting to the hourglass shape predicted by the theoretical models. However, we lack the sensitivity to clearly identify the hourglass shape. The lack of such a pinching on the other side of the core may be linked to the presence of an extended structure linked to the central core. The polarisation vectors detected in this filament show a significantly different direction of the magnetic field, perhaps hinting that material is being dragged towards the central object and that the infalling gas is perturbing the magnetic field morphology. Further line observations, which allow to infer the velocity field in the source, will help us enlighten this point.
Using the modified CF method we estimate the magnetic field strength on the plane of sky to be B pos = 78 µG, accurate within a factor of two. This value is in the lower end of the range found on similar spatial scales performed in star forming regions for example with the JCMT telescope (see e.g. Crutcher et al. 2004;Kwon et al. 2018;Soam et al. 2018;Liu et al. 2019), where B-strengths of 80 − 5000 µG have been reported.
The mass-to-flux ratio that we derive (λ = 0.95) is close to the transcritical regime (λ = 1). Usually cores are found to be super-critical (see e.g. Troland & Crutcher 2008), in contrast with molecular clouds which are often subcritical: for instance, Franco & Alves (2015) found λ = 0.027−0.057 in Lupus I. IRAS15398 hosts a protostellar object, which indicates that indeed gravitational collapse has happened and that, on some scales, the source must be supercritical. Our observations are likely missing the angular resolution to probe those scales, and therefore we are seeing an intermediate state due to the subscritical surrounding medium. On the other hand, the uniform-toturbulent ratio is smaller than one (B t /B 0 = 0.267), suggesting a strongly magnetised core. This scenario is also supported by the lack of a large Keplerian disk in the source (Yen et al. 2017 found an upper limit to the disk size of 30 AU). In fact, magnetic braking is an efficient way to remove angular momentum from the infalling and rotating material and prevent the formation of large disks (Li et al. 2014).
Overall, our data suggest that IRAS15398 evolved in an highly magnetised environment, and that the ordered magnetic field was preserved from cloud scales down to core scales. Future polarimetric observations with ALMA at high resolution will allow us to investigate wether this uniform magnetic field is preserved down to the disk/envelope scales.