© The Authors 2024

1. Introduction

Young stars are surrounded by circumstellar disks in which planets form (e.g., Williams & Cieza 2011). These disks dissipate in a few million years (e.g., Fedele et al. 2010). The processes governing their evolution and dispersal, and their impact on planet formation remain open questions (e.g., Manara et al. 2023, for a review). The first proposed mechanism involves the gradual spread of the disk due to viscosity, resulting in most of the material accreting onto the central star, with a small fraction of material carrying the angular momentum outward (Lynden-Bell & Pringle 1974; Hartmann et al. 1998). An alternative scenario sees the evolution dominated by the ejection of angular momentum from the disk in magnetohydrodynamic (MHD) winds (Blandford & Payne 1982; Lesur 2021b,a; Tabone et al. 2022). At the same time, high-energy radiation from the star, such as far-ultraviolet (FUV), extreme ultraviolet (EUV), or X-ray radiation, removes material from the inner disk in the form of photoevaporative winds. This process may potentially create a gap in the inner disk and ultimately lead to the dissipation of the disk (e.g., Clarke et al. 2001; Alexander et al. 2014; Ercolano & Pascucci 2017). Understanding the impact of viscosity, MHD winds, and internal photoevaporation is critical for our understanding planet formation. Indeed, while these processes are ongoing, the dust grains present in the disks are expected to grow from small grains into planetesimals (e.g., Johansen et al. 2007; Booth & Clarke 2016; Carrera et al. 2017, 2021). These planetesimals can further grow into planetary systems (e.g., Alibert et al. 2005; Liu et al. 2019; Lyra et al. 2023), or may be revealed in the form of collisionally active debris disks once the gas in the disk has dissipated (e.g., Wyatt 2008; Hughes et al. 2018; Lovell et al. 2021).

Each of these disk evolution mechanisms makes specific predictions about the time evolution of accretion rates through the disk, and thus on the central star. Determining the mass accretion rate of young stars is therefore a key parameter to shed light on whether viscous or MHD-wind driven evolution dominates disk evolution and which role photevaporation plays in dissipating the disk (e.g., Mulders et al. 2017; Lodato et al. 2017; Somigliana et al. 2020, 2023).

Numerous studies attempted to constrain disk evolution by analyzing theoretically expected and observed relationships between the mass accretion rate and stellar or disk properties, such as the disk and stellar mass (e.g., Dullemond et al. 2006; Hartmann et al. 2006; Alexander et al. 2006; Ercolano et al. 2014; Mulders et al. 2017; Rosotti et al. 2017; Manara et al. 2020). Currently no definitive conclusions have been drawn in favor of either model.

Recently, Alexander et al. (2023) showed that the distribution of ≳300 homogeneous accretion rate measurements of young stars in the 0.5 − 1.0 M_⊙ range should suffice to distinguish between a viscous+photoevaporation and a MHD-wind driven disk evolution scenario. Currently, only ∼100 accretion rate measurements are available in this mass range, and this sample lacks homogeneous upper limits for the non-detections. Providing such upper limits could enable even stronger statistical approaches to further constrain model parameters. Ercolano et al. (2023) used the accretion rate distribution of a sample of potential extreme low accretors presented by Thanathibodee et al. (2022, 2023) to discern between different types of photoevaporation that could dissipate the disk at the end of a viscously driven evolution scenario. This sample was limited as it only contained 24 sources. A thorough characterization of more low accretors can therefore provide valuable constraints on disk evolution processes.

Tracers of the accretion and outflow processes can be found in the spectra of young stars. The current paradigm for accretion in young, low-mass (≲1 M_⊙) stars is known as magnetospheric accretion (Hartmann et al. 2016, for a review). In this scenario, accreting material is funneled by the star’s magnetosphere onto its surface, creating hotspots that emit excess continuum radiation in the ultraviolet and optical spectra (Calvet & Gullbring 1998). The hot gas associated with the accretion streamers leads to the appearance of various prominent emission lines (e.g., Hα, Hβ Ca K, and Paβ) in the spectra (e.g., Muzerolle et al. 2001; Campbell-White et al. 2021).

The spectroscopic study of young stars can, therefore, provide vital constraints on the physics governing disk evolution. Accurately representing the emission from the central star is crucial in these investigations. Several methods are employed to estimate accretion rates in young stars. As a first option one can fit the line profile using magnetospheric accretion flow models (Muzerolle et al. 2001; Thanathibodee et al. 2023). The downside of these models is that they assume the stellar magnetic field to be both dipolar, aligned with the stellar rotation axis and disk, which is not necessarily accurate (e.g., Donati et al. 2008; Alencar et al. 2012; Singh et al. 2024). In a second method, the line luminosity of various emission lines is converted into an accretion luminosity using empirical relations (e.g., Herczeg & Hillenbrand 2008; Alcalá et al. 2014, 2017). These relations are calibrated using the third method which is the only one that directly probes the energy released at the accretion shock. This method consists of measuring the UV continuum excess emission (e.g., Valenti et al. 1993; Calvet et al. 2004; Herczeg & Hillenbrand 2008; Rigliaco et al. 2012; Ingleby et al. 2013; Pittman et al. 2022; Robinson et al. 2022). In order to self-consistently determine both the stellar and accretion parameters, the UV-excess and the spectral features in the optical part of medium-resolution flux-calibrated spectra must be fit simultaneously (Manara et al. 2013a). For this aim, both the photospheric and chromospheric emission in the UV part of the spectrum must be correctly accounted for (Ingleby et al. 2011a; Manara et al. 2013b).

Available synthetic spectra (e.g., Allard et al. 2011) do not fully reproduce the observed spectra of chromospherically active low gravity objects, such as pre-main sequence low-mass stars. In particular, none of the current models contain emission originating in the chromospheres of young stars, since only the stellar photospheres are modeled, which dominate only at wavelengths longer than the Balmer Jump. At shorter wavelengths, bound free and free free emission originating in the chromosphere starts to dominate the continuum (Houdebine et al. 1996; Franchini et al. 1998). The chromosphere also emits in several spectral lines that are used to constrain accretion properties of accreting young stellar objects (YSOs, Stelzer et al. 2013). The spectra of field dwarfs will also poorly represent those of young stars as the former have a significantly higher log g. The best solution appears to use spectra of non-accreting young stars. Class III stars (Greene et al. 1994), which are defined according to their infrared classification where they display a lack of infrared excess emission dlog(νF_ν)/dlog(ν) <  − 1.6 (Williams & Cieza 2011), make ideal candidates for such templates, as this category often overlaps with Weak lined T Tauri Stars (WTTS), which present much fainter emission lines than CTTS.

A grid of class III templates with broad wavelength coverage medium-resolution spectra obtained with the X-shooter spectrograph on the ESO Very Large Telescope (VLT) was previously provided to the community by Manara et al. (2013b) (hereafter MTR13) and further expanded by Manara et al. (2017a) (hereafter MFA17). This grid includes 41 spectra and contains spectral types ranging from G4 to M8. Both MTR13 and MFA17 used these spectra to estimate the contribution of chromospheric emission to the emission lines that are commonly used for the determination of mass accretion rates. To study the contribution of the chromosphere to the UV continuum emission, Ingleby et al. (2011a) compared the spectrum of the M0 class III RECX-1 to the photospheric spectrum of a standard dwarf star. Here it was found that attributing this excess to accretion would result in an estimated accretion luminosity of log(L_acc/L_⋆) =  − 1.3. The chromospheres of young stars can therefore significantly affect measured accretion rates in low accretors. Moreover, the Balmer continuum excess emission and the Balmer Jump are more difficult to detect in the spectra of early-type (< K3) YSOs than in the later types due to the lower contrast between photospheric emission and accretion induced continuum excess emission (Herczeg & Hillenbrand 2008). Constraints on the influence of the chromospheric emission in the UV on measurements of (low) accretion rates is still lacking. A better understanding of the chromospheres of young stars is therefore needed to characterize the lowest accretors to better understand the late phases of disk dispersal and hence, provide constraints on disk evolution models.

Here we expand the library of X-shooter class III templates previously provided by MTR13 and MFA17 by an additional 18 templates and present a method for interpolating between them. We also aim to constrain the influence of the UV continuum chromospheric emission on determination of mass accretion rates. The paper is structured as follows. The sample selection, observations and data reduction are described in Sect. 2. The analysis of the stellar properties of our sample is described in Sect. 3. In Sect. 4 the final grid and our method for interpolating it are discussed. Sect. 5 discusses an application of this interpolated grid using a self consistent method to derive the mass accretion rates and validate this method on a set spectra of accreting young stars in the Chamaeleon I star-forming region. In Sect. 6.1 we obtain a lower limit on the mass accretion rates that we can measure from the UV excess, taking into account uncertainties on the chromospheric emission and discuss the implications for studies of mass accretion rate (Ṁ_acc) in CTTS. Finally, we summarize our conclusions in Sect. 7.

2. Sample, observations, and data reduction

The new targets considered in this work come mainly from the sample of Lovell et al. (2021), who studied 30 class III stars (age ≲ 10 Myr) using the Atacama Large Millimeter Array (ALMA). Lovell et al. (2021) reaffirmed that these sources are class III YSO through a comparison of the K-band (2MASS) and either 12 μm (WISE) or 24 μm (Spitzer) fluxes. In addition to this Lovell et al. (2021) fitted the SEDs of these targets and found no significant NIR excess emission. These stars are likely members of the Lupus star forming region, although their membership still needs to be confirmed (Lovell et al. 2021; Michel et al. 2021). A total of 19 objects of this sample were recently observed using X-shooter (Pr.ID. 109.23D4.001, 110.23P2.001, PI Manara), a broad-band, medium-resolution, high-sensitivity spectrograph mounted on the ESO/VLT. The two components of a binary system, THA15-36A and THA15-36B, were observed simultaneously in the slit. Three additional archival X-shooter spectra were available for the targets MT Lup and NO Lup (Pr.ID. 093.C-0506(A), PI Caceres) and MV Lup (Pr.ID. 085.C-0764(A), PI Günther). Other targets in the sample of Lovell et al. (2021) were already included in the grid of MTR13 and MFA17. In addition to this, we use in this work two targets observed with VLT/X-shooter in the PENELLOPE Large Program (Pr.ID. 106.20Z8, Manara et al. 2021), namely RECX-6 and RXJ0438.6+1546. RECX-6 was identified as a class III YSO by Sicilia-Aguilar et al. (2009) based on a fit of its SED. RXJ0438.6+1546 was already included in the grid of MFA17, who selected their targets based on available Spitzer data. The latter was already observed with X-shooter and included in the sample of MFA17, but the new observations are used here. In total the targets considered in the analysis in this work are 24, observed in 23 observations.

