© The Authors 2026

1. Introduction

Star formation occurs in dense cores within molecular clouds. These clouds are magnetised, and the cores, called pre-stellar cores (Ward-Thompson et al. 1994), collapse under gravity, eventually giving birth to protostars. During this early evolution, magnetic fields are thought to play a significant role in several processes, such as accretion, ejection, and interaction with the circumstellar environment (Lada 1987). Through their role in these processes, which significantly impact stellar rotation and evolution, magnetic fields shape the close environment of the protostar (Bouvier et al. 2007; Vaytet et al. 2018).

Regulation of angular momentum is made through magnetospheric accretion and ejection processes during the pre-main sequence (PMS) phase (Bouvier et al. 2014; Lee et al. 2017). In classical T Tauri stars (CTTSs), the presence of large-scale magnetic fields is well established (Donati et al. 2011). These fields, typically of kilogauss strength, are responsible for truncating the inner disk and channelling the accreted material onto the stellar surface (Bouvier et al. 2007, 2014; Hartmann et al. 2016). The magnetic properties of CTTSs (between 0.5 and 1.3 M_⊙) are characterised by a dominant poloidal field, with an average strength between 1 and 2 kilogauss (Zaire et al. 2024).

While magnetism and its impact on the star–disk interaction have been intensely explored during the PMS phase (Romanova et al. 2015), our current knowledge of the magnetic properties and of the star–disk interaction in earlier phases remains limited. During those phases, the protostar experiences intense accretion, often through bursts, and grows in mass until it approaches its final mass (Baraffe et al. 2009; Fischer et al. 2023). Such major internal structural changes may have a strong impact on the stellar magnetic field, which in turn may influence the connection between the star and its accretion disk (Kunitomo et al. 2017). The progenitors of the CTTSs (Class II) are the Class I protostars. These objects are still embedded within their envelopes of gas and dust (Lada 1987). The protostars transiting from Class I to Class II are observed as flat-spectrum (FS) sources (Greene et al. 1994). Due to their dense and obscuring environments, the direct observation of Class I protostars and the detection of photospheric features is challenging (Connelley & Greene 2010; Doppmann et al. 2005), making it difficult to study their surfaces, and therefore their magnetic fields.

For several reasons we expect Class I protostars to be strongly magnetised, as in CTTSs. They have evolved enough to set off convection in their full interior (Palla 1998; Baraffe et al. 2017), and are thus capable of hosting strong magnetic fields generated by dynamo processes. In other types of fully convective stars, such as the M dwarfs on the main sequence, strong dynamo magnetic fields are routinely observed (Donati et al. 2008; Morin et al. 2010; Donati et al. 2023). Signs of magnetospheric accretion, as in the CTTSs, have been detected in Class I sources (Doppmann et al. 2005; Fiorellino et al. 2023), suggesting the presence of strong large-scale magnetic fields at the surface of the embedded star. Powerful flares have been detected in X-ray observations of Class 0 and Class I protostars (Grosso et al. 1997, 2020; Montmerle et al. 2000), similar to those present in the Sun or more evolved stars, but much more energetic (Feigelson et al. 2007). These emissions are clues of intense magnetic activity. These studies therefore suggest that strong magnetic fields may already be in place during the Class I phase and may regulate mass inflow and angular momentum loss early on.

A recent study by Flores et al. (2024) has investigated a sample of 42 Class I and FS protostars, and measured their total surface magnetic field from the Zeeman broadening of spectral lines. They found a median magnetic field strength of 1.9 kG, which is consistent with the value found in Class II sources and with earlier studies of Class I protostars (Johns-Krull et al. 2009), indicating that Class I and FS sources are indeed strongly magnetic. However, this study cannot inform us on the large-scale components of their magnetic fields, which interact with the disk. To understand whether magnetospheric accretion and ejection processes are in place at that age, and are potentially at the origin of the jets observed in almost all Class I sources, we need to characterise their large-scale magnetic field.

Furthermore, dynamo processes are not well understood, especially in fully convective stars. The imprint of the dynamo processes in stellar evolution is still an open question (Moss 2003; Donati et al. 2008). In particular, simulations have shown that the initial conditions at the launch of the dynamo may impact the future properties of dynamo magnetic fields, and may explain the bistability of magnetic fields observed in fully convective M dwarfs (Gastine et al. 2012, 2013; Raynaud et al. 2015). While the magnetic energy produced by dynamo processes is an important characteristic to constrain the models, another factor is the ratio of poloidal to toroidal fields, which is only accessible by measuring the magnetic geometry.

The topology and strength of the large-scale magnetic field can be obtained with the Zeeman Doppler imaging (ZDI) technique (Semel 1989). This technique uses time series of circularly polarised and intensity spectra well sampled over the rotation period of the star using a high spectral resolution spectropolarimeter. By assuming that the magnetic field components can be modelled as spherical harmonics, and by using an inversion technique, ZDI provides us with the maps of the radial, azimuthal, and meridional components of the field (Donati et al. 2006; Kochukhov 2016). This technique has been largely employed to map the magnetic field of various types of stars, including CTTSs (e.g. Donati & Gregory 2009; Zaire et al. 2024). To date, ZDI has not been largely performed in Class I and FS protostars.

Recent advances in instrumentation have opened new possibilities for investigating these elusive protostars. Among these instruments is the SPectropolarimètre InfraRouge (SPIRou), a spectropolarimeter installed on the 3.6-meter Canada-France-Hawaii Telescope (CFHT). SPIRou (Donati et al. 2020) operates in the near-infrared (nIR) domain, which allows us to pierce through their dusty envelope to probe the photosphere of Class I protostars. SPIRou’s high-resolution, high-sensitivity measurements in the YJHK bands are particularly well suited to study Class I protostars, which are largely inaccessible to optical spectropolarimeters.

To investigate the magnetic properties of Class I protostars, and their accretion–ejection mechanisms, we have built the Protostellar Magnetism: Heritage and Evolution¹ (PROMETHEE) project. This project will determine for the first time the large-scale magnetic properties of protostars, and the structure of their magnetospheres, and will build unique MHD models of magnetic protostars from their birth to the T Tauri phase. This paper is the first in a series that describes the first SPIRou analysis of a sample of Class I and FS protostars. These data have also been used to study the accretion and ejection diagnostics available in the SPIRou spectra (e.g. HeI 1083 nm, Paβ or Brγ emission lines, or continuum veiling). These results will be submitted in an associated paper (Sousa et al., in prep.). In the meantime, within the PROMETHEE project, we are developing dynamo models adapted to young fully convective protostars, and interacting with their disk (Guseva et al., in prep.). In order to understand the magnetic and stellar properties at the start of the protostellar phase, and how they can subsequently affect the dynamo processes, we are also developing models of molecular cloud collapses until the development of the second Larson core (Ahmad et al. 2025a,b; Ahmad et al., in prep.).

Before performing a monitoring for each source in our sample of Class I and FS protostars, to map their magnetic fields, we performed a snapshot survey with SPIRou. Our aim was to first detect their Zeeman signature in circularly polarised (Stokes V) spectra. In future work, we will monitor with SPIRou only those with a Zeeman detection to perform ZDI.

In this study, we present the analysis of our snapshot sample. It is structured as follows. We describe the observations, data reduction, and spectral and polarimetric analyses in Sect. 2. Our results are presented in Sect. 3 and discussed in Sect. 4. We summarise our conclusions in Sect. 5.

2. Observations

2.1. Sample selection

Our goal is to characterise the magnetic fields of a sample of Class I and FS objects. To this aim, we need to detect the circular polarisation signal induced by the Zeeman effect in the photospheric lines. To build our sample, we first listed the objects in which photospheric lines had been detected in previous surveys using high-resolution optical or nIR spectrographs (White & Hillenbrand 2004; Doppmann et al. 2005; Evans et al. 2009; Connelley & Greene 2010; Rebull et al. 2010; Viana Almeida et al. 2012; Dunham et al. 2015). Then we first kept the sources north to a declination of –40°, to ensure that they were observable from the CFHT. In a second step, we selected sources brighter than 10 magnitudes in the H band so that the total exposure time needed to detect a typical dipolar magnetic field of ∼0.5 kilogauss is within 1.5 h, assuming a typical veiling in the H band of about 2 (Doppmann et al. 2005). In the star-forming regions that we explored, Class I and FS sources display H band magnitudes typically ranging from 8 to 12, depending on how deeply embedded they are or how inclined their bipolar cavity is with respect to our line of sight. As a consequence, the magnitude criterion may introduce a bias towards the less embedded objects, i.e. more evolved Class I objects, or towards low inclination. Two of our targets have magnitudes H greater than 10. We had an opportunity to obtain additional observations, so we included them in our sample.

Our initial sample is composed of 18 Class I (6) and FS (12) objects located in different nearby star-forming regions: Perseus, Taurus, rho Ophiuchi, Serpens, and Corona Australis. When known, the spectral types of the observed sources range from G7 to M6, with the majority falling between K1 and M4. The projected rotational velocities (v sin i), when known, range from 5 to 52 km/s. These properties are representative of the typical characteristics observed in embedded protostellar objects (White & Hillenbrand 2004; Doppmann et al. 2005).

