GIARPS High-resolution Observations of T Tauri stars (GHOsT). IV. Accretion properties of the Taurus-Auriga young association

In the framework of the GIARPS@TNG High-resolution Observations of T Tauri stars (GHOsT) project, we study the accretion properties of 37 Classical T Tauri Stars of the Taurus-Auriga star forming region (SFR) with the aim of characterizing their relation with the properties of the central star, of jets and disk winds, and of the global disk structure, in synergy with complementary ALMA millimiter observations. We derive stellar parameters, optical veiling, accretion luminosity ($\rm L_{acc}$) and mass accretion rate ($\rm \dot M_{acc}$) in a homogeneous and self-consistent way using high-resolution spectra acquired at the Telescopio Nazionale Galileo with the HARPS-N and GIANO spectrographs, and flux-calibrated based on contemporaneous low-resolution spectroscopic and photometric ancillary observations. The $\rm L_{acc}$-$\rm L_{\star}$, $\rm \dot{M}_{acc}$-$\rm M_{\star}$ and $\rm \dot{M}_{acc}$-$\rm M_{disk}$ relationships of the Taurus sample are provided and compared with those of the coeval SFRs of Lupus and Chamaeleon I. Our results demonstrate the potential of contemporaneous optical and near-infrared high-resolution spectroscopy to simultaneously provide precise measurements of stellar and accretion/wind properties of young stars.


