© ESO 2021

1 Introduction

The way in which circumstellar disks evolve and form proto-planets is deeply influenced by the processes of mass accretion onto the star, ejection of outflows, and photo-evaporation of disk material through winds (Hartmann & Calvet 2016; Ercolano & Pascucci 2017, and references therein). In order to understand planet formation it is necessary to explain how optically thick accretion disks surrounding the youngest low-mass (M_⋆ ≲ 2.0M_⊙) stars evolve into optically thin debris disks (Morbidelli & Raymond 2016). In this framework, the mass accretion rate, Ṁ_acc, is a fundamental parameter for the evolution of accretion disks around young low-mass stellar objects (YSOs). Ṁ_acc measurements set important constraints for disk evolution models (Hartmann et al. 1998; Hartmann & Calvet 2016) and disk clearing mechanisms (Alexander et al. 2014; Ercolano & Pascucci 2017, and references therein).

In the current magnetospheric accretion paradigm for classical T Tauri (CTT) stars, the strong stellar magnetic fields truncate the inner disk at a few stellar radii (Donati & Landstreet 2009; Johns-Krull et al. 2013). Gas flows from this truncation radius onto the star along the stellar magnetic field lines, crashing onto the star and forming an accretion shock. The ~10⁴ K optically thick post-shock gas and optically thin pre-shock gas produce emission in the Balmer and Paschen continua and in many lines, including the Balmer and Paschen series and the Ca II IR triplet (Hartmann et al. 1998; Hartmann & Calvet 2016) observed in the optical spectra of CTTs.

The mass accretion rate can be derived from the energy released in the accretion shock (accretion luminosity L_acc; see Gullbring et al. 1998; Hartmann 1998) given the stellar properties. Observationally, this requires measurements of excess flux in continuum and lines with respect to similar nonaccreting template stars. Such measurements are best performed at ultraviolet (UV) wavelengths (λ ≲ 4000Å) with the Balmer continuum excess emission and the Balmer jump (see Herczeg & Hillenbrand 2008; Ingleby et al. 2013; Alcalá et al. 2014, 2017; Manara et al. 2017a, and references therein). In the past, L_acc has been calculated using veiling measurements in high-resolution optical spectra (e.g., Hartigan et al. 1991; Hartigan & Kenyon 2003; White & Hillenbrand 2004, and references therein). Also, it is well known that L_acc, and therefore Ṁ_acc, is correlated with the line luminosity, L_line, of H I, He I, and Ca II lines (e.g., Muzerolle et al. 1998; Calvet et al. 2004; Herczeg & Hillenbrand 2008; Rigliaco et al. 2012; Alcalá et al. 2014, 2017, and references therein). These latter works provide L_acc –L_line correlations simultaneously and homogeneously derived from the UV to the near-infrared (NIR), underlying the importance of these emission features as accretion diagnostics. These accretion tracers are key diagnostics with which to estimate L_acc via the correlations mentioned above when flux-calibrated spectra below λ ~ 3700 Å are missing.

On the other hand, accretion is a highly variable process (Basri & Batalha 1990; Jayawardhana et al. 2006; Cody & Hillenbrand 2010; Venuti et al. 2014), which leads to a range of Ṁ_acc values for a given object when measured at different epochs (see Costigan et al. 2012, 2014; Biazzo et al. 2012). Variability in YSOs induces dispersion in the observed Ṁ_acc–M_⋆ and Ṁ_acc–M_disk scaling relationships, but cannot explain the large scatter of more than 2 dex in log Ṁ_acc at a given YSO mass. Such scaling relationships are predicted by the theory of viscous disk evolution (Lynden-Bell & Pringle 1974; Hartmann & Calvet 2016; Rosotti et al. 2017, and references therein), but the Ṁ_acc–M_disk relationship has been confirmed observationally only recently by spectroscopic surveys in strong synergy with ALMA surveys of disks in star forming regions (Ansdell et al. 2016; Manara et al. 2016; Pascucci et al. 2016; Mulders et al. 2017).

As concluded in previous works (e.g., Rigliaco et al. 2012; Alcalá et al. 2014), the average L_acc and Ṁ_acc derived from several diagnostics, measured simultaneously, has a significantly reduced error. This suggests the need to use spectroscopic observations simultaneously performed from visible to NIR wavelengths and with the highest possible resolution to overcome problems due to line blending. These requirements can be achieved with the GIARPS (GIAno and haRPS, Claudi et al. 2017) high-resolution spectrograph at the Telescopio Nazionale Galileo (TNG) in the Canary Islands, Spain. The spectral coverage and resolution of GIARPS allows one to probe the properties of the accretion columns, hot spots, the inner gaseous disk, the stellar and disk winds, and the collimated jets (see Gangi et al. 2020; Giannini et al. 2019), making TNG/GIARPS a powerful instrument with which to investigate accretion in YSOs of the northern hemisphere, such as the Taurus-Auriga star forming region (SFR).

The Taurus SFR contains a rich population of pre-main sequence (PMS) stars, most of them T Tauri stars (see Kenyon et al. 2008, and references therein). Several works have addressed the problem of measuring Ṁ_acc in the Taurus population using the different methodologies mentioned above (see Herczeg & Hillenbrand 2008, 2014, and references therein). Also, a large number of CTTs have already been observed with ALMA (Andrews et al. 2018; Long et al. 2019), probing the outer regions of their disk, and allowing the scaling relationships predicted by the viscous disk evolution theory to be investigated further. However, despite being a general benchmark for star formation studies, a comprehensive and systematic analysis of the Taurus population, where the relevant stellar and accretion parameters are derived simultaneously with sufficient accuracy, is still missing. As a first step in filling this gap, here we present a pilot study of the accretion of seven YSOs in Taurus observed with TNG/GIARPS as part of the GHOsT (GIARPS High-resolution Observations of T Tauri stars) project. GHOsT is a survey of a flux-limited complete sample of T Tauri stars in the Taurus star forming region that is to be used to derive, in a homogeneous fashion (thus avoiding systematic errors due to the use of different sets of nonsimultaneous observations), the stellar, accretion, and outflow parameters, and to constrain the properties of both the inner disk and the associated winds and jets. The ultimate goal of GHOsT is to provide reliable measures of the mass-accretion and mass-loss rates of the Taurus population and to put them in relation with the properties of the central star and its disk, in synergy with the complementary ALMA observations (Andrews et al. 2018; Long et al. 2019).

A comprehensive study of the jet line emission of the targets investigated in this paper was published by Giannini et al. (2019, henceforth Paper I) and an investigation of the link between atomic and molecular winds was published in Gangi et al. (2020, henceforth Paper II), while the results of the complete GHOsT survey will be presented in forthcoming papers. One of the main goals of the present paper is the definition and assessment of the methodologies for the determination of the stellar and accretion properties, in a homogeneous and self-consistent way, to be adopted for the analysis of the complete set of GHOsT data.

The paper is organized as follows. In Sect. 2, we present the target selection, the observations, and the data processing. In Sect. 3, we show the procedures used for the determination of the stellar parameters, and veiling estimates both in the optical and in the NIR. In Sect. 4, the methodologies to derive the accretion parameters are presented, while in Sect. 5 we discuss results on the stellar and accretion properties of the studied sample, comparing with those of the 90% complete sample of Lupus YSOs. The results of the veiling measurements are also reported and discussed in the same section. Conclusions are outlined in Sect. 6.

2 Target selection, observations, and data processing

The targets,observations, and data reduction are extensively described in Paper I and Paper II, and are briefly summarized in this section.

2.1 The targets

For the pilot study, we selected six well-known CTTs in Taurus-Auriga with clear signatures of accretion. The sources are listed in Table 1 together with relevant properties from the literature (see also Paper I). Noteworthy is the range of values for both thestellar physical parameters and mass accretion rate, which is likely due to variability, but also to the distinct methodologies used to derive them in previous works. The extreme CTT RW Aur A (Alencar et al. 2005) was included in the sample as a case of highly veiled object spectrum (Herczeg & Hillenbrand 2014) to test the limits of our procedures. In addition, we also study the well-known intermediate-mass CTT CQ Tau. This object was recently investigated in detail with ALMA (Ubeira Gabellini et al. 2019) and has a low level of accretion (e.g., Mendigutía et al. 2011); therefore, it has been included to test the ability of our methods to trace very low accretion rates. However, this object displays a remarkable UX Ori-type variability (Shakhovskoj et al. 2005), which complicates the analysis of the data (Dodin & Suslina 2021). Several of the selected CTTs have been included in the recent ALMA survey by Long et al. (2019). Hereafter, we use the terms CTT and YSO interchangeably to refer to the objects in our sample.

2.2 Observations and data reduction