Most of the data have been obtained with PI programs allocated over different semesters (PI: E. Alecian). We have also included data from a few Class I sources (V347 Aur, HL Tau, the first observation of V512 Per 2020, GY92 214 and GY92 378) obtained during the SPIRou Legacy Survey (SLS) large program (PI: J.-F. Donati, Donati et al. 2020).

The complete log of the observations is provided in Table D.1. The data of V347 Aur, which have already been published in Donati et al. (2024a), do not appear in that table.

2.2. SPIRou observations and data reduction

The primary goal of SPIRou is to detect and characterise exoplanets around M dwarfs, but also to study magnetic fields in low-mass mature and PMS (Donati et al. 2020). SPIRou covers a spectral range from 950 nm to 2500 nm, with a resolution of 70 000 and a velocity bin per pixel of 2.3 km s⁻¹ (Donati et al. 2020). It consists of an achromatic polarimeter, which splits light into two beams associated with orthogonal states of the selected polarisation, transmitted via optical fibres to an ultra-stable thermally controlled cryogenic spectrograph. The spectrograph disperses light using an echelle grating and a double-pass cross-dispersing prism train, recording the spectra on an H4RG detector.

The polarimeter consists of Fresnel rhombs and a Wollaston prism, enabling simultaneous measurements of both orthogonal states of the selected polarisation (circular polarisation in the present case). The operational process includes stabilising the stellar image on the entrance aperture, injecting it into optical fibres, and using calibration sources (a Fabry-Pérot interferometer, in our case) for precise wavelength and radial velocity calibration.

The camera captures raw spectral data, with echelle orders spread across a 4K × 4K pixel detector. These data include wavelength and light intensity information, enabling the extraction of parameters such as radial velocities and magnetic signatures of stars. The spectrum is divided into 49 orders, each of which requires individual normalisation to maintain the uniformity and precision of the analysed data. Each order comprises two spectra, corresponding to left and right circular polarisations. For each stellar observation, four sub-exposures are taken at successive time intervals with a 90° rotation of the polariser, which helps mitigate instrumental biases during the extraction of the polarisation signal.

We used V0.7.292 version of the dedicated APERO² pipeline (Cook et al. 2022) to reduce the observations. Compared to version V0.7.255 described in Cook et al. (2022), in the APERO versions there have been changes to the background subtraction. The thermal background spectrum is obtained by averaging the spectra acquired from telescope dark frames (i.e. with no sky scene injected). The thermal spectrum is then scaled to match the thermal emission of saturated absorption lines in the K band. An issue was noted, as persistence structures in the J band led to spurious subtraction, being mistaken for low-level thermal emission. Since one does not expect a measurable level of thermal emission below 1.9 μm (considering the rapid decline of a near-room-temperature blackbody bluewards of this wavelength), the thermal emission spectrum was set to zero below 1.9 μm in the newest versions of APERO. It should be noted that, in practise, the thermal background is only subtracted within the K band; thus, this zeroing of the thermal emission at the blue end does not affect the overall thermal background correction. Finally, all individual spectra were corrected for telluric absorption (including OH lines), as described in Cook et al. (2022).

2.3. Continuum normalisation of individual spectra

To analyse the magnetic signals within the photospheric lines, a normalisation of the spectra to their continuum is needed. Automatic continuum normalisation can be performed during the extraction of the polarised spectra. However, the spectra of Class I sources have several variable emission lines that need to be taken into account on a case-by-case basis. While normalisation tools already exist, they lack specific options we aimed to incorporate (necessary for emission line spectra such as those of Class I sources). We developed a normalisation tool for SPIRou data files called Fantasio³. The code aims to reproduce the principles of the IRAF normalisation tool called continuum. The selection of the continuum points is performed by removing non-continuum points using a sigma-clipping procedure over several iterations. The normalisation is performed order by order when a 2D SPIRou spectrum is given as input. The developed interface allows users to modify normalisation parameters interactively for each order, including selecting the polynomial order of the spline and the number of spline interpolations. The sigma-clipping parameters, such as the upper and lower thresholds, and the number of iterations, can also be adjusted to optimise the fit. We allow different threshold values for points above and below the continuum to provide flexibility when both photospheric and emission lines are present within one order. Another interactive feature enables the removal of photospheric and emission lines to exclude them from the fitting process by selecting the regions to reject. This is particularly useful when dealing with spectra containing numerous strong emission lines, as in Class I sources. Normalisation can be performed manually by selecting parameters for each order, or it can be run automatically for multiple observations using pre-defined parameters.

We performed a manual normalisation on both A- and B-spectra, corresponding to the two orthogonal polarisations, for one observation of each star. We then applied the automatic normalisation to all other observations of the same star, using the same fitting parameters per order as those chosen for the fit of the first observation. In general, the order of the spline fit ranges from 1 to 3, as the number of knots of the spline. The sigma clipping can vary from 1 to 3 sigma, and 3 iterations are used. An example of the normalisation process is shown in Fig. 1.

Fig. 1. Upper panel: Order 32 of the spectrum of VV CrA SW without normalisation (solid black line). In blue are the data points selected to fit the continuum using a sigma-clipping of 3 and 2 iterations. The spline fit of order 3 is plotted in red with a number of knots of 2. The two regions with emission lines were manually removed using the interactive tool. Lower panel: Same spectrum after normalisation. The dashed red line shows the continuum normalised to 1.

2.4. Computation of the polarised spectra

The spirou-polarimetry package of Martioli (2020) was used to extract the polarised spectra from our normalised individual intensity spectra in both polarimetric channels. The Stokes V spectrum is obtained by using the ratio method, when the channels A and B spectra are recombined (Martioli 2020). The unpolarised intensity spectrum Stokes I is obtained by adding all channels and sub-exposure spectra. Finally, the null polarisation spectrum (N) is obtained by combining the A and B channels of the four sub-exposures, in a way that cancels the input circular polarisation. It is used to diagnose potential spurious polarised signal (due to instrumental, data reduction issues, or stellar pulsations).

The polarimetry pipeline computes the Stokes spectra by identifying all sequences consisting of four sub-exposures, ensuring that the polarimeter is positioned at the correct angle for each, and calculating the corresponding polarimetric spectra for each sequence, including their associated uncertainties, as described in Bagnulo et al. (2009). The Barycentric Earth Radial Velocity (BERV) shift is corrected before calculating the polarimetric spectra. This pipeline was originally developed to use non-blazed corrected, and non-continuum normalised individual spectra. We have therefore adapted the code to allow the use of our normalised individual spectra. This has also impacted the computation of the uncertainties, which needed to be divided by the blaze and normalised by the continuum. All spectra are extracted with respect to the intensity continuum (I_c), namely I/I_c, V/I_c, and N/I_c. In the following, we adopt the simplified writing, I, V, and N as their normalised value.

3. Analysis and results

In this section we present the results of our analysis performed on the photospheric lines of the spectra, both in intensity and polarisation, with the final objective being to detect the Zeeman signature in the Stokes V spectra, and to determine either the magnetic strength or an upper limit on the dipolar magnetic field when no magnetic field is detected. We detected photospheric lines in the spectra of all our sources, except for one: the FS protostar IRAS 03413+3202. Furthermore, the signal-to-noise ratio (S/N) values obtained for two Class I sources with H magnitude greater than 10 (GY92 214 and GY92 378) are too low to constrain either their effective temperatures or their magnetic fields. In the following, we therefore exclude those three objects. In addition, the dataset of V347 Aur has already been analysed and published in Donati et al. (2024a). The methods used in this paper are similar to those we have used on our data. We have therefore not re-done the analysis for this object, but we kept the source in our sample for the statistical analysis we performed and in the discussion. We present below the photospheric line analysis we performed on the 14 remaining sources.

3.1. Effective temperature determination

We employed two complementary approaches to estimate the effective temperatures (T_eff) of our targets. Our primary method consisted of fitting the observed spectra with synthetic spectra computed using the MARCS atmospheric models (Gustafsson et al. 2008; Plez 2012) and the spectrum-synthesis code ZeeTurbo, following the procedure described in Cristofari et al. (2023).

We assumed solar metallicity and a fixed microturbulent velocity of 1 km s⁻¹ and broadened with a Gaussian profile of full width at half maximum (FWHM) of 4.3 km s⁻¹ to account for the instrumental width. In our fitting process, we fixed the macroturbulence to 0 km s⁻¹ and allowed the projected rotational velocity (v sin i) and veiling in the H and K bands to vary. First, we also allowed the magnetic field to vary. However, the strong veiling and limited S/N of the individual spectra do not allow us to constrain the magnetic field strengths. We therefore choose to perform the fit by assuming null magnetic fields. We indeed observe that the derived effective temperatures are not significantly affected by this choice. Appendix B provides an example of the fits obtained with this procedure.

