© ESO 2019

1. Introduction

With an ever improving understanding of stellar magnetism, we now realise that magnetic fields have a paramount impact throughout the entire life of a star. Particularly on the pre main-sequence (PMS), where stars are relatively cool and magnetic fields seem ubiquitous, magnetism has a significant impact on stars themselves, their formation, accretion properties, rotation rate, flares, and wind characteristics among others. It also influences potential orbiting exoplanets, especially their chemical evolution and habitability.

The origin of stellar magnetic fields and their evolution along the PMS, from stars that evolve from a fully-convective regime to a partially radiative and then eventually to a fully radiative regime, is not fully understood. While cool PMS T Tauri stars exhibit ubiquitous magnetic fields, it seems that most stars, that evolve into PMS Herbig Ae/Be stars and eventually to A/B type main-sequence stars, lose their fields at some point (Alecian et al. 2013; Sikora et al. 2019). The 5–10% of stars remaining magnetic on the main sequence display simple, fossil fields that are not maintained by an active dynamo. Studying the evolution of magnetic fields in the T Tauri star regime in more detail can allow us to understand when and how fast this transition occurs, and perhaps identify different populations with distinct magnetic properties. This could improve our understanding of the mechanisms and timescales at play regarding the evolution and survival of the magnetic fields in the later evolutionary stages (Alecian et al. 2019).

There are two direct methods widely used to measure stellar magnetic fields that both rely on the Zeeman effect. The first is the measurement of Zeeman broadening, or splitting of spectral lines in intensity spectra, and the second is the analysis of the time series of polarised spectra using the Zeeman Doppler Imaging (ZDI) method. These two methods complement each other well. On the one hand, the first method provides a robust estimate of the magnetic field strength integrated over the stellar surface, and even takes small-scale magnetic structures into account. However, it provides almost no information on magnetic field topology. On the other hand, ZDI can reveal the topology of the large-scale magnetic field component. Yet, this method tends to miss some small-scale field structures due to the cancellation of opposite polarities, particularly when only circular polarisation is used. This leads to the severe underestimation of the surface field strength. The two methods have usually not been applied to the same datasets as spectropolarimeters mostly operate in the optical regime (up to 1050 nm), whereas the Zeeman broadening effect increases with wavelength and therefore is best studied at longer wavelengths (for example, in the near-infrared at ≈1500 nm or ≈2200 nm). In a few cases it was possible to obtain both types of measurements from the same data (Kochukhov & Lavail 2017; Hill et al. 2019; Kochukhov & Shulyak 2019) by capitalising on the wide wavelength coverage of optical spectropolarimeters, such as the ESPaDOnS (Echelle SpectroPolarimetric Device for the Observation of Stars) instrument (Donati 2003) at the Canada–France–Hawaii Telescope. The possibilities for self-consistent magnetic analyses will be greatly improved with the advent of a new generation of near-infrared spectropolarimeters such as SPIRou (SpectroPolarimètre Infra-Rouge) at the Canada–France–Hawaii Telescope and CRIRES+ (the upgraded CRyogenic high-resolution InfraRed Echelle Spectrograph) at the European Southern Observatory (ESO) Very Large Telescope (Artigau et al. 2014; Dorn et al. 2016).

A substantial sample of T Tauri stars has been studied using Zeeman broadening (Basri et al. 1992; Guenther et al. 1999; Yang et al. 2005, 2008; Johns-Krull 2007; Yang & Johns-Krull 2011; Lavail et al. 2017; Sokal et al. 2018). A once smaller, but rapidly expanding sample of PMS stars has had its range of large-scale fields characterised with spectropolarimetry and mapped using ZDI (Villebrun et al. 2019; Nicholson et al. 2018; Hill et al. 2019, and references therein). It is worthwhile to note that all these ZDI studies only used Stokes IV parameters.

In this work we aim to expand the sample of T Tauri stars with a magnetic field characterised from high-resolution near-infrared spectroscopy. Our sample includes stars that have been studied with the ZDI technique, allowing us to compare, in detail, the results obtained with the two methods, which are likely to be sensitive to magnetic fields at different spatial scales. Additionally, all but one star were observed repeatedly, enabling us to characterise rotational variability of the magnetic field responsible for Zeeman broadening. We also advance the methodology of magnetic measurements by testing different prescriptions of the distribution of magnetic field strengths over the stellar surface, including a data-driven approach which lets the information contained in the observed spectra determine the appropriate parameterisation of the magnetic field strength distribution.