The GIARPS observing mode combines the HARPS-N (Pepe et al. 2002; Cosentino et al. 2012) and GIANO-B (Oliva et al. 2012; Origlia 2014) high-resolution (resolving power of 115 000 and 50 000, respectively) spectrographs, simultaneously covering a wide spectral range of 390–690 nm for HARPS-N, and 940–2420 nm for GIANO-B. In order to perform a flux calibration as accurately as possible, avoiding additional uncertainties due to variability, the GHOsT survey is complemented with simultaneous/contemporaneous (within less than 2 days of the GIARPS observations) low-resolution spectroscopy in the optical, as well as with optical and NIR photometry (see Paper I).

The TNG/GIARPS observations were performed in 2017 during two nights, one on October 29 and the other on November 13. The journal of the GIARPS observations is reported in Paper I. The reduction steps of the GIARPS spectra are thoroughly described in Paper I and Paper II, but a summary is provided here.

The HARPS-N spectra were reduced using the latest version (November 2013) of the HARPS-N instrument Data Reduction Software and applying the appropriate mask depending on the spectral type of the object (Pepe et al. 2002). The basic processing steps for the data reduction consist in bias and dark subtraction, flat fielding, wavelength calibration, spectrum extraction, and cross-correlation computation. For the removal of the spurious features caused by the telluric lines we first used the molecfit tool (Smette et al. 2015; Kausch et al. 2015) to produce a synthetic telluric spectrum and then the package telluric in IRAF¹ to remove the telluric lines.

In order to flux-calibrate the HARPS-N spectra, the contemporaneous absolute flux-calibrated spectroscopy was acquired within two nights of our GHOsT runs with the 1.22 m telescope at the Asiago observatory, Italy. The spectra cover the wavelength range 330–790 nm and were fully reduced and flux-calibrated against a spectrophotometric standard observed during the same night. Their flux zero-point was also checked with BV R_CI_C photometricmeasurements collected during the same nights with the ANS Collaboration telescopes (see Munari et al. 2012). The photometry is reported in Paper I. Given the short temporal distance between the two datasets we assume that the continuum shape did not change significantly between the Asiago and TNG observations. Thus, for each source we first fitted the continuum of the Asiago spectrum and then multiplied it for the continuum-normalized HARPS-N spectrum.

The GIANO-B spectra were extracted following the data-reduction prescriptions described in Carleo et al. (2018). Halogen lamp exposures were employed to map the order geometry and for flat-field correction, while wavelength calibration was based on lines from a U-Ne lamp acquired at the end of each observing night. The molecfit tool (Smette et al. 2015) was used for removing the telluric contribution in the near infrared.

To flux-calibrate the GIANO-B spectra, NIR photometry in the JHK_S bands was acquired with the REMIR instrument on the Rapid-Eye-Mount (REM) telescope (Vitali et al. 2003), at the La Silla Observatory, Chile, during the night of November 11, 2017. We assume that the NIR magnitudes did not change significantly between TNG, Asiago, and REM observations. We then performed an interpolation between the considered flux measurements using a spline function to derive a smooth continuum function in the wavelength range 940–2420 nm, that is a response function, which was then applied to flux-calibrate the various (continuum-normalized) segments of the GIANO-B spectrum. Finally, the flux-calibrated spectral segments most relevant for this pilot study of accretion were selected.

We estimate that the flux calibration of the GIARPS spectra is precise within 20%. The final products do not contain the region of the Balmer Jump, but include 17 well-resolved accretion diagnostics (12 in the optical and 5 in the NIR) that form the basis for our measurements of accretion in the sample (see Sect. 4).

3 Stellar parameters

A first step in any study of accretion is the determination of the YSO physical properties. Estimates of the physical parameters, derived using a number of methodologies, exist in the literature for all the YSOs in our sample. In order to minimize uncertainties due to the different procedures used in the literature, and in view of the forthcoming GHOsT data, we need to adopt a methodology that allows us to determine the YSO properties in a self-consistent and homogeneous way.

3.1 Effective temperature and optical veiling

We used the ROTFIT code to determine effective temperature, T_eff, surface gravity, log g, radial velocity, RV, projected rotational velocity, v sin i, and veiling as a function of wavelength, r_λ. The code has been successfully applied on YSOs optical spectra at different resolutions (Frasca et al. 2006, 2015), and is particularly well suited for the HARPS-N spectra. In short, the code finds the best photospheric rotationally broadened template spectrum that reproduces the target spectrum by minimizing the χ² of the difference between the observed and template spectra in specific spectral segments. The spectral segments (~ 100Å each) are normalized to unity, and therefore a measure of extinction is not provided by the code. We adopted as templates a grid of high-resolution spectra of slowly rotating, low-activity stars with well-known atmospheric parameters, which were retrieved from the ELODIE Archive (Moultaka et al. 2004). The HARPS-N spectra were degraded to the ELODIE resolution (R = 42 000) before running the analysis code. We have chosen this template grid because real spectra are best suited for the determination of veiling and v sin i (Frasca et al. 2015, 2019). Moreover, they allow us to perform a careful subtraction of the photospheric spectrum (see Sect. 4.1). The spectral intervals analyzed with ROTFIT contain features that are sensitive to the effective temperature and/or log g.

The whole set of physical parameters derived with ROTFIT (T_eff, log g, RV, v sin i, and veiling) for the complete GHOsT sample will be presented in a forthcoming paper, while the effective temperature, essential for this work, is provided in Table 2. Optical veiling values, also crucial for this pilot study, are given in Table 3. The spectral types, determined using the T_eff from ROTFIT and the relation between spectral type versus T_eff by Herczeg & Hillenbrand (2014), are also reported in Col. 3 of Table 2.

For RW Aur A it was not possible to derive the parameters using the ROTFIT method, as the spectrum is totally crossed by emission lines that hide the photospheric lines at optical wavelengths. As discussed in Sect. 3.2, the templates best matching the NIR spectrum have a spectral type between K1 and K0, the latter being consistent with the determination by Herczeg & Hillenbrand (2014) for this star. We therefore adopted the T_eff value as determined from the K0 spectral type and the SpT versus T_eff relationship by these authors. For this star, Herczeg & Hillenbrand (2014) report a veiling of 0.5 at λ =750 nm. Based on a linear relationship between the veiling at 580 nm and 710 nm drawn from our previous studies of the Lupus population (Alcalá et al. 2014, 2017), we estimate that the veiling at λ = 600 nm is r600 ≈ 2.3 × r710, where r710 is the veiling at λ =710 nm. Within the errors, this is consistent with the measured veiling in Frasca et al. (2017) where r620 ≈ 2.0 × r710. For RW Aur A, we calculate a value of r600 ~ 1.2. We find a similar value of veiling considering the equivalent width of the lithium line at 6708 Å and applying the method described in Biazzo et al. (2011). However, we consider this a tentative value only, as the veiling of this CTT may be strongly variable (0.3 < r_λ < 6.1, for λ in the range 500–650 nm, Stout-Batalha et al. 2000).

The corresponding T_eff values derived from ROTFIT are in good agreement, within the errors, with those drawn from the Herczeg & Hillenbrand (2014) spectral types. As a logical consequence, the spectral types are also in good agreement, within a spectral subclass. We note that the values for CQ Tau, which have not been included in the Herczeg & Hillenbrand (2014) study, are in agreement with those published inMendigutía et al. (2011) and Meeus et al. (2012). Also, the negligible veiling at optical wavelengths for this star is in agreement with the recent results by Dodin & Suslina (2021).

3.2 Veiling and v sin i estimates in the NIR

The veiling of IR photospheric lines is usually attributed to the emission of the inner edge of the dusty disk, where the dust is heated by stellar and accretion radiation. As such, it can be used to infer properties of the emission in excess of the stellar photosphereoriginating in the inner disk and complementary to the UV and optical excess. In addition, veiling may have an important impact on extinction (see Fischer et al. 2011; Herczeg & Hillenbrand 2014, and references therein).

The procedure we used to measure veiling of the photospheric lines in the NIR follows the prescription thoroughly described in Antoniucci et al. (2017b). Briefly, we use template spectra acquired with GIANO-B and with similar spectral types to the targets. The templates are chosen to have very low v sin i values. Each template spectrum is broadened by convolution with a rotational profile (Gray 2005) by increasing v sin i and it is artificially veiled by adding a continuumexcess, with both vsini and r_λ as free parameters (see Eq. (1) in Antoniucci et al. 2017b), until the photospheric features match those of the target and minimum residuals between the target spectrum and the rotationally broadened and veiled template are obtained. Figure 1 shows an example of the result of the fitting procedure for HN Tau A in two spectral intervals. The r_λ values are assumed to be constant within each spectral order. Values of r_λ derived from this procedure are given in Table 3. The procedure also yields an estimate of v sin i, which is also given in Table 3. The errors in veiling are on the order of 20%.