Each spectrum was visually inspected, and only the fits obtained from the highest quality data were retained. When this method did not yield a reliable solution, due to excessive line broadening, low S/N, or strong veiling, we adopted another method. We have computed synthetic spectra using the MARCS atmospheric models and SYNTH3 radiative transfer code. We used the Vienna Atomic Line Database⁴ (VALD) to compute spectra over the SPIRou wavelength range. We have adopted microturbulent and macroturbulent velocities of 2 km s⁻¹, solar abundances, and a surface gravity of log g = 3.5 (cgs). We have computed a grid of models of effective temperatures between 3500 and 6000 K with a step of 250 K, and of v sin i values from 5 to 80 km s⁻¹ with a step of 5 km s⁻¹. For each observation we have selected one or two spectral regions in the H band, which gather spectral lines of enough S/N, and sensitive to the effective temperature. We used the independently determined veiling and radial velocity values, and then computed a χ² value for each model. The effective temperature was selected based on the χ² minimisation criterion. This method has the advantage to not fit too many parameters that cannot be constrained with the quality of our data, but to constrain the temperature with uncertainties around 250 K, which is enough for the computation of the least-squares deconvolution (LSD) profiles (see Sect. 3.2). A variation of ±500 K of the mask does not significantly change the LSD profile shapes, the B_ℓ values, and the corresponding uncertainties. We emphasise that the quality of our data does not allow us to constrain the microturbulent velocity, and we therefore choose to fix it to usual values for young cool stars, around 1 or 2 km s⁻¹ (Doppmann et al. 2005). As the broadening of the spectral lines is dominated by v sin i (see Table E.1), the choice of the microturbulent velocity value has a negligible effect on the temperature fit. Similarly, the macroturbulent velocity cannot be constrained with our data. It affects only the wings of the spectral lines, which are usually lost in the noise. Our choices of 0–2 km s⁻¹ for the macroturbulent velocity have a negligible impact on the effective temperature determination.

The veiling in the H- and K- bands ⁵ (r_H and r_K) has been measured using our spectra by comparing photospheric spectral lines with those of a cool star template of the same effective temperature, which is evolved enough so that circumstellar material is absent and cannot contaminate the spectrum any further. The method is described in Sousa et al. (2023), and the measured veiling values are summarised in Table E.1.

Determining the effective temperature of VV CrA SW is challenging due to the weakness of its photospheric lines. According to Prato et al. (2003), the SW component does not differ significantly from the NE component. However, the colour-colour diagram (see Fig. 3 in Prato et al. 2003) suggests a slightly lower temperature. VV CrA SW is therefore expected to be in the range of 3500 and 4000 K. We adopt an intermediate value of 3750 K for the effective temperature. This choice has only a minor impact on the LSD profiles and remains consistent with the presence of photospheric features observed in the LSD profiles (see Appendix E), indicating that the corresponding temperature is appropriate.

Our analyses are strongly impacted by the large broadening and strong veiling for several stars. To mitigate the impact on our temperature estimates, we inflated the error bars to reflect the dispersion of the results obtained over several nights, when possible. The final temperatures were rounded to the nearest 250 K, matching the temperature grid of the line masks used for the LSD analysis. The adopted effective temperatures range from 3250 K to 5000 K with a median of ∼4000 K, and are summarised in Table 1.

3.2. Least-squares deconvolved profiles

The LSD is a multi-line technique for computing a pseudo-average line profile as a function of velocity, from a spectrum (Donati et al. 1997). By combining the signal contained in all the spectral lines simultaneously, this method allows us to significantly increase the S/N of our data. The circular polarisation induced by the Zeeman effect at the surface of protostars is expected to be of the order of a few 10⁻⁴ times the continuum intensity, for a typical average line-of-sight magnetic strength of a few hundred gauss, in a moderate rotator (v sin i ∼ 20 km s⁻¹). This signal cannot be detected in individual lines of spectra, even with S/N of 100–200. The LSD method allows an increase of a factor of a few tens in the S/N (depending on the spectral type, and on the level of circumstellar contamination in the spectra), enabling the detection of stellar magnetic fields.

The fundamental principle of the LSD method is based on the hypothesis that the observed spectrum is the result of the convolution of an intrinsic line profile, inherent to the star, with a mask that indicates the wavelength, depth, and magnetic sensitivity of all the lines expected to form in the stellar spectrum. The method works by essentially reversing this process, i.e. deconvolving the observed spectrum by the line mask.

The line masks are generated using VALD (Piskunov et al. 1995), which selects the lines that can be formed in the photosphere of a star for a given effective temperature, surface gravity, microturbulence, and chemical composition. We used the effective temperatures tabulated in Table 1. We assumed a log g of 3.5, expected for Class I sources (Flores et al. 2024), solar abundances, and a microturbulent velocity of 2 km s⁻¹ that is typical in stars with convective envelopes (Doppmann et al. 2005).

In the following, we used the tools provided by the SpecpolFlow package⁶ developed and maintained by Folsom et al. (2025). Before proceeding with LSD, the interactive mask cleaning tool cleanMaskUI was used to clean and tweak the line mask of each source. In the cleaning process, we removed the lines from the original mask that do not comply with the LSD assumption that all spectral lines have the same shape (e.g. emission lines or heavily damped lines like H lines). Additionally, we removed regions strongly contaminated by emission lines, the CO region (CO rovibrational lines from 2.29 μm and above), and the Y and J bands because they were too noisy in our sources. The tweaking process consists of adjusting the relative strengths of the lines of the mask to those observed, within one order. Careful cleaning and tweaking of the mask for each star significantly improve the LSD profile. Finally, to compute the profiles, we used the LSDpy package of SpecpolFlow.

Each LSD profile is computed within ±200 km s⁻¹ around the central velocity of the individual lines, with a velocity pixel size of 2.4 km s⁻¹ to match the mean spacing of the spectral pixels, or 4.8 km s⁻¹ for the fastest rotators (v sin i > 35 km s⁻¹). For all sources, we used a magnetic Landé factor of 1.2, an average depth of 0.2 and a mean wavelength of 1660 nm, as normalisation factors.

The S/N of the LSD profiles is significantly higher than the S/N of the individual spectra, with a gain ranging from 20 to 75 (see Table D.1). The LSD profiles for the I, V, and N spectra are shown in Fig. E.1.

3.3. Magnetic field detection

A Zeeman magnetic signature manifests as a coherent non-zero signal spanning the full width of the photospheric line profile in the circular polarisation profile (V profile) with no counterpart in the null profile (N profile) as shown in Fig. 2. This ensures that the detected signal results from an astrophysical magnetic field, ruling out other phenomena that might influence the observations.

Fig. 2. LSD profiles of V806 Tau. Different colours are used for the different observations. The Y-axes are I/Ic, N/Ic, and V/Ic, where Ic is the intensity of the stellar continuum. The I profile represents the total intensity. Stokes V shows a clear Zeeman signature, which corresponds to an actual magnetic detection because the signal is significantly non-null compared to the surrounding noise and the N profile remains flat in the velocity domain.

We can observe in Fig. E.1 that the N profiles are consistent with zero, confirming that the polarimeter and data reduction perform as expected. If the V profile shows a significant non-zero signal, while the N profile remains null, we are confident that this signature is of stellar origin.

To assess the significance of the non-null Zeeman signature in the Stokes V profile, we can derive the false alarm probability (FAP). It quantifies the likelihood that the observed signal is due to random fluctuations rather than a real physical phenomenon. To ensure robust results, the FAP should be as low as possible. Following Donati et al. (1997), the criteria for interpreting the FAP and evaluating magnetic detections are:

• Definite Detection: the FAP must be below 10⁻⁵, indicating a very low probability that the signal is due to noise, which means the magnetic field is detected.

• Marginal Detection: the FAP should be between 10⁻⁵ and 10⁻³, suggesting that the detection is probable but requires further verification or more observations to confirm it.

• No Detection: the FAP must be greater than 10⁻³, indicating no magnetic field is detected.

We computed the FAP as in Donati et al. (1997) using SpecpolFlow, by computing the reduced χ² of the profile compared to a null model, and applied the χ² statistics to convert them to detection probabilities (P_det). The FAP is then computed as 1 − P_det. This value depends on the uncertainties of the LSD V profile. The cleaning and tweaking of the mask are therefore important. The FAP values are calculated within the photospheric LSD profiles. We therefore used a different velocity range for each star (see Table 3). Those ranges are defined as the regions where the LSD profile intensity is below the continuum in Stokes I. We have also computed the FAP in the N profiles using the same velocity range, as a sanity check. If a detection is obtained in N therefore the V profile is not reliable, and no conclusion can be made. However, it is not the case in our sample. The FAP values are listed in Table E.1

We considered a magnetic field to be definitively detected in a star if we have at least one definite detection in one V observation, with no detection in the associated N observation. According to this detection criterion, our sample contains five magnetic stars (V806 Tau, MHO 3, MHO 2, IRAS 16191-1936, L1689 SNO2), while nine stars do not show magnetic detection. Magnetic field appears to be detected only in FS protostars, this point is discussed in Sect. 4.1.

3.4. Mean longitudinal magnetic field (B_ℓ)

The mean longitudinal magnetic field component (B_ℓ) represents the projection of the magnetic field along the observer’s line of sight, integrated over the observed stellar surface. The mean longitudinal magnetic field B_ℓ is determined, in gauss (G), using the equation of Wade et al. (2000):