The wavelength range covered with X-shooter is divided into three arms, the UVB arm (λλ ∼ 300 − 550 nm), the VIS arm (λλ ∼ 500 − 1050 nm), and the NIR arm (λλ ∼ 1000 − 2500 nm). Different slit widths were used for the different arms and for fainter or brighter targets. For brighter sources, we used $1\stackrel{″}{.}0$, $0\stackrel{″}{.}4$, $0\stackrel{″}{.}4$ wide slits in the UVB, VIS and NIR arm respectively. These slit widths provide typical spectral resolutions of R ∼ 5400, 18 400 and 11 600 in the three arms. Fainter, later spectral type objects were observed using $1\stackrel{″}{.}0$, $0\stackrel{″}{.}9$, $0\stackrel{″}{.}9$ in the three arms, providing resolutions of R ∼ 5400, 8900 and 5600. For MV lup, RXJ0438.6+1546, MT Lup and NO Lup, the $0\stackrel{″}{.}5$ wide slits were used in the UVB arm, resulting in spectra with resolution R ∼ 9700. The observations with the slits with the aforementioned widths were preceded by short exposures with a slit with the significantly larger width of $5\stackrel{″}{.}0$ in all arms to obtain spectra not affected by slit losses to be used for absolute flux calibration. The only exception of the use of the wide slits were MV Lup, MT Lup, and NO Lup.

Appendix A contains a log of the observations presented here. In the observing log, several targets are highlighted as spatially unresolved binaries. We exclude these spectra from our sample since they can not be used as templates to represent the stellar emission of individual targets and an analysis of these spectra is outside the scope of this paper. In particular, we exclude CD-39 10292, THA15-38, CD-35 10498, V1097 Sco and NN Lup. NN Lup was discovered to be SB2 binaries in our analysis with the ROTFIT tool (see Sect. 3.2). CD-39 10292 is an SB2 binary first identified by Melo et al. (2001). V1097 Sco and CD-35 10498 were identified as binaries by Zurlo et al. (2021) and were not spatially resolved in the X-shooter slit. THA15-38 was resolved as a binary in the acquisition image of ESPRESSO observations that are part of the PENELLOPE program. This was possible due to the excellent seeing conditions ($\lesssim 0\stackrel{″}{.}6$) during the night of the ESPRESSO observations.

The data reduction was performed using the ESO X-shooter pipeline v.4.2.2 (Modigliani et al. 2010) in the Reflex workflow (Freudling et al. 2013). The pipeline executes the standard reduction steps: flat fielding, bias subtraction, extraction and combination of orders, rectification, wavelength calibration, flux calibration using a standard star observed during the same night and the final extraction of the 1D spectrum. The telluric lines were removed from the narrow slit observations using the molecfit tool (Smette et al. 2015). This was done by fitting the telluric features on the spectra themselves, rather than by using a telluric standard star. As a final step, the narrow slit observations were rescaled to the continuum flux of the wide slit spectra to get flux-calibrated spectra, using the procedure developed for the PENELLOPE program (see Manara et al. 2021). For the UVB arm this was done using a correction factor constant with wavelength. The VIS and NIR arms were rescaled using a factor with a linear dependency on wavelength. To test the flux calibration of our spectra they were compared to archival photometry. The overall agreement between the spectra and photometry is excellent (Δmag < 0.2 mag) with the exception of 2MASSJ16090850−3903430 and THA15-36. This is however within the typical variation range at optical wavelength which is mostly due to starspots rotational modulation. A description of how we flux calibrated targets with Δmag > 0.2 mag can be found in Appendix A. Appendix A also contains a description of how we obtained an accurate flux calibration for the apparent visual binary THA15-36, which was spatially resolved in the X-shooter slit. Zurlo et al. (2021) resolved this system using VLT/NACO, and argued that both components are unbound given their different Gaia distances (146.5 ± 0.8 pc for the primary and 154.3 ± 1.8 pc for the secondary).

3. Stellar parameters of the new class III spectra

In this section, we derive the stellar parameters for all the 19 resolved stars whose spectra are presented here for the first time. This was done by first determining the spectral type from atomic and molecular features present in the spectra, then comparing these estimates with results from the fitting of individual absorption lines, and finally determining the stellar luminosity to place the targets on the Hertzsprung–Russel Diagram (HRD) to derive their stellar masses.

3.1. Spectral type and extinction determination

Determining spectral types (SpT) for (non-accreting) young stellar objects is a complex process that is best performed by comparison with other targets of well known SpT. In this case, we consider the grid of non-accreting class III young stellar object of MTR13 and MFA17 as a starting point, and later also check and refine their grids using the new spectra presented here.

The first estimate of the SpT of our targets was obtained applying to the new spectra the same spectral indices used by MTR13 and MFA17, namely the spectral indices presented by Riddick et al. (2007), the TiO index of Jeffries et al. (2007) and the indices presented by Herczeg & Hillenbrand (2014, hereafter HH14). The indices of Riddick et al. (2007) are accurate for spectral types later than M3. The indices by Jeffries et al. (2007) hold for spectral types later than K6, since the used TiO feature disappears for earlier SpTs. The indices of HH14 are accurate for SpTs from early K to late M stars. The results for these different methods are listed in Table 1. Spectral indices are less reliable estimates of the SpT in cases in which the extinction is substantial since they are based on ratios of features at different wavelength ranges. In general, extinction is low for our targets, but it can still introduce additional uncertainty to the SpT estimates obtained from the spectral indices.

Therefore, to obtain a simultaneous extinction and SpT estimate, we performed a comparison of the spectra in our sample with the spectra of MTR13 and MFA17, which have negligible extinction. We obtained the first estimates of SpT and A_V using a simple χ²-like comparison. The final values of both the spectral type and extinction are derived through a visual comparison between the spectra presented here and those presented by MTR13 and MFA17. The main indicators used in this comparison are the depths of various molecular bands. The K6 to M9.5 spectra include a variety of molecular features whose depth increases almost monotonically for later SpTs in the spectral region between 580 and 900 nm. This includes various absorption bands from TiO (λλ 584.7–605.8, 608–639, 655.1–685.2, 705.3–727, 765–785, 820.6–856.9, 885.9–895 nm), CaH (λλ 675–705 nm) and VO (λλ 735–755, 785–795, 850–865 nm). Figure 1 displays spectra with SpTs ranging from M1 to K7.5. Similar figures for the other spectra with SpT later than K5 can be found on ZENODO. The previously mentioned molecular features can be seen to increase in depth for later spectral types.

Fig. 1. X-shooter spectra of class III YSOs with spectral type ranging from M1 to K7. All the spectra are normalized at 731 nm and offset in the vertical direction for clarity. The spectra are also smoothed to the resolution of 2500 at 750 nm to make easier the identification of the molecular features. The black colors indicate spectra presented by MTR13 and MFA17. The red color is used for spectra presented for the first time here.

For SpT earlier than K6 a MgH and Mg b absorption feature (λλ 505–515 nm) was used to assess the spectral type. The use of this feature was first discussed by HH14 who present it in the context of a spectral index. The depth of this absorption feature is best determined by comparing the flux at 510 nm to the flux expected at the same wavelength from a linear fit between the median flux of the spectral regions of λλ 460–470 nm and λλ540–550 nm. Figure 2 displays this feature for spectral types from K1 to K6. Here the linear fit between λλ 460–470 nm and λλ 540–550 nm is also indicated. Additionally, a linear fit for the regions of λλ 488–492 nm and λλ 514–516 nm was included in this figure to highlight the increasing depth of this feature with SpT. We perform this comparison at different values of extinction to estimate it simultaneously. The results are listed in Table 1, where we have rounded the SpT to half a subclass.

Fig. 2. X-shooter spectra of class III YSOs with a SpT ranging from K6 to K2. The spectra have been normalized at 460 nm and an offset has been added for clarity. The spectra have been smoothed to have a resolution of ∼2500 at 750 nm. The observations presented here are displayed in red, while the spectra of MTR13 and MFA17 are indicated in black. We highlight the R510 spectral index feature using the solid blue line. The dashed blue line indicates the slope of the surrounding continuum.

Our estimates are typically consistent within at most 1.5 subclass of the literature values. The TiO spectral indices agrees within 1 subclass for targets later than K5. The spectral indices of HH14 agree within a subclass at all spectral types except for the K5 to M0 range. Finally the indices of Riddick et al. (2007) have a similar agreement for the targets later than M3. We note that the use of features in the NIR region (λλ > 1000 nm) can give rise to significantly different spectral types. On top of the fact that the NIR indices are based on features covering a larger wavelength range and thus are more affected by extinction and by imperfect telluric removal and or flux calibration in that wavelength range (e.g., MTR13), the likely explanation for this deviation is the presence of cold spots on the stellar surface (e.g., Stauffer et al. 2003; Vacca & Sandell 2011; Pecaut 2016; Gully-Santiago et al. 2017; Gangi et al. 2022).

We adopted similar uncertainties on the SpT as MFA17. For objects later than K6 we estimate the uncertainties to be 0.5 subclass. For objects earlier than K6 we estimate the uncertainties to be 1 subclass. For the binary components of THA15-36 we assumed larger uncertainties of 1 subclass due to the less certain flux calibration of the spectra. We estimate the uncertainty on the extinction to be 0.2 mag. A description of how we confirmed the uncertainties on the SpT and estimated those on the extinction follows in Sect. 4 and Appendix B.

3.2. Photospheric properties from ROTFIT

A large number of absorption lines are resolved in the spectra, allowing us to obtain photospheric properties. We analyzed the absorption lines in the VIS arm of the X-shooter spectra using the ROTFIT tool (e.g., Frasca et al. 2015, 2017) to derive the effective temperature (T_eff), surface gravity (log g), projected rotational velocity (v sin i), radial velocities (RV), and veiling at three wavelengths (λ 620 nm, λ 710 nm and λ 970 nm). The veiling is given by r = EW(template)/EW(observation) − 1, where EW is the equivalent width of absorption lines near the wavelength of interest. ROTFIT searches, within a grid of templates, for the spectrum that minimizes the χ² to the target spectrum in different spectral regions. The templates consist of a grid of BT-Settl model spectra (Allard et al. 2011) of solar metallicity, log g ranging from 0.5 to 5.5 dex and an effective temperature ranging from 2000 to 6000 K. ROTFIT achieves this by first deriving the radial velocity of the star cross-correlating the target and template spectrum and uses this information to shift the observed spectrum to the rest frame. ROTFIT then convolves the templates with both a Gaussian to match the X-shooter resolution. Then the templates are iteratively broadened by convolution with a rotational profile with increasing v sin i until a minimum χ² is attained. ROTFIT analyzes spectral intervals that contain features sensitive to log g and/or the T_eff such as the KI doublet at λ ≈ 766 − 770 nm and the NaI doublet at λ ≈ 819 nm. An additional continuum emission component is added to the template in order to estimate the veiling. As for the v sin i, per each template, the veiling is also a free parameter in the fit.