As mentioned in the previous section, the optical HARPS-N spectrum of RW Aur A is heavily veiled and crossed by numerous emission lines. Yet its NIR GIANO-B spectrum shows some photospheric lines making it possible to estimate veiling and v sin i. During this process, we were also able to determine that the spectral type of the template best matching the photospheric lines corresponds to a type between K1 and K0. We stress that in the case of DG Tau, the r2130 and r2255 estimates are rather uncertain, and so we were only able to determine a range of values (see Table 3).

Fig. 1 Two portions of the continuum-normalized spectrum of HN Tau A (black solid line) in the HK bands, with the rotationally broadened and veiled spectrum of the stellar template overplotted (red dashed line). Absorption features, mainly of the iron-peak group, are indicated with small vertical lines. We refer to Rayner et al. (2009) for the details on the absorption features. Emission lines such as the H2 line, indicated with the blue dotted lines, are excluded from the analysis. The value of veiling, which minimizes residuals with respectof the rotationally broadened and veiled template, is indicated in the lower left.

3.3 Extinction

As a first step to derive the visual extinction, A_V, we adopt the methods described in our previous works (Alcalá et al. 2014, 2017; Manara et al. 2017a). As input data we use our primary flux calibrated spectra, that is the Asiago spectra, which cover a sufficient wavelength range (330–790 nm) for the purpose. To derive A_V for a given YSO, its spectrum was compared with the spectra of nonaccreting YSOs best matching in spectral type. All our nonaccreting templates were taken from Manara et al. (2013, 2017b) and have very low or negligible extinction (A_V < 0.5 mag; see Manara et al. 2013, 2017b). The templates were rebinned to the Asiago spectra and then artificially reddened by A_V in the range 0.0–4.0 mag in steps of 0.25 mag until the best match to the continuum slope of the YSO spectrum was found. To redden the spectra we used the extinction law by Weingartner & Draine (2001) for R_V = 5.5, which has been found to be particularly well suited for YSOs in general (Evans et al. 2009, and references therein). The A_V values derived in this way are listed in Col. 4 of Table 2.

The main sources of uncertainty on A_V are the errors in spectral type when associating a template to a given YSO, the error in the extinction of the templates, and the errors in flux calibration. The combined effect leads to an error of ~ 0.35 mag, except for the case of RW Aur A for which we estimate a larger error of about 0.5 mag. Adopting extinction laws by different authors yields results consistent, within the errors of about 0.35 mag, with those of Herczeg & Hillenbrand (2014). For instance, the latter authors find that for a star with their measured A_V = 1.0 mag, adopting R_V = 5 from Cardelli et al. (1989) or R_V = 5.1 from Weingartner & Draine (2001) would yield A_V = 1.2 and 1.15 mag, respectively.

As mentioned in Sect. 3.2, another source of uncertainty is strong veiling, which makes the YSO spectra intrinsically bluer than the templates. The extinction as derived above is not severely affected by the veiling, as soon as the analysis is restricted to wavelengths longer than about 450 nm (see analysis in Appendix B of Alcalá et al. 2014), but also shorter than about 800 nm to avoid the effects of veiling in the NIR (Fischer et al. 2011). Thus, to minimize the impact of veiling, the extinction was derived only from the red portion of the Asiago spectra, starting at 500 nm. In this case, the effects of NIR veiling are automatically excluded as the Asiago spectra do not cover wavelengths longer than 790 nm.

We also applied an alternative method to derive A_V using the prescription explained in Fischer et al. (2011, see their Eq. (4)), originally proposed in Gullbring et al. (1998, their Eq. (5)). This requires measurements of the veiling, r_λ, of the observed flux of the spectral template, $F_{\lambda }^T$, and of the observed flux of the object, $F_{\lambda }^O$, as a function of wavelength. The method is based on the fact that the quantity ${\Gamma }_{\lambda }=2.5\cdot \log [(1+r_{\lambda })\cdot F_{\lambda }^T/F_{\lambda }^O]$ is a linear function of the extinction curve, A_λ∕ A_V, and is equal to A_λ ∕A_V ⋅ (A_V^O-A_V^T)-2.5 ⋅ logC, where A_V ^O and A_V ^T are the visual extinction of the object and template, respectively, and C is a constant. The slope of a linear fit to the Γ_λ versus (A_λ ∕A_V) relationship yields the difference between the extinction of the object and the template, assuming the same extinction law for both. In our case, A_V ^T is always very close to zero and therefore the slope of the linear fit yields A_V ^O.

Figure 2 shows the plots of Γ_λ versus A_λ ∕A_V for the YSOs in the sample. For every YSO, we used the r_λ values of Table 3, and the corresponding continuum $F_{\lambda }^T$ and $F_{\lambda }^O$ values reported in Tables A.1–A.7. In these tables, we also report the r_λ values for convenience. The plots show the best linear fit and the corresponding slope A_V^O. The error on A_V indicated in each panel of Fig. 2 is only the error on the slope of the fit and does not represent the full error on A_V, which we estimate to be also on the order of 0.35 mag based mainly on flux calibration errors and mismatch in spectral type of the templates and the CTTs. The A_V results using this method are reported in Col. 5 of Table 2.

In the case of RW Aur A, where the fit is not well constrained because only one rough estimate of veiling in the optical was possible, the error may be as high as 0.5 mag. Nevertheless, the resulting extinction values as derived from both methods are consistent, within the errors. We also stress that the UX Ori-type variability of CQ Tau may lead to large uncertainties on the A_V value, as part of the extinction may be gray (Dodin & Suslina 2021).

The A_V derived from the spectral templates alone and those including the veiling measurements (respectively Cols. 4 and 5 in Table 2) are in very good agreement within the error of about 0.35 mag, although the veiling method provides systematically higher values by up to ~30%. As discussed in Fischer et al. (2011), the adoption of different values of the total-to-selective extinction, R_V, has little impact on the results. We adopt the A_V values derived from the “veiling" method for the subsequent analysis.

The veiling methodologies used in the NIR by Fischer et al. (2011) yielded A_V values systematically higher than those derived from other previous measurements. Yet, applying the same methods here, but extended to the optical veilings, provides A_V values that are significantly lower than those in Fischer et al. (2011) for the same stars. Herczeg & Hillenbrand (2014) also calculate much lower A_V values than Fischer et al. (2011) and interpret the higher values by the latter authors as possibly due to spectral mismatch between the templates and the CTTs and to the lower sensitivity of NIR spectra to extinction. Our extinction values are slightly higher than those of Herczeg & Hillenbrand (2014), but still consistent within the small amount of about 0.5 mag, and in some cases the agreement is very good. This is consistent with the fact that the veiling method gives A_V values that are systematically higher by ~30% with respectto the templates method alone.

3.4 Luminosity and mass

The stellar luminosity, L_⋆, was derived as follows. For a given object, we extracted first the BTsettl model (Allard et al. 2012) best matching the T_eff values derived with ROTFIT above (see Sect. 3.1). The model was then normalized to the extinction corrected flux, $F_{600}^{\text{corr}}$, of the Asiago spectrum at λ = 600 nm, where the veiling estimates are expected to be more reliable, and corrections subject to less uncertainty. To take into account the contribution of veiling, the extinction-corrected flux at 600 nm was multiplied by the factor 1∕(1 + r600), where r600 is the veiling at λ= 600 nm derived from ROTFIT and listed in Table 3 (see example for DL Tau in Fig. 3). This factor, which is always ≤1, effectively reduces the observed flux to that emitted by the stellar photosphere only, excluding the emission due to accretion.

We therefore assume that the BT-settl model normalized in this way best represents the spectral energy distribution (SED) of the object’s photosphere at the star distance. Integration of the normalized model at all wavelengths yields the bolometric flux corresponding to the object’s photosphere. The stellar luminosity was then calculated using this flux and adopting the distance reported in Table 1. The main sources of uncertainty on L_⋆ are the errors in flux calibration of the Asiago spectra and the error on the veiling correction. On this basis, we estimate an average uncertainty of about 0.2 dex in log L_⋆. We verified that using veiling values at other wavelengths yields consistent results of the veiling-corrected bolometric flux. For instance, in the case of DL Tau, which has the strongest veiling at 450 nm, we derive a corrected bolometric flux of 5 × 10⁻¹⁰ erg s⁻¹ cm⁻² and 4.55 × 10⁻¹⁰ erg s⁻¹ cm⁻² when using the r600 and r450 values, respectively. The stellar radius, R_⋆, was calculated from the effective temperature and stellar luminosity.