$$\begin{array}{c}\hfill B_{\ell }=-2.14\times {10}^{11}\frac{\int \frac{V(v)}{I_c}v\mathrm{d}v}{\lambda zc{\displaystyle \int [1-\frac{I(v)}{I_c}]}\mathrm{d}v}·\end{array}$$(1)

Here λ is the central wavelength of the spectral line in nm, z is the Landé factor, c is the speed of light in vacuum, v is the velocity, V and I are the Stokes parameters, and I_c is the intensity of the continuum. We apply this equation to LSD profiles. We have therefore used the normalisation values from the LSD computations: 1660 nm for λ and 1.2 for z. We used the same velocity range as used for the FAP determination (Table 3). The LSD profiles are not shifted to the stellar rest frame, resulting in asymmetric velocity windows. The uncertainties on B_ℓ are computed by propagating the errors of the I and V profiles in Eq. (1). The values of B_ℓ are listed in Table E.1. They range from ∼80 G to ∼200 G in absolute value for the five magnetic protostars.

4. Discussion

4.1. Study of observational bias

In our final sample, six protostars (including V347 Aur) exhibit clear detections of large-scale magnetic fields, through Zeeman signatures in their Stokes V profiles, confirmed with the FAP analysis. However, no magnetic field was detected in nine targets. This could indicate that these sources host significantly weaker magnetic fields, but it is first necessary to verify that this result is not driven by any observational or methodological bias. To better understand why some stars in our sample did not show detectable magnetic fields, we therefore investigated potential biases by comparing stars with definite magnetic detection (DD) to those with no detection (ND). We analysed how magnetic detectability correlates with various stellar and observational parameters, including the effective temperature, projected rotational velocity, veiling, signal-to-noise ratios, and the uncertainty on the longitudinal magnetic field measurements.

For each parameter, we constructed two cumulative distribution functions (CDFs), one for the DD sample, and one for the ND sample. We then performed a non-parametric Kolmogorov–Smirnov (KS) test to quantify the differences between the two distributions. The KS statistic s represents the maximum distance between the two CDFs, providing a measure of their discrepancy, while the p-value quantifies the statistical significance of this difference. A p-value of p > 0.05 indicates no significant difference, while p ≤ 0.05 suggests a significant difference between the distributions.

For some parameters such as T_eff, v sin i, and H mag, the tests are performed per source, while for the other parameters (e.g. veiling, σ(B_ℓ), S/N) they are performed per observation. The source sample contains the 14 protostars analysed in Sect. 3, to which we added V347 Aur, analysed in Donati et al. (2024a), for which we used the observation of November 22nd 2022, representative of the 80 observations obtained for that source. The observation sample contains all observations listed in Table D.1 plus the V347 Aur observation. The veiling is not available for VV CrA SW because of the weak photospheric lines. The results of the tests are summarised in Table 2.

On the birthline, at effective temperatures below ∼4750 K (Fig. 1 in Villebrun et al. 2019), stars are expected to be fully convective at solar metallicity. For a given mass, as the temperature (hence age) increases, the convective envelope shrinks, possibly resulting in weaker fields, as observed in Class II stars (Gregory et al. 2012). Furthermore, the higher the temperature, the lower the number of spectral lines contained in the mask. Thus, we could, expect a lower incidence of magnetic field detections in hotter stars. Nonetheless, we find that the effective temperature does not appear to significantly influence the occurrence of magnetic detections in our sample. In Fig. 3 we indeed observe that the detected magnetic stars are spread over our temperature range. However, this lack of correlation may not be significant, given the small sample size and the fact that our temperature range is relatively narrow.

Fig. 3. Distribution of projected rotational velocity for each star in the sample as a function of effective temperature. The magnetic stars are represented in red and the non-detected magnetic stars in black. For the magnetic stars, the sizes of the symbols are proportional to the maximum absolute values of the longitudinal magnetic field. For the non-detected magnetic stars, it is proportional to the minimum of the 3σ uncertainty on Bℓ.

Our sample exhibits relatively simple Zeeman signatures, suggesting that the magnetic geometry is dominated by low-order large-scale magnetic field. The rotational broadening of the spectral lines could also affect magnetic detectability. For a similar magnetic field, the Stokes V signature in a star with a low v sin i will have a larger amplitude than in a star with a large v sin i, as the same Zeeman signature will be spread over a smaller velocity range in the first case compared to the second one (presuming that the large-scale field is simple enough and regions of opposite polarities do not cancel out their respective signatures). If the S/N is the same in both cases, then the signal in the second case will be more diluted in the noise, and the FAP will be higher. When v sin i is larger, the magnetic detection probability is therefore lower.

Whatever the sample considered, we find that the null hypothesis cannot be rejected. Therefore, no significant difference is found between the DD and ND v sin i distributions. Remarkably, we are even able to detect magnetic fields in stars with v sin i as high as ∼55 km s⁻¹. At the same time, large-scale magnetic fields are not detected in many protostars with low v sin i (Fig. 3). Because, for almost all stars in our sample, the v sin i was unknown or poorly constrained at the time of the observations, we could not take this parameter into account in the computation of our exposure times. Our observational strategy therefore appears adequate as it allowed a potential v sin i observational bias to be avoided.

Another parameter that may influence magnetic detections is the veiling. Veiling refers to the dilution of photospheric absorption lines due to excess continuum emission, typically interpreted as a signature of active accretion in young stars (Connelley & Greene 2010). If a star exhibits strong veiling, this can affect the detectability of magnetic fields by reducing the line contrast with respect to the intensity continuum, and therefore increasing the relative noise in Stokes I and V profiles (see Appendix C). The KS tests reject the null hypothesis at a confidence level higher than 95% when comparing the veiling distributions of the DD and ND samples, in both H and K bands, and in both statistical samples. In Fig. 4, we indeed find that the higher the veiling, the lower the detection rate.

Fig. 4. Veiling in the K band as a function of the veiling in the H band for all observations of the sample. The size of the symbol is proportional to the uncertainty of the longitudinal magnetic field and each symbol corresponds to a different target. Magnetic stars are represented in red and stars without magnetic detections in black.

This indicates that the detectability of magnetic fields in our sample is influenced by the veiling. For some objects we have more observations than for others. We therefore verified that the veiling–detectability trend remains when these targets (HL Tau and IRAS 03247+3001) are excluded. If the veiling is not known before setting up the exposure time of the observation, and because the depths of the photospheric lines in a veiled spectrum are shallower than in an unveiled one, the veiled photospheric line profiles will be noisier than the unveiled ones, which lowers the contrast of the Zeeman signatures and may hamper Stokes V detections (see Appendix C). The fact that we do not detect magnetic fields in the most veiled observations could therefore be an observational bias. Moreover, veiling in Class I is higher than in Class II, which could also explain the lower detectability in our sample of Class I and FS protostars.

Figure 4 illustrates that magnetic detections are found only for the lowest values of r_H and r_K. The targets with the highest veiling are those for which the analysis of the photospheric lines is most challenging. A significant discrepancy between the adopted and the true effective temperature could also affect the derived veiling values. However, this would modify the veiling by at most a factor of two, which does not alter our conclusions, as the associated uncertainties remain smaller than the overall veiling amplitude observed across our sample. We also remark that magnetic detections are only observed in FS protostars, while none of the Class I protostars is detected. Class I sources are usually more veiled than FS sources, which could explain this result.

Finally, we may also wonder whether the quality of our data affects magnetic detection. We find no significant differences between the S/N distributions of the DD and ND samples, considering the S/N in either the spectra (S/N_s) or the LSD profiles (S/N_LSD).

Our magnetic detection rate is not affected by the S/N of our data but is affected by veiling. This means that we should observe a correlation between the magnetic field detections and the uncertainties of the longitudinal field measurements (see Appendix C). We find that the distributions of σ(B_ℓ) for the DD and ND samples are significantly different (Table 2). Furthermore, Fig. 4 shows that the uncertainties on B_ℓ are systematically lower for the least veiled spectra. The veiling phenomenon could therefore be the main source of difficulty for detecting magnetic fields in Class I and FS sources.

4.2. Upper limits on the dipole magnetic field

Stars in which we did not detect a magnetic field may still host one at their surface, but it may be too weak or too complex to be detected with our current data quality. In this section we aim to estimate the maximal magnetic field strength a star can have that would still remain hidden in the noise of our data. In the case of a complex magnetic field, this is not easy because such fields can exhibit a large variety of topologies, and a large amplitude of strength, which could remain hidden within a given noise level. On the other hand, magnetic fields observed in fully convective stars (i.e. in CTTSs or M dwarfs) are mainly dipolar (Donati et al. 2011; Gregory et al. 2012; Morin et al. 2008, 2010). Furthermore, our sample mostly contains simple Zeeman signatures, most likely produced by low-order magnetic fields. We therefore only consider a dipole to compute the upper limit.

To compute the upper limits on the magnetic field strength at the pole for each of these stars, we assume an oblique rotator (Stibbs 1950), i.e. a star hosting a dipole of strength B_pol at the pole, inclined by a magnetic obliquity β with respect to the rotation axis Ω, which is itself inclined by an angle i with respect to the line of sight. For a given star with a given v sin i, and a given depth d_I of the LSD I profile, we can predict the shape of the Stokes V profile at a given rotation phase, compare its amplitude to the noise of our observations, and use the FAP method to determine whether it would be detectable or not.