This paper is structured as follows. The observations and the data reduction process are described in Sect. 2, Sect. 3 explains the methods of our analysis, we present our results in Sect. 4 and finally discuss them in Sect. 5

2. Observations and data reduction

We acquired high-resolution spectra with the CRIRES spectrograph (Käufl et al. 2004) mounted at one of the Nasmyth foci of the 8.2 m Very Large Telescope Unit Telescope 1 at the ESO Paranal observatory. Our observations were obtained between mid-2008 and mid-2009, and were performed in nodding mode in order to improve the sky subtraction. The log of our observations is presented in Table 1. A 0.2″-wide slit yielded a resolving power of R ≡ λ/Δλ ≈ 10⁵. The spectra were recorded on a mosaic of four Aladdin III InSb detectors amounting to an effective 4096 × 512 pixels focal plane detector array. The data were reduced with the standard ESO CRIRES pipeline. We then used the MOLECFIT software (Smette et al. 2015) to remove the telluric lines from the reduced spectra, and simultaneously improve the wavelength solution obtained by the CRIRES pipeline. MOLECFIT retrieves the atmospheric and telescope conditions at the time and location of the observations, and fits the telluric lines in the observed spectra by performing a high-resolution radiative transfer modelling of the terrestrial atmosphere along the line of sight. An example of telluric removal is illustrated in Fig. A.1. Recent work by Ulmer-Moll et al. (2019) demonstrated that MOLECFIT is an excellent tool for removing telluric lines from CRIRES spectra and that it performs better than the classical approach of using a telluric standard star. In our correction, each of the CRIRES detectors was treated independently. In the final reduction step, the spectra from each detector were normalised using a custom routine, which fitted continuum level with a third order polynomial.

The stars in our sample were observed in either one or two CRIRES standard wavelength settings. All targets were observed in the first wavelength setting 25/−1/i at the reference wavelength 2242.8 nm. This setting contains the strong, magnetically sensitive Ti I lines at λ 2223.3, 2227.4, and 2231.1 nm, as well as the strong Na I line at λ 2208.4 nm. We used this set of lines to determine the mean magnetic field modulus ⟨B⟩. All targets but TW Hya were observed in the second wavelength setting 24/−1/n at the reference wavelength 2329.3 nm. This setting covers several lines from the CO band-head, which are magnetically insensitive and therefore are useful for constraining different sources of non-magnetic broadening such as v sin i. The spectra yielding the highest signal-to-noise ratio (S/N) of each star acquired in the wavelength setting 25/−1/i are displayed in Fig. 1. In the particular case of the resolved binary TWA 9 AB, the slit was aligned along the two components, separated by 6″, and the spectrum of the primary TWA 9A was extracted individually.

Fig. 1. Highest S/N normalised spectra in wavelength setting 25/−1/i for each target in our sample. Left panel: strong Na I line present on the first CRIRES detector. Right panel: entire second CRIRES detector covering several magnetically sensitive Ti I lines. The spectral lines used in this analysis and their respective effective Landé factors geff are indicated at the top.

Except for V2247 Oph, stars in our samples were observed repeatedly (2–5 times), in order to investigate potential rotational variability. Some observations were acquired during consecutive nights while others were months apart. The observing dates are indicated in Table 1.

3. Magnetic spectrum synthesis

Our analysis relies on a direct comparison between observations and theoretical synthetic spectra. These spectra were calculated using the Synmast code (Kochukhov 2007; Kochukhov et al. 2010), which numerically solves the polarised radiative transfer equation for a given stellar model atmosphere, spectral line list and magnetic field vector under the assumption of local thermodynamical equilibrium (LTE). Synmast incorporates a modern molecular equilibrium solver (Piskunov & Valenti 2017), allowing to treat both hot and cool star atmospheric conditions. The plane-parallel hydrostatic 1D stellar atmospheres employed in our study were taken from the MARCS model atmospheres grid (Gustafsson et al. 2008)¹. The MARCS atmospheres are not magnetic and hence might not be fully realistic models of magnetised T Tauri stars atmospheres where the single surface temperature assumption might be invalid. However, as we lack appropriate theoretical models (such as 3D magnetoconvection atmospheric simulations for classical T Tauri stars), it is extremely challenging to model the impact of the magnetic field in this context. Nevertheless, the impact of multiple temperatures on the stellar surface would be rather mild on our final results as our analysis is quite insensitive to changes in effective temperature (as detailed in Sect. 4).