In some cases (e.g., DG Tau) the L_⋆ values from the oldest literature may be overestimated most likely because such determinations did not consider thecontribution of veiling and/or overestimate extinction. Our stellar luminosities are well consistent with those derived in Herczeg & Hillenbrand (2014) within the error of about 0.2 dex in log L_⋆. It is worth noting that in the case of RY Tau our T_eff, A_V, and L_⋆ values are in good agreement, within the errors, with those reported in the recent and thorough study by Garufi et al. (2019). These values are also fairly consistent with those reported in the paper by Calvet et al. (2004).

We note that HN Tau A is subluminous on the HR diagram with both the value of Herczeg & Hillenbrand (2014) for L_⋆ and the one derived here. This is likely due to obscuration of the stellar photosphere by the highly inclined disk (i ≈70°; Long et al. 2019). The L_⋆, R_⋆, and M_⋆ values for this star may be underestimated by a factor of ~17, ~ 4.3, and ~ 2, respectively (see Sect. 4.4). Likewise, our derived luminosity for CQ Tau places the star slightly below the main sequence, most likely as a consequence of the fact that during the periods of dimming, the extinction may be gray, and therefore the luminosity is underestimated. The GHOsT observations of this star were performed when it was in a faint stage, that is, B = 10.89 mag as compared with B = 9 mag in its bright phase (Grinin et al. 2008), and therefore its luminosity may be underestimated by a factor of about six (see also Sect. 4.4).

Finally, the mass, M_⋆, of the sevenCTTs was estimated by comparison of the position of the objects on the HR diagram with the theoretical PMS evolutionary tracks by Siess et al. (2000). The uncertainties on L_⋆ and T_eff lead to a typical uncertainty of ~0.15 dex in M_⋆. The derived L_⋆, R_⋆, and M_⋆ values for the sample are given in Table 2. The mass reported in this table for the subluminous objects HN Tau A and CQ Tau corresponds to that of the evolutionary track closest to these stars on the HR diagram, but corrected values are estimated in Sect. 4.4.

Fig. 2 $2.5\times \log [(1+r_{\lambda })\cdot F_{\lambda }^T/F_{\lambda }^O]$ vs. Aλ ∕AV for the YSOs inthe sample (blue squares). The best linear fit for each object is shown as a black line. The value of the slope AV is shown in each panel. In the case of CQ Tau the r1565 value is shown as an upper limit.

Fig. 3 Example of BT-settl model normalization. The extinction-corrected Asiago spectrum of DL Tau is shown with the red line, while the BT-settl model normalized to DL Tau at λ = 600 nm is plotted with the blue line. We note the increase of excess emission of DL Tau at wavelengths shorter than about 450 nm. Theveiling-corrected BT-settl model is plotted with the black line. See text for details

4 Accretion diagnostics

The GIARPS spectra include several emission lines that are well correlated with L_acc (Alcalá et al. 2014, 2017) and will be the basis of our measurements of accretion in the sample. In particular, we extracted 17 spectral portions with well-resolved and flux-calibrated accretion diagnostics, namely, eight hydrogen recombination lines (Hα, Hβ, Hγ, Hδ, Paβ, Paγ, Paδ, Pa8), eight helium lines (He I4026, He I4471, He I4713, He I4922, He I5016, He I5876, He I6678, He I10830), and the Ca II3934 line.

Figure 4 shows examples of GIARPS spectra for several emission lines of every CTT star in our sample. For comparison, the spectra are normalized to the local continuum and shifted in velocity to the rest wavelength. The latter was determined from the profiles of the Li I λ6708Å photospheric line, assuming weighted λ = 670.7876 nm, and in the NIR from the Al Ilines at λ = 2019.884, 2116.958, 2121.396 nm. All the spectra are shown in the same velocity range for comparison.

There isa variety of line profiles with a range of widths and intensities, typical of accreting YSOs. The widest lines are from RY Tau, RW Aur A, and HN Tau A, while DG Tau and DL Tau show lines with intermediate width. The narrowest lines are from DO Tau, yet the width of the Hα line at 10%of the line peak that we measure for the latter is more than 400 km s⁻¹, pointing to significant accretion activity (White & Basri 2003; Natta et al. 2006). Within the qualitative classification scheme of the Balmer lines proposed by Antoniucci et al. (2017a), RY Tau, RW Aur A, and HN Tau A display wide and multi-peak profiles, DG Tau andDL Tau show multi-peak profiles, while DO Tau shows narrow and almost symmetric profiles. All these morphologies are displayed by the Lupus sample studied by Antoniucci et al. (2017a). Two exceptions in our sample are RY Tau and CQ Tau, where many lines are dominated by the photospheric absorption and therefore the subtraction of the photospheric spectrum is needed in order to measure the residual emission (see Sect. 4.1).

4.1 Flux and luminosity of lines

The flux at Earth in permitted lines was computed by directly integrating the GIARPS flux-calibrated spectra using the splot package under IRAF. Three independent measurements per line were made, considering the lowest, highest, and the middle position of the local continuum, depending on the local noise level of the spectra. The flux and its error were then computed as the average and standard deviation of the three independent measurements, respectively.

In the case of RY Tau and CQ Tau, we performed the photospheric subtraction in the following way: the spectral templates used to derive the stellar parameters and NIR veiling were artificially broadened for v sin i and veiled at the same values derived from ROTFIT and the NIR veiling procedure. The rotated and veiled templates were then subtracted from the YSO spectra. In this procedure, both the YSO and template spectra were normalized to the continuum around the selected lines to be analyzed. Three measurements of the equivalent widths of the corrected lines were made and the average was computed. The error was estimated by the standard deviation of the observed fluxes on the difference spectra in two spectral regions near the line. Finally, the flux of the corrected lines was computed as the product of the equivalent width times the absolute continuum flux adjacent to the lines. The procedure, illustrated in Fig. 5 for the Hα and Hβ lines and in Fig. 6 for the Paβ line, was applied to lines strongly affected by the photospheric contribution namely the Balmer, Paschen, and the He I λ10 830Å lines. We note that in some cases a residual contribution of nonphotospheric self-absorptions remained, in particular in the case of some Paschen and the He I λ10 830Å lines (see example in Fig. 6). Lines with strong self-absorption may underestimate L_acc and therefore the measurement was not considered when such a contribution was more than ~ 30% of the flux. The procedure adopted here to measure the integrated flux is the same adopted in Alcalá et al. (2017), where the L_line versus L_acc calibrations were derived.

Not all the 17 lines in the GIARPS spectral range were detected or could be measured in every CTT star. A summary of the number of lines with flux measurements in every CTT in our sample is provided in the first column of Table 4, and every panel of Fig. 7 shows the lines detected and measured. The absence of a point in these plots means that the line was not detected or could not be extracted with the photospheric subtraction, mainly because of large residuals of telluric lines and/or a strong contribution of the nonphotospheric self-absorption component. DL Tau is the only star for which the complete set of lines could be detected and measured, while CQ Tau has the lowest number of lines with measured fluxes. For most stars, fluxes have been measured for more than 12 lines. In the case of RW Aur A, the Balmer lines higher than H4 (Hγ, Hδ) are dominated by a nonphotospheric self-absorption component. Therefore, we did not attempt a measurement of the flux of those lines. Also, the Hα and Hβ lines in this star may be affected by a similar self-absorption, but we measured their flux with the awareness that it may be underestimated. A contribution of absorption components in the emission lines of DG Tau, DL Tau, HN Tau A, and DO Tau is not significant, and therefore we use the flux measured directly by the integration of the lines. The observed fluxes, equivalent widths, and their errors are reported in several tables provided in Appendix B (Tables B.1–B.5)².

The luminosity of the different emission lines was computed as L_line = 4πd² ⋅ f_line, where d is the YSO distance listed in Table 1 and f_line is the extinction-corrected flux of the lines.

Fig. 4 Examples of GIARPS spectra showing the line profiles of CTTs in our sample.

Fig. 5 Examples of photospheric line subtraction for the Hα and Hβ lines in CQ Tau and RY Tau. The black, red, and blue lines represent the observed, photospheric, and photospheric-subtracted spectra, respectively. See text for details.

Fig. 6 Examples of photospheric line subtraction for the Paβ line for CQ Tau (upper panels) and RY Tau (lower panels). The red lines represent the photospheric templates. The left and right panels show, in black lines, the observed and photospheric-subtracted spectra, respectively. The green vertical lines mark the intervals for the equivalent width measurement. The equivalent widths of the emission line before and after the subtraction are also indicated.

4.2 Accretion luminosity

We derived L_acc via empirical relationships between accretion luminosity and the luminosity of permitted emission lines. Such relationships have been derived by several authors (e.g., Herczeg & Hillenbrand 2008; Rigliaco et al. 2012; Alcalá et al. 2014, 2017). Here, we used the most recent ones by Alcalá et al. (2017), where the relationships are simultaneously derived for lines from the UV to the NIR for a more than 90% complete sample of Lupus YSOs.