Introduction
Pre-main-sequence star evolution and planet formation are connected to the interplay of mass accretion onto the star, ejection of outflows, and photoevaporative disk winds (Hartmann et al. 2016;Ercolano & Pascucci 2017). The mass loss from outflows and slow disk winds leads to the extraction of angular momentum, which in turn enables the accretion of matter onto the star and alters the disk density, which affects the formation and migration of protoplanets.
In this context, the mass accretion rateṀ acc (e.g., Hartmann 1998) is a key parameter for the study of young stellar objects (YSOs) during their first several million years of evolution toward the main sequence. Reliable measures ofṀ acc for YSOs located in different star-forming regions (SFR) and at different evolutionary phases are necessary to set important constraints on disk evolutionary models and disk clearing mechanisms (e.g., Manara et al. 2022;Pascucci et al. 2022).
In the current magnetospheric accretion (MA) paradigm (e.g., Hartmann et al. 2016), strong large-scale stellar surface magnetic fields of a few hundred to a few kilogauss truncate the inner disk at a few stellar radii. Driven by the stellar magnetic field lines, gas flows from this truncation radius onto the star and produces localized shocks. The characteristic emission line spectrum of classical T Tauri Stars (CTTs), including the Balmer and Paschen series, is then partly formed in the aforementioned accretion funnel flows (Hartmann et al. 1994). Ac-Article number, page 1 of 30 arXiv:2208.14895v2 [astro-ph.SR] 11 Oct 2022 A&A proofs: manuscript no. 44042corr cretion shocks are also responsible for continuum excess emission mainly from the UV to the optical wavelength range. This continuum emission is overimposed on the photospheric spectrum and thus weakens the stellar photospheric lines. This effect is usually referred to as veiling , and it is particularly remarkable in the UV Balmer continuum region at λ ≤ 400 nm. A proper modeling of the continuum excess emission is necessary to derive the energy that is released in the accretion shocks, the so-called accretion luminosity L acc , which in turn allows measuring the mass accretion rate,Ṁ acc , given the stellar mass and radius ).
On the other hand, L acc and henceṀ acc can also be measured from the well-known correlations between L acc and the luminosity of specific lines in emission (L line ) of H i, He i, and Ca ii, which have been derived for a wide spectral range from the UV to the near-infrared (NIR) (e.g., Calvet et al. 2004;Herczeg & Hillenbrand 2008;Rigliaco et al. 2012;Alcalá et al. 2014Alcalá et al. , 2017. The simultaneous use of many diagnostic lines in a wide spectral region turned out to be of enormous importance in reducing the uncertainties in the L acc andṀ acc measurements at typical values of 0.2-0.3 dex (e.g., Rigliaco et al. 2012;Alcalá et al. 2014).
Accurate determination of L acc andṀ acc is essential to derive the dependence of accretion on the properties of the central star and its disk.Ṁ acc and L acc are usually compared to the stellar and disk properties, and empiricalṀ acc -M ,Ṁ acc -M disk , and L acc -L scaling relations are used as benchmark tests for theoretical predictions of disk evolution (e.g., Manara et al. 2016;Rosotti et al. 2017;Lodato et al. 2017;Mulders et al. 2017;Tabone et al. 2022). Previous surveys have shown that these scaling relations present large scatters (more than 2 dex in logṀ acc or log L acc at a given YSO mass and luminosity) that cannot be explained alone in terms of the high variability characterizing the accretion processes (e.g., Costigan et al. 2012Costigan et al. , 2014Biazzo et al. 2012). One of the current challenges is then to understand the nature of this dispersion. Even though purely viscous evolution disk models predict the above correlations (e.g., Lynden-Bell & Pringle 1974), the observed high dispersion suggests a more complex scenario for disk evolution, in which magnetohydrodynamics (MHD) winds, external photoevaporation, and local dust processing may play a crucial role (e.g., Manara et al. 2022, and references therein). In addition, recent works have pointed out that stellar multiplicity can be another ingredient at the basis of the high accretion observed in some sources (Zagaria et al. 2022).
From an observational point of view, this type of investigation requires simultaneously and in a self-consistent fashion deriving all the relevant parameters of as many as possible statistically complete samples of YSO in SFRs, avoiding methodological systematics and biases due to temporal variability. This was attempted in several SFRs such as Lupus (Alcalá et al. 2014, Chamaeleon I , η−Chamaeleon (Rugel et al. 2018), TW Hydra (Venuti et al. 2019), and Upper Scorpius (Manara et al. 2020) using the capabilities of the X-Shooter instrument, which is the mediumresolution spectrograph (R=10000-20000; Vernet et al. 2011) covering the wide spectral range between 0.3 and 2.5 µm at the ESO Very Large Telescope (VLT). However, a similar survey is lacking for Taurus-Auriga, although it is a paradigmatic SFR.
To cover this gap, we started the GIARPS High-resolution Observations of T Tauri stars (GHOsT) project, a flux-limited survey of T Tauri Stars in the Taurus-Auriga SFR based on data obtained with the GIARPS instrument. This high-resolution spectrograph simultaneously covers the optical and NIR wave-length ranges and is located at the Italian Telescopio Nazionale Galileo (TNG).
The GHOsT project aims at providing reliable measurements of stellar and accretion parameters of Taurus-Auriga members and to place them in relation with the disk properties, in synergy with complementary ALMA observations that are available for the majority of the targets (Andrews et al. 2018;Long et al. 2019). In addition, the high sensitivity and spectral resolution of GIARPS also allows us to investigate the sub-au disk and jets environments by studying specific diagnostic lines. A first study of the GHOsT program was focused on the jet line emission for a subsample of bright sources (Giannini et al. 2019, henceforth Paper I), and then was followed by an investigation of the link between atomic and molecular disk winds (Gangi et al. 2020, henceforth Paper II). The definition and assessment of the methods for determining stellar and accretion properties was then presented in Alcalá et al. (2021, henceforth Paper III).
In this fourth work of the GHOsT series, we present the results of the accretion measurements for the complete sample. The paper is organized as follows. In Sect. 2 we present the source selection and the sample properties, while observations and data processing are reported in Sects. 3 and 4, respectively. In Sect. 5 we describe the method we used to derive the stellar and accretion properties, while results are shown in Sect. 6. Discussion and conclusions are finally presented in Sects. 7 and 8, respectively.

Sample
The GHOsT original sample was selected based on the most recent Taurus-Auriga young population census at the time of the beginning of the GHOsT project, that is, the one presented by Esplin et al. (2014). The initial selection was driven by the GI-ARPS sensitivity. We considered sources with J < 11 mag and R < 13 mag, which reduced the total sample to about 70 objects with spectral type (SpT) earlier than M5. We then started the observational campaign giving priority to the sources with known jets and with outer disks that were already characterized by ALMA observations. It was possible to observe 46 of the 70 objects with the telescope time available to GHOsT. Their masses are in the range between ∼ 0.2 and ∼ 2.2 M , the spectral types are between G1 and M3, and the luminosities are between ∼ 0.05 and ∼ 8.9 L . Although the observed sample is not complete, it is representative of the distribution in the complete sample with SpT earlier than M5. We observed 9 class III sources in the sample as classified in the new census by Esplin & Luhman (2019), with the purpose of refining the empirical relations between the equivalent widths (EWs) and fluxes of the Hα and Paβ lines and the effective temperature (T eff ) for a preliminary selection between class II and class III sources to be observed in the NIR.
Information such as disk inclination, dust mass, effective disk radius, and disk structure is available for the majority of the observed sources. In Table 1 we report the respective data collected from the most recent ALMA literature. Disk classification was made preferentially on the basis of ALMA images. In addition to the full (i.e., no substructures in the dust density profile detected at ALMA resolution) and transitional (TD; i.e., with an inner cavity around the star) disk categories, we indicate as substructured the class of disks that are characterized by the presence of gaps and rings. Finally, ∼ 35% of the sample consists of multiple systems. The type of multiplicity, separation, and the relevant references are reported in columns 4 and 5 of Table 2. Article number, page 2 of 30 M. Gangi et al.: GIARPS High-resolution Observations of T Tauri stars (GHOsT) Notes.
(a) Disk inclinations are measured from ALMA observations and refer to the outer disk. An exception is made for GG Tau, where the inclination of the resolved inner disk is reported.

Observations
Observations were carried out from October 2017 to January 2020. They consist of high-resolution optical and IR spectra obtained with the GIARPS instrument that were flux-calibrated through ancillary low-resolution spectroscopy and photometric data. The logbook of observations is given in Table 2.

GIARPS observations
The sample was observed in four distinct runs, hereafter run I in October -November 2017, run II in December 2018, run III in November 2019 -January 2020, and run IV in October -December 2020. The GIARPS observing mode consists in the simultaneous use of the HARPS-N (Pepe et al. 2002;Cosentino et al. 2012, resolving power R = 115000) and GIANO-B (Oliva et al. 2012;Origlia et al. 2014, R = 50000) spectrographs. HARPS-N is a fiber-fed echelle with a FoV = 1 and covers the spectral range between 390 and 690 nm while GIANO-B is a nearinfrared cross-dispersed echelle with a slit on-sky with dimensions of 6 × 0.5 and a spectral range between 940 and 2400 nm.

Ancillary observations
To accurately flux-calibrate the GIARPS spectra we carried out contemporaneous ancillary observations consisting of absolute flux-calibrated low-resolution spectroscopy and photometry. In particular we acquired optical low-resolution spectra (R=2400, 330-790 nm) with the 1.22 m telescope of the University of Padova at the Asiago observatory (Italy). Spectra were reduced and flux-calibrated against spectrophotometric standards observed on the same night. We also checked and refined the flux zeropoint with BVR c I c photometry acquired with the ANS Collaboration telescopes (Munari et al. 2012) in runs I-III, while in run IV, griz photometry was taken with the ROS2 instrument of the REM telescope (Molinari et al. 2014). In the NIR we obtained JHK photometry that was acquired in run I, II, III-2019 with the REMIR instrument at the REM telescope (Vitali et al. 2003) and in run III-2020 and IV with the NICS camera (Baffa et al. 2001) at the TNG telescope. In addition, we also acquired low-resolution (R∼ 50) NIR spectra in run IV using the NICS Amici prism.
Details of the observation date and instruments we used are listed in Table B.1. The photometry of runs I-III is reported in Gangi et al. (2020), while that of run IV is reported in Table B.2 for completeness.

Data reduction
The data reduction processes are described in detail in Papers I-III, but are briefly summarized here for completeness. The HARPS-N spectra were reduced according to the standard procedures with the HARPS-N data reduction software pipeline (Pepe et al. 2002). Spectra were corrected for heliocentric and radial velocity (RV) using the Li i photospheric profile, assuming the weighted λ air = 670.7876 nm. We estimate a final accuracy on the wavelength calibration of about 2 km s −1 . We then extracted and normalized to the continuum the 16 spectral portions, 50-100 Å wide, that contained the diagnostics that we used for the accretion measurement, namely the Balmer recombination lines H3, H4, H5, H6, H7, H8, the helium lines He ii469, He i403, He i447, He i471, He i492, He i502, He i588, He i668, and the Ca ii H and K doublet.
Article number, page 3 of 30 A&A proofs: manuscript no. 44042corr  The GIANO-B spectra were reduced following the prescriptions given in Carleo et al. (2018). We then considered the five spectral segments containing the accretion diagnostics, namely the hydrogen recombination lines of the Paschen series (i.e. Pa5, Pa6, Pa7 and Pa8) and the He i1083 nm line. They were continuum-normalized and corrected for heliocentric and RV. For this latter task, we used the three NIR Al i lines at λ vac = 2109.884, 2116.958, 2121.396 nm as reference after verifying that they agreed well with the velocity measured from the optical Li i line. We obtained a typical final accuracy on the wavelength calibration of about 1-2 km s −1 . Finally, the five spectral segments were corrected for telluric features. For this purpose, we used the molecfit tool (Smette et al. 2015) to compute the synthetic telluric spectrum and the IRAF task telluric to properly correct the observed spectral segments for telluric contribution.
We flux-calibrated the GIARPS spectra on the basis of the ancillary data. They were acquired close in time with the GI-ARPS observations, with a maximum temporal distance of ten days. On these timescales, the typical flux variability of CTTs can reach values of 10% (e.g., Venuti et al. 2021). We then conservatively considered this latter as the typical uncertainty to the flux calibration due to the source variability.
For each source, we computed a polynomial curve representative of the continuum flux within the GIARPS spectral range on the basis of the available ancillary data. In particular, for runs I-III, we built this curve through a polynomial fit of the continuum of the Asiago spectrum and an interpolation of the I JHK photometric points. In run IV, we fit the total spectrum resulting from the sum of the Asiago and Amici spectra after refining the continuum level with the photometric points.
We then multiplied the 21 continuum-normalized spectral segments for the computed curve. Considering both the source variability and the errors on the photometric points, we estimate an accuracy of 20% on the final flux calibration.

Methods for the derivation of stellar and accretion properties
In this section we present the methods we adopted for measuring the stellar and accretion parameters. The general method has been extensively discussed in Paper III, to which we refer. Here, we briefly present it and complement it with a detailed discussion of the case of heavily spotted stars.

General method
Photospheric and kinematical parameters (i.e. effective temperature T eff , surface gravity logg, RV, and projected rotational velocity v sin i) and veiling as a function of wavelength r λ were determined with the ROTFIT code (Frasca et al. 2003). In short, the code performs a χ 2 minimization between observed and template spectra in specific optical spectral segments and has been successfully applied both on medium-and high-resolution data (e.g., Frasca et al. 2015Frasca et al. , 2017. We used the HARPS-N spectra and a grid of templates retrieved from the ELODIE archive (Moultaka et al. 2004) of real spectra of slowly rotating and lowactivity stars with well-known atmospheric parameters. T eff and veiling at 600 nm, r 600 , which are the only two parameters used for our analysis, are reported in Table 2, while the complete set of photospheric parameters derived with ROT-FIT will be presented and discussed in a future dedicated work. Spectral types were determined using the SpT-T eff conversion of Herczeg & Hillenbrand (2014), and they are also reported in Table 2.
As discussed in Paper III, the visual extinction A v was estimated by computing the ratios of the flux-calibrated lowresolution Asiago spectra and artificially reddened templates of the same spectral type. We adopted the grid of nonaccreting YSO templates with negligible extinction of Manara et al. (2013Manara et al. ( , 2017 reddened by A v in the range 0.0-5.0 mag in steps of 0.10 mag using the extinction law by Weingartner & Draine (2001) with R v = 5.5. The A v was then estimated as that of the reddened template that minimized the slope of the corresponding ratio. This analysis was limited to the portion of the spectra between 550 and 800 nm, which is least affected by both the optical and infrared veiling (e.g., Fischer et al. 2011). Values derived in this way are reported in column 7 of Table 2.
To compute the stellar luminosity L , we first considered the BTSettl model (Allard et al. 2012) for each source that best matched the photospheric parameters derived with ROTFIT and normalized it to the extinction-corrected Asiago spectrum at λ = 600 nm. We then integrated the BTSettl spectrum over all wavelengths and multiplied it by the factor 1/(1+r 600 ) to take the veiling contribution into account, thus obtaining the bolometric flux of the object photosphere. The stellar luminosity was then computed as L = 4πd 2 F, where F is the bolometric flux and d is the Gaia EDR3 distance reported in Table 2. As reported in Paper III, we estimated an average uncertainty of 0.2 in log L on the basis of the typical errors in flux calibration of the Asiago spectra and in veiling correction. The stellar radius, R , was computed from L and T eff , and an average uncertainty of 0.3 in log R is estimated.
Finally, the mass, M , was computed from the pre-mainsequence evolutionary tracks by Siess et al. (2000). From the typical uncertainties of L and T eff , we estimated an average uncertainty of about 0.15 in log M . Values of L , M , and R are also listed in Table 2.

Stellar parameters in spotted sources
Magnetic activity in young stars can generate starspots over a large area of the stellar surface that can significantly contribute to the total stellar flux, with a spectrum that is distinct from the ambient photosphere. Observationally, this can lead to a strong spread in the effective temperature and luminosity derived by different methods and wavelength ranges. It has been often found that effective temperatures measured in young stars with NIR spectra are systematically lower than those derived at optical wavelengths (e.g., Cottaar et al. 2014;Flores et al. 2021). This offset is particularly remarkable in the temperature range between ∼ 3800 and ∼ 4500 K where it can reach average values of 500 K. This discrepancy can also introduce a significant uncertainty in the mass and age values estimated from the position of the star in the HR diagram (e.g., Gully-Santiago et al. 2017). Flores et al. (2021) compared the masses estimated from optical and NIR temperatures with the dynamical masses measured from ALMA for a sample of YSOs, finding that neither the optical nor the NIR temperatures reproduce the stellar dynamical masses.
A detailed analysis of the influence of starspots on the measurement of photospheric properties and a discussion of whether the NIR-optical T eff offset can be fully explained in terms of starspots are beyond the scope of this paper. Here we limit our analysis to heavily spotted stars in our sample and evaluate the best strategy to properly estimate their luminosity and mass be-  cause these latter are the parameters that can directly affect the L acc -L andṀ acc -M distributions. Fig. 1 shows the offset between the T eff we measure in the HARPS-N range with ROTFIT and the T eff values reported in the literature measured from an analysis mostly based on TiO bands (Herczeg & Hillenbrand 2014) and from NIR spectral features (Nofi et al. 2021;López-Valdivia et al. 2021). Similarly to what was found in other YSOs samples, we confirm that the main T eff discrepancies lie in the temperature range between ∼ 4000 and ∼ 4500 K. In particular, the detection of optical TiO bands (which are typically sensitive to T eff ≤ 3800 K) in this temperature range strongly suggests the presence of regions colder than 4000 K.

Class III sources
To evaluate the impact of multiple temperature components on the broadband emission of the source, we first focused on the subsample of class III sources, whose analysis is greatly simplified because they are not affected by significant veiling. Fig. 2 (left panel) shows the extinction-corrected low-resolution Asiago spectra and the contemporaneous broadband photometry compared with the BTSettl model corresponding to the T eff measured from the analysis of many optical features and scaled to the Asiago flux, as explained in the previous section. For the six sources lying in the temperature range between ∼ 4000 and ∼ 4500 K (namely DI Tau, LkCa 4, IW Tau, V1070 Tau, V827 Tau, and V819 Tau), the synthetic models significantly underestimate the NIR fluxes.
This discrepancy is not attributable to the variability because class III sources are known to be quite stable with small changes of amplitude of about a few tenths of magnitude (e.g., Grankin et al. 2008;Lanza et al. 2016) while the mismatch we found reaches values of more than 1 mag in the NIR bands. An additional source of bias could be the assumption of an incorrect extinction value. In principle, we can still reproduce the broad-band spectrum with a single model corresponding to the optical T eff by significantly increasing the A v value with respect to that determined with the method described in the previous subsection. However, as reported in Sec. 5.2, for class II sources, we can adopt a independent way, free from the choice of the SpT, to compute the extinction based on the minimization of the spread around the L acc -L line relations. The results of this latter approach are always consistent within the errors with that obtained with the former method, confirming the robustness of our measures. On the other hand, the synthetic models corresponding to the T eff measured from the NIR are not able to reproduce the optical region of the spectra, showing deeper TiO molecular bands and an incorrect extinction (Fig. 2, right panel). In conclusion, our evidence indicates that a single T eff model is not sufficient to reproduce the total flux of a subsample of class III sources and suggests the need of multiple-temperature modeling.
In these sources, large spotted regions have been revealed on V819 Tau, V827 Tau, and V1070 Tau from photometric longterm time-series analysis (Grankin et al. 2008). The models show that the observed light-curve amplitudes are compatible with spotted regions covering from 17% to 73% of the visible stellar surface with mean temperatures 500-1400 K lower than the ambient photosphere. In addition, spectral features of LkCa4 have been interpreted in terms of 2 T eff components by Gully-Santiago et al. (2017), who found a hot photosphere of ∼ 4100 K and a cool contribution of ∼ 2700 − 3000 K due to starspots covering a surface area of ∼ 80%.
For the six class III sources with evidence of starspots we have followed the method described in Gully-Santiago et al. (2017) which is similar to that proposed by Frasca et al. (2005). In short, we assumed that the stellar photosphere can be represented by two components, one with a hot temperature T hot and one with a cool temperature T cool . Each of these components contributes to the total flux with appropriate weighting factors, w hot and w cool . The effective temperature best representing the total stellar flux was computed as while the luminosity was computed from the integral of the total flux in the whole spectral range. We used BTSettl models by fixing T hot to the closest value found by ROTFIT, and we varied the T cool models from 2500 to 4000 K in steps of 100 K. The weighting factors were left as free parameters of the fit. The adopted method (hereafter broadband fit) allowed us to reproduce the broadband spectrum of the six class III sources in the T eff range between 4000 and 4500 K, and hence to properly estimate the corresponding luminosities. Fig. 3 (top panel) shows an example of the fitting procedure for IW Tau, and the whole subsample of class III sources and the fitting results (i.e., temperatures and weighting factors) are reported in Fig. B.1 and Table  B.3, respectively.
We found cool components with temperatures between ∼ 3100 and ∼ 3500 K and weighting factors between ∼ 27% and ∼ 80%, leading to a decrease in the effective temperature between ∼ 5% and ∼ 18%. LkCa4 has the largest contribution (∼ 80%) from a cool component (T cool ∼ 3100 K), which perfectly agrees with the results of Gully-Santiago et al. (2017).
However, we stress that the physical meaning of the T cool and T hot components and their relative weighting factors cannot be directly linked to starspots, ambient photosphere, and filling factor. Localized hotspots such as plages can additionally contribute to the total spectrum, which is indistinguishable from Left: Extinction-corrected Asiago spectra (black) with contemporaneous photometry (orange points). BTSettl models corresponding to the optical T eff are shown in blue. The filled light blue regions highlight the NIR flux that is underestimated by the models. Right: As in the left panel but assuming Av=0 and BTSettl models corresponding to the T eff measured from NIR high-resolution spectra . Source names and T eff are labeled for each spectrum.
the ambient photosphere with the current technique. More advanced methods such as the Zeeman Doppler Imaging (Semel 1989;Brown et al. 1991) or line-bisector analysis (e.g., Prato et al. 2008) are necessary for this purpose. Fig. 4 (left panel) shows the position of these sources in the HR diagram with an effective temperature and luminosity based on the 2 T eff broadband fit compared with that obtained assuming a single temperature. The 2 T eff analysis leads to a systematic shift toward younger ages and lower masses.

Class II sources
The broadband fitting procedure cannot be directly applied to class II sources because of the accretion excess, which can contribute significantly to the observed flux. Our strategy for these sources was to perform a 2 T eff spectral fitting of the lowresolution Asiago spectra in the spectral range between 550 and 800 nm. This latter is in fact less affected by the veiling and contains the optical TiO and VO molecular bands which are particularly sensitive to the presence of a cool T eff component. The adopted procedure is the same as for the broadband fit with the difference that we now consider specific spectral features. For this reason the spectral resolution of the models was degraded to that of the Asiago spectra and the portion of the spectrum containing the Hα line was excluded with a specific mask.
We first applied the spectral fitting procedure to the subsample of class III sources (an example is shown in the bottom panel of Fig. 3) and compared the results with those obtained with the broadband fit. The results shown in Figure B.3 show that differences between the two methods are less than ∼ 10% for both luminosities and temperatures, thus indicating that the spectral fitting procedure unambiguously allows a proper estimation of the stellar luminosity and can also be applied in the case of accreting sources.
We then applied this method to the whole class II sample by again fixing T hot to the closest value found by ROTFIT. The second T cool component was considered non-negligible only when it improved the χ 2 of the spectral fit by more than 20% of its minimum value. Following this criterion, we found that a second T cool component was necessary for ten class II sources with optical T eff ≥ 4000 K, which we therefore identified as heavily spotted sources. Independent evidence of large spotted regions in these sources has already been found for DQ Tau (Kóspál et al. 2018) and GI Tau (Guo et al. 2018) from the analysis of multiband light curves, for DN Tau from RV variations (Prato et al. 2008), and for V836 Tau from both multiband light curves (Grankin et al. 2008) and RV variations (Prato et al. 2008). The results of our modeling are listed in Table B.3 and shown in Fig.  B.2. They show that cool regions with mean temperatures 900-1300 K lower than the optical T eff contribute as much as 24% to 90% to the total flux.  BTSettl -3200 K -4100 K Fig. 3. Example of the 2 T eff fitting procedures for the class III source IW Tau. The extinction-corrected Asiago spectrum is shown in black, and contemporaneous photometry is shown as filled orange squares. The best-fit spectrum is shown in purple, and the corresponding T hot =4100 K and the T cool =3200 K BTSettl models are reported in blue and red, respectively. Top: Broadband fitting. Bottom: Spectral fitting. Flux units are 10 −13 ergs −1 cm −2 nm −1 . Fits to the complete subsample of heavily spotted sources are reported in Appendix B.2 and B.1.

Accretion properties
The accretion luminosity of the 37 class II sources was estimated as described in Paper III, that is, by using the empirical relations between L acc and the luminosity of emission lines (L line ) given in Alcalá et al. (2017) and determined from X-shooter observations of a sample of class II sources in Lupus. The L line was computed as L line = 4πd 2 F line , where d is the Gaia EDR3 distance reported in Table 2 and F line is the extinction-corrected line flux. This latter was in turn determined through the integration of the line profile with the IRAF task splot. For each line we performed three independent measurements at the lowest, highest and middle position of the local continuum, to take the uncertainties introduced by the local noise into account. F line and its error were then computed as the average and standard deviation of these three measurements, respectively. The correction for the extinction was performed using the A v value, which was determined as described in the previous section.
For sources with spectral types earlier than K0 the F line measurement is strongly influenced by photospheric absorptions. For these sources, we therefore performed the photospheric subtraction using appropriate spectral templates as described in Paper III. In addition to CQ Tau and RY Tau, which were analyzed in Paper III, we successfully performed this subtraction in another four sources, namely HQ Tau, MWC480, UX Tau and SU Aur. Details are reported in Appendix A. We also note that regardless of the level of accretion, the strong photospheric continuum expected for these sources can limit the number of detected accretion diagnostics.  In all sources the Ca ii H and H7 lines were found to be blended. We attempted to separate the two contributions following the Gaussian decomposition procedure of Gangi et al. (2020). We obtained reliable results for nine sources (i.e., CI Tau, CY Tau, DE Tau, DN Tau, DQ Tau, GH Tau, IP Tau, UY Aur, and V836 Tau) for which the limited broadening of the profiles allowed us to unambiguously disentangle the two lines at the HARPS-N resolution.
When no diagnostic line was detected, we estimated a 3σ upper limit as 3× rms ×∆λ, with rms the local flux noise and ∆λ the expected line width. The latter was estimated from the other detected lines and assumes typical values between 0.1 and 0.2 nm. Fig. 5 shows an example of the derived L acc values plotted as a function of the line diagnostics, while the plots for the complete sample are reported in Fig. B.4. We note that the A v calculated in the previous section is always consistent with that obtained by minimizing the spread of L acc , confirming the consistency of the adopted methods. For each source, we assumed the median value L acc and its standard deviation as the accretion luminosity and error, respectively.
The mass accretion rateṀ acc was then computed aṡ with R star and R in the stellar and the inner disk radius, respectively. We assumed R star /R in ≈ 1/5 as in Gullbring et al. (1998) and Hartmann (1998), and used the stellar mass M and radius R derived in the previous section and reported in Table  2. As discussed in Paper III, we estimate an average uncertainty of about 0.4 in logṀ acc , which derive from the uncertainties on L acc , M , R , and Gaia EDR3 distances and from the differences in the adopted evolutionary tracks. The obtained values of L acc andṀ acc are reported in Table 3, together with the number of lines used for the determination of L acc .
Finally, we analyzed the EWs and fluxes of Hα and Paβ lines as a function of effective temperature for the complete sample and provided a new empirical criterion to distinguish class III from class II sources based on the EW Paβ -T eff and F Paβ -T eff planes. Details of the adopted method and results are reported in Appendix A.  Inspecting the disk classification, we note that sources with full and substructured disks show the largest spread both in L acc and L . In contrast, sources with transitional disks fill the area between the L acc = 0.1L and L acc = 0.01L boundaries. No In the top panels the symbols are distinguished according to the disk structure and in the lower panels the symbols are allotted according to the multiplicity. The bottom panels provide a comparison with Lupus (blue) and Chamaeleon I (magenta). Solid lines are the best fit to data of the same color.
Article number, page 11 of 30 A&A proofs: manuscript no. 44042corr Table 3. Accretion properties of the selected accreting CTTs derived in this work and number of diagnostics used to compute the L acc . clear trend between different source multiplicities is found with the current statistics.
Although the completeness level of our sample is lower than that of Chamaeleon I and Lupus, there is an overall similarity in the distribution of the three SFRs in the slope and in the spread (Fig. 10). We performed a linear fit of the data with the idl linmix_err tool, which performs a linear regression using the Bayesian method of Kelly (2007) and allows us to take the uncertainties on both axes into account. We adopted the median values and standard deviations of the chains as best-fit coefficients and errors, respectively. The fitting procedure, applied to the same L range as was probed with our data, confirms this agreement: Taurus : log L acc = 1.69(±0.93) log L − 0.74(±0.29), (3) Lupus : log L acc = 1.73(±0.27) log L − 1.33(±0.16), (4) Cham : log L acc = 2.11(±0.23) log L − 1.10(±0.22).
The fitting coefficients found here for the Lupus and Chamaeleon I distributions are compatible within the errors with those obtained using Gaia DR2 distances by Alcalá et al. (2017Alcalá et al. ( , 2019 and Manara et al. (2017).

Mass accretion rate versus stellar mass
The dependence ofṀ acc on M in the Taurus sample is shown in the right panels of Fig. 6. Again, there is a general similarity between the Taurus, Lupus, and Chamaeleon I distributions. For the latter two, M , and thusṀ acc , was determined basing on the nonmagnetic models of Baraffe et al. (2015) for targets with T eff ≤ 3900 K and of Feiden (2016) for hotter stars (see Manara et al. 2022, for details).
The GHOsT points do not fill the region at masses lower than ∼ 0.2 M , and hence it is not possible to investigate whether a double power trend, as theoretically predicted by Vorobyov & Basu (2009) We therefore find a lower slope than was proposed in the past (∼ 1.6-2.2; e.g., Antoniucci et al. 2014;Alcalá et al. 2014;Manara et al. 2016;Venuti et al. 2019), which is probably due to the small M range we probed. Our result is comparable within the errors with the values obtained for the distributions of Lupus and Chamaeleon I when the analysis is restricted to the highmass (i.e., M ≥ 0. Again, the fitting results for the Lupus and Chamaeleon I distributions agree well with those reported in Alcalá et al. (2017Alcalá et al. ( , 2019 and Manara et al. (2017).
We also note that except for a few cases, the Taurus distribution falls below the theoretical prediction of Vorobyov & Basu (2009). This is in line with the Lupus distribution, while the Chamaeleon I dataset presents a substantial fraction of points placed above this locus. Finally, transitional disks tend to be placed on the lower side of the distribution, while no particular trend is seen between the full and substructured disks, nor for the source multiplicity.

Mass accretion rate versus disk mass
The dependence ofṀ acc on the disk dust mass (M disk,dust ) is shown in the left column of Fig. 7, where again the disk structure and the source multiplicity is highlighted and the Taurus distribution is compared with the distributions of the Lupus and Chamaeleon I populations. M disk,dust were retrieved from the most recent literature based on ALMA 3.1 mm flux measurements (see Table 1). When only the flux measurement has been reported, we consistently computed the M disk,dust as described by Hildebrand (1983), with F ν the integrated 3.1 millimeter flux, d the Gaia EDR3 distance, k ν the dust opacity, and B ν (T dust ) the blackbody emission at temperature T dust = 20 K. We adopted a power-law opacity k ν = 2.3(ν/230 GHz) 0.4 cm 2 g −1 as in Long et al. (2019).
We find a moderate correlation between M disk,dust anḋ M acc , which is similar to what is observed in the Lupus and Chamaeleon I samples. As pointed out by Manara et al. (2016), for example, a large range of M disk,dust needs to be probed to find a stronger correlation. Considering the disk classification, we note that transitional and substructured disks tend to have higher values of M disk,dust , while sources with full disks homogeneously fill theṀ acc -M disk,dust plane.
In the right column of Fig. 7, we consider the surface mass density defined as Σ disk,dust = M disk,dust /4πR 2 eff instead of the disk mass, with R eff the disk radius enclosing 68% of the 3.1 mm flux. In this case, the distribution of transitional and substructured disks becomes similar to that of full disks. We discuss this different behavior in Section 7.
The distribution of sources with multiplicity is well distinct from that of single sources in theṀ acc -M disk,dust plane. This trend can be interpreted both in terms of lower M disk,dust at a givenṀ acc or higherṀ acc at a given M disk,dust in multiple stars compared with single stars. In any case, stellar multiplicity provides a substantial contribution to the observed spread in theṀ acc -M disk,dust relation.

Accretion luminosity versus wind and jet properties
The physical properties of protoplanetary-disk (PPD) winds and jets can be investigated through the emission of atomic or weakly ionized forbidden lines in the optical and infrared spectral ranges (e.g., Giannini et al. 2019). These lines usually display a composite profile, with a high-velocity component (HVC, |v| > 40 km s −1 ), associated with extended collimated jets, and a lowvelocity component (LVC, |v| < 40 km s −1 ), attributed to 0.5-10 AU PPD winds (Hartigan et al. 1995). In some cases, the LVC shows a composite profile where a broad low-velocity component (BLVC; FWHM 40 km s −1 ) and a narrow low-velocity component (NLVC; FWHM 40 km s −1 ) can be distinguished (Rigliaco et al. 2013;Simon et al. 2016;Banzatti et al. 2019;Gangi et al. 2020  630 nm line is the brightest and is most frequently studied to infer the properties of PPD winds. Molecular winds have also been observed, and from the ground, they are usually traced through high-resolution spectroscopy from the UV to the IR by studying specific bands such as those of H 2 , CO, H 2 O, and OH (e.g., Najita et al. 2007;France et al. 2012;Banzatti et al. 2022). The band that falls within the GIARPS spectral range is the ro-vibrational 1-0 S(1) transition of molecular hydrogen at 2.12 µm.
The properties of the PPD atomic and molecular winds and jets have previously been analyzed for some subsamples of the GHOsT program in Paper I and II, while those of the complete sample will be presented in a forthcoming work (Nisini et al., in preparation). Here we investigate the connection with the accretion properties. The [O i] line luminosity (L line ) was found to correlate with L acc and was used to measure the mass accretion rate (e.g., Herczeg & Hillenbrand 2008). This correlation is still valid when the LVC and HVC luminosities are considered separately, as found by Rigliaco et al. (2013) and Simon et al. (2016) in a sample of Taurus sources with an inhomogeneous set of data, and by Nisini et al. (2018) for a large sample of CTTs from the Lupus, Chamaeleon, and σ-Orionis SFRs. This supports the scenario in which the physical origin of the two line components is related to the accretion mechanism. In Fig. 8 we report the line luminosities of the [O i] NLVC (which is the most frequently detected LVC, see Paper II) and HVC and H 2 as a function of the accretion luminosity. Line luminosities were computed from the deconvolved fluxes of Paper II, corrected for the values of A v reported in Table 2 and rescaled to the new Gaia EDR3 distances. In addition, as suggested by Mendigutía et al. (2015), both quantities were normalized at the stellar luminosity to take off the well-known L acc − L and L line − L correlations. Despite the low statistics and the small L acc range that was probed, we can confirm the correlation between the luminosity of the For both components we then find a slope that agrees with the results of Nisini et al. (2018), namely 0.4 for the LVC and 0.7 for the HVC, respectively, confirming that the HVC line luminosity decreases more steeply for low accretion luminosities. This behavior was also pointed out in Rigliaco et al. (2013). Regarding the hydrogen molecule component, we find no correlation between the H 2 line luminosity and L acc . As we pointed out in paper II, the detection of H 2 2.12 µm depends on the degree of H 2 dissociation in the wind, a process that might depend on several causes (i.e., stellar luminosity or radial wind extension) that are not necessarily connected with accretion.
In addition to the L line −L acc correlation, Banzatti et al. (2019, henceforth B19) found a correlation between the deprojected HVC peak velocity (HVC vel ) and L acc . In Fig. 9 we plot HVC vel as a function of L acc for our sample and for the B19 sample. The latter is based on the Taurus, Lupus, Ophiucus, Corona Australis, and TW Hya star-forming regions, and their deprojected HVC vel values were updated here considering the latest measurements of disk inclinations available in the literature. Our data confirm the trend found in B19. The best-fit linear regression gives Moreover, we find a smaller scatter (∼ 40 km s −1 ) than B19, which could be explained by the homogeneity and simultaneous measurements of L acc and dHVC vel in our sample.
Finally, in Fig. 9 we also indicate with a cross the L acc for the sources for which the [O i] HVC was not detected. The distribution of these points within the wide range of L acc supports the idea, already suggested in Nisini et al. (2018), that the presence of a jet only marginally depends on the accretion level of the source.

Discussion
Our survey of the accretion properties of YSOs in the Taurus-Auriga population confirms the results found in other starforming regions with similar age. Although the incompleteness of our sample prevents us from deriving global results for the whole Taurus population, we can analyze some general properties of the accretion process at this stage of stellar evolution, taking advantage of (i) the self-consistency of the method we adopted to derive both stellar and accretion properties, (ii) the detailed knowledge of the disk structures from the complementary ALMA observations that are available for the majority of our sample, and (iii) the knowledge of the jet and disk wind properties derived in Paper I and II using the same set of data.
7.1. Dispersion on the L acc − L andṀ acc − M planes We find a dispersion of ∼ 0.6 dex for the L acc −L andṀ acc −M relations. The dispersion is computed as the 1σ of the Gaussian fit to the distribution of the residuals, which are defined as the distances of the individual L acc andṀ acc values from the respective best-fit relations. Fig. 10 shows the normalized distributions of the residuals obtained for the Taurus sample and compared with those of the Lupus and Chamaeleon I. For these latter we limited the computation of the residuals to the same L and M ranges as were probed with the GHOsT data to avoid possible biases due to the adoption of different abscissa ranges. While the dispersion on the L acc − L relation is similar for the three samples, the dispersion around theṀ acc − M relation in our sample appears to be lower than those measured in the Lupus and Chamaeleon I datasets. We emphasize that the stellar and accretion parameters were measured in a self-consistent way in all three samples, so that this discrepancy cannot be an effect of the inhomogeneity of the data. On the other hand, for Taurus, we adopted a proper modeling to derive the stellar properties in case of heavily spotted sources by performing a 2 T eff component fit for the photospheric emission. In Fig. 11 and 12 we compare the L acc − L andṀ acc − M distributions obtained assuming both a single and a 2 T eff modeling for the same subsample of sources with a strong evidence of spots. The L acc − L relation is not significantly influenced by the choice of the method because this latter only affects the L measurements. As already noted in Sect. 5.1.2 the stellar luminosities derived with a 2 T eff modeling are systematically higher than those measured from a single T eff , but the differences are always within the typical errors. In contrast, the dispersion in theṀ acc − M plane seems to be affected in a significant way: the adoption of a single T eff modeling for the determination of stellar parameters in spotted sources tends to increase the dispersion in theṀ acc − M plane by ∼ 0.25 dex (Fig. 12) in such a way that we obtain a dispersion value that is more in line with those measured in the Lupus and Chamaeleon I samples. This suggests that the lower dispersion found in thė M acc − M relation for the Taurus sample might be the result of our accurate determination of the stellar parameters in sources with a strong presence of spots.
Even considering possible biases in theṀ acc measurements, it is well established that large dispersions are a common characteristics of different SFR populations, with typical values of ∼ 0.75 dex and a total range in logṀ acc and log L acc of about two orders of magnitude at a given stellar mass and luminosity (e.g., Hartmann et al. 2016;Testi et al. 2022). All of the observational efforts currently conclude that much of this dispersion is real, that is, not due to observational and methodological biases, and it cannot be entirely attributed to the scatter due to the vari-ability of the accretion processes (e.g., Biazzo et al. 2012;Venuti et al. 2014).
Several processes are usually invoked to explain the observed spreads. In the framework of a pure viscous evolution of the disk, it has been suggested that the mass accretion rate should decline with age and would be prevented by the stellar magnetosphere at ages 10 Myr . A tentative decrease in mass accretion rates with age has been found in samples from SFR with different ages (e.g., Sicilia-Aguilar et al. 2010;Antoniucci et al. 2014;Testi et al. 2022), but the uncertainties in measuring reliable individual ages make it very difficult to test any dependence of the mass accretion rates on stellar age in the individual SFR.
A different way to indirectly investigate the effect of disk evolution on the accretion process is to search for any trend with the disk dust structure. In the current idealized picture of the protoplanetary disk evolution, the disk dispersal is predominantly guided by the interplay between accretion and mass loss through disk winds. While the disk evolution is initially driven by accretion, mass-loss through low-velocity winds starts to significantly alter the disk density and dominates the accretion on a timescale of a few million years. As a consequence the disk gas is rapidly cleared inside-out leading to the formation of transitional and substructured disks (e.g., Ercolano & Pascucci 2017). In this framework, at a given stellar mass, it is expected that transitional disks and/or disks with substructures have relatively lower accretion rates than full disks. In our sample, with the limited statistics we have on hand, we note that TD sources tend to have lower values of accretion of L acc (they fill the region of the plane where L acc = 0.1L and L acc = 0.01L ) andṀ acc on average, in line with this scenario and with other previous observations (e.g., Espaillat et al. 2014, and references therein). However, this does not effectively contribute to an increase in spread in the relations because sources with full disks present the largest spread and homogeneously fill the L acc − L andṀ acc − M planes. This indicates that the complexity in disk structure is not the main tracer of the accretion level.

Accretion and disk mass
Viscous evolution models predict a correlation between disk mass and mass accretion rate (e.g., Lynden-Bell & Pringle 1974;Hartmann et al. 1998;Rosotti et al. 2017) that can be expressed as with t disk the disk lifetime and t v the viscous time. The M disk /Ṁ acc ratio is then driven by the disk age if t t v or by the viscous time if t t v (e.g., Lodato et al. 2017).
From an observational point of view, theṀ acc − M disk correlation was confirmed for the 1-3 Myr old SFRs of Lupus  and Chamaeleon I (Mulders et al. 2017). However, these studies found large spreads of ∼ 1 dex around theṀ acc − M disk relation, which could be compatible with a purely viscous model only assuming large spreads on the viscous timescales and typical values of about the stellar population age. As a consequence, the spread around theṀ acc − M disk relation in the old populations of SFRs is expected to be reduced to the typical uncertainties of theṀ acc measurements. However, this was not the case for the > 5Myr old SFR of Upper Scorpius, where a large spread > 0.9 dex and median high values ofṀ acc were found by Manara et al. (2020). This evidence has driven the development of different models for disk evolution, alternative or in parallel to the simple viscous-driven model. For example, recent models predict that disk winds might have a large influence in driving accretion by removing the disk angular momentum excess (e.g., Bai 2016). These wind-driven models moreover appear to be able to reproduce the large spreads observed in thė M acc − M disk correlation. In addition, both internal and external photoevaporation are expected to play an important role in the evolution of the M disk /Ṁ acc ratio with the stellar age and might contribute to increase the scatter of the correlation (Jones et al. 2012;Rosotti et al. 2017;Somigliana et al. 2020). Our finding that TD shows lowerṀ acc values than sources with the same M disk supports the hypothesis that photoevaporation can play an important role. In addition, Zagaria et al. (2022) recently found that binaries or high-order multiple star systems present higher values ofṀ acc than single stars at a given M disk for the SFR populations of Lupus, Chamaeleon I and Upper Scorpius. This can be explained by the fact that disks in multiple stars are subject to tidal truncation, which increases the grain radial drift and reduces the disk lifetime. In our sample we have found two distinct distributions in theṀ acc −M dust,disk plane between single and multiple sources, which is evidence that can be compatible with the scenario described in Zagaria et al. (2022).
Finally, an additional cause for the high dispersion observed in theṀ acc − M disk plane can be found in the lower spread in theṀ acc − Σ disk,dust correlation with respect to theṀ acc − M disk (Fig. 7). Despite the lower statistics of our sample and the small M disk range we probed, it seems that plotting the mass accretion rate as a function of the disk surface density significantly reduces the spread in the correlation. We speculate that the reason might be that millimeter emission is not totally optically thin, which can cause M dust,disk to be not properly estimated. A number of works have already pointed out that M dust,disk might be underestimated when measured from millimeter flux density alone (e.g., Ballering & Eisner 2019;Zhu et al. 2019;Ribas et al. 2020). Millimeter emissions partially or totally dominated by optically thick structures have been invoked to explain the size-luminosity relation found in 50 nearby PPDs (Tripathi et al. 2017) as well as in multiband ALMA observations for a sample of 26 bright PPDs in the Lupus SFR (Tazzari et al. 2021). Ueda et al. (2022) found that in the case of CW Tau, the disk mass can be underestimated due to opacity effects if the measurement is carried out from ALMA bands 6-8 (1.3-0.75 mm), highlighting again the importance of multiband observations and a proper modeling of the dust opacity.
In conclusion, the large spread observed in theṀ acc − M dust,disk plane can be partially due to source multiplicity, as found in Zagaria et al. (2022). However, it might also be due to systematics in the computation of M dust,disk stemming from the adoption of the same dust opacity for all three disk categories. In this case, considering the surface density instead of the disk mass would mitigate this systematics.

Accretion, wind and jets
We have shown the clear correlation between the luminosity of the [O i] LVC and HVC and L acc for the Taurus sample with the use of self-consistent and homogeneous data. We have found a relation similar to that found by Nisini et al. (2018) for the sources in Lupus, Chamaeleon, and σ-Ori in terms of slope and spread, suggesting a common behavior among different SFRs. As pointed out by Nisini et al. (2018), this correlation, which persists even when both luminosities are normalized to the stellar luminosity, suggests a common mechanism connected to the accretion process for these two components. However, we note that the slight difference in the slope between the L [O i],LVC -L acc and L [O i],HVC -L acc relations seems to be real and might support a scenario in which the two components originate from distinct mechanisms, but both are related to the accretion. In particular, while MHD jet-formation models predict a direct correlation between L [O i],HVC and L acc driven by theṀ jet /Ṁ acc relation, the LVC could be dominated by photoevaporation, as suggested by Weber et al. (2020). Models that simultaneously address the physics of PPD winds and jets are necessary to confirm this scenario.
Regarding the kinematics of the jets, the correlation between the [O i]630 nm HVC velocity and accretion luminosity confirms the finding of Banzatti et al. (2019), that is, faster jets are driven by stronger accretion. However, we note that the opposite is not always valid because we found strong accretors with a lack of HVC in the [O i] profile. Two possible interpretation can be at the basis of the observed correlation.
In MHD jets, the mass ejection rate,Ṁ jet , is directly related to the mass accretion rate and thus to the accretion luminosity. M jet can be expressed asṀ jet = M gas v jet /l jet , where M gas is the mass of the gas, and v jet and l jet are the velocity and length of the jet, respectively. Thus, a direct correlation between v jet and L acc can be found provided that the mass of the gas in the considered jet section is constant. This latter condition is hardly met, however. Studies of the physical conditions in some of the considered sources (e.g., Paper I) show that they can vary significantly from one source to the next, with the total density ranging between 10 5 cm −3 and more than 10 7 cm −3 , making a direct proportionality between L acc and v jet difficult.
On the other hand, in the framework of the standard MHD model, the poloidal jet velocity v jet is proportional to (r A /r 0 )v k , with r A the Alvén radius, r 0 the footprint radius (i.e., the radius of the disk where the jet originates), and v k the disk Keplerian velocity at distance r 0 . The latter can be expressed as v K = √ GM /r 0 , thus making v jet proportional to √ M . In this  Figure 11, but for mass accretion rate and stellar mass. Differences between the two methods affect theṀ acc and M directions. Fits to the black and red distributions and black and blue distributions are shown as solid red and blue lines, respectively. Right: Distribution of the dispersion from the individual best-fitting M acc − M relation considering the total sample with the 2 T eff (red) and the single T eff (blue) modeling for the spotted source. The dispersion of the distributions is labeled in the panel.
case, the observed HVC velocity and L acc dependence would be a consequence of the fact that both the v jet andṀ acc correlate with M . In support of this, Fig. 13 shows the deprojected HVC peak velocity as a function of √ M . Except for a few exceptions, we find a moderate correlation between the two quantities, which might be compatible with this scenario. We note that RY Tau, which presents the largest discrepancy with respect to the supposed correlation, is the only source with a transitional disk showing an HVC in the considered sample. In this case, the lower deprojected HVC velocity could be the result of a larger footprint radius, which is expected if the very inner disk region is devoid of both dust and gas. Other cases that deviate from the suggested correlation include four sources (i.e. DQ Tau, DK Tau, GI Tau and IQ Tau) for which the deprojected HVC peak velocities are particularly low (lower then ∼ -80 km/s), thus making the interpretation of these line components as tracers of jets doubtful. On the other hand, three of them were identified in this work as heavily spotted sources. In this case, we might also speculate that the MHD jet acceleration could be significantly reduced due to the higher-order complexity of the stellar magnetic field expected in these sources. A higher statistics and a detailed knowledge of the magnetic field topology are necessary to draw any firm conclusion on this latter point.

Conclusions
In the framework of the GHOsT project we have presented a study of the accretion properties of 37 CTTs of the Taurus-Auriga population with masses > 0.2 M . We used highresolution spectroscopic observations taken with the GIARPS instrument in a wide range from the optical to the NIR to measure the stellar and accretion parameters in a homogeneous and self-consistent fashion. The complementary ALMA data also enabled us to compared these parameters with the physical properties of the disks. The main results of our study are summarized below.
-The standard approach of inferring effective temperatures solely from high-resolution optical or infrared spectra could be inappropriate in case of heavily spotted sources because it can lead to effective temperatures that are systematically too hot or too cold when the measurement is performed from the visible and IR spectral range, respectively. As a consequence, the determination of the stellar luminosity, mass, and age will be strongly affected by this systematics. We showed that the adoption of a 2-T eff modeling of the photospheric emission can mitigate this discrepancy and provides reliable estimates of stellar parameters for stars in which large areas are affected by spots.
-The L acc − L andṀ acc − M relations derived for the Taurus sample are consistent with those measured in SFRs of similar age such as Lupus and Chamaeleon I but with a lower spread in theṀ acc − M relation. This latter could be a consequence of properly addressing the stellar properties of heavily spotted sources. -Transitional disk sources are placed on the lower side of the distributions of accretion luminosity and mass accretion rate when compared to the bulk of the CTTSs with full disks. This is in line with the currently accepted scenario for the protoplanetary disk evolution. However, the complexity of the disk structure is not the main tracer of the accretion level because the largest spread around the L acc −L andṀ acc −M relations is dominated by sources with a full disk. -Multiple star systems display higher values ofṀ acc than single stars at a given M dust,disk , which partially explains the ob-servedṀ acc − M dust,disk dispersion. The high dispersion can also be partially due to a bias in computing the M dust,disk with the same disk opacity for all sources. This latter can be mitigated inspecting theṀ acc − Σ dust,disk dependence, with Σ dust,disk the surface dust mass density of the disk, which is less sensitive to the assumptions on the optical thickness of the emitting dust.    GI_Tau -01.26.2020 Log<L acc /L ô > = -1.37 ± 0.20 A V = 1.9  LkCa_15 -12.14.2020