For the purpose of modelling the Zeeman broadening, we divided the visible stellar hemisphere in seven annular regions of equal area and computed the local intensity spectra for each region. These local spectra were then broadened, disk-integrated, and normalised to the continuum. Here we generated a grid of synthetic spectra for a purely radial, homogeneous magnetic field, ranging from 0 to 14 kG with 1 kG steps. These synthetic spectra were later linearly combined to simulate more complex surface magnetic field strength distributions, following different prescriptions discussed in Sect. 3.2. Calculations adopted solar abundances according to Asplund et al. (2009).

The line list used in this study was extracted from the VALD3 database (Ryabchikova et al. 2015)². We corrected this list in the following two occurrences where data were missing or produced a bad fit to the observed spectra. One of the four spectral lines used to determine the magnetic field modulus, the Na I line at λ 2208.4 nm, has no van der Waals broadening coefficient in the VALD3 database. In this case Synmast defaults to the van der Waals broadening treatment according to Gray (1992), which appears to underestimate the broadening. We therefore adopted an extended broadening coefficient of 1300.210 calculated by Barklem (priv. comm.) using the ABO theory (Anstee & O’Mara 1995; Barklem & O’Mara 1997). This changes significantly the shape of the Na I line and produces substantially broader wings. Furthemore, we noticed that the synthetic spectra showed systematically stronger blending lines in the blue wing of the Ti I line at λ 2227.4 nm when compared to observations. This was due to two OH transitions around λ 2227.2 nm. We reduced the oscillator strengths of these two transitions by 7 and 5 dex respectively to reach a good fit to observations.

3.1. Stellar parameters

The stellar parameters that are necessary inputs for our analysis are the effective temperature T_eff, the surface gravity log g, the projected rotational velocity v sin i, and both the micro- and macrotubulent velocities, v_mic and v_mac. Effective temperatures were adopted from the literature according to the information in Table 2. Surface gravities (rounded to 0.5 dex) were determined from the position of stars on the HR diagram, as illustrated in Fig. B.1. We then adopted a MARCS model atmosphere with the closest T_eff and log g values. For stars that we observed in the wavelength setting 24/−1/n, we determined v sin i through Markov chain Monte Carlo (MCMC) fitting of the magnetically insensitive CO lines present in that wavelength region. The free parameters that were allowed to vary were v sin i and a scaling factor k. The observed spectra and synthetic calculations with the best-fit v sin i values are shown in Fig. 2. The typical errors on v sin i introduced by changes in the stellar parameters are: ±0.1 km s⁻¹ for T_eff changes of ±200 K, and ±0.2 km s⁻¹ for log g variations of ±0.5 dex. For TW Hya, we adopted a v sin i value from Donati et al. (2011a). We assumed fixed values of v_mic = 2 km s⁻¹ and v_mac = 0 km s⁻¹ following Lavail et al. (2017). The stellar parameters for all targets are listed in Table 2.

Fig. 2. Normalised CRIRES spectra (black histogram) in the wavelength setting 24/−1/n covering the magnetically insensitive CO lines used to measure v sin i (marked with red arrows). The corresponding synthetic spectra computed with the best-fitting v sin i values are overplotted in orange.

3.2. Magnetic field analysis

To determine the magnetic field modulus ⟨B⟩, we fitted magnetically sensitive lines in the observed spectra with magnetic synthetic spectra. These synthetic spectra were generated using different prescriptions for the distribution of surface magnetic field strength. We started by adopting the same three prescriptions as in Yang et al. (2008), also following their model nomenclature (Models 1, 2, and 3).

In Model 1, a uniform magnetic field B covers a fraction of the stellar surface. The rest of the surface is field-free. The free parameters of this model include the magnetic field strength B in the grid from 0 to 14 kG and the magnetic filling factor f. The magnetic field modulus is given by ⟨B⟩ = Bf. In Model 2, the stellar surface is split into three regions: one non-magnetic region and two magnetic regions with variable field strengths B₁ and B₂ (within the 0–14 kG grid) and corresponding filling factors f₁ and f₂. The free parameters of this model are the two magnetic filling factors and the two field strength values. The magnetic field modulus is calculated as ⟨B⟩ = B₁f₁ + B₂f₂. Model 3 postulates three magnetic regions on the stellar surface with fixed magnetic field strengths B₁, B₂, B₃ set respectively to 2, 4, and 6 kG. The remainder of the surface is non-magnetic. The free parameters of this model are the three magnetic filling factors f₁, f₂, f₃. The resulting magnetic field modulus is given by $⟨B⟩=\sum _{i=1}^3{B}_i{f}_i$.