In Fig. 7, the derived L_acc values are plotted as a function of the line diagnostics for the seven targets. The error bars include the errors in L_line, as well as the errors in each L_acc–L_line calibration. The individual L_acc values corresponding to each diagnostic are reported in the various tables in Appendix B (Tables B.1–B.5) for every CTT in the sample. An average accretion luminosity, ⟨L_acc⟩, was then calculated for each target. These values are reported in Table 4 together with the corresponding ⟨L_acc⟩/L_⋆ ratio.

The typical standard deviation of ~0.25 dex on log ⟨L_acc⟩ is within the expected error estimated from the fit of the UV continuum excess emission in other samples using slab models (see Alcalá et al. 2014, 2017; Manara et al. 2017a), although the error on L_acc for the individual diagnostics is larger than this value. This confirms that the average L_acc, derived from several diagnostics measured simultaneously, has a significantly reduced error.

4.3 Mass accretion rate

The average accretion luminosity of each target (see Table 4) was converted into mass accretion rate, Ṁ_acc, using the relation $${\dot{M}}_{\text{acc}}={\left(1-\frac{R_{\star }}{R_{\text{in}}}\right)}^{-1} \frac{\langle L_{\text{acc}}\rangle R_{\star }}{G{M}_{\star }}\approx 1.25 \frac{\langle L_{\text{acc}}\rangle R_{\star }}{G{M}_{\star }},$$(1)

assuming $\frac{R_{\star }}{R_{\text{in}}}=\frac{1}{5}$, where R_⋆ and R_in are the stellar radius and inner-disk radius, respectively (see Gullbring et al. 1998; Hartmann 1998), and using the stellar parameters reported in Table 2. The resulting Ṁ_acc values are reportedin Table 4.

The sources of error in Ṁ_acc are the uncertainties on L_acc, stellar mass, and radius. Propagating these, we estimate an average error of ~ 0.35 dex in Ṁ_acc (see Appendix A of Alcalá et al. 2017). However, additional errors on these quantities come from the uncertainty in distance, as well as from differences in the adopted evolutionary tracks. The largest uncertainty on the YSOs distance from the Gaia EDR3 is estimated to be less than ~ 10%, yielding relative uncertainty of ~0.2 dex in the mass accretion rate³. On the other hand, it has been shown (see Appendix A in Alcalá et al. 2017) that using different sets of PMS evolutionary tracks leads to uncertainties of 0.04 dex to 0.3 dex in Ṁ_acc. We therefore estimate that the cumulative relative uncertainty in Ṁ_acc is about 0.4 dex.

We note that our L_acc = 0.08L_⊙ for CQ Tau is more than an order of magnitude lower than the value reported in Donehew & Brittain (2011, L_acc = 3.8 L_⊙; see also Tables 1 and 4 for the Ṁ_acc estimates). Yet, using the Brγ luminosity reported bythese authors and the Alcalá et al. (2017) relationships we derive a L_acc = 0.15 L_⊙, which is only a factor two greater than our value and within the limits of long-timescale variable accretion (< 0.4 dex, see Costigan et al. 2012, 2014). A correction estimate due to UX Ori-type variability yields a L_acc ~ 0.38L_⊙ for CQ Tau (see Sect. 4.4), which is still lower by an order of magnitude than the Donehew & Brittain (2011) determination. It is worth noting that our estimate of log Ṁ_acc = −8.68 for this star is consistent with the upper limit derived by Mendigutía et al. (2011, log Ṁ_acc < −8.3), suggesting that these latter authors also observed the star in its faint phase.

We also stress that the L_acc values derived here for RY Tau, DL Tau, HN Tau A, and DO Tau are consistent, within the errors, with those by Ingleby et al (2009) and Ingleby et al. (2012), although in the case of RY Tau our value is a factor of about four lower. We think this may be due to variable accretion. It is also worth mentioning that our Ṁ_acc estimates for RY Tau and DG Tau are in good agreement with the range of values derived by Frasca et al. (2018) based on the Hα and He I6678 emission lines.

Fig. 7 Lacc as a function of the different accretion diagnostic as labeled for the seven targets (blue dots). The vertical error bars consider the error in Lline as well as the errors in the Lacc–Lline relationships. The horizontal red lines in each panel represent the average Lacc.

4.4 The cases of HN Tau A and CQ Tau

The highly inclined disk of HN Tau A may occult, at least partially, the emission from the accretion flows and from the shock onto the stellar surface, and therefore L_acc for this star will be underestimated in a similar way to the L_⋆ value (see Sect. 3.4). Nevertheless, the L_acc/L_⋆ ratio for this star is at the level of the highest accretors, possibly showing that both L_acc and L_⋆ of the star are obscured in the same manner. This interpretation has proven to be correct for subluminous objects in Lupus (see Sect. 7.4 in Alcalá et al. 2014).

Following the same reasoning as in Alcalá et al. (2020, their Appendix C), we can use the luminosity of the [O I ] λ6300 line to estimate a correction factor on L_acc and L_⋆. The [O I] line is found to originate relatively far from the star and therefore should not be significantly affected by obscuration effects from the inner disk, and the line luminosity is also correlated with L_acc and L_⋆ (Natta et al. 2014; Nisini et al. 2018). From the HARPS-N spectrum of HN Tau A, we measure a line flux F_[O I] = 9 × 10⁻¹⁴ erg s⁻¹ cm⁻² and derive the extinction corrected (A_V = 1.53 mag) flux F${}_{[\text{ion}Oi]}^{\text{corr}}=3\times {10}^{-13}$ erg s⁻¹ cm⁻². From Paper I we estimate that about one-third of the flux comes from the low-velocity component (LVC), and therefore we derive a line luminosity, in solar units, of $\log L_{[\text{ion}Oi]}^{LVC}=-4.2$. Using the L_acc-L_line and L_acc-L_⋆ relationships for the LVC⁴ derived by Nisini et al. (2018) we estimate the log L_acc and log L_⋆ values given in Table 5, which are a factor of approximately 7 and 17 greater than the subluminous values, respectively, meaning that both L_acc and L_⋆ are affected by a similar obscuration factor, which is consistent with the interpretation for subluminous objects (Alcalá et al. 2014). For comparison, we also include the measured values in Table 5. The corrected stellar luminosity implies a stellar radius a factor about 4.3 larger (i.e., 2.6 R_⊙), and places the star on the HR diagram in a position consistent with the other CTTs.

Using the corrected L_⋆ value and the Siess et al. (2000) evolutionary tracks, we estimate a mass of 1.6 M_⊙ for HN Tau A, which is a factor of about two higher than the value provided in Table 2. The log g = 3.8 calculated from the corrected radius and mass is consistent with the typical gravity for YSOs, while the subluminous values yield a much higher gravity of logg = 4.7. The corrected L_acc, R_⋆, and M_⋆ values and Eq. (1) would imply a log Ṁ_acc ≈ −7.3, that is, HN Tau A would be among the strongest solar-mass accretors in Taurus.

The work by Nisini et al. (2018) also provides L_line–M_⋆ and L_line–Ṁ_acc relationships allowing us to estimate M_⋆ and Ṁ_acc indirectly from the $\log L_{[\text{O\hspace{0.17em}I}]}^{LVC}$ value. The results are also provided in Table 5 and are in agreement with the values derived from the evolutionarytracks and Eq. (1).

We warn the reader about the uncertainties on extinction in high-inclination objects. The above calculations for HN Tau A in this section are based on the assumption that the visual extinction of the LVC of the [O I ] line is the same as measured for the star, i.e., A_V = 1.53 mag, which is not necessarily true. However, we note that adopting A_V = 0 mag yields a corrected L_⋆ a factor of about three lower than when assuming A_V = 1.53 mag, which is still underluminous on the HR diagram. On the other hand, a much higher value of A_V would make the star unreliably luminous. We therefore conclude that a reasonable value for the visual extinction of the LVC of the [O I ] line should be in the range from about 0.7 mag to 2 mag. As explained above, adopting A_V = 1.53 mag yields a corrected L_⋆ = 2.6 L_⊙, which leads to consistent results on the stellar parameters. We therefore used this value, but warn the reader that the genuine L_⋆ value might be in the range between ~ 1.5 L_⊙ and ~ 3.5 L_⊙.