We did not set this continuum component and therefore the veiling to 0 despite most of our targets being non-accretors. We preferred to let the veiling vary for our class III stars as a check. As expected, we always found zero or low veiling values, which could be the result of small residuals in the scattered light subtraction or the effect of the interplay between different parameters. The largest value of r = 0.4 is found at 620 nm for three K-type stars. A possible explanation could be the effect of starspots on the spectra. Another possible explanation is that of an undetected companion. Apart from these extreme values, we do not consider veiling values as small as 0.2 to be significant. The photospheric parameters, v sin i, veiling measurements, and RV derived with ROTFIT can be found in Table 2. As stated in Sect. 2, during this analysis, we found two targets to be double lined spectroscopic binaries, namely NN Lup and CD-39 10292. We excluded these from our analysis. Figure 3 shows the comparison between the T_eff derived using ROTFIT and obtained using different SpT–T_eff conversions. Here it can be seen that the effective temperature obtained from the SpT agrees best with ROTFIT when using the SpT–T_eff relation by HH14. The low log g found are mostly compatible with that of young objects (≲10 Myr).

Fig. 3. Comparison between the temperatures obtained with ROTFIT and from the spectral type for the targets presented here. The relation of Herczeg & Hillenbrand (2014) is labeled as HH14, the one by Pecaut & Mamajek (2013) as PM13, that of Kenyon & Hartmann (1995) as KH95 and the relationship of Luhman et al. (2003) as Luh03. The relations are only applied to spectral type ranges where they are valid. The solid line represents the one-to-one relation.

3.3. Luminosity determination

The stellar luminosities were obtained using a method similar to that used by MTR13. In this method, the dereddened spectra are extrapolated to wavelengths not covered by the spectra and integrated to compute the bolometric flux.

The X-shooter spectra were extrapolated using a BT-Settl (Allard et al. 2011) synthetic spectrum appropriate for the target. We used the effective temperature obtained from the SpT–T_eff relation of HH14. We prefer this T_eff, since in Sect. 5.1 we will use the obtained SpT to compute the luminosity. For the metallicity of the synthetic spectra we assumed solar values, typical for targets in the Lupus region (Biazzo et al. 2017), and for the surface gravity we chose log g = 4.0 in all cases. This parameter has very limited impact on the global luminosity estimate, since most of the stellar emission is covered by the X-shooter spectra and the synthetic spectra are used to measure a minor fraction of the global emission. The synthetic spectra were then matched to the spectra at 400.5 nm and 2300 nm and used to represent the flux at wavelengths shorter and longer than these respective values. The wavelengths chosen here differ from those used by MTR13. We chose these wavelengths in order to avoid the high noise level at short wavelengths and poor flux calibration at the end of the K-band present in some spectra. Due to this difference, we also reapplied this method to the sample of MTR13 and MFA17. We also performed a linear interpolation across the strong telluric absorption features between the J, H, and K band (λλ 1330–1550 nm, λλ 1780–2080 nm) in order to account for the stellar flux in these regions. An example of a prolonged spectrum created in this way is shown in Appendix C. The bolometric flux was then computed by integrating over the prolonged spectrum.

We used the photogeometric distances of Bailer-Jones et al. (2021) to convert the fluxes to luminosities for all targets with exception of TWA 26. For this star we used the geometric distance of the same authors since a photogeometric distance is unavailable. The distances provided by Bailer-Jones et al. (2021) are based on the Gaia DR3 (Gaia Collaboration 2016, 2023) parallaxes but include a photometric and geometric prior to improve the distance estimate. The inclusion of these priors significantly improves the uncertainties for the most distant targets (namely those in σ Orionis). For more nearby targets the distance estimate and its uncertainty do not change significantly from the inverse Gaia parallaxes.

The uncertainties on these estimates are computed by propagating the uncertainties of matching the BT-Settl models to the observations, the distance, the extinction, and the photometric flux. We adopt an extinction uncertainty of σ_AV = 0.2 mag. For the targets calibrated using wide slit spectra, we assumed the uncertainties on the photometric flux 5% of the total flux, which is slightly more conservative than the 4% of Rugel et al. (2018). For the spectra calibrated using the available photometry we assumed an uncertainty of 0.2 dex on the luminosity, the same as that of MTR13.

We also computed the luminosities of the grid using the bolometric correction of HH14. This correction uses the flux measured at 751 nm to estimate the total flux of the star. We noticed a small but systematic discrepancy between the luminosities derived from the two methods at effective temperatures lower than 4500 K. We therefore propose a correction to the relationship provided by HH14 which is further discussed in Appendix C. This new relation is assumed in this work.

Figure 4 shows our sample in a Hertzsprung-Russell diagram together with the samples of MTR13 and MFA17. Here, we used the luminosities obtained by integrating the extended spectrum and the effective temperature is obtained from the spectral type using the relationship of HH14. The isochrones and evolutionary tracks of Baraffe et al. (2015a) are also plotted. The isochronal age of the targets in our sample appears to be dependent on the stellar mass. Targets with M_⋆ ≳ 0.7 M_⊙ tend to be ∼10–30 Myr old according to the isochrones, whereas those with lower mass have ages of ∼1–10 Myr, despite all except one target being members of the Lupus star forming region. This trend is also present in the samples of MTR13 and MFA17 as well as other star samples such as those of Bell et al. (2014), Herczeg & Hillenbrand (2015), and Pecaut (2016). Stelzer et al. (2013) see a similar trend when comparing their sample to isochrones in the log g − T_eff diagram. It is unclear what causes this discrepancy but possible origins include the presence of starspots affecting the position in the HRD (e.g., Gangi et al. 2022; Pérez Paolino et al. 2024) or affecting overall evolution (Somers & Pinsonneault 2015), the effect of accretion during earlier stages of the star’s evolution (Baraffe et al. 2017) or the effect of the stellar magnetic field on the subphotospheric convective motions in the stars (Feiden 2016). The evolutionary tracks by Feiden (2016) include the later effect and are similar to the models of Baraffe et al. (2015a) in the nonmagnetic case. The SPOTS models of Somers et al. (2020) includes the influence of starspots at different filling factors on the isochrones. An analysis of the spot coverage of our sample and comparison with the SPOTS isochrones is beyond the scope of this work and deferred to a future work.

Fig. 4. HR diagram of the objects analyzed here (highlighted with larger black outlined star symbols) and those analyzed by MTR13 and MFA17. The model isochrones and evolutionary tracks by Baraffe et al. (2015a) are also shown. The isochrones are the 1.2, 3, 10 and 30 Myr ones.

4. New combined grid

The first goal of this work is to improve upon the grid of template spectra of YSOs presented by MTR13 and MFA17 by adding new templates to the grid and interpolating the spectra to generate a continuous grid. This is discussed in this section.

4.1. Description of the new grid

We combined the observations presented here with those of MTR13 and MFA17 to create an enhanced library of spectral templates that is made available to the community. We add one additional spectrum presented by Manara et al. (2016), that of the K0 star HBC407. The sources selected by MTR13 and MFA17 were classified to be class III YSOs and using Spitzer data (e.g., Evans et al. 2009). Manara et al. (2014) presented additional X-shooter observations of class III targets, namely IC348-127, T21, CrA75. These targets were highly extincted, with A_V ∼ 6 mag, A_V = 3.2 mag and A_V = 1.5 mag, respectively. This high extinction causes additional uncertainties in the dereddened spectra due to the uncertainties in the assumed extinction law. Therefore we do not include these spectra in our grid.

There are 6 stars from our new observations that are not included in the grid. The 5 unresolved binaries previously mentioned in Sect. 2 are excluded from our grid since they cannot be used to represent the stellar emission of individual stars. To avoid potential accretors in our sample we we used the Hα equivalent width criterion of White & Hillenbrand (2004) and Barrado y Navascués & Martín (2003), among others. The spectrum of 2MASSJ16075888−3924347 shows evidence of ongoing accretion, therefore we excluded it from our grid. This is further discussed in Appendix E.

From the grid of MTR13, Sz121 and Sz122 were excluded since they are likely spectroscopic binaries or ultrafast rotators (MTR13, Stelzer et al. 2013). For RXJ0438.6+1546 a spectrum was presented by MFA17. This target was also observed as a part of the PENELLOPE VLT Large Program (Manara et al. 2021). We only use the more recent PENELLOPE spectrum since the observations were performed with better sky transparency conditions.

We confirmed the spectral types of all except for one target of MTR13 and MFA17 through visual inter-comparison. The TiO absorption bands of CD 36-7429A appear to best match that of the K7 templates. We reassigned the spectral type of CD 36-7429A from K5 (MTR13) to K7. This spectral type was further confirmed by the TiO index of Jeffries et al. (2007) which yields a spectral type of K7.0 and the analysis of Fang et al. (2017, 2021) who also adopted a SpT of K7 for this target. Pecaut & Mamajek (2013) also adopt a SpT of K7 for this source. Interestingly, the estimate of T_eff from the ROTFIT analysis of the absorption lines lead to a result more consistent with a SpT of K5 (Stelzer et al. 2013). It is possible that spots covering the stellar surface lead to this discrepancy (e.g., Stauffer et al. 2003; Vacca & Sandell 2011; Pecaut 2016; Gully-Santiago et al. 2017; Gangi et al. 2022). We however assumed K7 as the SpT for this target based on the molecular features. The spectra first presented in this work have been dereddened using the values listed in Table 1 to provide an extinction-less grid. This grid can be found on GitHub.

The uncertainties on both A_V and SpT were estimated by fitting the spectra of targets within 1 spectral type subclass of other targets in the combined grid. In this procedure, we search for the best fitting template spectrum in the remainder of the grid at different values of extinction. The best fit template and extinction are found by searching for the minimum of a χ²-like metric that includes the spectral features discussed in Sect. 3.1.