Furthemore, we found that, given the high quality of our observations, proper fitting of spectral lines of some of the stars in our sample required fields stronger than 6 kG. Consequently, we introduced a new prescription (Model 4), representing a generalisation of Model 3 to include stronger field components. Model 4 allows up to seven magnetic regions on the stellar surface with field strengths B₁, B₂, B₃, B₄, B₅, B₆, and B₇ set to 2, 4, 6, 8, 10, 12, and 14 kG while the rest of the surface is non-magnetic. The number of magnetic components is selected via an information criterion within the framework described in Sect. 3.4. The free parameters of this model is the set of magnetic filling factors f_i; the magnetic field modulus is given by $⟨B⟩=\sum _{i=1}^7{B}_i{f}_i$.

The spectra of T Tauri stars are known to be affected by veiling, which often significantly reduces the residual line depths. In addition, some of the stars might exhibit deviations from the solar chemical composition assumed in our calculations. It is impossible to disentangle these effects based on the short wavelength intervals observed with CRIRES. For this reason, we introduced a single parameter k (scaling factor) that accounts for both the true veiling and other deviations in residual line depths (due to e.g. chemical composition, continuum placement). The scaling factor k modifies the synthetic spectra s according to the relation s_scaled = s ⋅ k + 1 − k. This parameter was optimised together with the magnetic parameters for each of the four field parameterisation models described above.

In practice, our analysis consisted of: firstly calculating a library of synthetic spectra using a given set of stellar parameters and a grid of magnetic field strengths with the Synmast code, secondly convolving these spectra for appropriate non-magnetic broadening (v sin i, instrumental broadening corresponding to R = 10⁵) using the s3div code³, thirdly combining the individual magnetic synthetic spectra into a composite synthetic spectrum following one of the four models described above, fourthly fitting the observed spectra by performing a full-grid search (Models 1–2) or a MCMC search (Models 3–4) in the free-parameter space. The scripts and example data used to perform the analysis are available on GitHub⁴.

3.3. Models 1–2: full-grid search

We generated a grid of synthetic magnetic field spectra for the set of stellar parameters (T_eff, v sin i, and log g) of our sample, with the magnetic field strength B ranging from 0 to 14 kG in increments of 1 kG. For each of the two magnetic field prescriptions, we then conducted a search for a set of free parameters which minimised the ${\chi }_r^2$ of the fit between observed and synthetic spectra. The scaling factor k was allowed to range from 0 to 2 in increments of 10⁻², the magnetic filling factors varied from 0 to 1 in increments of 2 × 10⁻² and 4 × 10⁻² for Models 1 and 2, respectively. The non-magnetic filling factor was not a free parameter, as it is the remainder of the surface that is not covered with a magnetic field. This translated to an additional constraint that the sum of all filling factors must be equal to unity.

3.4. Model 3–4: MCMC optimisation

Our MCMC analysis made use of the Goodman & Weare’s affine invariant ensemble (AIE) sampler (Goodman & Weare 2010) in its implementation in the Python package emcee (Foreman-Mackey et al. 2013). The free parameters of Models 3 and 4 were the magnetic filling factors f_i of the magnetic regions and the scaling factor k. We adopted a uniform prior U(0, 1) for each of the individual magnetic filling factors f_i and an additional prior U(0, 1) for the sum $\sum _{i=1}^N{f}_i$. This prior means that the fraction of the stellar surface covered by any magnetic field must be less than or equal to unity, the remainder was set non-magnetic. A uniform prior U(0, 2) was adopted for the scaling factor. We generated a set of randomised filling-factor combinations as an initial value for each Markov chain. The initial value of the scaling factor was set to k = 0.8 with an added random noise. A simple Gaussian likelihood was used to measure the deviation of the model to the observed data adopting a standard deviation σ = S/N⁻¹.

The MCMC AIE sampler was run with 60 chains, each initially 1000 steps long. We discarded the first 800 “burn-in” samples after ensuring that each chain had reached a stable state well within the burn-in length. Formally, the MCMC was iteratively run with 1000 additional steps until the chains for each parameter reached a minimum of 2500 effective sample size (ESS) defined as