For CQ Tau, the [O I] line is barely seen, and therefore we cannot use the above methodologies to correct the stellar and accretion parameters for obscuration effects. However, we can estimate upper limits based on the photometric variations and assuming gray extinction. As pointed out by Dodin & Suslina (2021), CQ Tau did not show variations before 1940 and was approximately constant at B ~ 9 mag (Grinin et al. 2008). Adopting this value as the unobscured magnitude of the star and based on our B = 10.9 mag measurement during the GHOsT observations, we estimate a correction factor on the bolometric flux of the star of ~ 5.7. This correction provides a maximum bolometric flux, yielding upper limits of 15 L_⊙, 2 M_⊙ and 2.8 R_⊙ for L_⋆, M_⋆, and R_⋆, respectively. These corrected values of R_⋆ and M_⋆ yield a log g = 3.9, which is more consistent with the gravity of a YSO than the logg = 4.5 derived from the observed R_⋆ and M_⋆ values in Table 2. Assuming that L_acc is affected by the same obscuration factor as L_⋆ (see Sect. 7.4 in Alcalá et al. 2014), we derive an upper limit of − 7.7 for log Ṁ_acc. The estimated upper limits are listed in Table 5.

5 Results and discussion

The results of the previous sections show that the studied CTTs in Taurus are highly accreting objects. The novel science products and aspects of this pilot study are as follows.

• Contemporaneous low-resolution spectroscopic and photometric observations, allowing an accurate flux calibration of the high-resolution spectroscopy;

• simultaneity of the high-resolution, wide-band spectroscopic observations, from the optical to the NIR;

• simultaneous use of veiling measurements, both in the optical and NIR, to determine A_V;

• use of more than ten line diagnostics to estimate accretion luminosity, yielding a much reduced error in L_acc as compared with determinations using single diagnostics.

All of the above was achieved using well-defined and assessed procedures for deriving the stellar physical and accretion parameters in a self-consistent and homogeneous way. Therefore, the properties derived here can be considered as more robust and reliable than in previous studies. Some limitations to the application of our procedures are related to the extremely veiled CTTs like RW Aur A, but this type of object is not common. We are therefore confident that the same procedures can be applied to the entire GHOsT data sample, which will be presented in forthcoming papers (Gangi et al., in prep.). In the following, we discuss a few aspects of the stellar and accretion properties of our sample, as well as of the continuum excess emission in the NIR of these CTTs.

5.1 Stellar and accretion properties

This GIARPS/TNG pilot study confirms the high levels of accretion of the selected CTTs. It is worth noting from Table 4 that the L_acc/L_⋆ ratio for every CTT in this sample is well above the level of chromospheric noise emission in YSOs (max (L_acc/L_⋆)_noise ≈ 0.01 Manara et al. 2017b). We point out that in three sources, namely, DG Tau, DL Tau, and RW Aur A, the total luminosity is accretion dominated (i.e., L_acc/L_⋆ > 1), which is more typical of the class I stage of evolution (e.g., Antoniucci et al. 2008, Fiorellino et al. 2021).

The most actively accreting object in this sample is RW Aur A. Noteworthy, the recent study by Takami et al. (2020) shows that the star was in its bright stage during the GIARPS observations (13 November 2017). Moreover, from Paper I we know that the magnitude of the star was 10.44 mag in the V -band and this is consistent with the high-accretion activity scenario by Takami et al. (2020).

The apparently least active object in the sample is CQ Tau. This might be mostly related to the much higher photospheric flux in comparison with the other objects, as CQ Tau is an intermediate-mass YSO with T_eff ~ 6800 K. Nevertheless, its Ṁ_acc is comparable with that of many actively accreting stars of similar mass in other star forming regions. Confirming the upper limits derived in Sect. 4.4 as true values would suggest CQ Tau has a similar level of accretion activity to RY Tau. On the other hand, if the determination of L_acc based on the [O I] λ6300 line for HN Tau A is correct (see Sect. 4.4), this would also be one of the most actively accreting objects in Taurus.

To investigate the levels of accretion in the Taurus subsample in more detail, we use in the following the 90% complete Lupus sample studied in Alcalá et al. (2017) for comparison purposes, but using the rescaled values reported in Appendix A of Alcalá et al. (2019). Given the very small number statistics of this pilot study, we cannot provide any result on the sample as a whole, but only at the level of the individual objects.

Fig. 8 Accretion luminosity as a function of stellar luminosity for the stars in Lupus (black symbols) and the seven Taurus CTTs studied here (blue circles). The latter are labeled with their names. The Lupus transitional disks are shown with crossed squares. The continuous lines represent the three Lacc vs. L⋆ relations as labeled. The long-dashed line represents the shift of HN Tau A on the diagram when correcting its Lacc and L⋆ values for obscuration effects by the disk. The corrected values are shown with the red dot. The leftward and downward red arrows represent the upper limits on L⋆ and Lacc for CQ Tau, respectively. The average errors for the Taurus sample are shown in the upper left. Figure adapted from Alcalá et al. (2017).

Fig. 9 Mass accretion rate as a function of stellar mass for the seven Taurus CTTs studied (blue circles) compared with stars in Lupus (black symbols). The Lupus transitional disks are shown with crossed squares. Lupus objects classified as weak or negligible accretors are plotted with downward arrows. The long-dashed black line represents the shift of HN Tau A on the diagram when correcting its Ṁacc and M⋆ values for obscuration effects by the disk. The corrected values are shown with the red dot. The leftward and downward red arrows represent the upper limits on M⋆ and Ṁacc for CQ Tau, respectively. The black dashed line shows the double power law theoretically predicted by Vorobyov & Basu (2009), and the continuous magenta lines represent the fits to the data as in Eqs. (4) and (5) of Alcalá et al. (2017). The average errors for the Taurus sample are shown in the upper left. Figure adapted from Alcalá et al. (2019).

5.1.1 Accretion luminosity versus stellar luminosity

Figure 8 shows the accretion luminosity as a function of stellar luminosity for the Lupus sample and the seven Taurus CTTs. While most Lupus objects lie at L_acc/L_⋆ values between 0.01 and 0.1, five of the Taurus CTTs have higher values, between 0.3 and 1.5 (see Table 4). The two CTTs with values compatible with most of the Lupus YSOs are RY Tau and CQ Tau. We note that recent ALMA data have shown that the dusty disk of the latter has a large cavity (R_cavity =53 au, Ubeira Gabellini et al. 2019), likely suggesting a transitional disk. Also, in the case of RY Tau, submillimeter data show evidence for a protoplanetary disk with structures such as rings and an inner cavity (Isella et al. 2010; Pinilla et al. 2018; Long et al. 2018). The L_acc/L_⋆ ratio for both CQ Tau and RY Tau is compatible with those of other transitional disks in Lupus.

Tilling et al. (2008) presented simplified stellar evolution calculations for stars subject to time-dependent accretion history, and derived evolutionary tracks on the L_acc–L_⋆ diagram for a variety of fractional disk masses, f_disk ≡ M_disk/M_⋆, and YSO masses. Alcalá et al. (2014) showed that the relation of the Lupus data in Fig. 8 is steeper than the L_acc/L_⋆ = constant lines, more or less following the slope of the Tilling et al. (2008) tracks. These latter authors also concluded that the disks of the Lupus objects with the lowest masses should have masses lower than 0.014× M_⋆. The five high accretors in Taurus (DG, DL, HN, DO, and RW) fall, instead, in the region of the diagram consistent with the Tilling et al. (2008) tracks for YSOs of one solar mass and f_disk = 0.2. These five stars have indeed masses on the order of 1M_⊙ hence, according to such tracks one would expect their disk mass to be on the order of 0.2 ×M_⋆. Although a correlation between the total mass of gas+dust in the disk and the stellar mass has not yet been confirmed, correlations between the dust mass in the disk and the stellar mass of YSOs in the Lupus and Chameleon star forming regions have been found (Ansdell et al. 2016; Pascucci et al. 2016).

In this scenario, the most massive disk would be RY Tau. The presence of the important jet in RY Tau was recently explained in terms of a very massive disk around this intermediate-mass T Tauri star (Garufi et al. 2019), which might still be fed by the interstellar matter in which the object is embedded. For DG Tau, Podio et al. (2013) estimate a total disk mass in the range 0.015–0.1 M_⊙, depending on the assumed dust size distribution. The upper limit would be consistent with the results from the Tilling et al. (2008) tracks. We note that in the case of HN Tau A, the L_acc and L_⋆ values corrected for obscuration effects would still yield similar results for the estimated fractional disk mass, while confirmation of the upper limits for CQ Tau would imply a disk mass similar to that of RY Tau.

5.1.2 Mass accretion rate versus stellar mass

The relationship between mass accretion rate and stellar mass is a fundamental aspect of the study of disk evolution in YSOs. During the CTT phase, that is, after the protostar has almost entirely dispersed its envelope but is still actively accreting from the optically thick accretion disk, the stellar mass undergoes negligible changes. Therefore, the Ṁ_acc versus M_⋆ relation represents a diagnostic tool for the evolution of Ṁ_acc (Clarke & Pringle 2006) and for the process driving disk evolution (Ercolano et al. 2017).