For both the difference in spectral type and extinction we find a distribution with a median value of 0 and standard deviations of ∼0.5 subclasses and ∼0.25 mag respectively. We note that the standard deviation is larger for spectra earlier than K6, in part due to the lack of spectra of similar spectral type. We therefore estimated the uncertainties on SpT’s earlier than K6 to be ∼1 subclasses, those at later SpT to be ∼0.5 subclasses. For the extinction we adopt an uncertainty of 0.2 mag. THA15-36B appears as a strong outlier with ΔSpT = −1.5 and ΔA_V = 0.8. This supports our decision to provide higher uncertainties on both THA15-36A and THA15-36B. Because of this we also excluded THA15-36A and THA15-36B when constructing our interpolated grid in Sect. 4.2.

The final grid of templates is reported in Table F.1. Figure 5 shows the spectral type distribution of all the spectra in our final grid. Our final grid of templates includes 57 targets and spans a range of spectral types from G5 to M9.5. The spectral types later than K6 appear to be well represented. The earlier spectral types appear less well sampled, with the most prominent gaps between G5 to G8 and K4 to K5.5.

Fig. 5. Histogram of SpT for the grid of class III templates.

4.2. Interpolation

The grid presented in Sect. 4.1 is limited in two key ways. First, the grid is sparse and possesses gaps at a number of spectral types. Secondly, individual spectra have intrinsic uncertainties. Directly applying one of these spectra to represent the photosphere of a class II star or to analyze other class III targets may therefore bias the results. In order to mitigate these downsides we interpolated the grid of template spectra for the characterization of young stellar photospheres.

We first started by normalizing the dereddened spectra to the flux measured at 731 nm in order to remove the dependency of individual target’s distance and luminosity. We chose 731 nm as our reference region, since it is free of strong telluric features, avoids features dependent on log g and the photospheric emission dominates at this wavelength. Median normalized fluxes were then computed in a wavelength interval of choice for each of the spectra. To compute the uncertainty on this flux, we propagated the uncertainties on the flux in the normalization range (∼731 nm) and the flux in the wavelength range. We also propagated the uncertainties on the extinction. The spectra of some of the later spectral type stars in our sample have a signal to noise of ∼0 at the shortest wavelengths in the UVB arm. In this case, we excluded these normalized fluxes from the procedure. This procedure results in a set containing a normalized flux f_λ/f_731 nm and associated uncertainty for each template in our library.

We used a non-parametric local polynomial fit (Cleveland 1979) to interpolate between these fluxes as a function of SpT. In this type of interpolation, the value at a given position (spectral type in our case) is obtained by fitting a polynomial to a weighted version of the data points, this is expressed in Eq. (1). The polynomial at a given x, representing SpT subclass, is defined by the parameters ${\widehat{a}}_1(x),\dots ,{\widehat{a}}_n(x)$ as

$$\begin{array}{c}\hfill {\widehat{a}}_1(x),\dots ,{\widehat{a}}_n(x)={\mathrm{argmin}}_{a_0,\dots ,a_n}{\displaystyle \sum _i}K\left(\frac{x-x_i}{h}\right){(y_i-{\displaystyle \sum _{k=0}^n}{a}_n·x_i^k)}^2.\end{array}$$(1)

Here x_i represents the SpTs of our set of median fluxes in a given wavelength range and y_i is the values of the normalized fluxes themselves. K is a chosen kernel used for weighting the values that are fitted and h is the selected bandwidth. We used a polynomial of degree 2, a bandwidth of 2.5 SpT subclasses, and a Gaussian kernel.

To obtain the value of the local polynomial fit (y) at spectral type position (x), the expression

$$\begin{array}{c}\hfill y={\displaystyle \sum _{k=0}^n}{\widehat{a}}_n·x\end{array}$$(2)

is evaluated. This procedure is repeated at several equally spaced x values representing spectral types from G8 to M9.5, allowing us to retrieve the interpolated value at any SpT within this range. We limited our interpolation to range from G8 to M9.5 because of the large gap between our only G5 spectrum and the rest of our grid. We made use of the implementation provided by the localreg PYTHON package¹.

To account for the heteroscedastic uncertainties we used a Monte-Carlo simulation of 1000 iterations. In each iteration, we resampled the normalized flux and SpT of each data point by adding a value sampled from a Gaussian distribution with a standard deviation equal to the respective error and a mean of zero. The non-parametric fit was computed for each iteration. At each of the spectral type points we adoped the median values of these fits as the final interpolated model spectrum and for the error we adopt the 1-σ  interval around this median. Figure 6 shows an example of this procedure for one wavelength range.

Fig. 6. Example of a local polynomial fit to the normalized fluxes (black points) extracted in the wavelength range of 399–402 nm. The red point have been excluded from the fit due to the low S/R of this spectrum. The blue line indicates the median fit resulting from the Monte Carlo simulation. The transparent blue region indicates the 1-σ  uncertainty interval. The residuals are computed using only the uncertainties on the non-parametric fit.

We repeated this procedure for multiple wavelength ranges. An interpolated spectrum is then obtained by evaluating the values of these multiple non-parametric fits at a given spectral type. We applied it to wavelength ranges of 1 nm width over the entire UVB and VIS arms. A comparison between our interpolated spectra and class III templates of the same spectral type is shown in Fig. 7. Here it can be seen that individual templates can have a flux that is lower (PZ99J160843.4−260216) or higher (Par Lup 2) than that of the interpolated templates at short wavelengths. This effect could be a consequence of the uncertainties on the SpT and/or A_V of the spectra and/or differing levels of chromospheric activity. We find that roughly half of the class III spectra have a Balmer continuum flux higher than that of the interpolated templates while the other half have a lower flux. Figure 7 also shows that the interpolation can deviate significantly from the observed spectra in the region around 950 nm. This is a consequence of a poor telluric correction in several template spectra. This extinction-less interpolated grid can be found on GitHub along with the Python script used to generate it and obtain the interpolated templates.

Fig. 7. Comparison between individual class III templates and the interpolated spectra of the same spectral type. The templates are plotted in black and the interpolated spectra in blue and orange. The transparent blue/orange regions represent the 1-σ uncertainties in the interpolated spectra. The blue and orange colors correspond to the wavelength ranges of the X-shooter VIS and UVB arms respectively. The templates have been convolved with a Gaussian kernel for clarity.

5. Interpolated spectra for spectral fitting

In this section, we apply the new class III interpolated grid to the classic problem of deriving in a self-consistent way the stellar (SpT, A_V, L_⋆) and accretion properties (L_acc) of an actively accreting TTS. In particular, we compare our results to those obtained using a method in which the individual template spectra (those listed in Table F.1) are used to account for the stellar contribution.

Determining accretion and stellar properties in YSOs is not trivial. Excess continuum emission due to accretion and extinction alters the observed spectrum in opposite ways, although with a different wavelength dependence. A self consistent method for deriving stellar properties, accretion luminosity, and extinction simultaneously was previously presented by Manara et al. (2013a). The method presented here is a further development of this prescription, but with a different approach to the way the stellar photosphere is accounted for when modeling an observed spectrum. This also implies changes in how the stellar luminosity is calculated. We refer to this method as FRAPPE: FitteR for Accretion ProPErties. Here, we apply this method to the 14 accreting T Tauri stars in the Chamaeleon I star-forming region observed with X-shooter during the VLT/PENELLOPE Large Program (Manara et al. 2021).

5.1. Fitting the UV excess

The fitting procedure used here combines three components to fit an observed spectrum: (1) A continuum slab model is used to represent the emission from the accretion shock; (2) a reddening law to account for extinction; (3) our interpolated templates are used to represent the stellar contribution. Our fitting procedure consists of a Python code that searches in a grid of parameters the combination that best reproduces the observed spectrum of the input accreting YSO.

For the continuum excess emission we used the grid of isothermal hydrogen slab models developed by Rigliaco et al. (2012) and Manara et al. (2013a), the details of which details are described in Manara (2014). Prior to these works, such an approach was used by Valenti et al. (1993) and Herczeg & Hillenbrand (2008). Here we describe the key parts of the model. The models assume local thermodynamic equilibrium (LTE) and include emission from both H and H⁻. A given models are described by an electron density (n_e), electron temperature (T_e⁻), and the optical depth at 300 nm (τ). The grid contains slab models for n_e ranging from 10⁻¹¹ cm⁻³ to 10⁻¹⁶ cm⁻³, T_e⁻ ranging from 5000 to 11 000 K and τ from 0.01 to 5. The slab models are matched to the observed, dereddened spectrum of an accreting YSO using a scaling factor (H_slab). The models extend below the minimum wavelength of X-shooter (≲330 nm), allowing us to compute L_acc from the total flux of the slab model when scaled to match the observed spectrum during the fitting procedure.

We used the Cardelli extinction law (Cardelli et al. 1989) to account for extinction in our procedure. We fixed the total-to-selective extinction value R_V to the average interstellar value of R_V = 3.1 for all of the fits performed in this work. To find the best fitting A_V value we let the extinction free to vary from 0.0 to 2.0 mag in steps of 0.1 mag. When an extinction value is found at the upper edge of this range, we reran the procedure for larger A_V.

The main difference with respect to the method presented by Manara et al. (2013a) is the use of an interpolated grid of templates created with the method described in Sect. 4.2., while Manara et al. (2013a) used the individual spectra in the grid of class III templates from MTR13 (and later MFA17) to represent the photospheric and chromospheric contribution of the central star. During the fitting procedure, this grid is sampled at specific SpT values. We ran our procedure for SpTs at steps of 0.5 subclasses. The interpolated spectra were scaled to the observed, dereddened spectrum of the accreting YSO using a scaling factor (K_cl3) which sets the stellar luminosity.

The best fit is found by minimizing the likelihood function given by

$$\begin{array}{c}\hfill {\chi }_{\mathrm{like}}^2={\displaystyle \sum _{\mathrm{W}.\mathrm{L}.\mathrm{ranges}}}\frac{{(f_{\mathrm{mod}.}-f_{\mathrm{obs},\mathrm{dered}})}^2}{{\sigma }_{\mathrm{obs},\mathrm{dered}}^2+{\sigma }_{\mathrm{mod}.}^2}.\end{array}$$(3)

Here f_mod. is the flux density of the model spectrum, which comprises the sum of the scaled slab model and scaled interpolated class III template, within one of the selected wavelength range. f_obs, dered is the flux of the observed spectrum in the corresponding range after dereddening, σ_obs is the noise of the dereddened observation in the same range. This expression is different than that from Manara et al. (2013a) in that it includes σ_mod., the uncertainties on the interpolated grid, taking into account the scaling applied to match the observations (i.e., σ_mod. = K_cl3 ⋅ σ_interp.). At the edges of the interpolations SpT range its uncertainties become larger in all wavelength ranges. This will bias the results toward the edges of our spectral type range. To prevent this, we set the errors at SpT’s earlier than K0 equal to those obtained at K0.