$$\begin{array}{c}\hfill \text{ESS}=\frac{N}{2{\tau }_{\text{int}}},\end{array}$$

where N is the total number of samples in the chains, and τ_int is the integrated autocorrelation (Sharma 2017, Sect. 3.7.1). The resulting jump acceptance fractions were comprised between 44 and 65% (median: 53%).

For each observed spectrum, we ran the MCMC analysis as described above multiple times using field strength parameterisation models of increasing complexity. We started with the simplest model comprised of only two components: 0 and 2 kG. We iteratively added more components: 4, 6, 8, 10, 12, and ultimately 14 kG. More complex models generally provide a better fit to observations, but can also overestimate the magnetic field modulus by introducing a small fraction of very strong field that effectively fits the noise or systematics in the spectra. To avoid this issue, we selected the model complexity that is favoured by the data by choosing the model that yields the smallest Bayesian Information Criterion (BIC). The BIC is a model-selection tool that estimates the relative quality of a model applied for fitting a given dataset with regards to the number of model parameters. The criterion is defined as

$$\begin{array}{c}\hfill \mathrm{BIC}\equiv -2\text{ln}{\mathcal{L}}_{\text{max}}+2n\text{ln}N,\end{array}$$

where ℒ_max is the maximum likelihood value, n is the number of model parameters, and N is the number of data points (Sharma 2017).

4. Results

We obtained a set of four magnetic field modulus ⟨B⟩ values for each observation, corresponding to each of the four magnetic models. These results are presented in Table 3, which summarises magnetic measurements together with the corresponding reduced ${\chi }_r^2$ values. The simplest Model 1 fails to reproduce the observed spectra satisfactorily in many cases yielding high ${\chi }_r^2$ values. For stars with moderate field strengths (≲3 kG), Models 2 and 3 produce an acceptable fit to the data, with ${\chi }_r^2$ values very close to those obtained with Model 4. However, these parameterisations lead to larger ${\chi }_r^2$ for stars with stronger fields. Model 4 appears to provide a reliable fit to all spectra in our sample, up to high magnetic field strengths (> 4 kG), without overfitting spectra of stars with weak apparent magnetic field (e.g. V2058 Oph). For this reason, we judge Model 4 to be the most reliable and will use results obtained with this model in the remainder of the paper.

A typical parameter inference result from the application of Model 4 is displayed in Fig. 3, where the top subplot shows a corner plot with the posterior probability distributions for the model parameters, and the bottom subplot shows the posterior probability distribution of the surface magnetic field ⟨B⟩. Typical 99.7% confidence intervals for ⟨B⟩ obtained from the MCMC analysis have a width of 0.3 kG. However, these statistical uncertainties can be smaller than the systematic uncertainties that would arise from errors of the adopted atomic data, such as the Van der Waals broadening parameters. The typical errors on the magnetic field modulus introduced by errors on stellar parameters are: ±0.02 kG for v sin i variations of ±0.5 km s⁻¹, ±0.1 kG for T_eff changes of ±200 K, and ±0.3 kG for surface gravity errors of ±0.5 dex.

Fig. 3. Illustration of the MCMC analysis for the spectrum of TW Hya obtained on 2009-02-12 using Model 4. The corner plot at the top shows the posterior distribution of the model parameters (magnetic filling factors fi and the scaling factor k). The bottom histogram displays the resulting posterior distribution of the magnetic field modulus ⟨B⟩, overplotted with its median value (yellow thick line) and the limits of the 3σ confidence interval (red dashed lines).

The best fits obtained with Model 4 to the highest-S/N observations of each target are illustrated in Fig. 4. For comparison, this figure also shows non-magnetic synthetic spectra calculated with the same stellar parameters. It is evident that, in all cases, magnetic field has a large impact on the studied spectral lines.

Fig. 4. Highest-S/N observation of each target around the four magnetically-sensitive spectral lines used in this analysis (black histogram), overplotted with the best non-magnetic fit (grey dashed line) and the best magnetic fit with Model 4 (orange line).

By obtaining repeated observations and magnetic field measurements for all stars (with the exception of V2247 Oph), we could assess the degree of rotational variability of the magnetic field modulus ⟨B⟩. This analysis establishes a relatively low level of night-to-night changes of the mean field strength. Figure 5 displays individual ⟨B⟩ measurements and the average value for each star. The maximum deviation of an individual measurement from the average is 0.4 kG (TWA 9A). The mean deviation is 0.1 kG. The maximum peak-to-peak variation of individual ⟨B⟩ measurements for a given star is 0.7 kG (for TWA 9A, which amounts to a 21% variation relative to the mean value of ⟨B⟩ = 3.4 kG). The mean peak-to-peak variation is 0.3 kG. The absence of dramatic magnetic field changes with the stellar rotation indicates that magnetic regions responsible for the Zeeman broadening signal are not concentrated in discrete patches but distributed approximately homogeneously over the stellar surface.