Figure 9 shows the accretion rate as a function of stellar mass for the seven CTTs in Taurus as compared with the Lupus sample. The Taurus CTTs populate the upper right part of the diagram. Most of the Lupus YSOs fall well below the theoretical prediction by Vorobyov & Basu (2009, short-dashed black line), but the latter is relatively consistent with the upper envelope of the Lupus data points distribution. The few Lupus YSOs on the upper envelope of the Ṁ_acc–M_⋆ relationship are the strongest Lupus accretors at a given mass and are also among the most luminous on the HR diagram. Interestingly, the five Taurus CTTs with the highest L_acc/L_⋆ ratios tend to follow the upper envelope⁵, which qualitatively is well fitted by the theoretical prediction. These results demonstrate that the level of accretion of these five CTTs is as high as that of the strongest accretors in Lupus (Alcalá et al. 2017).

On the other hand, CQ Tau falls in the lower envelope of the Ṁ_acc–M_⋆ relationship, with its accretion properties more closely resembling those of transitional disks than those of the full disks, while RY Tau follows the Ṁ_acc–M_⋆ trend for the mostmassive Lupus stars with full disks. Should the upper limits for CQ Tau be confirmed, its accretion properties will be similar to those of RY Tau.

5.2 The continuum NIR excess emission

Clear evidence of continuum excess emission or veiling increasing with wavelength, from the optical to the NIR, was found in CTTs over adecade ago (Fischer et al. 2011, and references therein). As mentioned in Sect. 3.2, the behavior of IR veiling with wavelength can provide information on the physical properties of the dust in the inner edge of the disk.

Figure 10 shows the NIR veiling as a function of wavelength for the seven CTTs in the sample. The increase of veiling with wavelength in the NIR is evident in all the objects. Adopting the same strategies as in Antoniucci et al. (2017b), a fit of r_λ versus λ with a black-body power law yields the temperature, $T_{\text{eff}}^{\text{BB}}$, for each CTT indicated in the corresponding panel of Fig. 10 and in Table 6. We also include in this table the results for DR Tau and XZ Tau from Antoniucci et al. (2017b). Typical errors on the $T_{\text{eff}}^{\text{BB}}$ are on the order of 100–150 K. The fitting procedure requires another parameter, $F_{factor}^{BB}$, namely the factor by which the black-body must be multiplied to fit the data. Such a factor is based on the ratio of the areas of the stellar disk and the emitting region producing the NIR veiling. The results using this factor are also given in Table 6. We warn the reader that the fit for DG Tau must be taken with care because only a range of values for the two reddest points could be determined (see Sect. 3.2).

While the resulting temperatures for RY Tau, DG Tau, DL Tau, and within errors DO Tau, are consistent with the origin of the NIR continuum excess emission being the inner rim of the dusty disk, the derived $T_{\text{eff}}^{\text{BB}}$ value for the other two CTTs is significantly higher than the dust sublimation temperature (~ 1500 K) and is more difficult to interpret in the same way, although temperatures as high as 2000 K may be expected for the sublimation of silicate dust (Pollack et al. 1994).

In the simple model by Dullemond et al. (2001), where the inner disk is directly irradiated by the central star, the inner disk edge is located at a radius, R_sublim, which is given by the following equation: $$R_{\text{sublim}}=\sqrt{(1+f)\left(\frac{L_{\star }+L_{\text{acc}}}{4\pi \sigma T_{\text{sublim}}^4}\right)},$$(2)

where f is the ratio of the inner edge height to its radius, and is estimated to be 0.1 for T Tauri stars. Assuming this model and using the L_⋆ and L_acc results of the previous sections, and adopting $T_{\text{eff}}^{\text{BB}}$ in Table 6 as sublimation temperature, we computed the R_sublim values listed in Table 6 for our sample of CTTs in Taurus. The error in R_sublim was calculated by error propagation in Eq. (2) and adopting errors of 0.2 dex, 0.25 dex, and 150 K on log L_⋆ logL_acc, and $T_{\text{eff}}^{\text{BB}}$, respectively. We can compare these values with the inner radius, $R_{\text{in}}^{\text{cont}}$, derived from interferometric observations of the continuum NIR emission by Eisner et al. (2010) for a few of the CTTs in common. These values are listed in the last column of Table 6. We note that except for RY Tau, where R_sublim and $R_{in}^{\text{cont}}$ are in very good agreement, the R_sublim values for the other four stars with interferometric observations are lower than the $R_{\text{in}}^{\text{cont}}$ measurements. Interestingly, the upper limits on L_⋆ and L_acc for CQ Tau would yield a R_sublim value more similar to those typically found using interferometric observations (see Table 6).

Previous works (i.e., Eisner et al. 2007; Anthonioz et al. 2015) have shown that the irradiated disk model for CTTs predicts radii for the inner rim of the dusty disk that are underestimated with respect to interferometric measurements. We therefore expect the R_sublim values in Table 6 for DL Tau, HN Tau A, DO Tau, CQ Tau, and XZ Tau to be underestimated by a factor of between about two and three. However, one possibility is that the $T_{\text{eff}}^{\text{BB}}$ we measure is not the sublimation temperature, but being instead related to the excess emission in the NIR, may be the temperature of hot gas inside the sublimation radius.

Several authors (e.g., Fischer et al. 2011) discussed possible scenarios for the NIR continuum emission. Among the suggested contributions that may enhance the NIR continuum excess at the shorter IY J wavelengths there are warm annuli around accretion hot spots, hot gas inside the dust sublimation radius, and hot gas in the accretion flows and/or winds. Fischer et al. (2011) modeled and discussed the presence of multiple temperature components in the shocked photosphere. Also, detailed multiple-component modeling of the accretion emission has been performed by Ingleby et al. (2013).

The results in Paper I showed that DG Tau and HN Tau A share similar physical conditions in their jets. The temperature and ionization gradients of the jets in these two objects would favor a magneto-hydrodynamical shock heating in which the warm and ionized streamlines originate in the internal and mainly gaseous disk, while the low-velocity and almost neutral stream lines come from the dusty regions of the outer disk. Therefore, the hypothesis of a gaseous hot disk inside the dust sublimation radius in HN Tau A and the similarity of the physical conditions of its jet with those in the DG Tau jet would also favor a much higher accretion rate in HN Tau A than what is measured without obscuration correction. The possibility for the formation of the jets in DO Tau and RW Aur A in a gaseous inner disk was also explored in Paper I. In this context, jets signify the presence of dense hot regions inside the gaseous inner disk.

In conclusion, it may be possible that some of the NIR continuum excess emission of the CTTs, in particular in the case of HN Tau A, RW Aur A, and DR Tau, originates in a thick gaseous disk inside the dust sublimation radius, as already suggested by Fischer et al. (2011) and Antoniucci et al. (2017b) for other similar cases. In fact, the size of an emitting region estimated as $R_{\text{emiss}}^{\text{cont}}~$ R_⋆ $\cdot \sqrt{F_{\text{factor}}^{BB}}$ is always smaller than the typical innerrim radius of ~0.1 au in CTTs measured from interferometric observations in the NIR. This is consistent with the results by Koutoulaki et al. (2019) who detected ro-vibrational emission of CO in the NIR X-Shooter spectrum of RW Aur A. The modeling presented by these latter authors of five band-heads shows that the CO emission comes from a region at a distance of ~ 0.06–0.08 au from the star.

Fig. 10 NIR veiling as a function of wavelength for the stars in our sample (blue dots). The red dashed lines show the black body fits to the data with $T_{\text{eff}}^{\text{BB}}$ as labeled. See text for details.

6 Conclusions

In this pilot study, we report GIARPS@TNG high-resolution observations of seven CTTs in the Taurus-Auriga star forming region, namely RY Tau, DG Tau, DL Tau, HN Tau A, DO Tau, RW Aur A, and CQ Tau. The spectra simultaneously cover a wide spectral range from the optical to the NIR. Contemporaneous spectrophotometric and photometric observations were also performed, allowing us to flux-calibrate the high-resolution spectra with an estimated accuracy of < 20%. We show that the GIARPS@TNG data, together with the ancillary observations, allows the stellar and accretion parameters of YSOs to be derived in a self-consistent and homogeneous way.

High-resolution spectroscopy was used to derive the veiling throughout the wide spectral range. The impact of veiling on the estimates of extinction was accounted for. Deriving extinction on the basis of well-flux-calibrated spectral templates alone may lead to underestimation of A_V by up to ~30%, if veiling is neglected. Best matching the spectral type of templates and YSOs reduces errors in the A_V estimates. Simultaneously deriving veiling, stellar parameters, and vsini avoids degeneracy, which in turn allows more accurate stellar parameters to be obtained.