The weak field approximation is assumed to compute the local Stokes V profiles (Landi Degl’Innocenti & Landolfi 2004). We can then run Monte-Carlo simulations assuming random i, β, and rotation phase for one value of B_pol, and count the number of detections. We repeat this process for several values of B_pol, and trace the detection probability curve (e.g. Fig. A.1). To be consistent with the literature, we chose the upper limit B_pol at a detection probability of 90%. The full description of the methodology is detailed in Alecian et al. (2016), and recalled in Appendix A.

We first computed the upper limits for the non-detected stars using marginal and definite FAP criteria. The resulting values are summarised in Table 3. These upper limits are relatively high, ranging from approximately 500 G to more than 5 kG, suggesting that relatively strong dipoles might indeed be present in these stars, but are below the sensitivity of our observations.

To assess the reliability of these upper limit values, we also computed them for stars with detected magnetic fields. We find values between approximately 2 kG and more than 5 kG. We observe that the highest upper limits can be found independently in detected and non-detected magnetic protostars. The projected rotational velocity affects the upper limit on the dipole field strength. A higher v sin i dilutes the Stokes V signature, implying that stronger magnetic fields are required for a detection in stars with large v sin i, compared to stars with low v sin i. We find that all sources with a limit above 5 kilogauss have a v sin i larger than 35 km s⁻¹. The fact that such large limits are found in DD and ND protostars means that large v sin i is not the main parameter that hampers the detection of magnetic field. This is in line with the similar v sin i distributions that we find between the DD and ND samples (Sect. 4.1).

It is interesting to note that both ND and DD samples have similar upper limits. This suggests that if our non-detected sample hosts dipolar magnetic fields of comparable strength to those of the detected stars, they would likely have been detected. The non-detected stars might then host weaker or more complex magnetic fields, if any.

4.3. Impact of our magnetic measurements on the origin and evolution of stellar magnetic fields

We conducted a spectropolarimetric analysis of a sample of Class I and FS protostars, aiming to investigate the presence of large-scale magnetic fields. All young, low-mass observed protostars are expected to be fully convective or to have a large convective envelope (Palla 1998; Baraffe et al. 2017), we expect all of them to host detectable magnetic fields, generated by convective dynamos, as in low-mass main-sequence stars. However, magnetic fields were detected in only 6 out of the 15 stars of our final sample, corresponding to 40% of the sample. This relatively low detection rate raises several questions regarding the magnetic origin of these young stars, their internal processes, and their magnetospheric environments. In particular, the absence or apparent non-detection of large-scale magnetic fields in a substantial fraction of Class I and FS objects challenges the assumption that such fields are ubiquitous at these evolutionary stages.

A central question that emerges is whether large-scale magnetic fields, commonly observed in more evolved systems such as CTTSs, are truly universal across all largely convective stellar objects. Assuming that our protostars lie on the birthline, we expect them to be fully convective up to an effective temperature of 4750 K, and have a convective envelope of mass larger than 40% of the stellar mass up to 5250 K (Villebrun et al. 2019). Under this assumption, ∼80% of our sample is fully convective, while the remainder has a radiative core plus a large convective envelope. We can also estimate the range of mass probed by our sample using the effective temperatures derived in Sect. 3.1. We find a range of 0.2–2.5 M_⊙. This means that their direct descendants are the CTTSs of low- and intermediate-mass (Herbig 1962; Villebrun et al. 2019). The magnetic properties of these objects have been intensely characterised thanks to dense spectropolarimetric monitoring, and ZDI analysis (e.g. Hussain et al. 2009; Donati et al. 2013; Pouilly et al. 2020; Nowacki et al. 2023; Zaire et al. 2024). Magnetic fields are detected in fully convective stars, and almost all of them host a strong (> 1 kG) magnetic field, dominated by a dipolar component, and sometimes a multipolar magnetic field. When published, the maximum absolute B_ℓ values are larger than 400 G and can reach 900 G. Our snapshot observations do not allow us to perform ZDI to retrieve the topology of the magnetic field we have detected in our sample. However, we do not detect strong magnetic fields in all of them: the maximum absolute B_ℓ value is 200 G. A single measurement of B_ℓ is strongly dependent on the rotation phase. Our phase sampling is less complete than that of the ZDI studies, and therefore we may have missed the rotation phases when B_ℓ is maximum. However, out of the 16 observations we obtained for the 5 magnetic protostars, none are higher than 200 G, which are at most similar, but mainly weaker than in fully convective CTTSs (Donati et al. 2024b,c; Nowacki et al. 2023).

On one hand, Class I and FS protostars have a larger radius than CTTSs as they are less evolved. Class I and FS protostars are younger than 1 Myr, but probably older than 0.5 Myr (Fiorellino et al. 2023). CTTSs are 1–3 Myr old (e.g. Zaire et al. 2024). According to stellar evolutionary models (e.g. Baraffe et al. 2015), from 0.5 Myr to 3 Myr, the radius decreases by a factor up to 2. If we assume a constant magnetic flux as the protostars evolve towards the pre-main sequence, then the surface magnetic field can increase by a factor up to 4, which could explain why the B_ℓ in our sample are in general weaker than at later stages. On the other hand, it is not clear whether magnetic flux conservation can be invoked when magnetic fields are of dynamo origin. The major structural difference between Class I and fully convective Class II stars is the surface gravity, hence the gravity gradient within the convection zone. Magneto-hydro-dynamical models of convective stars of various stratifications are required to understand the role of stratification in the efficiency of dynamo processes. Within the PROMETHEE project, we are developing models for different stratification levels that will help us to interpret this result (Guseva et al., in prep.).

Finally, it is interesting to mention that in fully convective main-sequence M dwarfs, a phenomenon known as dynamo bistability has been observed (Morin et al. 2010). In this scenario, stars with identical global properties (mass, rotation rate, convection regime) can exhibit either strong predominantly dipolar magnetic fields or much weaker, complex, and multipolar fields. It is plausible that a similar mechanism is at play in our sample, and could explain part of the non-detections.

4.4. Comparison with the total magnetic fields

A study of 32 Class I and FS protostars was conducted by Flores et al. (2024), focussing on measuring the average strength of the total magnetic field. They used K-band observations with the high-resolution near-infrared echelle spectrograph iSHELL on the IRTF (Hawaii), achieving a spectral resolution of 50 000. Their targets share similar star-forming regions, magnitudes, and spectral indices with ours. They used the radiative transfer code Moog Stokes (Deen 2013) to derive stellar parameters (effective temperature, surface gravity, mean magnetic field strength ⟨B⟩, and projected rotational velocity v sin i), as well as stellar masses and ages, based on magnetic stellar evolutionary models from Feiden (2016).

As their data are of similar quality to ours, and they used a model similar to ours to fit the spectra, it is first interesting to compare their results with ours. Both studies cover similar effective temperatures, and v sin i ranges (see Fig. 5). Flores et al. (2024) reports a mean effective temperature of 3450 K and an average v sin i of 26 km s⁻¹, while our sample shows a higher average effective temperature of 4000 K and an average v sin i of 36 km s⁻¹. The distributions of the two samples are therefore slightly different. In proportion, we have more stars with v sin i larger than 40 km s⁻¹, and with effective temperatures higher than 4000 K. In addition, we find that we have four targets in common with their sample: MHO 2, V806 Tau, L1689 SNO2, and V347 Aur, which enable a direct comparison. We find similar projected rotational velocities, v sin i, although our derived effective temperatures are slightly higher, but still consistent within the uncertainties.

Fig. 5. Projected rotational velocities v sin i as a function of effective temperatures measured in this work (blue) and in the work of Flores et al. (2024, red). The points connected by the dashed line indicate the common stars in the two samples.

Out of their 32 sources, Flores et al. (2024) was able to measure non-null ⟨B⟩ in 27 of them, and report upper limits for the 5 remaining sources between 0.8 and 1.2 kG. Among the 27 sources with a measured magnetic field, they report an average magnetic field strength of approximately 2 kG, consistent with the values typically found in more evolved Class II sources (Johns-Krull 2007). These measurements, based on Zeeman broadening, are sensitive to the total magnetic field, including both large and small-scale components. In contrast, our spectropolarimetric method is sensitive only to large-scale magnetic fields. These two approaches probe different components. While Zeeman broadening captures the cumulative surface magnetic field regardless of orientation, spectropolarimetry detects only the longitudinal component of the field. This means that regions of opposite polarity on the stellar surface can cancel each other out. If the stellar magnetic field is only organised on small scales, or is a combination of large- and small-scales fields, the spectropolarimetric method leads to an underestimation of the total magnetic field (See et al. 2019; Lehmann et al. 2019).