Fig. 5. Individual ⟨B⟩ measurements for each target from Model 4 (filled orange circles with 3σ error bars) and average ⟨B⟩ value for each star (black diamonds).

5. Discussion

In order to put our results in a wider context, we present our mean magnetic field measurement for each star in the sample on the Hertzsprung–Russell (H–R) diagram alongside literature measurements of Zeeman broadening from near-infrared spectra of other low mass and intermediate mass T Tauri stars, as shown in Fig. 6. The mass tracks, the zero-age main sequence (ZAMS), and the fully-convective transition boundary displayed in the H–R diagram are all extracted from the Yale–Potsdam Stellar Isochrones grid (Spada et al. 2017), which is available online⁵. It is worthwile to note that this stellar evolution model does not include the effect of stellar magnetic fields, which are making a PMS star appear cooler and therefore potentially younger. Incorporating magnetic fields into stellar evolution models is an ongoing effort in order to improve their accuracy (e.g. Feiden 2016). To place the stars on the H-R diagram, we used effective temperatures from the literature, and derived consistent stellar luminosities following a methodology adapted from Villebrun et al. (2019). We obtained J-band magnitudes m_J from 2MASS (Cutri et al. 2003), V-band magnitudes m_V from Bouvier & Appenzeller (1992) and Torres et al. (2006, and references therein), and parallaxes π from Gaia DR2 (Gaia Collaboration 2016, 2018). Then, we obtained the stellar luminosities through the following set of equations:

$$\begin{array}{cc}& E(V-J)=(m_V-m_J)-{(V-J)}_0,\hfill \end{array}$$

$$\begin{array}{cc}& A_J=R_J\times E(V-J),\hfill \end{array}$$

$$\begin{array}{cc}& M_J=m_J+5+5{\text{log}}_{10}(\pi )-A_J,\hfill \end{array}$$

$$\begin{array}{cc}& M_{\mathrm{bol}}=M_j+{\mathrm{BC}}_J,\hfill \end{array}$$

$$\begin{array}{cc}& \frac{L}{L_{\odot }}={10}^{-(M_{\mathrm{bol}}-M_{\mathrm{bol},\odot })/2.5},\hfill \end{array}$$

Fig. 6. HR diagram summarising near-infrared spectroscopic magnetic field measurements of T Tauri stars. The filled grey circles represent the intermediate-mass T Tauri stars sample from Lavail et al. (2017). The two shades of grey represent results obtained for two different sets of spectral lines. The open dark circles show the samples from Johns-Krull (2007), Yang et al. (2008), and Yang & Johns-Krull (2011), and the open red circles display the average measurement for the sample presented in this work. The radius of the circles is proportional to ⟨B⟩2. The ZAMS is plotted with dashed blue line, the limit between fully and partially convective regimes is shown with the dash-dotted yellow curve, and evolutionary tracks for different masses are plotted with thin grey lines.

where (V − J)₀ is the empirical intrinsic colour and E(V − J) the colour excess, A_J the extinction in the J-band, R_J the total to selective extinction set to 0.437 following Villebrun et al. (2019), π the parallax in arcseconds, and BC_J the bolometric correction in the J-band. The empirical intrinsic colours and bolometric corrections for 5–30 Myr PMS stars were computed according to Pecaut & Mamajek (2013), and we adopted the associated value of M_bol, ⊙ = 4.755. The determined luminosities and intermediate values from the computation are listed in Table 4. It is important to note that the stellar luminosities of these objects can be quite uncertain, due to e.g. their variability, extinction, and the non-contemporaneous V- and J-band measurements.

From examination of Fig. 6 it appears that the magnetic field modulus is not a function of the position in the H-R diagram, at least in the fully convective regime where we have a substantial number of observations showing stars in the same location having quite different magnetic field modulus ⟨B⟩. Across the fully convective boundary the sample of stars is small. The two stars from our sample in this regime (albeit quite close to the boundary) are TWA 9A and TW Hya with mean field modulus of ⟨B⟩ = 3.4 and ⟨B⟩ = 3.5 kG respectively. This is stronger compared to the field modulus reported for the intermediate mass T Tauri stars sample of Lavail et al. (2017), in which four stars are partially convective.