We selected a limited number of wavelength ranges to create an interpolated grid of templates using the method presented in Sect. 4.2. These wavelength ranges were selected starting from those used by Manara et al. (2013a), but with some additions, as they carry information about the SpT and about the excess emission due to accretion. Wavelength ranges not considered by Manara et al. (2013a) are highlighted in Table 3. We note that the way Manara et al. (2013a) used the information in these ranges differs from ours. In the best fit determination of Manara et al. (2013a) the slopes of the observed Balmer and Paschen continuum and the flux ratio at both sides of the Balmer jump are matched to the model. We, however, chose to directly match the flux within these wavelength ranges. We also modified some of the regions adopted from Manara et al. (2013a) to prevent an overlap in wavelength. We included a number of TiO bandheads around ∼710 nm to constrain the spectral type of targets ranging from K6 to M9.5. To constrain targets with SpT from G5 to K6 we included wavelength ranges around ∼465, ∼510, and ∼540 nm. These are based on the ones used in the R515 index presented by HH14 with slight adjustments being made to avoid the sharp edges of features which are present in the spectra of later type stars.

To constrain the accretion spectrum we included two wavelength ranges that carry information about the size of the Balmer jump to be constrained (∼361 nm and ∼400.5 nm), three additional ranges to further constrain the Balmer continuum (∼335 nm, ∼340 nm and ∼355 nm) and three more to constrain the Paschen continuum (∼414.5 nm, ∼450 nm, and ∼475 nm). At longer wavelengths, the inner disk of class II YSO can contribute to the observed flux (e.g., Pittman et al. 2022). Since our models do not include such a component, we omitted the use of wavelength ranges in the NIR arm. Suboptimal telluric correction of both the templates and spectra being analyzed can affect the analysis when included. Wavelength regions that include telluric lines were therefore also omitted. The exact wavelength ranges can be found in Table 3.

The accretion emission is known to veil photospheric absorption features in the observed YSO spectrum. Therefore we compared the Ca I line at ∼420 nm in the observed spectrum with that of a veiled class III template to verify the goodness of our fits. We used a random class III template with the SpT nearest to the best fit solution and the best fitting slab model for this test.

Figure 8 shows how the normalized TiO 710 nm features vary as a function of SpT in our interpolated grid. Here it can be seen that the different wavelength rages disperse at late spectral types, showing the usefulness of this feature to classify M stars.

Fig. 8. Example of several interpolated normalized fluxes and their associated uncertainties. The values of the fluxes in the TiO 710 nm absorption band can be seen to diverge at late spectral types, indicating its usefulness to constrain the spectral type of these later than K5. The shaded areas indicate the respective 1 sigma uncertainties.

The procedure works as follows. For each spectral type we dereddened the observed spectrum for each value of A_V. Then for each of these combinations, we considered all of the slab models. The scaling factor for each combination was obtained by having the combined model match the dereddened observation at ∼355 nm and ∼731 nm. We note that the latter wavelength range is the same one as used to normalize our grid of interpolated spectra. K_clIII therefore corresponds to the stellar flux at 731 nm. In other words, we solved the following system of equations numerically for K_cl3 and H_slab.

$$\begin{array}{c}\hfill \{\begin{array}{c}{f}_{355}=K_{\mathrm{cl}3}·f_{\mathrm{clIII},355}+H_{\mathrm{slab}}·f_{\mathrm{slab},355}\hfill \\ f_{731}=K_{\mathrm{cl}3}·f_{\mathrm{clIII},731}+H_{\mathrm{slab}}·f_{\mathrm{slab},731}\hfill \\ =K_{\mathrm{cl}3}·1+H_{\mathrm{slab}}·f_{\mathrm{slab},731}\hfill \end{array}\end{array}$$(4)

We then used the parameters to calculate ${\chi }_{\mathrm{like}}^2$. This was done for every point in the grid. When this is done for the entire grid of parameters we select the model with the lowest ${\chi }_{\mathrm{like}}^2$ as the best fit. This optimization procedure is the same as that used by Manara et al. (2013a).

From this procedure, we directly obtained SpT and A_V. The accretion luminosity was calculated by first integrating the scaled slab model over its entire wavelength range (λλ= 50 nm to 2500 nm). This results in the total accretion flux F_acc which can then be converted to L_acc using L_acc = 4πd²F_acc, with d being the distance to the target. T_eff was obtained from the SpT using the relationship by HH14.

Due to our use of interpolated spectra, our stellar luminosity calculation differs from that of Manara et al. (2013a). To obtain the stellar luminosity we first computed the stellar flux at 751 nm, the reference wavelength of our bolometric correction. We obtained f₇₅₁ by first computing the flux at this wavelength in the dereddened observed spectrum and then subtracting the scaled slab model flux at the same wavelength. For sources with T_eff < 3500 K we applied the same interpolation over the VO feature as discussed in Appendix C. The bolometric flux F_bol. is then computed using our adjusted relationship of HH14 and converted to L_⋆ using L_⋆ = 4πd²F_bol.

The stellar radius was computed from T_eff and L_⋆ using R_⋆ = L_⋆/(4πσ_BoltzT_eff⁴). The stellar mass (M_⋆) and age are derived by placing the target on the HRD and interpolating evolutionary tracks at its position. We preferred to use to more recent models of Baraffe et al. (2015b). However, for 5 targets the T_eff and/or L_⋆ range of the Baraffe et al. (2015b) models did not extend to the found values, in these cases we used the tracks of Siess et al. (2000) instead. Finally, we can obtain the mass accretion rate Ṁ_acc using

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{acc}}=\frac{L_{\mathrm{acc}}{R}_{\star }}{G{M}_{\star }}(1-\frac{R_{\star }}{R_{\mathrm{in}}}{)}^{-1},\end{array}$$(5)

where we use the typical assumption for the inner disk radius of R_in = 5R_⋆ (Gullbring et al. 1998).

5.2. Accretion rates in accreting young stars in the Chamaeleon I region

To test our method we applied it to the sample of 14 Chamaeleon I targets observed with the VLT/X-shooter spectrograph as part of the VLT/PENELLOPE Large Programme (see Manara et al. 2021, for details on the survey). The results are listed in Table 4, and are presented here for the first time. A comparisons between the dereddened observations and best fit models can be found on ZENODO. For the purpose of comparison, we applied the procedure of Manara et al. (2013a) to the same sample and list the results in Table 5. In Fig. 9 the results of both methods are plotted on a L_acc–SpT diagram. Figure 10 displays both the obtained stellar luminosity and accretion luminosity. Results obtained for the same target are connected with a blue line to make the comparison easier. The values for the nearly complete sample of accreting young stars in the Chamaeleon I star-forming region presented by Manara et al. (2017a), using the values from Manara et al. (2023), is also shown with blue circles.

Fig. 9. Accretion luminosities vs. spectral type as obtained using FRAPPE (gray squares) and obtained using the method of Manara et al. (2013b) (red triangles). Results obtained for the same target with the two methods are connected with a blue line. The green triangle indicates the M5 solution of WZ Cha. The light blue circles are the nearly complete Chamaeleon I sample presented by Manara et al. (2013a). The black errorbars indicate the uncertainties on the accretion luminosity. The gray arrows indicate sources listed as upper limits by Manara et al. (2023).

Fig. 10. Accretion luminosities vs. stellar luminosity as obtained using FRAPPE (red triangles) and obtained using the method of Manara et al. (2013b) (gray squares). Results obtained for the same target with the two methods are connected with a blue line. The green triangle indicates the M5 solution of WZ Cha. The light blue dots is the nearly complete Chamaeleon I sample presented by Manara et al. (2013a). The black errorbars indicate the uncertainties on the stellar luminosity and accretion luminosity. The gray arrows indicate sources listed as upper limits by Manara et al. (2023).

It can be seen that the results obtained here are generally in agreement within the typical uncertainties for the fitting procedure of Manara et al. (2013a) (σ_AV ≈ 0.2 mag, σ_Lacc ≈ 0.25 dex and ${\sigma }_{{\dot{M}}_{\mathrm{acc}}}\approx 0.35\mathrm{dex}$, Manara et al. 2023). For the SpT we expect typical uncertainties of half a subclass for SpTs later than K6, for earlier SpTs we expect uncertainties of about 1 subclass. For 3 targets we find differences larger than the typical uncertainties. These targets are WZ Cha, CV Cha and VZ Cha. We describe the results for these targets in Appendix D.

The derived accretion luminosities of all targets agree within the aforementioned uncertainties. The obtained stellar luminosities are also in good agreement.

Figure 11 displays the resulting Ṁ_acc–M_⋆ diagram. For WZ Cha we obtained two degenerate solutions using the method of Manara et al. (2013a) one at SpT M3 (indicated in red) and one at SpT M5 (indicated in green). We prefer the M5 solution as the class III template to which the M3 solution was fitted appears to be an outlier (see Appendix D). Here it can be seen that, with exception of the M3 solution of WZ Cha, our results are within the typical uncertainties.

Fig. 11. Mass accretion rate vs. stellar mass as obtained using FRAPPE (red triangles) and obtained using the method of Manara et al. (2013b) (gray squares). Results obtained for the same target with the two methods are connected with a blue line. The green triangle indicates the M5 solution of WZ Cha. The light blue dots are the nearly complete Chamaeleon I sample presented by Manara et al. (2013a). The blue errorbars indicate the uncertainties on the stellar mass and mass accretion rate. The gray arrows indicate sources listed as upper limits by Manara et al. (2023).

The result obtained for WZ Cha and VZ Cha illustrates the benefits of our method (see Appendix D). For WZ Cha, FRAPPE breaks a degeneracy between spectral types and in the case of VZ Cha the Paschen continuum appears better reproduced. The analysis for CV Cha highlights the need for improved constraints when fitting late G to early K stars. Such constraints can not only be obtained from the inclusion of class III spectra around these spectral types in our grid but also from additional information on the emission and absorption lines present in the spectra of the class II YSOs.

6. Limits on the derived accretion properties

6.1. Limits on accretion luminosity measurements from the UV excess