A large number of emission line diagnostics were used to calculate the accretion luminosity in the seven CTTs via accretion luminosity versus line luminosity relationships. We confirmed that the average L_acc derived from several diagnostics measured simultaneously has a significantly reduced error.

We therefore conclude that the GHOsT data sets and the procedures adopted here yield more robust results on the stellar and accretion parameters than those in previous studies of Taurus CTTs. However, in the case of extremely veiled objects, our procedures may fail or provide uncertain results, although this type of object is expected to be rare.

Assuming magnetospheric accretion, we calculated the mass accretion rate of the CTTs in the sample, confirming the high-levels of accretion in these objects. Being in its bright stage during the GIARPS observations, RW Aur A is the most actively accreting object in this sample, which is consistent with the high-accretion activity scenario by Takami et al. (2020). The apparently least active objects in the sample are CQ Tau and HN Tau. We identify these two objects as subluminous on the HR diagram, the former because of its UX Ori-type variability and the latter because of the high inclination of its disk. Correction for disk obscuration makes HN Tau one of the most actively accreting objects in the sample, at a level close to RW Aur A.

A comparison of the accretion properties of the Taurus CTTs with those of Lupus YSOs yields the following results: the Taurus CTTs have values of L_acc/L_⋆ of between 0.3 and 1.5, which are normally higher than those of the Lupus YSOs. The two Taurus objects classified as transitional disks, namely RY Tau and CQ Tau, have values similar to those of the Lupus transitional disks. The five CTTs with the highest L_acc/L_⋆ ratios, namely DG Tau, DL Tau, HN Tau A, DO Tau, and RW Aur A, tend to follow the upper envelope of the Ṁ_acc-M_⋆ relationship for the Lupus population, and have accretion rates comparable to those of the strongest accretors in Lupus.

The NIR veiling increases with wavelength in all the studied CTTs. The analysis of this behavior shows that these CTTs display a significant continuum excess emission in the NIR. In some cases, such excess can be ascribed to the thermal emission from the inner rim of the dusty disk, while in others may be more compatible with emission from a thick gaseous disk inside the dust sublimation radius. The origin of the jets studied in Paper I is compatible with the latter possibility.

Acknowledgements

We very much thank the anonymous referee for her/his comments and suggestions. This work has been supported by PRIN-INAF-MAIN-STREAM 2017 “Protoplanetary disks seen through the eyes of new-generation instruments” and by PRIN-INAF 2019 “Spectroscopically tracing the disk dispersal evolution (STRADE)”. We warmly thank the GAPS team for sharing the solar spectrum acquired with GIARPS, and the TNG personnel for their help during the observations. We dedicate this work to Mr. G.Atusino, whose premature loss deeply saddened all of us; for his immense presence, unconditional support and long-life great friendship. This research made use of the SIMBAD database, operated at the CDS (Strasbourg, France). 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.

Appendix A Measurements of continuum fluxes for the CTTs and templates

In Tables A.1–A.7, we report the observed continuum fluxes for the CTTs and the spectral templates adopted for the determination of the extinction according to the methods described in Fischer et al. (2011). We define the quantity: $${\Gamma }_{\lambda }\equiv 2.5\cdot \log [(1+r_{\lambda })\cdot \frac{F_{\lambda }^{\text{T}}}{F_{\lambda }^{\text{O}}}],$$(A.1)

where r_λ, $F_{\lambda }^{\text{T}}$, and $F_{\lambda }^{\text{O}}$ are the veiling, the observed continuum flux of the spectral template, and the observed continuum flux of the object, as a function of wavelength, respectively. As mentioned in Sect. 3.3, the method is based on the fact that $${\Gamma }_{\lambda }=(A_{\text{v}}^{\text{O}}-A_{\text{v}}^{\text{T}})\cdot \frac{A_{\lambda }}{A_{\text{v}}}-2.5\cdot \log C,$$(A.2)

so that Γ_λ is a linear function of $\frac{A_{\lambda }}{A_{\text{v}}}$, where $A_{\text{v}}^{\text{O}}$ and $A_{\text{v}}^{\text{T}}$ are the visual extinction of the object and template, respectively, and C is a constant. The tables also provide Γ_λ as a function of wavelength and the r_λ values are also included for convenience.

Appendix B Individual fluxes, equivalent widths and L_acc estimates

Tables B.1 to B.5 report the observed fluxes, equivalent widths, for every CTTs in the sample, as well as the corresponding L_acc values derived from the individual accretion diagnostics and using the L_acc–L_line relationships by Alcalá et al. (2017).

References

All Tables

All Figures

Fig. 1 Two portions of the continuum-normalized spectrum of HN Tau A (black solid line) in the HK bands, with the rotationally broadened and veiled spectrum of the stellar template overplotted (red dashed line). Absorption features, mainly of the iron-peak group, are indicated with small vertical lines. We refer to Rayner et al. (2009) for the details on the absorption features. Emission lines such as the H2 line, indicated with the blue dotted lines, are excluded from the analysis. The value of veiling, which minimizes residuals with respectof the rotationally broadened and veiled template, is indicated in the lower left.

In the text

	Fig. 2 $2.5\times \log [(1+r_{\lambda })\cdot F_{\lambda }^T/F_{\lambda }^O]$ vs. Aλ ∕AV for the YSOs inthe sample (blue squares). The best linear fit for each object is shown as a black line. The value of the slope AV is shown in each panel. In the case of CQ Tau the r1565 value is shown as an upper limit.
In the text

	Fig. 3 Example of BT-settl model normalization. The extinction-corrected Asiago spectrum of DL Tau is shown with the red line, while the BT-settl model normalized to DL Tau at λ = 600 nm is plotted with the blue line. We note the increase of excess emission of DL Tau at wavelengths shorter than about 450 nm. Theveiling-corrected BT-settl model is plotted with the black line. See text for details
In the text

	Fig. 4 Examples of GIARPS spectra showing the line profiles of CTTs in our sample.
In the text

	Fig. 5 Examples of photospheric line subtraction for the Hα and Hβ lines in CQ Tau and RY Tau. The black, red, and blue lines represent the observed, photospheric, and photospheric-subtracted spectra, respectively. See text for details.
In the text

Fig. 6 Examples of photospheric line subtraction for the Paβ line for CQ Tau (upper panels) and RY Tau (lower panels). The red lines represent the photospheric templates. The left and right panels show, in black lines, the observed and photospheric-subtracted spectra, respectively. The green vertical lines mark the intervals for the equivalent width measurement. The equivalent widths of the emission line before and after the subtraction are also indicated.

In the text

	Fig. 7 Lacc as a function of the different accretion diagnostic as labeled for the seven targets (blue dots). The vertical error bars consider the error in Lline as well as the errors in the Lacc–Lline relationships. The horizontal red lines in each panel represent the average Lacc.
In the text

Fig. 8 Accretion luminosity as a function of stellar luminosity for the stars in Lupus (black symbols) and the seven Taurus CTTs studied here (blue circles). The latter are labeled with their names. The Lupus transitional disks are shown with crossed squares. The continuous lines represent the three Lacc vs. L⋆ relations as labeled. The long-dashed line represents the shift of HN Tau A on the diagram when correcting its Lacc and L⋆ values for obscuration effects by the disk. The corrected values are shown with the red dot. The leftward and downward red arrows represent the upper limits on L⋆ and Lacc for CQ Tau, respectively. The average errors for the Taurus sample are shown in the upper left. Figure adapted from Alcalá et al. (2017).

In the text

Fig. 9 Mass accretion rate as a function of stellar mass for the seven Taurus CTTs studied (blue circles) compared with stars in Lupus (black symbols). The Lupus transitional disks are shown with crossed squares. Lupus objects classified as weak or negligible accretors are plotted with downward arrows. The long-dashed black line represents the shift of HN Tau A on the diagram when correcting its Ṁacc and M⋆ values for obscuration effects by the disk. The corrected values are shown with the red dot. The leftward and downward red arrows represent the upper limits on M⋆ and Ṁacc for CQ Tau, respectively. The black dashed line shows the double power law theoretically predicted by Vorobyov & Basu (2009), and the continuous magenta lines represent the fits to the data as in Eqs. (4) and (5) of Alcalá et al. (2017). The average errors for the Taurus sample are shown in the upper left. Figure adapted from Alcalá et al. (2019).

In the text

	Fig. 10 NIR veiling as a function of wavelength for the stars in our sample (blue dots). The red dashed lines show the black body fits to the data with $T_{\text{eff}}^{\text{BB}}$ as labeled. See text for details.
In the text