As already mentioned, our snapshot sample cannot be used to determine the topology and strength of the large scale magnetic field, and therefore estimate its contribution to the total magnetic field measured through Zeeman broadening. However, if we assume that our magnetic sample contains only dipolar fields, then we can use the statistical relationship between the cumulative longitudinal magnetic field (⟨B_ℓ⟩_rms) of a magnetic sample, and its average total magnetic strength ⟨B⟩_av. Kochukhov (2024) determined a statistical relationship between the two parameters by running many simulations of dipoles at the surface of an A0 main-sequence star for random magnetic obliquities and inclinations, and of a linear limb-darkening parameter of 0.5. The average value of limb-darkening parameter of our sample is 0.4 (Claret & Bloemen 2011). When correcting for our different limb-darkening value, the relationship can be written as

$$\begin{array}{c}\hfill {⟨B⟩}_{\mathrm{av}}=4.{04}_{-0.5}^{+1.1}{⟨B_{\ell }⟩}_{\mathrm{rms}},\end{array}$$(2)

where ⟨B_ℓ⟩_rms is the root-mean-square of our measurements of B_ℓ in our magnetic sample.

Based on Eq. (2), the mean surface field modulus for our magnetic sample would be approximately 450 ± 200 G, assuming it hosts only dipolar fields.

This value is well below the average value obtained from Zeeman broadening by Flores et al. (2024), as in most low-mass stars, i.e. consistent with dynamo fields. This suggests that a significant portion of the magnetic energy in these protostars is stored in small-scale magnetic structures, likely concentrated in active regions or spots, rather than in a globally organised dipole.

Complementary observations, including the monitoring of some targets, will allow us to better measure and characterise the full magnetic topology through ZDI. This method will provide a more detailed view of the large-scale field and a more complete picture of the magnetic environments in Class I and FS protostars.

When they are strong enough, large-scale fields are responsible for funnelling accretion flows, generating hot spots, and regulating angular momentum transport through magnetospheric interaction (Bouvier et al. 2014; Romanova & Owocki 2015). Meanwhile, small-scale magnetic activity contributes to coronal heating, X-ray emission, and potentially to disk ionisation and chemistry, impacting planet formation conditions from the earliest evolutionary stage of solar-type stars (Grosso et al. 2020). Therefore, characterising both the small and large-scale magnetic fields is essential to build a coherent picture of the magnetic environment in protostars. Such an analysis will be presented for a few stars in a following paper.

5. Conclusions

We have carried out a detailed analysis of 14 Class I and FS protostars using SPIRou, aiming at characterising their stellar and large-scale magnetic properties. When including V347 Aur, the final statistical sample comprises 15 objects. When possible, we derived effective temperatures, and using LSD profiles, we measured the detection probability of a Zeeman signature in the Stokes V profiles, the longitudinal magnetic fields, the projected rotational velocities, and the radial velocities. Our results can be summarised as follows:

1. We detected large-scale magnetic fields in 6 out of 15 protostars, with absolute longitudinal field strengths ranging from approximately 80 G to 200 G. In contrast, nine objects show no clear magnetic signature.

2. The effective temperatures for our sample of Class I and FS sources range from 3250 K to 5000 K, with a median value of approximately 4000 K.

3. We find no significant correlation between magnetic field detection and stellar parameters, such as effective temperature or projected rotational velocity.

4. We find a correlation between magnetic field detection and veiling, suggesting that veiling influences the detectability of magnetic fields by diluting the signal in the Stokes V profile.

5. We derived upper limits on the dipolar field strength, with a median value below 3 kG, indicating that strong magnetic fields may still be present, but are below our detection threshold when no magnetic field is detected.

6. The longitudinal magnetic field values of our magnetic sample do not seem to reach the higher levels observed in their evolutionary descendants, i.e. the fully convective CTTSs.

7. Our statistical estimate of the average total field strength of our magnetic sample is significantly lower than those obtained via Zeeman broadening studies (e.g. Flores et al. 2024). This may suggest that a significant portion of the magnetic energy in these protostars is stored in small-scale magnetic structures, likely concentrated in active regions or spots, rather than in a globally organised dipole.

Our analysis highlights the challenges of detecting large-scale magnetic fields in Class I and FS protostars, especially when the fields are weak, complex, or below our detection thresholds. This naturally raises the question of how early in the stellar formation process magnetic fields become detectable and dynamically important.

In CTTSs, the magnetic star–disk interaction is well documented (Bouvier et al. 2007; Hartmann et al. 2016). The strong dipolar component of the stellar magnetic field truncates the inner disk at a few stellar radii, and funnels accreting material along magnetic field lines, forming accretion columns. In a following paper, we will investigate whether there is any correlation between magnetic field detections and the equivalent widths of Brγ and Paβ emission lines, which are commonly used as accretion tracers in embedded protostars (Antoniucci et al. 2008; Fiorellino et al. 2021). We cannot rule out that magnetospheric accretion processes are at play in the all sources that we observed. If magnetospheric accretion occurs in Class I and FS objects, strong kilogauss large-scale magnetic fields are likely required, given the higher mass accretion rates of such objects compared to Class II (Fiorellino et al. 2023).

The comparison with Zeeman broadening measurements, which are sensitive to the total surface magnetic field, suggests that a significant fraction of the magnetic energy in these young stars may reside in small-scale field structures, undetectable through circular polarisation due to spatial cancellation, or that the magnetic topologies may vary widely, possibly due to dynamo bistability, as seen in fully convective M dwarfs (Morin et al. 2010). This underlines the importance of combining different observational diagnostics to obtain a comprehensive view of stellar magnetism.

Importantly, recent studies have pushed the observational frontier to even earlier stages. Although Class 0 protostars are even more deeply embedded and challenging to observe, recent studies have begun to detect photospheric absorption lines in these objects (Greene et al. 2018; Le Gouellec et al. 2025). These detections open promising prospects for future spectropolarimetric studies of Class 0 sources, offering the potential to provide critical constraints on the initial conditions of dynamo generation and accretion physics, and trace magnetic field evolution from the very first stages of stellar formation. Class 0 objects have accretion diagnostics that are similar to those in the more evolved Class I and II protostars (Laos et al. 2021; Le Gouellec et al. 2024), perhaps suggesting that strong and organised stellar magnetic fields are present.

This study demonstrates the capabilities of SPIRou to probe the large-scale magnetic fields of embedded protostars, making it a powerful tool to characterise magnetism in the early stages of star formation. The next step is to perform ZDI on our magnetic sample. This will unveil the topology of the magnetic field in Class I and FS objects, which has not been done to date. These topologies will be a key element in understanding the evolution of magnetic fields in protostars and of the star–disk interactions.

Acknowledgments

We want to thank the referee for helpful comments which significantly improved the clarity of the manuscript. This work has received funding from the French Agence Nationale de la Recherche (ANR) through the project PROMETHEE (ANR-22-CE31-0020). This research has made use of the SIMBAD database, the VizieR catalogue access tool, and the Aladin sky atlas, CDS, Strasbourg Astronomical Observatory, France. OK acknowledges support by the Swedish Research Council (grant agreement no. 2023-03667) and the Swedish National Space Agency. AK acknowledges support by the NKFIH NKKP grant ADVANCED 149943 and the NKFIH excellence grant TKP2021-NKTA-64. Project no. 149943 has been implemented with the support provided by the Ministry of Culture and Innovation of Hungary from the National Research, Development and Innovation Fund, financed under the NKKP ADVANCED funding scheme. CPF acknowledges funding from the European Union’s Horizon Europe research and innovation programme under grant agreement No. 101079231 (EXOHOST) and from UK Research and Innovation (UKRI) under the UK government’s Horizon Europe funding guarantee (grant number 10051045). H.N. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 101019653). V.J.M.L.G. acknowledges support by the Spanish program Unidad de Excelencia María de Maeztu CEX2020-001058-M, financed by MCIN/AEI/10.13039/501100011033, and by the MaX-CSIC Excellence Award MaX4-SOMMA-ICE. V.J.M.L.G. acknowledges support by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No. 101098309 – PEBBLES).

References

Appendix A: The dipole strength upper limit method

We implemented a method to derive the projected rotational velocity, line depth, and radial velocity from observed stellar absorption line profiles by fitting synthetic profiles computed from a physical stellar model. The model computes spectral line profiles based on the three key parameters. These parameters are iteratively updated in a configuration file before each model run. The synthetic profiles are compared to observed data by computing the total χ² statistic, accounting for measurement uncertainties. A Nelder–Mead simplex algorithm is employed to minimise χ² by adjusting the parameters simultaneously, searching best-fit values for line depth, v sin i and v_rad. Confidence intervals for each parameter are derived through one-dimensional χ² scans. For each parameter, values around the optimum are scanned while re-optimising the others, and uncertainties correspond to Δχ² = 1, representing 68 per cent confidence intervals.