The study of the class III spectra is important not just because they are used as templates of the photospheres of accreting stars, but also because they provide essential information on the limitations introduced by the stellar chromospheric activity to the measurements of low accretion luminosities. MTR13 and MFA17 discussed how chromospheric activity in the form of line emission defines the lower boundaries to L_acc derived from the well-established correlation between line emission and accretion luminosity, and named this L_acc, noise. Here we examined how the noise in the Balmer continuum emission in class III objects sets an upper boundary to the values of L_acc derived from the Balmer excess.

A variety of factors such as the uncertainties on the spectral type, extinction, and levels of chromospheric activity may contribute to the uncertainties on our interpolated model spectra. While such uncertainties may only cause minor contributions to the accretion luminosity measured in strong accretors, in weakly accreting objects such noise may have an important impact. Therefore, to characterize the limitations of our method we derived the accretion luminosities that correspond to the 1σ uncertainties on our interpolated spectra, and we refer to this value as L_acc, noise. To do this, we create a set of artificial spectra that include a UV excess equivalent to the 1σ uncertainties and fit these with FRAPPE.

Templates later than SpT M5.5 have a S/N < 1 in the Balmer and Paschen continuum. Therefore, including these templates in our analysis would strongly affect the uncertainty of the continuum UV flux. We therefore limited the analysis presented here to spectral types ranging from G8 to M5.5. We generate another interpolated grid based on our templates but excluding templates later than M5.5. We used this grid to create a set of blue enhanced artificial spectra. For wavelength ranges in the Balmer and Paschen continuum (< 500 nm) we set the value of the artificial spectrum equal to the interpolation spectra plus 1σ of the uncertainty thereon. The 1σ uncertainty is obtained from the Monte Carlo simulation discussed in Sect. 4.2. For longer wavelengths, we simply adopt the values of the interpolated spectrum. The artificial spectra are unitless since our interpolation is normalized to have the flux at 731 nm equal to 1. We therefore multiply the artificial spectra with a random value on the order of 0.1 × 10⁻¹⁴ to 100 × 10⁻¹⁴ erg nm⁻¹ s⁻¹ cm⁻² so that they have a “realistic” flux. We created such artificial spectra for each half subclass ranging from G8 to M5.5.

During this fitting procedure, we kept the SpT fixed to the value used for constructing the spectrum. The extinction was fixed to A_V = 0 mag. In doing so we obtained the accretion luminosities corresponding to spectra whose Balmer and Paschen continuum are 1σ interval above the interpolation. Here we assumed that the uncertainties in the Balmer and Paschen continuum are fully correlated.

The L_acc, noise and L_⋆ obtained with this method are arbitrary as they scale with the multiplicative factor used to create the artificial spectra. Therefore we do not report the obtained values of L_acc, noise and L_⋆. Instead, we report the ratio L_acc, noise/L_⋆ which is independent of this scaling and representative of the contrast down to which we can accurately measure a UV excess. The results are listed in Table 6. The dependence of L_acc, noise/L_⋆ as a function of SpT is shown in Fig. 12. Assuming the SpT–T_eff relation of HH14, we also plot L_acc, noise/L_⋆ as a function of T_eff in Fig. 12. Here it can be seen that the trend as a function of T_eff can be well approximated as a simple linear relation. This relationship is be expressed as

$$\begin{array}{c}\hfill log(L_{\mathrm{acc},\mathrm{noise}}/L_{\star })=(5.9\pm 0.2)·log({\mathrm{T}}_{\mathrm{eff}}[\mathrm{K}])-(23.3\pm 0.7).\end{array}$$(6)

The relationships in Eq. (6) and Table 6 can be seen as a criterion for the lowest measurable accretion luminosities. If the log(L_acc, noise/L_⋆) measured for a CTTS is of a similar order of magnitude or falls below it, then we can not conclude if a target is accreting based on the UV excess measurement. Several relations between the effective temperature and spectral type have been presented in the literature (e.g., HH14, Pecaut & Mamajek 2013; Kenyon & Hartmann 1995; Luhman et al. 2003). Due to this lack of uniformity, we recommend to directly use the values of L_acc, noise at each SpT reported in Table 6.

Fig. 12. SpT and Teff depedency of the accretion noise log(Lacc, noise/L⋆) values obtained from interpreting the uncertainties in the Balmer and Paschen continuum (< 500 nm) on our model class III spectra as accretion (green dots). The blue dots illustrate the log(Lacc, noise/L⋆) values obtained when fitting the class III spectra with an accretion slab model. The results obtained for the PENELLOPE Chamaeleon I sample is indicated with orange dots. The black line shows the Lacc, noise relationship measured from the emission lines by MFA17 and the red lines on the left show the uncertainty thereon. On the right, the red line illustrates the best fit to the log(Lacc, noise/L⋆) values.

Figure 12 also includes results obtained by fitting the class III templates themselves with FRAPPE. Here we once again kept the SpT fixed to the SpT of the respective class III and kept A_V = 0 mag. In this way, we obtain a similar L_acc, noise measurements. Only half of the class III templates provide a L_acc, noise measurement since half have a Balmer continuum below the interpolation. The L_acc, noise values obtained from the class III templates display a large scatter, this is a natural consequence of the scatter in the residuals between the templates and the interpolation in the Balmer and Paschen continuum. These results demonstrate how different the Balmer continua, and therefore the inferred UV-excess, of the various templates are. Using individual templates one may derive accretion luminosity that can change by fractions of the stellar luminosity of up to 10% relative to that found by using the interpolated spectra. The most extreme case is for late G SpT targets. This uncertainty is significantly weaker at later SpTs and in most cases it is only on the order of 1% of the stellar luminosity or lower.

We also plot the results for the class II targets presented in Sect. 5.2. One target, Sz19, clearly falls below the relationship in Eq. (6). Its measured mass accretion rate can therefore be highly affected by uncertainties on the chromospheric emission. Finally, for comparison, we also plot a similar relationship derived by MFA17 for the lowest accretion luminosities that one can measure by studying the emission lines. MFA17 measured the luminosities of several chromospheric emission lines in class III stars and converted them to accretion luminosities using the relations of Alcalá et al. (2014). We converted their results obtained as a function of spectral type to effective temperature using the SpT–T_eff of HH14 for consistency. It can be seen that the 1σ threshold that we derive for the UV continuum excess is higher in the entire SpT/T_eff range than the threshold obtained from the line luminosities. This implies that an analysis of accretion tracing emission lines allows us to measure lower accretion luminosities than the UV excess. However, the empirical relationships used to convert line luminosities to accretion luminosities (Alcalá et al. 2014, 2017) have been calibrated using measurements of the UV excess. It is therefore uncertain whether these relationships can be extrapolated to accretion luminosities below the threshold on the measurable UV excess. A more detailed modeling of the Hα emission lines, such as performed by Thanathibodee et al. (2023) may yet allow for the characterization of accretion properties of stars with low accretion rates, even well below the threshold of MFA17. Here we need to add the caveat that new work is currently being done on the influence of chromospheric emission lines on accretion rate measurements (Stelzer et al. in prep.), which will update the limit of MFA17 using the relationship between line and accretion luminosity of Alcalá et al. (2017).

We note that our criterion for measurable accretion luminosities and that of MFA17 are different in nature. The criterion of MFA17 represents typical chromospheric line luminosities, whereas ours represents the 1 sigma uncertainty on the stellar and chromospheric UV continuum emission in our interpolated class III model spectrum. The interpretation of both therefore differs. If the line luminosity measured in a CTTS is of a similar order of magnitude to the criterion of MFA17, then the chomospheric emission likely provides a dominant contribution to the measured L_acc. On the other hand, if a L_acc obtained from the UV excess is smaller than the limit presented here, then the uncertainties on this L_acc is dominated by the uncertainties on the stellar emission.

6.2. Implications for mass accretion rate estimates.

The limitations on measuring L_acc through the UV excess in turn reflects on the Ṁ_acc values that we can reliably measure. It is of particular interest to investigate how these limitations affect the Ṁ_acc − M_⋆ relation. We used the values of L_acc, noise from Eq. (6) to compute the lower limit of Ṁ_acc as a function of stellar mass using the nonmagnetic isochrones of Feiden (2016) at 1, 2, 3, and 6 Myr. The results are shown in Fig. 13. We use these isochrones since they extend to higher masses than those of Baraffe et al. (2015b) and therefore better cover the spectral types in Table 6. The isochrones of Feiden (2016) are mostly consistent with those of Baraffe et al. (2015b).

Fig. 13. Mass accretion rate vs. stellar mass for the Chamaeleon I sample presented by Manara et al. (2023) (blue dots) and the Chamaeleon I sample analyzed in this work. Objects that have an accretion luminosity lower than the criterion given in Eq. (6) are marked with a red dot. The limits we derived based on the uncertainties on the chromospheric emission are shown with the blue, orange, and green lines for 1.2, 3, and 6 Myr old objects respectively. The red errorbars indicate the uncertainties on the stellar mass and mass accretion rate. The gray arrows indicate sources listed as upper limits by Manara et al. (2023).

We also applied the relationship of Eq. (6) directly to the measured L_acc and L_⋆ of the Chamaeleon I sample from Manara et al. (2023) and the sample presented in Sect. 5.2. The results are also plotted in Fig. 13. The majority of class II targets fall well above the UV excess measurement limit. However, a considerable fraction appears to fall below the limit. The measured accretion rates of these targets are therefore dominated by the uncertainties on the stellar UV-continuum emission. More detailed analyzes are needed to disentangle the possible low accretion rate in these targets from the stronger chromospheric emission.

Manara et al. (2017b) used the Chamaeleon I sample shown here in combination with a sample of targets in Lupus (from Alcalá et al. 2017) to discuss whether the correlation between Ṁ_acc and M_⋆ is driven by the detection limits and if the observed scatter in Ṁ_acc fills the observable range. Only 3 sources from the sample of Manara et al. (2023) that fall below our detection limit were previously not indicated as upper limits. These targets all fall in a locus were upper limits were already present. Our limit on the measurable mass accretion rate therefore does not change the conclusion reached by Manara et al. (2017b). Namely, that this correlation is real and not a consequence of detection limits.

7. Summary and conclusions

We presented the analysis of 24 new VLT/X-shooter spectra of class III young stars for their use as photospheric templates to complement the previous samples of MTR13 and MFA17. All spectra are available in reduced and dereddened form on GitHub. We obtained the spectral type through the use of spectral indices and by comparing the relative strength of photospheric features to the samples of MTR13 and MFA17.