We can also compare the magnetic field modulus ⟨B⟩ from this work, with the average large-scale magnetic field ⟨B_V⟩ obtained from Stokes V ZDI studies. Values for ⟨B_V⟩ are 113 G for TWA 9A (Nicholson et al. 2018), 1628 G for TW Hya (Donati et al. 2011a), ≈300 G for V2247 Oph (Donati et al. 2010), 1424 G for TWA 8A (Hill et al. 2019), 56 G for TWA 7 (Nicholson et al., in prep.), and we computed ⟨B_V⟩≈1020 G from the magnetic field configuration described in Donati et al. (2011b) for V2129 Oph.

On the one hand, a relatively large fraction (around, respectively, 42% and 25% of the mean field strength) of the fields of TW Hya and V2129 Oph are recovered from Stokes V. The ZDI mapping indicates that these fields are relatively strong and simple, being mostly poloidal, axisymmetric, and dominated by octupole and dipole components (Donati et al. 2011a,b). A large fraction (34%) of the field of TWA 8A is also recovered from Stokes V, with the ZDI mapping also indicating a mostly poloidal field, with most of the poloidal component being axisymmetric (Hill et al. 2019). These values are significantly higher than those reported for M dwarfs: less than 25% and typically in the 5–20% range (Reiners & Basri 2009; Reiners 2012; Kochukhov & Shulyak 2019). On the other hand, a smaller fraction of the field is recovered by the Stokes V analyses for V2247 Oph (11%), whose field is weaker, more complex, and is comprised mostly of a toroidal and a poloidal non-axisymmetric components (Donati et al. 2010), for TWA 7 (2%) which exhibits a predominantly toroidal and non-axisymmetric field (Nicholson et al., in prep.), and for TWA 9A (3%), which also exhibits a weak and mostly non-axisymmetric large-scale field, but with a dominant poloidal component (Nicholson et al. 2018). This assessment seems to corroborate the hypothesis that simple strong axisymmetric fields are better recovered by Stokes V ZDI mapping, but complex fields with structures at smaller scales tend

to cancel out and are missed by ZDI, at least when only Stokes IV parameters are used. This is illustrated in Fig. 7. Naturally, this comparison is somewhat limited by the fact that spectropolarimetric observations used for the tomographic mapping, and spectroscopic observations for this work have been acquired with different instruments and at different epochs. A similar exercise was recently performed by See et al. (2019) who collected ⟨B⟩ (noted ⟨B_I⟩ in their publication) and ⟨B_V⟩ values for 83 F-G-K-M stars from multiple literature sources, and derived empirical formulae linking ⟨B⟩ and ⟨B_V⟩ through power laws. Our ⟨B⟩ measurements do not seem to follow the trends with ⟨B_V⟩ predicted by either of the two empirical power laws obtained by See et al. (2019) (which generally predict too strong fields with ⟨B⟩ values as strong as 10 or 14 kG for our sample). Although the discrepancy could possibly arise from the variety of methods used to determine ⟨B⟩ values in the literature and different constraints in the ZDI mapping used to recover ⟨B_V⟩, it is not at all guaranteed that T Tauri stars should follow the same trend as established using main sequence stars.

Fig. 7. Fraction of the mean magnetic field recovered by ZDI ⟨BV⟩ over ⟨B⟩ as a function of the axisymmetric large-scale magnetic field component recovered with ZDI (magnetic field components in spherical harmonic modes with m <  l/2). The radius of the symbols is proportional to ⟨B⟩2.

It will be worthwhile to repeat simultaneous determination of ⟨B⟩ and ⟨B_V⟩ with data from upcoming instruments such as SPIRou or CRIRES+ which will allow one to use both magnetic analysis methods simultaneously on the same dataset. The large spectral grasp of these instruments will also provide the opportunity to analyse molecular features that predominantly form at the cooler spot temperatures so the same dataset can be used to measure both the spot filling factors (Gully-Santiago et al. 2017) and the magnetic filling factors, another clue towards a better understanding of how the extremely strong field is concentrated at the smallest scales. This will enable us to obtain a complete picture of young-star magnetism reconciling our understanding of both small-scale and large-scale magnetic components.

Acknowledgments