The Nelder-Mead algorithm is particularly well suited for minimising functions that are difficult to differentiate, such as a χ² function. To better constrain the initial parameters, we estimated them with a Gaussian fit, which provides a reliable estimate of the most probable initial values. To estimate the uncertainties associated with the optimised parameters, we performed a detailed χ² analysis. Starting from the optimal parameter values, each parameter was varied individually by 20% around its best-fit value while keeping the other parameters fixed. For each parameter, the χ² function was evaluated at 20 evenly spaced points within this range. The minimum total χ² value is computed, the upper and lower bounds of the parameter were found by locating the points where ${\chi }^2={\chi }_{\text{min}}^2+1$. These bounds define the asymmetric uncertainties (σ₋ and σ₊) for each parameter, reflecting the sensitivity of the χ² function to changes in the parameter values. To further validate the reliability of the uncertainties, we compared the results with other methods, such as Monte Carlo simulations and Probability Density Function analysis. Additionally to the one-dimensional analysis, we performed a two-dimensional χ² analysis to examine correlations between pairs of parameters. In this case, the criterion for the 68% confidence region corresponds to Δχ² = 2.4, as there are two free parameters. These complementary approaches confirmed the robustness of the uncertainty estimates, ensuring that the reported confidence intervals are both accurate and reliable.

The detection limit of the magnetic field, denoted B_lim, represents the value below which the magnetic field signal can remain below the noise level (Alecian et al. 2009, 2016). To estimate the maximum strength, we use the oblique rotator model (Stibbs 1950). We assume a dipolar configuration, which is the simplest magnetic model (see Fig. 9 in Alecian et al. 2016), and fit the corresponding local profiles to our LSD I observations. The three free parameters in the fitting procedure are the I line depth, the v sin i, and the radial velocity. We fitted each observation independently, which results in slightly different values of these parameters for different observations of the same star. However, the measurements remain consistent within the error bars. The variations observed in v sin i can be explained by line profile variability in the Stokes I profiles (see Fig E.1).

Next, a Monte-Carlo method is used to simulate a grid of 1000 synthetic V profiles for each observation, where random values are assigned to the model parameters. Each model uses a random inclination angle i, obliquity angle β, and rotational phase Ω, with various values of the polar magnetic fields B_pol using the same mean Landé factor (1.2) and wavelength (1660 nm). The local line profiles are computed using an instrumental broadening corresponding to a resolving power of 70 000 and a macroturbulent velocity of 5 km s⁻¹. The shapes of the wings of the LSD I profiles are not compatible with a null macroturbulent velocity. Furthermore, if we let it be free in the fit, we obtain values that are unphysical. We therefore made a compromise and chose to fix it to a value that reproduces well the wings of the LSD profiles, but is not too high. We should caution that the computation of the LSD profiles themselves may produce artificially broadened wings, for example because of line blending, or may produce depressed continuum. It can therefore be higher than the value needed to fit the individual spectral lines for the temperature determination in Sect. 3.1. This value should therefore be interpreted with caution. The surface integration of the local profiles is done using a linear limb-darkening coefficient. These values have been extracted from the Claret & Bloemen (2011) database for the corresponding effective temperatures of our stars (Sect. 3.1), assuming a surface gravity of 3.5, and solar abundances. Because our LSD masks are dominated by H-band lines, we choose the coefficients computed in this band, and we used the mean values of both methods (least-squares or flux conservation) used by Claret & Bloemen (2011) to compute them.

Gaussian noise is added to each profile. This process is repeated for all the observations, with the 1000 models of V profiles enabling the calculation of the detection probability for a range of magnetic field strengths. For each set of models, the probability of detecting a magnetic field in each simulated V profile is computed. This allows us to determine whether a magnetic field is detected or not, using the FAP. We ensured that the FAP does not exceed a prescribed value, using FAP = 10⁻³ for marginal detection (Donati et al. 1997). We chose the marginal criterion to be consistent with previous publications, and be able to compare with other samples. We also compute the upper limit using definite criterion. However, the choice of this criterion is arbitrary, as all values will be affected in a similar way. For each value of the field strength, we calculated the number of detections across several models for each observation and plotted the detection rates as a function of field strength. A detection rate of 90 per cent is required to statistically consider that the field would be detected in one observation. This enables us to determine the upper limit of the field strength, below which we are unable to detect it in our observations. Assuming the observations are all independent, we can compute the combined probability following the method outlined by Neiner et al. (2015), for various polar field strengths, and apply the 90 per cent detection threshold to derive the upper limit. By analysing this curve, we can determine a value of B_lim. This corresponds either to a profile containing only noise or to one in which the signature of a weak magnetic field is hidden within the spectral noise in Stokes V (see Fig A.1).

Fig. A.1. Example of dipolar magnetic detection probability curve as a function of the dipole strength. The coloured curves represent the probability for each individual observation, while the black curve represents the combined probability for all observations. The 90 per cent threshold is used to determine the field strength upper limit.

Appendix B: Fitting example with ZeeTurbo

The best fit for one observation of IRAS 04113+2758 N, showing all the lines used in the temperature fitting.

Fig. B.1. ZeeTurbo fitting process for the first observation of IRAS 04113+2758 N. The observed spectrum is shown in black, and the best-fit model in red. The blue shading indicates the regions used for parameter fitting, and the corresponding absorption lines are labelled.

Appendix C: The effect of veiling on the uncertainties of B_ℓ

As our sources can be strongly affected by veiling, i.e. a continuum excess that adds to the stellar continuum, we explore its effect on the uncertainties on the B_ℓ measurement. We can re-write Eq. (1) as

$$\begin{array}{cccccc}& B_{\ell }=\alpha \frac{X}{Y}\hfill & \hfill & \text{where}\hfill & & X={\displaystyle \int \frac{V(v)}{I_c}}vdv\\ \hfill & \hfill & \hfill & Y={\displaystyle \int [1-\frac{I(v)}{I_c}]}dv\hfill \end{array}$$

and α = −2.14 × 10¹¹/(λzc)≈ − 358 G km⁻¹ s when using the normalisation values of λ and z of our LSD profiles. The numerical computation of X and Y, and of their uncertainties σ_X and σ_Y can be written

$$\begin{array}{cccc}& X={\displaystyle \sum _{i=1}^N}{\left(\frac{V}{I_c}\right)}_i{v}_i\mathrm{\Delta }{v}_i\hfill & & {\sigma }_X=\sqrt{\sum _{i=1}^N{\left[{\left(\frac{{\sigma }_V}{I_{\mathrm{c}}}\right)}_i{v}_i\mathrm{\Delta }{v}_i\right]}^2}\hfill \\ \hfill & Y={\displaystyle \sum _{i=1}^N}[1-{\left(\frac{I}{I_c}\right)}_i]\mathrm{\Delta }{v}_i\hfill & & {\sigma }_Y=\sqrt{{\displaystyle \sum _{i=1}^N}{\left[{\left(\frac{{\sigma }_I}{I_{\mathrm{c}}}\right)}_i\mathrm{\Delta }{v}_i\right]}^2}\hfill \end{array}$$

where σ_I and σ_V are the uncertainties on the unnormalised I and V spectra, and N is the number of pixels within the photospheric profile. The uncertainty on B_ℓ is therefore

$$\begin{array}{cc}& \frac{{\sigma }_B}{B_{\ell }}=\sqrt{{\left(\frac{{\sigma }_X}{X}\right)}^2+{\left(\frac{{\sigma }_Y}{Y}\right)}^2}\hfill \end{array}$$

We consider a star that we observe twice. During the first observation the star is not affected by veiling, and emits a spectrum I₁(v) = I₀(v) with a continuum I_1c = I_★. During the second observation the star is contaminated with an additional excess continuum I_ex(v) = rI_★(v), where r is defined as the veiling parameter. We assume that r does not depend on the wavelength, and that the region emitting the veiling does not emit circular polarisation. During the second observations we therefore record the spectrum I₂(v) = I₀(v)+I_ex(v), with a continuum I_2c = I_★ + I_ex. Both polarised spectra are assumed identical, and an associated polarised spectrum V₂(v) = V₁(v) = V(v). X and Y for both observations are therefore

$$\begin{array}{cccccc}& X_1\hfill & \hfill ={\displaystyle \sum _{i=1}^N}{\left(\frac{V}{I_{★}}\right)}_i{v}_i\mathrm{\Delta }{v}_i& & X_2& ={\displaystyle \sum _{i=1}^N}{\left(\frac{V}{I_{★}+I_{\mathrm{ex}}}\right)}_i{v}_i\mathrm{\Delta }{v}_i\\ \hfill & Y_1\hfill & \hfill ={\displaystyle \sum _{i=1}^N}[1-{\left(\frac{I}{I_{★}}\right)}_i]\mathrm{\Delta }{v}_i& & Y_2& ={\displaystyle \sum _{i=1}^N}[1-{\left(\frac{I+I_{\mathrm{ex}}}{I_{★}+I_{\mathrm{ex}}}\right)}_i]\mathrm{\Delta }{v}_i\\ \hfill & ={\displaystyle \sum _{i=1}^N}{\left(\frac{I_{★}-I}{I_{★}}\right)}_i\mathrm{\Delta }{v}_i\hfill & & ={\displaystyle \sum _{i=1}^N}{\left(\frac{I_{★}-I}{I_{★}+I_{\mathrm{ex}}}\right)}_i\mathrm{\Delta }{v}_i\hfill \end{array}$$

As the continuum in the LSD profiles is not dependent on v, when dividing X by Y in both observations, the numerators cancel, and X₁/Y₁ = X₂/Y₂, hence