We employed ROTFIT to derive the photospheric properties (log g, T_eff, v sin i, RV) of the sample. The results of both methods are generally in agreement with each other and the photospheric properties confirm the young nature of these targets. As previously noted in the literature (e.g., Bell et al. 2014; Herczeg & Hillenbrand 2015; Manara et al. 2017a; Pecaut 2016), we see a mass dependent trend of the derived isochronal ages for targets in the same star forming region. We compared the stellar luminosities obtained through extending and integrating the X-shooter spectra in our sample with luminosities obtained by applying the bolometric correction of HH14. We found a small but systematic difference between the results of both methods for effective temperatures lower than about 4500 K, leading us to propose an adjustment to the bolometric correction of HH14.

To remove the nonuniform SpT sampling of the templates and mitigate the uncertainties associated with their use, we proposed a method of interpolating between the templates. We applied this type of interpolation to the entire wavelength range covered by the spectra in 1nm wide bins and made the results publicly available on GitHub⁶. We also applied the method to specific regions in the spectra that are useful to obtain the stellar and accretion properties of class II targets. This empirical method has been implemented in a code that self consistently measures the extinction, stellar, and accretion properties. We verified this framework by analysing the Chamaeleon I sample observed as part of the PENELLOPE program and compared the results to those obtained using the method presented by Manara et al. (2013a) on the same spectra. We found good agreement in the results obtained with the two methods. Several future improvements, such as additional class III templates at early SpTs or the inclusion of veiling measurements in the best fit determination, can still be made to FRAPPE. Such improvements are particularly useful to constrain the stellar properties stars with a spectral type earlier than K6 and will make FRAPPE more user friendly.

We measured the typical accretion to stellar luminosity ratios that would be obtained if the 1σ uncertainties in our model spectra were interpreted as excess emission due to accretion. This ratio is representative of the typical lower limit on the measurable mass accretion rates in young stars when using continuum UV excess. The value has a strong dependency on the spectral type. This trend is likely a consequence of the worse constraints on earlier SpT stellar emission within our sample of class III sources, as well as of a worse contrast with the UV excess due to earlier stars having the peak of the photospheric emission at bluer wavelengths, hampering the detection of an excess for contrast reasons. This limit is higher than the lower limit presented for accretion rates derived from emission lines presented by MFA17. We, therefore, conclude that an in-depth analysis of emission lines is needed to obtain accurate measurements of mass accretion rates in low accreting objects. We applied our lower limit to the Chamaeleon I sample of Manara et al. (2023). The majority of objects in this region were found to have accretion rates well above our lower limit.

This shows that the uncertainties associated with the use of class III templates do not significantly affect the observed correlations between mass accretion rate and disk mass or stellar mass (see Manara et al. 2023 for a review). Our improved methodology for self consistently deriving accretion properties will facilitate further studies into these relations.

Data availability

Additional data for this article are available at https://zenodo.org/doi/10.5281/zenodo.12686015 and at https://github.com/RikClaes/FRAPPE

Acknowledgments

We thank the anonymous referee for the useful comments to our submitted manuscript. R. Claes acknowledges the PhD fellowship of the International Max-Planck-Research School (IMPRS) funded by ESO This work was partly funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) in the framework of the YTTHACA Project 469334657 under the project code MA 8447/1-1. Funded by the European Union (ERC, WANDA, 101039452). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. J. B. Lovell acknowledges the Smithsonian Institute for funding via a Submillimeter Array (SMA) Fellowship. J.M. Alcalá acknowledges financial support from PRIN-MUR 2022 20228JPA3A “The path to star and planet formation in the JWST era (PATH)” funded by NextGeneration EU and by INAF-GoG 2022 “NIR-dark Accretion Outbursts in Massive Young stellar objects (NAOMY)”. J.M. Alcalá and A. Frasca acknowledge financial support from Large Grant INAF 2022 “YSOs Outflows, Disks and Accretion: towards a global framework for the evolution of planet forming systems (YODA)”. This work benefited from discussions with the ODYSSEUS team (HST AR-16129), https://sites.bu.edu/odysseus/. 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 the SIMBAD database, operated at CDS, Strasbourg, France.

References

Appendix A: Observations and data reduction

The sources MV Lup, NO Lup, and MT Lup were not observed using slits of 5.0″. The narrow slit spectra were therefore flux calibrated by using the available photometry (see Table A.2). This was done by converting the photometry in a given band to a corresponding flux density. The flux density in the same bands is also computed for the narrow-slit spectra by integrating them after being multiplied to the filter bandpass. The ratio of these fluxes allows us to rescale the spectra in absolute units correcting for slit losses. For MV Lup and NO Lup the spectrum in the UVB arm is matched to the photometric flux using a factor constant with wavelength. whereas VIS and NIR arm spectra are rescaled using a factor with a linear dependency on wavelength. This was done because variation of the seeing disk with wavelengths in the wide ranges of VIS and NIR arms and atmospheric refraction can produce a wavelength dependent light loss. For MT Lup a constant factor was used for the entire spectrum to preserve continuity between the UVB and VIS arms.

2MASSJ16090850-3903430 and THA15-36 were both observed consecutively during the night of 23/24-Jun 2022. Both are about ∼0.7 mag less bright than the archival photometry. The spectrum of 2MASSJ16090850-3903430 whas therefore rescaled with a constant factor to match the photometry (see Table A.2).

THA15-36 is a binary system with a separation of ∼2″ that was aligned-with and resolved-in the X-shooter slits. The traces were fitted with two Gaussians at a number of wavelengths. This customized spectra extraction allowed us to estimate the flux ratio between both components as a function of wavelength. The latter enabled us to calculate the total flux in the broad-slit spectra for both components. In the narrow slit observations components were extracted from the rectified, flux calibrated 2d spectra using an extraction window with a width of 0.8″ in all arms in order to avoid contamination from its companion. The extracted narrow slit spectra were then rescaled to their respective flux contributions to correct for wavelength dependent flux loss in the same way as mentioned before. The spectra were then multiplied by the same constant factor such that the sum of both components matches the available photometry (see Table A.2).

For RXJ1607.2-3839 we noticed noticed that the continuum in the NIR arm bent upward at shorter wavelengths and was less bright than the archival photometry. The bend in the continuum is likely caused by an issue with the used flux standard. The flux standard taken later that night was therefore used instead. This mitigated the upward bend of the VIS arm. After this correction, the overall flux density in the arm was still lower compared to the VIS arm and the NIR photometry. This was rectified by scaling the NIR to the VIS arm such that their overlapping regions have the same median flux density. This also improved the match between the NIR arm and the available 2MASS photometry.

The NIR arm of RECX-6 displayed a similar upward bend as RXJ1607.2-3839. To mitigate this issue the flux standard of the night before was used for the reduction of the NIR arm spectrum of RECX-6.

CD-31_12522 was observed in poor seeing conditions. This caused the slit losses in the UVB narrow slit observation to become wavelength depended. We, therefore, used a linear relation rather than the mean flux ratio to rescale the UVB narrow slit to the broad slit observation.

Appendix B: Uncertainties on extinction.

An independent confirmation of the extinction can be found through the ratio F_red = F(833.0nm)/F(634.8nm), as first suggested by HH14. This ratio is sensitive to both extinction and spectral type. Here we adopt the relationship of MFA17, who recalibrated the relations of HH14 using observations of class III YSO’s. MFA17 adjusted the relationships of HH14 to account for differences in observed M-type spectra and the corresponding synthetic and dwarf spectra on which (HH14) based their relations. Fig. B.1 shows these ratios for our observations after dereddening. Here it can be seen that our sample has a similar agreement with A_V = 0 as the samples of MTR13 and MFA17.

Fig. B.1. Fred = F(833nm)/F(634.8nm) ratio calculated for the class III templates analyzed here (orange stars) as well as for the samples presented by MTR13 and MFA17 (blue circles). The error bars represent the uncertainties due to the noise in the spectra and do not include the uncertainty on the extinction correction. The blue line shows the value of Fred for AV = 0 as a function of spectral type as derived by cite.manara17b. The green and orange lines are obtained by deredening and reddening the relationship of MFA17 by 1 mag respectively.

Appendix C: Luminosity calculation

Fig. C.1 shows an example of a extrapolated X-shooter spectrum used in our luminosity calculation. We notice a small but systematic discrepancy between the luminosities obtained by integrating our extrapolated X-shooter spectra and the values found by using the bolometric correction of HH14. Fig. C.2 shows that the ratio between the two luminosity estimates increases for later spectral types. The correction of HH14 was based on the same BT-Settl models that we used to extrapolate our X-shooter spectra. HH14 mention that for targets with T_eff < 3500 K they interpolate over a VO feature that is stronger in the BT-Settl models compared to their observations. We apply the same interpolation when using the bolometric correction.

The reason for the discrepancy appears to be due to additional flux in the NIR part of the observed spectra compared to the BT-Settl models. A likely cause for this could be the presence of cold spots on the stellar surface, which was discussed in Sect. 3.1. We therefore propose a correction to the relationship provided by HH14, which we expect to account for the typical contribution of stellar spots to the total luminosity. We compute this temperature dependent correction factor by using the non-parametric fit displayed on Fig. C.2. During this fitting procedure, we excluded Par Lup3 1, since it is an extreme outlier, with an extremely bright NIR flux for its SpT, which affects the resulting best fit. The adjusted bolometric correction can be found is illustrated in Fig. C.3 and its values are listed in Table C.1.

Fig. C.1. Example of a flux calibrated X-shooter spectrum extrapolated with a BT-Settl model spectrum. The model segments (dark green) are scaled at the blue and red edges of the X-shooter spectrum. The UVB, VIS, and NIR arms are plotted in blue orange and gray, respectively. The strong telluric bands in the NIR arm are replaced by a linear interpolation (red lines). In this figure, we show the example of RXJ1607.2-3839 which has been extended with a BT-Settl model of Teff = 4000 K. The BT-Settl model spectra have been convolved with a Gaussian kernel for clarity.

Fig. C.2. Ratio of the luminosities obtained from our spectral integration method and the value computed using the bolometric correction of (HH14). The blue dots indicate the values obtained for the spectra presented here aswell as in (MTR13) and (MFA17). The red line indicates the non-parametric fit, and the shaded area indicates the corresponding 1σ uncertainty interval obtained from bootstrapping the non-parametric fitting procedure. The orange point is Par-Lup3-1, which was excluded from this procedure.

Fig. C.3. The bolometric correction (F751/FBol) of HH14 and our suggested adjustment for the additional flux present in observations of class III YSOs.

Appendix D: Notable results for the Chameleon I sample

D.1. WZ Cha:

For WZ Cha our method obtains a SpT of M5 whereas the fit to the individual class III templates is degenerate with local solutions at M3 and M5. The local solution at SpT M3 appeared due to the CD -36 7429B template and since the solution appears to better reproduce the veiling in the Ca II line at ∼420 nm. CD 36 7429B does however appear as an outlier among the M3 templates in that its TiO features around 710 nm appear to be deeper than templates of the same SpT. When removing this template the degeneracy disappears. We therefore prefer the solution presented here.

D.2. CV Cha:

The best fit of CV Cha displays a strong degeneracy in terms of SpT, stretching from ∼K4 to ∼G9. The solution found here results in a SpT of K3. This degeneracy is a likely consequence of the poor constraint on early spectral types in our model spectra.

D.3. VZ Cha:

For VZ Cha we find an extinction that is ΔA_V = 0.5 mag lower than obtained using the method of Manara et al. (2013a). This differing result appears to be a consequence of the stricter constraint imposed on the Paschen continuum by the best fit metric used in our method, resulting in a better fit of this region. Despite the large difference in obtained A_V, the accretion and stellar properties still agree within errors.

D.4. Sz 19:

The luminosity of Sz19 appears large for its spectral type, resulting in a high mass estimate. Fig. D.1 illustrates this by showing Sz19 and our class III templates on an HRD. The spectra of this target agree with the available photometry within < 2 mag. The results presented here agree with those found by Manara et al. (2014) (although the slightly later SpT results in a very different stellar mass). An analysis of the emission lines results in accretion luminosities that are in agreement with those found from our analysis of the UV excess. A detailed investigation into why this source appears as an outlier is beyond the scope of this paper.

Fig. D.1. HRD containing the different results obtained for Sz19 for comparison the objects analyzed here (highlighted with larger black outlined markers) and those analyzed by MTR13 and MFA17 are displayed. The model isochrones and evolutionary tracks by Feiden (2016) are also shown. The isochrones are the 1.2,3,10 and 30 Myrs ones

Appendix E: 2MASSJ16075888-3924347

The spectrum of 2MASSJ16075888-3924347 appears to have spectral features typical of an accreting T Tauri star rather than a non-accreting class III object. For example, the Hα line luminosity appears much larger than in the other class III templates of similar spectral type (M5.5). We measure an equivalent width of EW = −3.914 ± 0.093 nm. This would classify the target as a CTTS rather than a WTTS according to the classification scheme of Barrado y Navascués & Martín (2003), White & Basri (2003). Similarly, the Ca II IR triplet is detected for this target but not in comparable class III templates. This result is puzzling considering that the SED of this target presented by Lovell et al. (2021), does not show any excess emission at longer wavelengths and is consistent with that of other class III objects. Moreover, Lovell et al. (2021) did not detect any flux at 856 μm, giving an upper limit of < 0.11 mJy. No previous Hα equivalent widths are available in the literature for this target. A detailed analysis of this star is beyond the scope of this work and will be presented by Stelzer et al. in prep..

Appendix F: The combined grid

All Tables

All Figures

Fig. 1. X-shooter spectra of class III YSOs with spectral type ranging from M1 to K7. All the spectra are normalized at 731 nm and offset in the vertical direction for clarity. The spectra are also smoothed to the resolution of 2500 at 750 nm to make easier the identification of the molecular features. The black colors indicate spectra presented by MTR13 and MFA17. The red color is used for spectra presented for the first time here.

In the text

Fig. 2. X-shooter spectra of class III YSOs with a SpT ranging from K6 to K2. The spectra have been normalized at 460 nm and an offset has been added for clarity. The spectra have been smoothed to have a resolution of ∼2500 at 750 nm. The observations presented here are displayed in red, while the spectra of MTR13 and MFA17 are indicated in black. We highlight the R510 spectral index feature using the solid blue line. The dashed blue line indicates the slope of the surrounding continuum.

In the text

Fig. 3. Comparison between the temperatures obtained with ROTFIT and from the spectral type for the targets presented here. The relation of Herczeg & Hillenbrand (2014) is labeled as HH14, the one by Pecaut & Mamajek (2013) as PM13, that of Kenyon & Hartmann (1995) as KH95 and the relationship of Luhman et al. (2003) as Luh03. The relations are only applied to spectral type ranges where they are valid. The solid line represents the one-to-one relation.

In the text

	Fig. 4. HR diagram of the objects analyzed here (highlighted with larger black outlined star symbols) and those analyzed by MTR13 and MFA17. The model isochrones and evolutionary tracks by Baraffe et al. (2015a) are also shown. The isochrones are the 1.2, 3, 10 and 30 Myr ones.
In the text

	Fig. 5. Histogram of SpT for the grid of class III templates.
In the text

Fig. 6. Example of a local polynomial fit to the normalized fluxes (black points) extracted in the wavelength range of 399–402 nm. The red point have been excluded from the fit due to the low S/R of this spectrum. The blue line indicates the median fit resulting from the Monte Carlo simulation. The transparent blue region indicates the 1-σ  uncertainty interval. The residuals are computed using only the uncertainties on the non-parametric fit.

In the text

Fig. 7. Comparison between individual class III templates and the interpolated spectra of the same spectral type. The templates are plotted in black and the interpolated spectra in blue and orange. The transparent blue/orange regions represent the 1-σ uncertainties in the interpolated spectra. The blue and orange colors correspond to the wavelength ranges of the X-shooter VIS and UVB arms respectively. The templates have been convolved with a Gaussian kernel for clarity.

In the text

	Fig. 8. Example of several interpolated normalized fluxes and their associated uncertainties. The values of the fluxes in the TiO 710 nm absorption band can be seen to diverge at late spectral types, indicating its usefulness to constrain the spectral type of these later than K5. The shaded areas indicate the respective 1 sigma uncertainties.
In the text

Fig. 9. Accretion luminosities vs. spectral type as obtained using FRAPPE (gray squares) and obtained using the method of Manara et al. (2013b) (red triangles). Results obtained for the same target with the two methods are connected with a blue line. The green triangle indicates the M5 solution of WZ Cha. The light blue circles are the nearly complete Chamaeleon I sample presented by Manara et al. (2013a). The black errorbars indicate the uncertainties on the accretion luminosity. The gray arrows indicate sources listed as upper limits by Manara et al. (2023).

In the text

Fig. 10. Accretion luminosities vs. stellar luminosity as obtained using FRAPPE (red triangles) and obtained using the method of Manara et al. (2013b) (gray squares). Results obtained for the same target with the two methods are connected with a blue line. The green triangle indicates the M5 solution of WZ Cha. The light blue dots is the nearly complete Chamaeleon I sample presented by Manara et al. (2013a). The black errorbars indicate the uncertainties on the stellar luminosity and accretion luminosity. The gray arrows indicate sources listed as upper limits by Manara et al. (2023).

In the text

Fig. 11. Mass accretion rate vs. stellar mass as obtained using FRAPPE (red triangles) and obtained using the method of Manara et al. (2013b) (gray squares). Results obtained for the same target with the two methods are connected with a blue line. The green triangle indicates the M5 solution of WZ Cha. The light blue dots are the nearly complete Chamaeleon I sample presented by Manara et al. (2013a). The blue errorbars indicate the uncertainties on the stellar mass and mass accretion rate. The gray arrows indicate sources listed as upper limits by Manara et al. (2023).

In the text

Fig. 12. SpT and Teff depedency of the accretion noise log(Lacc, noise/L⋆) values obtained from interpreting the uncertainties in the Balmer and Paschen continuum (< 500 nm) on our model class III spectra as accretion (green dots). The blue dots illustrate the log(Lacc, noise/L⋆) values obtained when fitting the class III spectra with an accretion slab model. The results obtained for the PENELLOPE Chamaeleon I sample is indicated with orange dots. The black line shows the Lacc, noise relationship measured from the emission lines by MFA17 and the red lines on the left show the uncertainty thereon. On the right, the red line illustrates the best fit to the log(Lacc, noise/L⋆) values.

In the text

Fig. 13. Mass accretion rate vs. stellar mass for the Chamaeleon I sample presented by Manara et al. (2023) (blue dots) and the Chamaeleon I sample analyzed in this work. Objects that have an accretion luminosity lower than the criterion given in Eq. (6) are marked with a red dot. The limits we derived based on the uncertainties on the chromospheric emission are shown with the blue, orange, and green lines for 1.2, 3, and 6 Myr old objects respectively. The red errorbars indicate the uncertainties on the stellar mass and mass accretion rate. The gray arrows indicate sources listed as upper limits by Manara et al. (2023).

In the text

Fig. B.1. Fred = F(833nm)/F(634.8nm) ratio calculated for the class III templates analyzed here (orange stars) as well as for the samples presented by MTR13 and MFA17 (blue circles). The error bars represent the uncertainties due to the noise in the spectra and do not include the uncertainty on the extinction correction. The blue line shows the value of Fred for AV = 0 as a function of spectral type as derived by cite.manara17b. The green and orange lines are obtained by deredening and reddening the relationship of MFA17 by 1 mag respectively.

In the text

Fig. C.1. Example of a flux calibrated X-shooter spectrum extrapolated with a BT-Settl model spectrum. The model segments (dark green) are scaled at the blue and red edges of the X-shooter spectrum. The UVB, VIS, and NIR arms are plotted in blue orange and gray, respectively. The strong telluric bands in the NIR arm are replaced by a linear interpolation (red lines). In this figure, we show the example of RXJ1607.2-3839 which has been extended with a BT-Settl model of Teff = 4000 K. The BT-Settl model spectra have been convolved with a Gaussian kernel for clarity.

In the text

Fig. C.2. Ratio of the luminosities obtained from our spectral integration method and the value computed using the bolometric correction of (HH14). The blue dots indicate the values obtained for the spectra presented here aswell as in (MTR13) and (MFA17). The red line indicates the non-parametric fit, and the shaded area indicates the corresponding 1σ uncertainty interval obtained from bootstrapping the non-parametric fitting procedure. The orange point is Par-Lup3-1, which was excluded from this procedure.

In the text

	Fig. C.3. The bolometric correction (F751/FBol) of HH14 and our suggested adjustment for the additional flux present in observations of class III YSOs.
In the text

	Fig. D.1. HRD containing the different results obtained for Sz19 for comparison the objects analyzed here (highlighted with larger black outlined markers) and those analyzed by MTR13 and MFA17 are displayed. The model isochrones and evolutionary tracks by Feiden (2016) are also shown. The isochrones are the 1.2,3,10 and 30 Myrs ones
In the text