A. L. would like to thank Colin Hill and Colin Folsom for partaking in useful discussions, and Martin Sahlén for helping with the MCMC analysis. A. L is also grateful to Peter K. G. Williams and the contributors to the “Recommendations for Reporting MCMC Analyses in Academic Literature” (https://github.com/pkgw/mcmc-reporting) which provided useful guidelines. O.K. acknowledges support by the Knut and Alice Wallenberg Foundation (project grant “The New Milky Way”), the Swedish Research Council (project 621-2014-5720), and the Swedish National Space Board (projects 185/14, 137/17). Finally, the authors would like to thank the anonymous referee for their report which improved the quality of this work. Based on observations collected at the European Southern Observatory under ESO programme 081.C-0555. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made use of NASA’s Astrophysics Data System Bibliographic Services. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.

References

Appendix A: Example of telluric removal with MOLECFIT

Fig. A.1. Example of telluric line removal performed with the MOLECFIT software on the spectrum of V2058 Oph acquired on 2008-07-10.

Appendix B: HR diagram with surface gravity tracks

Fig. B.1. HR diagram used for surface gravity determination. The red diamonds represent the sample of T Tauri stars presented in this paper. The ZAMS is plotted with dashed blue line, the limit between fully and partially convective regimes is shown with the dash-dotted yellow curve. Surface gravity tracks are traced with dotted lines in shades of blue.

Appendix C: Parameter inference with Model 4

All Tables

All Figures

	Fig. 1. Highest S/N normalised spectra in wavelength setting 25/−1/i for each target in our sample. Left panel: strong Na I line present on the first CRIRES detector. Right panel: entire second CRIRES detector covering several magnetically sensitive Ti I lines. The spectral lines used in this analysis and their respective effective Landé factors geff are indicated at the top.
In the text

	Fig. 2. Normalised CRIRES spectra (black histogram) in the wavelength setting 24/−1/n covering the magnetically insensitive CO lines used to measure v sin i (marked with red arrows). The corresponding synthetic spectra computed with the best-fitting v sin i values are overplotted in orange.
In the text

Fig. 3. Illustration of the MCMC analysis for the spectrum of TW Hya obtained on 2009-02-12 using Model 4. The corner plot at the top shows the posterior distribution of the model parameters (magnetic filling factors fi and the scaling factor k). The bottom histogram displays the resulting posterior distribution of the magnetic field modulus ⟨B⟩, overplotted with its median value (yellow thick line) and the limits of the 3σ confidence interval (red dashed lines).

In the text

	Fig. 4. Highest-S/N observation of each target around the four magnetically-sensitive spectral lines used in this analysis (black histogram), overplotted with the best non-magnetic fit (grey dashed line) and the best magnetic fit with Model 4 (orange line).
In the text

	Fig. 5. Individual ⟨B⟩ measurements for each target from Model 4 (filled orange circles with 3σ error bars) and average ⟨B⟩ value for each star (black diamonds).
In the text

Fig. 6. HR diagram summarising near-infrared spectroscopic magnetic field measurements of T Tauri stars. The filled grey circles represent the intermediate-mass T Tauri stars sample from Lavail et al. (2017). The two shades of grey represent results obtained for two different sets of spectral lines. The open dark circles show the samples from Johns-Krull (2007), Yang et al. (2008), and Yang & Johns-Krull (2011), and the open red circles display the average measurement for the sample presented in this work. The radius of the circles is proportional to ⟨B⟩2. The ZAMS is plotted with dashed blue line, the limit between fully and partially convective regimes is shown with the dash-dotted yellow curve, and evolutionary tracks for different masses are plotted with thin grey lines.

In the text

	Fig. 7. Fraction of the mean magnetic field recovered by ZDI ⟨BV⟩ over ⟨B⟩ as a function of the axisymmetric large-scale magnetic field component recovered with ZDI (magnetic field components in spherical harmonic modes with m <  l/2). The radius of the symbols is proportional to ⟨B⟩2.
In the text

	Fig. A.1. Example of telluric line removal performed with the MOLECFIT software on the spectrum of V2058 Oph acquired on 2008-07-10.
In the text

	Fig. B.1. HR diagram used for surface gravity determination. The red diamonds represent the sample of T Tauri stars presented in this paper. The ZAMS is plotted with dashed blue line, the limit between fully and partially convective regimes is shown with the dash-dotted yellow curve. Surface gravity tracks are traced with dotted lines in shades of blue.
In the text