$$\begin{array}{c}\hfill B_{\ell 1}=B_{\ell 2}=B_{\ell }.\end{array}$$

To explore how veiling affects the computation of the uncertainties, we assume first that both observations were obtained in identical weather conditions, with the same exposure time. The noise is therefore identical in both observations: σ_I1 = σ_I2 = σ_I, and σ_V1 = σ_V2 = σ_V

$$\begin{array}{ccccc}& {\sigma }_{X1}\hfill & \hfill =\sqrt{{\displaystyle \sum _{i=1}^N}{\left[{\left(\frac{{\sigma }_V}{I_{★}}\right)}_i{v}_i\mathrm{\Delta }{v}_i\right]}^2}& {\sigma }_{X2}\hfill & =\sqrt{{\displaystyle \sum _{i=1}^N}{\left[{\left(\frac{{\sigma }_V}{I_{★}+I_{\mathrm{ex}}}\right)}_i{v}_i\mathrm{\Delta }{v}_i\right]}^2}\\ \hfill & {\sigma }_{Y1}\hfill & \hfill =\sqrt{{\displaystyle \sum _{i=1}^N}{\left[{\left(\frac{{\sigma }_I}{I_{★}}\right)}_i\mathrm{\Delta }{v}_i\right]}^2}& {\sigma }_{Y2}\hfill & =\sqrt{{\displaystyle \sum _{i=1}^N}{\left[{\left(\frac{{\sigma }_I}{I_{★}+I_{\mathrm{ex}}}\right)}_i\mathrm{\Delta }{v}_i\right]}^2}\end{array}$$

When dividing σ_X by X and σ_Y by Y, the denominators of X and Y cancel, and therefore

$$\begin{array}{c}\hfill \frac{{\sigma }_{B1}}{B_{\ell 1}}=\frac{{\sigma }_{B2}}{B_{\ell 2}}=\frac{{\sigma }_B}{B_{\ell }}.\end{array}$$

We note that in that case, in the second observation the received signal is higher than in the first observation. Its S/N is therefore higher than that of the first one: ${(S/N)}_2={(S/N)}_1\sqrt{1+r}$. However, because we received the same flux from the star itself, it is normal that the uncertainty on B_ℓ is not affected.

We now assume a second case where we observe two identical stars of the same magnitude. One is veiled, and the other one is not. However, before scheduling the observation, we do not know that the second star is veiled. We therefore observe them with the same exposure time, in the same weather conditions. Because the magnitudes are the same, the recorded fluxes for both stars are the same, and therefore (S/N)₂ = (S/N)₁. The uncertainties on the flux from the stars, however, are not the same:

$$\begin{array}{c}\hfill {\sigma }_{I2}={\sigma }_{I1}\sqrt{1+r},{\sigma }_{V2}={\sigma }_{V1}\sqrt{1+r}.\end{array}$$

Substituting these expressions into the definitions of σ_X and σ_Y give

$$\begin{array}{ccccc}& {\sigma }_{X1}\hfill & \hfill =\sqrt{{\displaystyle \sum _{i=1}^N}{\left[{\left(\frac{{\sigma }_{V1}}{I_{★}}\right)}_i{v}_i\mathrm{\Delta }{v}_i\right]}^2}& {\sigma }_{X2}\hfill & =\sqrt{{\displaystyle \sum _{i=1}^N}{\left[\frac{\sqrt{1+r}}{(1+r)}{\left(\frac{{\sigma }_{V1}}{I_{★}}\right)}_i{v}_i\mathrm{\Delta }{v}_i\right]}^2}\\ \hfill & & & =\frac{{\sigma }_{X1}}{\sqrt{1+r}}\hfill \\ \hfill & {\sigma }_{Y1}\hfill & \hfill =\sqrt{\sum _{i=1}^N{\left[{\left(\frac{{\sigma }_{I1}}{I_{★}}\right)}_i\mathrm{\Delta }{v}_i\right]}^2}& {\sigma }_{Y2}\hfill & =\sqrt{\sum _{i=1}^N{\left[\frac{\sqrt{1+r}}{(1+r)}{\left(\frac{{\sigma }_{I1}}{I_{★}}\right)}_i\mathrm{\Delta }{v}_i\right]}^2}\\ \hfill & & & =\frac{{\sigma }_{Y1}}{\sqrt{1+r}}\hfill \end{array}$$

Hence, for the veiled observation we obtain

$$\begin{array}{c}\hfill \frac{{\sigma }_{X2}}{X_2}=\sqrt{1+r}\frac{{\sigma }_{X1}}{X_1},\frac{{\sigma }_{Y2}}{Y_2}=\sqrt{1+r}\frac{{\sigma }_{Y1}}{Y_1},\end{array}$$

and

$$\begin{array}{c}\hfill \frac{{\sigma }_{B2}}{B_{\ell 2}}=\sqrt{1+r}\frac{{\sigma }_{B1}}{B_{\ell 1}}.\end{array}$$

Therefore, the uncertainty on B_ℓ increases by a factor $\sqrt{1+r}$ compared to the unveiled case.

In summary, when we know the unveiled magnitude of a star, we can set up the exposure time to get sensitive measurements of the magnetic field, and in that case, the uncertainty on B_ℓ is unaffected by veiling if we don’t change the exposure time, whatever is the flux magnitude during the observation. On the contrary, if the H- or K-band magnitudes of a veiled star are used to set up the exposure time, then the uncertainty on B_ℓ is multiplied by a factor $\sqrt{1+r}$ compared to the unveiled star. Anyhow, the measured values of B_ℓ remain unchanged in both cases.

Appendix D: Log of observations

Table. D.1 summaries the log of observations obtained with SPIRou spectropolarimeter.

Appendix E: LSD Profiles and parameter values

The LSD profiles of each star are shown in Fig. E.1. Table. E.1 summaries for each observation, the corresponding values for the depth I, v sin i, v_rad, B_ℓ, and if the magnetic field is detected or not.

Fig. E.1. LSD profiles for each star. The different colours indicate the different observations. The Y-axes are I/Ic, N/Ic, and V/Ic, where Ic is the intensity of the stellar continuum. The I profile represents the total intensity, the Stokes V profile corresponds to the circular polarisation, and the Stokes N profile corresponds to the null spectra. A pixel size of 4.8 instead of 2.4 is applied for stars with v sin i > 35 km s−1.

Fig. E.2. Fig. E.1. continued.

Fig. E.3. Fig. E.1. continued.

All Tables

All Figures

Fig. 1. Upper panel: Order 32 of the spectrum of VV CrA SW without normalisation (solid black line). In blue are the data points selected to fit the continuum using a sigma-clipping of 3 and 2 iterations. The spline fit of order 3 is plotted in red with a number of knots of 2. The two regions with emission lines were manually removed using the interactive tool. Lower panel: Same spectrum after normalisation. The dashed red line shows the continuum normalised to 1.

In the text

Fig. 2. LSD profiles of V806 Tau. Different colours are used for the different observations. The Y-axes are I/Ic, N/Ic, and V/Ic, where Ic is the intensity of the stellar continuum. The I profile represents the total intensity. Stokes V shows a clear Zeeman signature, which corresponds to an actual magnetic detection because the signal is significantly non-null compared to the surrounding noise and the N profile remains flat in the velocity domain.

In the text

Fig. 3. Distribution of projected rotational velocity for each star in the sample as a function of effective temperature. The magnetic stars are represented in red and the non-detected magnetic stars in black. For the magnetic stars, the sizes of the symbols are proportional to the maximum absolute values of the longitudinal magnetic field. For the non-detected magnetic stars, it is proportional to the minimum of the 3σ uncertainty on Bℓ.

In the text

	Fig. 4. Veiling in the K band as a function of the veiling in the H band for all observations of the sample. The size of the symbol is proportional to the uncertainty of the longitudinal magnetic field and each symbol corresponds to a different target. Magnetic stars are represented in red and stars without magnetic detections in black.
In the text

	Fig. 5. Projected rotational velocities v sin i as a function of effective temperatures measured in this work (blue) and in the work of Flores et al. (2024, red). The points connected by the dashed line indicate the common stars in the two samples.
In the text

	Fig. A.1. Example of dipolar magnetic detection probability curve as a function of the dipole strength. The coloured curves represent the probability for each individual observation, while the black curve represents the combined probability for all observations. The 90 per cent threshold is used to determine the field strength upper limit.
In the text

	Fig. B.1. ZeeTurbo fitting process for the first observation of IRAS 04113+2758 N. The observed spectrum is shown in black, and the best-fit model in red. The blue shading indicates the regions used for parameter fitting, and the corresponding absorption lines are labelled.
In the text

Fig. E.1. LSD profiles for each star. The different colours indicate the different observations. The Y-axes are I/Ic, N/Ic, and V/Ic, where Ic is the intensity of the stellar continuum. The I profile represents the total intensity, the Stokes V profile corresponds to the circular polarisation, and the Stokes N profile corresponds to the null spectra. A pixel size of 4.8 instead of 2.4 is applied for stars with v sin i > 35 km s−1.

In the text

	Fig. E.2. Fig. E.1. continued.
In the text

	Fig. E.3. Fig. E.1. continued.
In the text