© The Authors 2024

1 Introduction

Mass outflows represent a fundamental ingredient for understanding the accretion process in young stellar objects (YSOs) and the evolution of proto-planetary disks (PPD). Jets and winds ejected from the disk at different spatial scales remove both mass and angular momentum, and thus they have a fundamental role in driving accretion and setting the timescale of disk dispersal (Pascucci et al. 2022).

Ejections of matter from the inner part of the disks in the form of extended collimated jets have been studied for decades. Their role in extracting angular momentum and thus driving accretion from the disk to the central star is now well recognised (Frank et al. 2014; Ferreira et al. 2006, and references therein). Recent attention has been given to the large influence that disk winds, launched from disk radial distances >0.5 au, might have on the evolution and dispersal of PPD. In particular, the classical paradigm of viscous accretion has been recently questioned due to the fact that both theoretical simulations and observations show that the level of disk turbulence is not high enough to support viscous disk evolution (e.g. Turner et al. 2014; Flaherty et al. 2020). This has led to the recognition that radially extended magnetised winds could be the dominant mechanism driving accretion in planet-forming disks by extracting angular momentum from the disk. Consequently, several theoretical works in recent years have developed wind-driven accretion and disk evolution models (Bai 2016; Suzuki et al. 2016; Lesur 2021; Tabone et al. 2022).

In addition to magnetised winds, the importance of photo-evaporative (PE) winds, i.e. winds thermally driven by stellar energetic photons, has also been recognised in late disk evolution. In particular, PE winds are so far the main mechanism able to explain the very rapid disk dispersal observed at later stages, when the wind mass-loss rates might become comparable to the disk accretion rates (e.g. Alexander et al. 2014; Ercolano & Pascucci 2017). Both type of winds, that is, magnetohydrody-namic (MHD) and PE winds, might compete with each other in driving the evolution and shaping the disk surface density, and thus also influencing the formation and migration of planets (e.g. Suzuki et al. 2016; Kunitomo et al. 2020). It is therefore clear that determining the properties of disk winds is a critical step towards understanding disk evolution and planetary systems’ formation.

From an observational point of view, atomic forbidden lines in a wide spectral range, from UV to IR wavelengths, represent the brightest and easiest way to access tracers of both collimated jets and radially extended disk winds. These lines have been widely used as diagnostic tools to determine the main physical parameters of collimated jets and infer the mechanism for their launching, as well as their role in mass and angular momentum extraction (e.g. Frank et al. 2014; Pascucci et al. 2022). From these previous studies it appears that magneto-centrifugal disk wind models are more suitable to explain the kinematical and dynamical properties of such jets (e.g. Ferreira et al. 2006), including the Ṁ_jet/Ṁ_acc ratio between ∼0.01-0.1 so far derived in large samples of YSOs (e.g. Nisini et al. 2018).

Low-velocity disk winds are more elusive, as they are not spatially resolved by the current instrumentation. Their presence and properties have thus far mainly investigated through high resolution spectroscopy, able to kinematically isolate, in the forbidden line emission profiles, the contribution from the slow winds from that of the high-velocity jets. Indeed, since the original spectroscopic works by Edwards et al. (2013) and Hartigan et al. (1995), it has been common practice to decompose the profile of bright lines, such as the [O I] at 630 nm, in a high-velocity component (HVC), which is usually blueshifted with peak velocities V_p larger than 30 km s⁻¹, and a low-velocity component (LVC), which presents V_p not exceeding a few km s⁻¹. As demonstrated by spatially resolved observations, the HVC is associated with the extended collimated jets, while the LVC appears spatially confined on scales of a few tens of au (Whelan et al. 2021). Spectroscopic observations at resolution −7 km s⁻¹ have further shown that the LVC itself presents a composite profile that can be assimilated to a superposition of two components, called, following Rigliaco et al. (2013), broad LVC (BLVC) and narrow LVC (NLVC) because of their different widths. Further studies conducted on larger samples of classical T Tauri stars (CTTs) suggest that the BLVC and NLVC originate from distinct regions in the disk and that the presence of both components is preferentially found in highly accreting sources, while more evolved disks (for example transitional disks, TD), usually present a single narrow component (Simon et al. 2016; Banzatti et al. 2019; McGinnis et al. 2018). Given that the BLVC is mainly observed in conjunction with the HVC, it has been suggested that it is probably associated with inner MHD winds that feed the high-velocity jets (Banzatti et al. 2019). Less clear is the origin of the NLVC, for which both photoevaporative and MHD winds have been proposed to explain its properties (Fang et al. 2018). Based on the above evidences, Banzatti et al. (2019) proposed an evolutionary trend in the line profiles, where the [O I] 630 nm emission region recedes to larger radii as the inner disk region becomes depleted.

The excitation and physical properties of the LVC have been studied in some details by Natta et al. (2014) and Fang et al. (2018) through an analysis of the [O I] 630 nm, [O I] 557 nm and [S II] 406 nm line ratios. They derive that the observed ratios are best explained by a thermally excited and neutral gas with electron density ~10⁷–10⁸ cm⁻³ and temperatures ~5000–10000 K. Fang et al. (2018) further show that the BLVC and NLVC broadly share similar conditions, although it is very difficult to tightly constrain the physical parameters. Even more difficult is to set constraints on the mass-loss rate, which requires the knowledge not only of the physical conditions in the wind but also of its geometry. Natta et al. (2014), adopting very crude approximations on the wind geometry, estimates Ṁ_wind ≈ 0.1–1 Ṁ_acc. Fang et al. (2018), thanks to the higher resolution of their spectra, provide constraints on the Ṁ_wind in both the NLV and BLV components. They conclude that the mass-loss rate in the BLVC dominates over that of the NLVC, and might be of the same order as the mass accretion rates, hence significantly higher than the HVC and thus jet rates.

In this paper, we aim at further investigating the kinematical and physical properties of the forbidden lines in PPD, and the LVC in particular, through observations at unprecedented spectral resolution (~2.6 km s⁻¹). For that, we use the database of high resolution optical spectra acquired with the GIARPS instrument at the Telescopio Nazionale Galileo in the framework of the GHOsT (GIARPS High-resolution Observations of T Tauri stars) project (e.g. Alcalá et al. 2021; Gangi et al. 2022, hereafter Papers III and IV), a flux-limited survey of T Tauri Stars in the Taurus-Auriga star forming region. We have already used sub-samples of data from the GHOsT project to address two specific topics related to the mass-loss in young stars. In a first paper (Giannini et al. 2019, Paper I), a pilot study was performed on five sources known to drive energetic jets in order to constrain the physical parameters of mass ejection as a function of velocity. A second study (Gangi et al. 2020, Paper II) was devoted to the comparison of the kinematical properties between the atomic NLVC traced by the [O I] 630 nm line, and the molecular winds traced by the H₂ 2.12 μm transitions. The present paper analyses a set of optical forbidden lines from oxygen, sulphur and nitrogen species in the total GHOsT sample of 36 sources. Moreover, this study differs from previous high resolution surveys performed on similar sources (e.g. Simon et al. 2016; Banzatti et al. 2019; McGinnis et al. 2018; Fang et al. 2018) in the following major ways: i) it is based on data with a spectral resolution higher by at least a factor of two with respect to previous studies, and thus allowing us to better disentangle the different kinematical components; ii) stellar and accretion parameters for the sample have been simultaneously and self-consistently determined from the same dataset; iii) for most of the considered sources their disks have been resolved by ALMA, providing the basic information for seeking possible links between the inner gaseous disk, traced by the forbidden lines, and the properties of the outer disks.

2 Sample

In this work, we analyse the forbidden line emission in the sample of YSOs of the GHOsT project. The sample has been observed with the GIARPS instrument that combines the HARPS-N (Cosentino et al. 2012) and GIANO-B (Origlia et al. 2014) high-resolution (115 000 and 50000, respectively) spec-trographs, simultaneously covering a wide spectral range of 390–690 nm for HARPS-N, and 940-2420 nm for GIANO-B. The GIARPS observations have been complemented with contemporaneous low-resolution spectroscopy and photometry, that have been used to flux calibrate the GIARPS spectra. All the details about the observing strategy and data reduction are fully described in Paper IV and will not be repeated here. The sample consists of 36 CTTs showing IR excess from the list of Taurus-Auriga members of Esplin et al. (2014). The source selection was driven by the sensitivity limits of the GIARPS instrument and resulted in a sample covering masses between 0.2 and 2.2 M_⊙ and luminosities between 0.05 and 8.9 L_⊙. The stellar (L_* M_*, R_*) and accretion (L_acc, Ṁ_acc) properties of the sources were consistently derived from the GHOsT flux-calibrated data (Alcalá et al. 2021; Gangi et al. 2022, Paper III and IV) and are listed in Table 1, together with other information taken from the literature. In particular, the disks of most of these sources have been spatially resolved by ALMA or other high resolution millimetre facilities, and Table 1 reports the inferred disk properties that will be used in the analysis. Based on the disk structure as observed in the ALMA images, disks have been divided into sub-structured (showing dust gaps and rings), transitional (TD, with an inner dust cavity) or full (uniform distribution at the resolution of the observations). We point out that this nomenclature refers only to the disk morphology as revealed at mm wavelengths, which is not necessarily associated with an evolutionary stage of the object. For example with this nomenclature we consider as TD also RY Tau, which is a highly accreting object with a jet, and GG Tau A, whose inner cavity is related to the binary nature of the source. The disk of a few of our sources have been also imaged with SPHERE in scattered light, (e.g. Ménard et al. 2020; Mesa et al. 2022; Garufi et al. 2019; Keppler et al. 2020), showing additional structures at large scale, such as spirals, streamers and flaring that we consider in the discussion of individual objects.

3 Spectra of forbidden lines

The GIARPS spectral range covers several forbidden lines from abundant atoms/ions (see e.g. Giannini et al. 2019). In this paper we consider only the five optical lines most relevant for the analysis and observed with the highest signal-to-noise ratio (S/N) in the HARPS-N spectra, which are listed in Table 2. In particular, GIANO IR forbidden lines with high enough S/N have been detected only in the five bright targets analysed by Giannini et al. (2019) and will not be considered here. We corrected the forbidden line profiles for the contribution of the absorption features due to the stellar photosphere, using HARPS spectra of main sequence stars of the same SpT as the sources (see Table C.1 of Manara et al. 2021). In particular, we divided the source spectrum around the line of interest by the photospheric template, corrected for both the vsin(i) broadening and veiling r_λ estimated in Paper IV.

Finally, to transform the spectra to the stellar velocity reference frame, we have corrected them by the stellar radial velocity (RV). This latter has been derived by cross-correlating each star optical spectrum with the synthetic BTSettl model (Allard et al. 2012) that best matched the T_eff derived in Paper IV, adopting the ROTFIT tool (Frasca et al. 2006). The cross-correlation is carried out separately for 14 spectral segments of 200 Å each from 4000 Å, to 6800 Å, after masking the emission lines and the regions strongly affected by telluric absorption lines. The model spectrum is degraded to the spectral resolution of HARPS-N and resampled to the spectral points of the target spectrum. To measure the centroid of the cross-correlation peak, we fitted it with a Gaussian by the procedure CURVEFIT (Bevington et al. 1993), taking the cross-correlation function (CCF) noise, σ_CCF, into account. The latter was evaluated as the standard deviation of the CCF values in two windows at the two sides of the peak. The RV error per each segment, σ_RV, was estimated as the error of the centre of the Gaussian fitted to the CCF (see also Frasca et al. 2017). Table 1 provides the stellar RV and associated uncertainty computed as the weighted mean and weighted standard deviation of the RVs obtained from individual segments, adopting $w=1/{\sigma }_{\text{RV}}^2$ as the weight. The RV uncertainties typically range between 0.2 and up to 2 km s⁻¹ for sources whose spectra are highly affected by emission lines and the depth of the photospheric absorption lines strongly reduced by the veiling. The method could not be applied to RW Aur as very few photospheric lines have been detected. For this source, we shifted the spectrum using the weighted λ_air = 670.7876 nm of the Li I photospheric line as reference (Campbell-White et al. 2023). The resulting corrected profiles are shown in Figs. A.1–A.6. In Table 2, we indicate the rate of detection of each line.

4 Gaussian profile decomposition

We adopt here the commonly used procedure to fit each spectral profile as the sum of different Gaussian functions. Models of emission lines in winds and disks show that the resulting profiles can be rarely described as a single Gaussian, being in most cases asymmetric or with double-peaked profiles (e.g. Weber et al. 2020). Nevertheless, Gaussian fitting helps in the identification of different kinematical components and allows us to compare our results with those of the recent literature. The limitations of the applicability of the Gaussian fitting will be discussed in the following sections.

The Gaussian fitting was already applied by Gangi et al. (2020) on the [O I]630 nm profile on a sub-sample of our sources. We have re-done the profile decomposition in all the sources adopting a slightly different fit procedure.

4.1 Decomposition procedure

The procedure consists of an IDL code based on the χ² minimisation with the MPFIT method and allows to perform multi-Gaussian fit of complex line morphology’s. The procedure provides, for each Gaussian component, width, peak velocities and peak intensity values from which flux and luminosity are computed. For a given profile, the total number of Gaussians is determined as the minimum number of Gaussians yielding a χ² stable within 20% of its minimum value, a criterion already adopted in previous similar studies (e.g. Banzatti et al. 2019). Errors on fit parameters were estimated by following a Gaussian randomisation approach. In short, we simulate for each line profile 10⁴ datasets having a random Gaussian distribution with central value the observed spectral points and standard deviation the local S/N. Each simulated profile is then decomposed and the errors are determined as the sigma of the fit parameter distributions. The error on the fitted peak velocity has been then added in quadrature to the uncertainty of the stellar RV, and reported in Table 1. Our procedure was already successfully applied on both absorption (Gangi et al. 2021) and emission line profiles (Gangi et al. 2023), from the UV to the NIR spectral range.

In this work, we firstly apply the procedure to the [O I] 630 nm line, which is the one observed at highest S/N. Next, we assume the solution of the [O I] 630 nm decomposition as the initial input for the fit of the other lines. Figure 1 shows the result of the procedure for some of the sources in which all the considered forbidden lines have been detected. The decomposition for all the other sources is given in Appendix A. Tables B.1 and B.2 give the result of the decomposition for the different line, reporting the kinematical information (V_p and FWHM, this latter deconvolved for the instrumental resolution) and the luminosity of each component. This latter has been determined from the extinction corrected integrated flux, using the distance and the A_V value given in Table 1 and applying the extinction law by Cardelli et al. (1989).

Fig. 1 Representative continuum-subtracted spectra of the forbidden lines observed in our survey. Flux units are erg s−1 cm−2 Å−1. Red lines shows the Gaussian profiles of the individual fitted components, while the blue line shows the combined fitted Gaussians.

4.2 Velocity components

For the different Gaussian components we adopt the name applied in previous similar studies, i.e. Narrow Low velocity Component (NLVC) if |V_p| < 30 km s⁻¹ and FWHM ≲ 50 km s⁻¹, Broad Low velocity Component (BLVC) if |V_p| < 30 km s⁻¹ and FWHM ≳ 50 km s⁻¹ and High Velocity Component if |V_p| > 30 km s⁻¹.

The separation between the HVC and LVC has a physical justification based on the evidence that the HVC is associated with spatially resolved high-velocity jets (e.g. Hirth et al. 1997) and it has excitation conditions different from the LVC, pointing to a physical distinct region. The distinction between the NLVC and BLVC is more subtle and the adopted width value that separates a narrow component from a broad component is somewhat arbitrary. This is even more evident when different lines of the same sources are considered, since a component which is classified as narrow in the [O I] 630 nm line can have a larger width in the [S II] 406 nm or [O I] 557 nm lines. In Table B.1, we adopt 55 km s⁻¹ as a separation value between a broad and a narrow LVC, which roughly corresponds to the maximum width displayed by [O I] 630 nm profiles when a single LVC component is detected. The only exceptions are DS Tau and IP Tau, where we detect a single line having a width of 90 and 160 km s⁻¹, respectively. Observations with MUSE in these two sources show indeed that the [O I] 630 nm emission is extended and can be associated with a large scale jet (Flores-Rivera et al. 2023).

Fig. 2 Examples of [O I] 630 nm profile types for the LVC. Type 1 is a single low velocity profile. Type 2 is a double low velocity profile, that could be fitted either with a double component (one narrow and one broad) or with a Lorentzian profile. In Type 3 profiles, the LVC is composed by more than one component, not necessarily one narrow and one broad, peaking at different velocities.

4.3 Profile types of the [O I]630 nm line

The above definitions, while having the advantage of giving a simplified interpretation of the profiles, might not reflect the complex dynamical processes of the different outflow manifestations. In more general terms, the observed [O I] 630 nm line profiles can be classified into three different types: Type 1 – a single component profile peaking at low velocity, such as those shown in Fig. 2, upper panel. This is observed in 10 sources (CY Tau, DN Tau, DS Tau, GG Tau, GM Aur, IP Tau, LkCa15, MWC480, V409 Tau and V836 Tau). Type 2 – triangular-like profiles whose LVC is decomposed with two Gaussians, one narrow and one broad, peaking at similar, close to zero, velocity. They could be alternatively fitted with a single component having a Lorentzian shape. They can have or not also a HVC as in the case of CW Tau shown in Fig. 2, middle panel. We remark that emission from a disk extending at small (<1 au) radii could in principle take into account this kind of profiles without invoking an additional component, as we show later. Type 3 – profiles where the LVC is composed by more than one component, not necessarily one narrow and one broad, peaking at different velocities. These represent the majority of the observed profiles. Here the distinction between a BLVC and a HVC is less compelling, as the separation between the LVC and HVC set at 30 km s⁻¹ is somewhat arbitrary and depends also on the velocity correction for the disk inclination. In addition, at the higher resolution of our observations, some of the profile components that were previously fitted with a single Gaussian appear even more structured, with more than one NLVC at increasing velocity shifts, as in the case of DF Tau, DG Tau, and HN Tau.

4.4 Comparison between different lines

The separation between BLVC and NLVC is often more difficult to be recovered in lines different from the [O I] 630 nm. In particular, sources having both low velocity components in the [O I] 630 nm line often have a single profile in the other lines. This is partially due to the lower S/N ratio that prevents to distinguish the two components. In general the LVC is detected less frequently in the [S II] lines, and in half of these cases only the BLVC is observed (see Table 2).

To better analyse whether the adopted decomposition retrieves the same NLVC components for the three brightest lines, i.e. [O I] 630 nm, [O I] 557 nm and [S II] 406 nm, Fig. 3 shows a comparison of the peak velocities and deconvolved FWHM of the [S II] 406 nm and [O I] 557 nm lines, vs. those of the [O I] 630 nm line. In these plots, we mark in blue the sources where the [O I] 630 nm line has been deconvolved in a NLV and a BLV component while the other lines were fitted with a single component. As indicated before, for 9 sources only a LVC has been identified in [S II] 406 nm. Looking at Fig. 4 we note that their centroid velocities tend to be more blue shifted than the [O I] 630 nm line, while the FWHM is wider than that of the [O I] 630 nm line in half of the sources. For the [O I] 557 nm line, the peak velocity is mostly consistent with that of [O I] 630 nm, with few exceptions where the [O I] 557 nm line is less blue shifted than [O I] 630 nm. This latter effect has been already noted (e.g. Hartigan et al. 1995; Giannini et al. 2019) and interpreted considering the different critical densities of the two lines (see Table 2): when the lines originate from an accelerating wind, [O I] 557 nm, which has a critical density an order of magnitude larger than that of [O I] 630 nm, is less blue shifted as it comes from gas at lower velocity in the inner and denser region of the wind. This interpretation cannot be however applied to the [S II] 406 nm line, which has a critical density similar to that of [O I] 630 nm. Figure 4 also shows that FWHM([O I] 557 nm) is compatible with FWHM([O I] 630 nm) except for the sources where it was not possible to separate the broad and narrow components from the [O I] 557 nm LVC.

Regarding the other two considered lines, [N II] 658 nm is detected only in the HVC, as a consequence of the higher ionisation fraction attained in the jets with respect to the neutral gas constituting the LVC (e.g. Giannini et al. 2019). Finally, The [S II] 673 nm line has a critical density an order of magnitude lower than the other lines, and thus intercepts the low density gas in the extended jet better than the LVC high density gas. Hence, its detection in the LVC is also limited (i.e. 14 sources in total).

5 Kinematical properties of the NLVC

We concentrate here on the kinematic properties of the NLVC, seeking any link between them and the properties of corresponding sources and their disk. We use the [O I] 630 nm for this analysis, as it has the largest number of detections at the highest S/N.

Figure 4 presents the deconvolved FWHM of the line versus its peak velocity. Each of the four panels highlights the distribution of sources with different properties. The two vertical dashed lines delimit the maximum uncertainty in velocity centroid attained by our observations. The figure shows that there is no correlation between the line width and line velocity, as already reported in other studies (e.g. Banzatti et al. 2019; Simon et al. 2016). It also shows that in about 15 of the sources, the NLVC V_p is compatible with zero velocity, i.e. for 40% of the sources of our sample there is no direct evidence that the NLVC originates in a wind. Inclination effects have little influence in this trend, as only four of these sources have disk inclined by more than 50 degrees, that would shift their peak velocity by a factor >1.5.

Upper left panel shows that full disks and disks with substructures are equally distributed in this plot, but noticeably all the known TD sources, with the exception of RY Tau and UZ Tau do not show any appreciable velocity shift. We point out that both RY Tau and UZ Tau are rather atypical TD source as, ad variance with the other TDs, they drive strong jets (Agra-Amboage et al. 2009). In addition, upper right panel shows that all the sources with a blue shifted NLVC, but MWC480 and GH Tau, are those with a high-velocity outflow. On the other side, sources with not appreciable velocity shift also comprise several objects with a prominent HVC, such as CW Tau. From bottom left panel, we also see that all the type 1 and 2 profiles, i.e. single LVC profiles, or profiles having a triangular-like shape, fall in the category of lines with no-velocity shifts. Finally, in bottom right panel we distinguish high and low accretors (i.e. sources with a mass accretion rate above or below the median value of the Ṁ_acc = log(Ṁ_acc) < 8.2 Ṁ_⊙ yr⁻¹). It can be noted that on average sources with low V_peak have a lower mass accretion rate with respect to the rest of the sample.

In conclusion, we derive that the NLVC with no-velocity shifts are also those having single peaked or Lorentzian profiles and lower mass accretion rates, and can be preferentially found among TD sources or sources with no high-velocity outflows. Full disks or disks with sub-structures do not show any clear segregation in the plots.

This result supports the view, proposed by Banzatti et al. (2019) and recently confirmed by Fang et al. (2023) on a survey of YSOs of the 5–10 Myr old Upper Sco population, that as disks evolve and accretion starts to settle down, the [O I] 630 nm profile tends to have a single narrow component while the peak velocities move from being mainly blue shifted to being closer to the stellocentric velocity.

The evidence that in more evolved disks the NLVC has peak centroids compatible with zero velocity implies that in these systems the emission mostly originates from gas bounded in the disk or from a very low velocity disk wind, poorly accelerated above the disk surface. Velocity shifts less than ~3 km s⁻¹ are accounted for by the X-ray PE wind models of Ercolano & Owen (2010) only for low inclination angles. More symmetric profiles with respect to zero velocity are predicted for disks with dust inner cavities such as in TD. In these case the wind is seen in both sides within the cavity as its red shifted emission is not blocked by the dusty disk. According to the Ercolano & Owen (2010) simulations however, this effect is appreciable only for sufficiently large inner holes (i.e. larger than about 8 au) and disk inclinations below about 50 degrees, which would not explain the symmetric profiles observed in TD with small inner dust holes or in full disks imaged at high resolution by ALMA (e.g. DR Tau or DN Tau, Long et al. 2019).

Ballabio et al. (2020) compared their analytical PE model with kinematical properties of the [O I] 630 nm emission concluding that observed low V centroids can be reproduced only assuming sound speeds around 3–5 km s⁻¹. These are however inconsistent with the observed FWHM which are better reproduced with larger sound speeds of the order of 10 km s⁻¹.

Finally, also in MHD winds the velocity centroid is significantly more blue shifted except for favourable wind inclinations (Weber et al. 2020).

Therefore wind models are not able to account for the large number of sources with peak velocities close to zero and different disk inclination angles. This might suggest that as the inner disk gas evolves, the contribution of the emission from bound gas becomes dominant with respect to wind emission. In such a case, the emission region should extend to large disk radii, to explain the lack of double peaks in our spectra. This will be further discussed in Sect. 7.2.

Fig. 3 Comparison of Vp and deconvolved FWHM between the [O I] 630 nm line and the [O I] 557 nm (upper panels) and [S II] 406 nm (bottom panels) lines for the NLVC. Red dots indicate sources where the NLVC has been clearly identified in both lines, while blue dots refer to sources where the LVC in the [O I] 630 nm line has been decomposed into a NLVC and a BLVC, while a single component has been identified in the other lines.

Fig. 4 Deconvolved FWHM vs. line peak velocity Vp for the NLVC. The four panels highlight the distribution of sources with different properties. Upper left: disk type, upper right: presence of a HVC, bottom left: types 1, 2 and 3 profiles as described in the caption of Fig. 2, bottom right: value of Ṁaccbelow (filled circles) or above (empty circles) the median value for our sample (log Ṁacc = −8.2).

6 Comparison with previous studies

Banzatti et al. (2019) analysed the kinematical properties of the [O I] 630 nm line in 22 out of the 36 sources of our sample, while Fang et al. (2018) conducted an analysis of the different components of 8 of our sources based on the [O I] 630 nm, [O I] 557 nm and [S II] 406 nm lines. With respect to these previous studies, our observations are conducted at higher resolution (i.e. 2.5 km s⁻¹ against 6 km s⁻¹) allowing one to better isolate the different velocity components and identify winds at lower velocities. In addition, through a comparison with the previous observations, we can analyse any kinematical variability on a temporal scale of a decade or more.

Figure 5 compares the profiles obtained in our observations with respect to those of Banzatti et al. (2019) for the 22 sources in common. The figure shows similarities in most profiles, where small differences are mainly due to different resolution/signal to noise ratio, but also remarkable differences attributable to variability. As also shown by Simon et al. (2016), who compared their observed line profiles with those presented by Hartigan et al. (1995) almost 20 yr before, the most important changes are observed in the HVC, while the profile of the LVC generally remains fairly constant. Sources with complex HVC profiles extending to very large velocities, such as DF Tau, DG Tau and DK Tau, are those with larger variability, both in the relative intensity of the components and in the terminal velocity. These can be partially due to different regions of the extended jets intercepted by the slit widths or fiber size of the instruments, but also to real variability of the ejection mechanism.

Indeed collimated jets consists of moving knots ejected from the star with different velocities and at irregular time intervals, that can account for the different observed peak intensity of the HVC in different epochs (e.g. Takami et al. 2022). DG Tau, in particular, is known to have episodic ejection phenomena on temporal scales of about 1.5 yr (e.g. White et al. 2014) that can be explained by the different observed peak intensity of the HVC. In addition, the terminal velocity of the DG Tau jet has considerably decreased in time. Indeed, the line profile shown by Hartigan et al. (1995), as well as the velocity resolved spectral images obtained more than 20 yr ago (Maurri et al. 2014; Laval-ley et al. 1997) have detected maximum velocities in excess of 400 km s⁻¹, compared with the less than 200 km s⁻¹ measured now.

Noticeable is also the case of DN Tau, where the [O I] 630 nm line remained undetected in the Banzatti et al. (2019) observations. The line was however detected by Hartigan et al. (1995) with a similar profile as the one we observe with GIARPS.

Figure 6 shows the difference in peak velocity of the NLVC between Banzatti et al. (2019) and this work, as a function of V_p. For this comparison, the V_p values reported in Banzatti et al. (2019) have been corrected by the slightly different adopted RV with respect to our work. Different colour symbols indicate both the disk type and the presence of a HVC. Differences in V_p are never larger than a few km s⁻¹ and could be due either to errors in wavelength calibrations or intrinsic variability. We note in particular that the larger differences arise in sources where we measure an appreciable blue or red shifted emission, while sources with a centroid close to zero are more consistent with the Banzatti et al. (2019) values. We also note that peak differences are more common in sources with a HVC. On this bases, we believe that this variability is real and indicates that the kinematics of low-velocity winds is also affected by variations in the jet.

Fig. 5 Comparison of continuum-subtracted and peak normalised [O I] line profiles observed in this work (black) and in Banzatti et al. (2019; red).

Fig. 6 Difference between the peak velocity measured in this work (Vp) for the NLVC and that measured by Banzatti et al. (2019; Vp[B19]), plotted against Vp. The Vp of Banzatti et al. (2019) have been corrected by the slightly different adopted RV with respect to our work. The two dashed lines at ± 1.5 km s−1 indicate the velocity uncertainty of our measurements. Symbols colours refer to the disk type (red: full disks, blue: disk with sub-structures, green: TD) while squares indicate sources with a HVC.

7 Correlations with disk structure and inclination

We discuss here whether there is any relationship among the properties of the forbidden line emission, and in particular the [O I] 630 nm line, and those of the gas and dust distribution in the circumstellar disks, which, for the majority of cases, have been resolved through ALMA observations. Table 1 reports the disk inclination and type, while Table 3 includes additional information relative to the TD. For these latter, it has been shown that inner and outer disks might be misaligned (e.g. Bohn et al. 2022; Francis & van der Marel 2020) therefore this table reports the properties of both inner and outer dusty disks, when available, as well as the size of the dust cavity. Measurements of the inner disk inclination, when available, have much larger uncertainties (up to 20–30 degrees) with respect to the outer disk inclinations estimated from ALMA resolved images. Therefore for the analysis we consider i_disk,in only when it differs from i_disk,out outside the errors, which basically is the case only for GG Tau and UX Tau.

In general, considering also the kinematical analysis described inSect. 5, we did not find any difference in the properties of the forbidden lines between full disk and disks with structures (gaps/rings). Both types of disks share the same kinematical properties in terms of LVC widths and velocity shifts. In addition, jets traced by the HVC are observed in sources with both type of disk structure. Therefore, as for the presence of outflows concerns, there is no evidence of any evolutionary trend between full and structured disks. This also confirms the result shown in Paper IV that there are no differences in accretion properties connected with the structure of the outer disk.

Hereafter, we focus on possible relationship between the [O I] 630 nm NLVC emission line region and the properties of the outer dusty disk.

7.1 [O I]630 nm Keplerian radius

If the NLVC comes from a slow wind or from gas bound in the disk, then it must be subject to broadening due to Keplerian rotation. In such an hypothesis the line width, corrected for the disk inclination, has been often used to have a rough estimate of the characteristic radius of the line emission (e.g. Simon et al. 2016; Banzatti et al. 2019; McGinnis et al. 2018).

To compare our results with those of these previous studies, we plot in Fig. 7 the HWHM divided by the square root of the stellar mass against the sine of the disk inclination, marking lines of constant Keplerian radii (i.e. R_kep = (sin(i)/HWHM)² G M_*) at 0.2, 0.5, 2 and 5 au. In the figure, blue symbols indicate sources having a peak velocity compatible with gas bounded to the disk within 2σ, while cyan symbols are sources with blue shifted velocity centroids, and thus having a clear wind component. The figure shows that sources with zero centroid velocity have Keplerian radii between 0.5 and 5 au, in line with previous studies (e.g. Simon et al. 2016; McGinnis et al. 2018). Sources with blue shifted peaks are more scattered extending further in the inner 0.2 au. We remind that our definition of the NLVC is broader than that of previous works, including sources with a single LVC with FWHM < 55 km s⁻¹ while in previous studies a limit of 40 km s⁻¹ was adopted. In the lower panel of the figure, colour symbols identify disks with different structure. In TD disks, the [O I] 630 nm emitting extends towards larger radii, with R_kep between 1 and 5 au, while there is no distinction for the other disks. Errors shown in Fig. 7 take into considerations only the errors on the fitting procedure. The final uncertainty on the R_kep values should also consider the uncertainty on the disk inclinations and source mass. The latter is of the order of 0.15 in log M_*(Gangi et al. 2022), while the error on i_disk is typically 5 degrees, if measured with ALMA. This results in a final relative uncertainty varying between ~30 to 60%.

Weber et al. (2020) pointed out that in line profiles resulting from wind models, the decomposition in a BLVC and a NLVC might not reflect the presence of distinct physical components originating in separate regions of the disk. This is because a single Gaussian component will contain emission from many different parts of the wind, and this effect is more prominent at large inclinations. However, Type 1 and 2 profiles should be less affected by this problem, as basically a single component is detected and thus contamination by emission at different velocities should be minimised.

Even in the case of emission from the disk, however, the emitting region will extend from an inner to an outer radius, and this latter should be rather large since we do not detect the double peaked profile typical of disk emission. It is therefore important, for a correct interpretation of the R_kep parameter, to evaluate to which extent it can be used as a proxy for the [O I] 630 nm emission region. We analyse this issue in details in the next section.

Fig. 7 Relation between the deconvolved HWHM normalised by the stellar mass and sin (idisk) for the NLVC. Black lines indicate line widths as a function of inclination at disk radii of 0.2, 0.5, 2 and 5 au, in the assumption of Keplerian motion. In the upper panel, blue and cyan symbols distinguish sources where the Vp of the line is consistent with the origin in the disk or consistent with a wind, respectively. In the bottom panel, the distribution with respect to the disk type is shown, with colours code as in Fig. 4.

7.2 Disk radial distribution of the [O I]630 nm NLVC and significance of the R_kep parameter

In this section we analyse whether a simple model of emission from bound gas in Keplerian rotation in the disk can account for the line profiles of LVCs having V_p consistent with zero velocity. This analysis will also allow us to understand to which extent the R_kep parameter can be considered an indicator of the size of the line emission region, as often assumed. To this aim, we modelled the observed line profiles assuming that the line intensity follows a power-law distribution as a function of the radial distance from the star R of the form I(R) = I(R_i) × (R/R_i)^−α, where I(R_i) is the intensity at the inner radius (Fedele et al. 2011).

The radial profile is converted into a velocity profile assuming that the gas is in Keplerian rotation. Finally, the resulting line profile is convolved with a velocity width $V\sqrt{V_{\text{in}}^2+V_{\text{th}}^2}$, where V_in is the instrumental broadening and V_th is the thermal broadening, computed assuming a gas temperature of 8000 K. For the Keplerian rotation and line convolution we used the IDL codes keprot and convolve described in Acke et al. (2005). We then used the IDL routines mpfitfun to find the best fit to the observed line profiles, adopting an uncertainty equal to the rms on the continuum around the line. The parameters that are varied in the fitting procedure are the inner (R_i) and outer (R_o) radii of the emitting gas and the power-law index α. The inner radius is constrained by the maximum velocity reached in the line wings. On the other hand, the outer radius regulates the extent of the double-peak profile characteristic of the disk emission: the larger is the extent of R_o, the smaller is the separation between the two peaks of the line. We do not detect any unambiguous double peak in our line profiles, which implies that R_o should be large or, alternatively, that an additional wind component at very low velocity is present. We applied the model to a sub-set of sources, observed at a sufficiently high S/N, where the NLVC is peaking at V_p ~ 0 and can be easily identified with respect to the other velocity components, and for which we know the inclination angle, needed to model the Keplerian rotation. With these criteria we selected 12 sources. For targets with velocity components in addition to the LVC, we subtracted the corresponding Gaussian fits to isolate the NLVC only. For the sources with Type 2 profiles, we apply the model to the entire LVC (thus to the NLV+BLV components), to analyse if a simple disk model can account for the triangular-like profile of the LVC. To perform the comparison with the model, a small velocity shift (up to a maximum of 3 km s⁻¹) has been applied to sources with V_p not strictly equal to zero in order to have a profile symmetric around zero velocity. The list of analysed sources is given in Table 4 together with the fitted parameters. Figure 8 shows the fitted profiles superimposed on the observed ones. The simple adopted model reproduces quite well the profiles. Looking at the fitted parameters, we see that the radial extent of the emission goes from r < 1 au to very large outer radii, up to more than 50 au, to reproduce the lack of double peak in the profiles. However, for some of the extreme cases of very large R_o, such as CQ Tau, GG Tau and UY Aur, the outer radius is poorly constrained as testified by the large errors provided by the fit. This reflects the fact that above a certain value of R_o the χ² of the fit remains rather constant. On the other side, DN Tau has a rather noisy flat peak that can be fitted with a double peak within the large errors, consistent with a R_o < 2 au. On the other side, the smallest R_i values are found in sources with the triangular-like profile and thus having a large full width at zero intensity (FWZI). The power-law index is always between 2 and 2.5, with the exception of a couple of cases with α < 2.

We point out that the fit provides unrealistically high R_o values for some of the sources, most notably CW Tau and UX Tau, as they are larger than the size of the corresponding disk as resolved by ALMA (e.g. Ueda et al. 2022; Bacciotti et al. 2018). Although the actual gaseous disk can be larger than the ALMA dusty disk, it remains to understand how the [O I] 630 nm line can be excited at such large radii. The presence of a flared disk that intercept the stellar photons at large distances could be considered, however high-contrast imaging SPHERE observations of some of these disks (e.g. Ménard et al. 2020; Keppler et al. 2020) do not provide evidences for them. Alternatively, we should consider that the absence of a double-peak in the [O I] 630 nm profile, at least in some of the sources, is not due to the emission extending at larger radii but rather to an additional wind component at very low radial velocity. For CW Tau and UX Tau, this hypothesis is supported by the rather high inclination angles of their inner disks, which are 59 and 74 degrees respectively.

Using the fitted disk surface brightness we calculate what is the percentage of flux arising within our empirically estimated R_kep radius. The total line flux within R_kep is given by $F\left(R_{\text{kep}}\right)={\displaystyle {\int }_{R_i}^{R_{\text{kep}}}I\left(R\right)\times \text{Ω}\left(R\right)=2\pi \text{\hspace{0.17em}cos\hspace{0.17em}}\left(i\right)/{\text{D}}^2}\times \int I\left(R\right)RdR$, where Ω(R) is the solid angle subtended by a disk ring at the radius R, and D is the distance of the source. This value is then compared with the total flux F_tot calculated between R_i and R_o. The corresponding F(R_kep)/F_tot ratio is reported in Table 4.

We see that, according to the model, a fraction between 40% and 80% of the line emission is confined within R_kep. Taking into consideration the uncertainties on the parameters, we can say that R_kep roughly indicates the size of the region within which, on average, about 50% of the emission arises.

Simon et al. (2016) applied a similar disk model to the [O I] 630 nm profiles of nine sources observed at relatively lower resolution with respect to our data. We find fairly similar results for the sources in common, the main differences being a larger Ro driven by the fact that we still do not see any double peak even at the high resolution of our spectra, and a smaller Ri in those cases, such as BP Tau and DR Tau, where we fit the NLVC and the BLVC as a single line profile.

7.3 Correlation of R_kep with disk properties

7.3.1 Properties of [O I] 630 nm in TD disks

In eight sources of our sample, high resolution sub-mm observations have resolved the structure of an inner dust cavity and in some cases also studied the gas distribution inside it. Table 3 summarises the main parameters derived for these disks. As we have pointed out in Sect. 2, we name all of them as TD, even if the presence of the cavity might not be necessarily related to an evolutionary effect.

We have already shown that the [O I] 630 nm emission in TD generally has a negligible velocity shift and a relatively large Keplerian radius R_kep, between 1 and 5 au. Figure 9 compares R_kep with the size of the dusty cavity R_cav. When two determinations for this latter are available in the literature, both of them are plotted in the figure. For all sources R_kep is smaller than R_cav. Considering that only about 50% of the emission arises from within R_kep, this result indicate that most of the cavity is not fully depleted of gas yet. The figure does not show any direct correlation between R_kep and the size of the dust cavity resolved by sub-mm observations. Banzatti et al. (2019) found a correlation between R_kep, and the spectral energy distribution spectral index between 13 and 30 μm for sources with a single component profile, and interpreted it as evidence that gas and dust in the inner disk are both being depleted in an inside-out fashion. A similar correlation has been also shown for other samples (McGinnis et al. 2018; Fang et al. 2023). While we confirm the general trend of a larger median R_kep radius for TD with respect to sources with full disks, we show that the [O I] 630 nm launching radius is maintained fairly constant among the TD sample and thus it is not directly related to the size of the resolved sub-mmm dust cavity. In this respect, it is noticeable that the only peculiar source in this plot is UZ Tau, which has R_kep < 1 au. We included this source among the TD as a small inner cavity of 9.5 au has been resolved by Long et al. (2019). This cavity however has only a factor of two of dust depletion, therefore this source is a TD less evolved than the others of our sample. This implies that the general trend delineated by the correlation of R_kep with the spectral index, which is a parameter related to long scale time evolution of the disk structure, is not preserved once the a well defined dust cavity is formed.

It should be also noted that the inner dusty cavity of TD disks is, at least partially, depleted of cold (<100 K) gas, as shown by ALMA observations, and sometimes also depleted of warm molecular gas (Carmona et al. 2017; van der Plas et al. 2009). For example, CO mm observations of CQ Tau indicates a gas depletion factor between 10⁻³ and 10⁻¹ inside a cavity of ~20 au (i.e. more than a factor of two lower than the dust cavity; Ubeira Gabellini et al. 2019). Also GM Aur shows an heavily depleted gas inner cavity of ~ 20 au (Bosman et al. 2021). In spite of that, for the TD of our sample the [O I] 630 nm luminosity scales with respect to the accretion luminosity in the same way as for disks without inner cavity (Gangi et al. 2022), indicating no significant differences in the temperature and density of the warm atomic gas on the surface layers of disk.

In fact, as seen in Sect. 8.2, the ratio [O I] 557 nm/[O I] 630 nm in TD sources is similar to those of full disks, suggesting similar physical conditions in the upper layers of the disk from where the lines originate. These evidences indicate that the innermost regions of the disk (i.e. < 5 au) are still not severely depleted by atomic gas, as also testified by the still relatively high Ṁ_acc measured in these and other TD sources (e.g. Gangi et al. 2022; Manara et al. 2023).

Fig. 8 Comparison between the [O I] 630 nm line profiles having Vp close to zero, and a model of emission from a Keplerian disk. In the observed profiles the NLVC (or the NLVC+BLVC in the cases of Type 2 profiles) has been isolated by subtracting the Gaussian fit derived for the other components at different velocities. The fitted parameters of the model are indicated in each panel.

Fig. 9 Keplerian radius of the [O I] 630 nm NLVC emission region for the TD of our sample is plotted against the size of the dusty disk inner cavity. For sources where multiple determinations of the inner cavity are given in literature, they are both displayed and connected with a line.

7.3.2 Correlation between R_kep and inner dusty disk

We here discuss how the [O I] 630 nm emitting region is located with respect to the inner dusty disk, by showing, in Fig. 10 the [O I] 630 nm Keplerian radius versus the dust inner radius. This latter is estimated through different assumptions. In particular, the colour filled symbols in the figure consider the dust sublimation radius, derived from the relationship $R_{\text{subl}}=\sqrt{\frac{L_{\star }+L_{\text{acc}}}{4\pi \sigma T_{\text{suble}}^4}}$ (Dullemond et al. 2001), where T_subl is the sublimation temperature that ranges between ~1500 and 2300 K, which we take here equal to 1700 K (e.g. Alcalá et al. 2021). All other parameters are taken from Paper IV and listed in Table 1. Empty squares refer instead to the magnetospheric truncation radius, i.e. the radius where the disk is truncated by the magnetic field lines and gas material is directly funnelled from here through the accretion columns (e.g. Koenigl 1991). This is estimated taking the stellar parameters from Table 1, and assuming a stellar magnetic field of 1 kG (e.g. Johns-Krull 2007).

Finally, with asterisks we consider the direct measurements of the inner dust radius provided by infrared interferometric facilities, taken from Akeson et al. (2005); GRAVITY Collaboration et al. (2021); Eisner et al. (2007). Interferometric measurements in T Tauri stars measure inner disk radii which are a few times larger than the sublimation radii predicted by disk models (e.g. Eisner et al. 2007; Alcalá et al. 2021). Eisner et al. (2007) have shown that in sources with relatively low mass accretion rates the derived large interferometric radii can be accounted for by magnetospheric truncation. In any case, the [O I] 630 nm Keplerian radii are always significantly larger than the other estimates, which give values typically below than 0.2 au.

We however remark that in Sect. 7.2 we estimated that the [O I] 630 nm emission can extend below the R_kep, down to inner radii ranging between 0.01 and 0.6 au, and thus might extend within the dust free inner disk region. Inner radii <0.1 au are found for sources where the triangular profile has been fitted with a single model that needs to take into account the highvelocity wings with gas at a relatively high Keplerian velocity. Similarly, most of the broad components of Type 3 profiles, if interpreted as Keplerian broadening, imply emission within 0.02 and 0.5 au (e.g. Banzatti et al. 2019; McGinnis et al. 2018) and thus might similarly extend within the inner dusty disk. These considerations indicate that the [O I] 630 nm emission covers a disk region starting from within the gaseous disk and extending well beyond the dust sublimation radius.

Except for sources resolved with interferometry, this plot clearly does not consider cases where dust has been already heavily depleted in the inner region, such as the TD sources, for which we have already shown in Fig. 9 that the [O I] 630 nm emitting region is always inside the dusty hole.

Fig. 10 Keplerian radius of the [O I] 630 nm NLVC emission region is plotted against the estimated and/or measured radius of the dusty disk inner rim, following different approaches. Color circles refer to the dust sublimation radius (colour code as in left upper panel of Fig. 4), squares refer to the magnetospheric radius, while asterisks are the dust inner rim radius measured through near-IR interferometric observations (see text for more details).

8 Line diagnostics in the LVC

Ratios of forbidden lines have been used in different studies to infer the physical conditions of the emitting gas in the LVC. In particular, Natta et al. (2014) considered the excitation conditions of the LVC concluding that the gas should be very dense (>10⁸ cm⁻³), warm (T = 5000–10000 K) and mostly neutral. The limited resolution of their observations prevented them to separately investigate the NLVC and the BLVC. Fang et al. (2018) were able to separately analyse the NLVC and the BLVC, finding similar conditions for both components in agreement with Natta et al. (2014). Both studies conclude that the observed line emission is consistent with thermal excitation, in contrast to a possible origin for the [O I] 630 nm emission from photodissociation of OH, as previously suggested for some of the sources (Gorti et al. 2011; Rigliaco et al. 2013). Giannini et al. (2019) further analysed the variations of physical conditions as a function of velocity in well known jet-driving sources, finding that gas at the lowest velocity is colder and denser than the high-velocity gas, and almost neutral.

Here we aim at addressing again the issue of the diagnostic potential of forbidden lines ratios, taking in particular advantage of the higher resolution of our observations that allows us to better isolate the NLVC with respect to the other components. In addition, we also consider in the analysis the observation of the [N II] 658 nm line, whose ratio with the [O I] can provide upper limits on the ionisation fraction (e.g. Giannini et al. 2019).

To this aim, our dereddened line luminosity ratios have been compared with a statistical equilibrium model for electron collisions of a gas at a single temperature and density (see also Natta et al. 2014; Giannini et al. 2019). The model assumes that the forbidden emission originates from the wind or from the uppermost layers of the disk where molecules are dissociated, and therefore considers that all oxygen is in atomic form and all sulphur singly ionised. These assumptions are considered valid for the LVC also on the bases of i) the non-detection of near-IR lines of [S I] in some of the sources (Giannini et al. 2019) and ii) the low ionisation fraction pertaining to the NLVC implying that the gas is mostly neutral (see Giannini et al. 2019, and next section). Solar abundance values from Asplund et al. (2005) have been assumed. Only collisions with electrons are assumed, taking into consideration that electronic collisional rates are more than three order of magnitudes larger than those for H atoms. Gorti et al. (2009) have for example shown that electrons always dominate collisions if the fractional ionisation x_e > 4 × 10⁻⁴.

We separately consider and discuss the three line ratios.

8.1 The [S II] 406 nm/[O I] 630 nm ratio

The [S II] 406 nm LVC is detected in 19 sources (Tables 2 and 6), and in only nine this detection can be attributed to a NLVC. In Fig. 11 we show the example of DO Tau and UX Tau, where bright [O I] lines are detected in the NLVC but not the sulphur lines. Figure 12 plots the expected [S II] 406 nm/[O I] 630 nm ratio as a function of the electron density n_e. Since [O I] 630 nm and [S II] 406 nm have a similar critical density and differ for ~ 10 000 K in excitation temperature, their ratio varies very little with density, and only a factor of three for temperatures in the range 5000–10000 K. In fact, Fig. 12 shows that the ratio is between 0.2 and 0.6 for T_e ≈ 5000–10000 K and n_e between 10⁵ and 10¹⁰ cm⁻³. We see from the figure that these values are consistent with the detected ratios and a large fraction of upper limits. Some of the upper limits are however below the values predicted by the thermal excitation model. We also point out that if an oxygen abundance lower than the Solar value is assumed, as in Fang et al. (2018), the discrepancy between observed ratios and the slab model would even be enhanced.

If we consider more physical models, such as the PE and MHD winds model by Weber et al. (2020), they predict a [S II] 406 nm/[O I] 630 nm and a [S II] 673 nm/[O I] 630 nm ratio larger than observed by at least a factor of 10 and also larger by a factor greater than five than those predicted by a simple statistical equilibrium model. Therefore we have in general evidence that sulphur lines are sub-luminous in the NLVC of some of the sources or, alternatively, that the [O I] 630 nm emission is enhanced with respect to both wind models and statistical equilibrium codes for collisional excitation.

Also in the Fang et al. (2018) study there are a few sources whose [S II] 406 nm/[O I] 630 nm is overpredicted by models. One of these sources is TW Hydrae, for which Gorti et al. (2009) suggested a non-thermal excitation for the [O I] lines due to OH dissociation. The detected [S II] 406 nm/[O I] 630 nm ratio in TW Hydrae is 0.06, which is in line with some of our observed upper limits and much below the ratio expected from thermally excited gas. Therefore we cannot ruled out that OH photo-dissociation partially contributes to the [O I] 630 nm emission, at least for some of the targets such as CY Tau, CW Tau and DG Tau (see Fig. 12).

Also, Nemer et al. (2020) investigated the possibility that the [O I] 630 nm emission is affected by a contribution due to fluorescence excitation by the stellar UV photons, which has been recently suggested by Gangi et al. (2023) by comparing the [O I] 630 nm and H₂ UV line emissions in young stars of Orion OB1 and s-Ori association, and σ Ori clouds.

Another possibility is that part of sulphur is depleted in the disk with respect to the solar value. Significant sulphur depletion in dense clouds with respect to the ISM is a well known observational result (e.g. Jenkins 2009). In addition, Spitzer observations suggest that in embedded protostars only 5–10% of total sulphur is in the gas phase (Anderson et al. 2013). Sources for which we find the largest discrepancy between the observed [S II] 406 nm/[O I] 630 nm ratio and the thermal excitation model are full disks (e.g. DO Tau and DG Tau), where the NLVC originates from the dusty disk region where depletion might be significant. A similar effect is observed with iron, whose abundance drops to less than 10% with respect to the solar value in the low velocity wind components (e.g. Agra-Amboage et al. 2011; Giannini et al. 2019; Nisini et al. 2016).

Fig. 11 Continuum-subtracted spectra of the [O I] and [S II] lines in DO Tau and UX Tau. The NLVC is well detected in the [O I] lines but not in the [S II] lines, with a [O I] 630 nm/[S II] 406 nm ratio >10.

Fig. 12 Slab model prediction for the [S II] 406 nm/[O I] 630 nm ratio for collisions with electrons. Solid lines refer to gas temperatures of 6000, 8000, and 10 000 K. The cyan area indicates the location of the range of values observed in the LVC (both NL and BL) for the sub-sample of sources where the [S II] 406 nm line has been detected. The grey area shows the location of the upper limits measured in all the other sources. Horizontal lines indicate the more stringent upper limits.

8.2 The [O I] 557 nm/[O I] 630 nm ratio

The profiles of the [O I] 557 nm and [O I] 630 nm lines in the LVC are very similar and their luminosity ratio span between 0.1 and 1. Figure 13 compares the observed [O I] 557 nm/[O I] 630 nm ratios against ne in the statistical equilibrium model. The observed range implies a gas at density >10⁶ cm⁻³ for a temperature in the range 5000–10 000 K, in line with previous studies.

In this plot, we also show the range of values for the [O I] 557 nm/[O I] 630 nm ratio detected in the HVC. This ratio is a few times lower than in the LVC, and implies electron densities ~ 10⁶ cm^–3 for electron temperatures >8000 K.

Given that there is an order of magnitude spread in the ratio for the NLVC, we investigated whether this spread is due to some intrinsic property of the source. We did not find any correlation between the observed [O I] 557 nm/[O I] 630 nm ratio and L_acc or Ṁ_acc. We see however an inverse correlation between the observed ratio and R_kep ([O I]) for those sources having a V_p close to zero, i.e. compatible with emission from the disk (Fig. 14, Pearson correlation coefficient -0.93). This trend was also found by Banzatti et al. (2019) and can be interpreted in terms of a decrease of both disk density and temperature as a function of radius. As already remarked in Sect. 7.3.1, we see that TD sources follow the same relationship as the other sources, indicating that the physical conditions of warm gas within the dusty-free cavity are similar to those of sources where dust is still present. A linear fit through the points in Fig. 14 suggests a relationship of the form: [O I] 557 nm/[O I] 630 nm = (−0.17 ± 0.05)R_kep − (0.20 ± 0.12).

Fig. 13 Slab model prediction for the [O I] 557 nm/[O I] 630 nm ratio for collisions with electrons. Solid lines refer to gas temperatures of 5, 8 and 10 × 103 K. The cyan area indicates the location of the range of values observed in the NLVC for the sub-sample of sources where the [O I] 557 nm line has been detected in this component. The grey area indicates instead the location of the detected ratios in the HVC.

Fig. 14 The [O I] 557 nm/[O I] 630 nm ratio observed in the NLVC with no velocity-shifted peak is plotted against the estimated [O I] 630 nm emission radius assuming Keplerian rotation. Green symbols indicate TD sources. The dashed line represents the linear fit through the points, whose analytical expression is indicated in the figure. The Pearson correlation coefficient is also displayed in the bottom right of the figure.

8.3 The [N II] 658 nm/[O I] 557 nm ratio

The [N II] 658 nm/[O I] 557 nm ratio can be used to constrain the degree of ionisation in the wind. This is a very important parameters both to derive the total gas density (n_H = n_e / x_e), and to constrain models that predict a different degree of ionisation in the wind. The [N II] 658 nm line is detected in our spectra only in the HVC, in agreement with the finding that the LVC is mostly neutral (e.g. Hartigan et al. 1995; Natta et al. 2014) and that the ionisation increases with the jet velocity (e.g. Giannini et al. 2019; Bacciotti et al. 2000). We can therefore use the upper limits on the [N II] 658 nm in the NLVC to derive limits on the ionisation degree in this component.

To this aim, we have used a ionisation model that consider as processes collisional ionisation, radiative and dielectronic recombination, and direct and inverse charge-exchange with hydrogen. Details are given in Giannini et al. (2015).

Figure 15 plots the ionisation equilibrium model predictions for the [N II] 658 nm/[O I] 630 nm ratio as a function of x_e, for different values of gas temperature and density. In this plot, horizontal lines with arrows indicate some of the upper limits derived from our data. They show that for the most stringent limits (e.g. in CW Tau), x_e is constrained to be <0.1 for a density of 10⁸ cm⁻³ and T = 6000 K or <0.01 for n = 10⁶ and T= 10000 K. These upper limits are consistent with the values of x_e between 0.05 and 0.1 derived in the LVC of DG Tau, HN Tau and DO Tau from the analysis of multiple optical/IR forbidden lines (Giannini et al. 2019). These values also imply total densities of the order of 10⁸–10⁹ cm⁻³. For the objects where the gas emission originates in the disk, this density corresponds to that predicted by thermo-chemical models in the inner disk region (r < a few au; e.g. Woitke et al. 2019).

The bottom panel of Fig. 15 shows that in the HVC, where the [N II] 658 nm line is often detected, the [N II] 658 nm/[O I]630 nm ratio implies x_e between ~0.1 and 0.5, for a density 10⁵–10⁶ cm⁻³ and T = 6000–10000 K, in line with the values estimated from optical observations of resolved jets, which are in the range 0.3–0.5 (Frank et al. 2014).

8.4 The [S II] 673 nm line

The [S II] 673 nm line has a critical density two order of magnitudes lower than the [O I] and [S II] lines, and therefore it is able to intercept low density wind regions located further out with respect to the other lines. This line has been detected more frequently in the HVC (36% of sources) while only five sources show the NLVC. In these cases, the [S II] 673 nm profile is similar to that of [S II] 406 nm both in peak velocity and FWHM. However, the observed [S II] 406 nm/[S II] 673 nm ratio in the range 1–10 is consistent with a density <10⁵–10⁶ cm⁻³ i.e. lower by at least 1–2 order of magnitudes with respect to what implied by the [O I] 557 nm/[O I] 630 nm ratio. Presumably [S II] 673 nm intercepts larger regions of the wind at a too low density to significantly excite [O I] 630 nm. This has been also demonstrated by spectro-astrometric observations of [O I] 630 nm and [S II] 673 nm performed on RU Lup (Whelan et al. 2021), showing that the [S II] 673 nm in the NLVC is extended in a wind region of 40 au, up to five times larger than [O I] 630 nm. Both the photoevaporative and MHD models by Weber et al. (2020) predict a [S II] 673 nm/[O I] 630 nm ratio close to 1, while in the five sources where [S II] 673 nm is detected the observed ratio is in the range 0.1–0.2. As for the case of the [S II] 406 nm line, the general lower than expected ratio of sulphur lines with [O I] 630 nm can be explained either assuming that sulphur is depleted by dust in the upper layer of the dusty disk atmosphere or by a non negligible contribution of UV pumping to the [O I] 630 nm excitation.

9 The mass-loss rate in winds and jets

The mass-loss rate is one of the fundamental parameters governing the evolution of proto-planetary disks as it sets the timescales for disk dispersal. Deriving the mass-loss rate for all the components gives us a measure of how much material is removed from the disk upper layers as a function of the radial location on the disk. The derivation of the mass-loss rate of the atomic outflow components from single on-source spectra is severely hampered by the lack of spatial information needed to resolve the jets/winds. Nevertheless, flux-calibrated high resolution spectroscopic data are very powerful to investigate the mass-loss from the different kinematic components of the wind in large sample of sources (e.g. Hartigan et al. 1995; Nisini et al. 2018; Fang et al. 2018).

Here we discuss the methods used to infer the mass-loss rate in both the low and high-velocity components from single source high resolution spectra, discussing their caveats and limitations and how they might affect the results.

Fig. 15 Ionization equilibrium model predictions for the [N II] 658 nm/[O I] 630 nm ratio for collisions with electrons, plotted as a function of the xe = ne/nH ratio. Different lines refer to different values of electron temperature and density as indicated in the legend. In the upper panel, horizontal lines indicate some of the more stringent upper limits of the [N II] 658 nm/[O I] 630 nm ratio in the NLVC, while in the lower panel the coloured region indicates the range of the [N II] 658 nm/[O I] 630 nm ratio observed in the HVC.

9.1 Ṁ_Wind in the LVC

The determination of the mass-loss rate in the wind traced by the LVC is subject to large uncertainties due both to the poorly determined physical conditions but also to the unknown physical extent and geometry of the wind.

Rough estimates of Ṁ_wind in the LVC have been provided by Natta et al. (2014), who computed the volume of the [O I] 630 nm line emission and considered the extreme cases of spherical and cylindrical geometry. A similar approach was used by Hasegawa et al. (2022). Fang et al. (2018) separately estimated the Ṁ_wind from the NLVC and the BLVC, assuming that the [O I] 630 nm line emission is in LTE and that the wind length is proportional to the Keplerian radius inferred from the line FWHM. They concluded that on average, the Ṁ_wind /Ṁ_acc ratio in the NLVC is around 0.01–0.1, if one assumes a wind length comparable with the wind radius at the base (i.e. around 1-2 au). This work also suggests that the Ṁ_wind in the BLVC is higher than that of the NLVC, basing this latter statement on the assumption that the broad component has the same or smaller wind length than the narrow component.

Here we revisit the Ṁ_wind estimate in the NLVC using a different approach, which is also subject to large caveats but allows us to analyse the issue from a different perspective.

The mass-loss rate can be written as the mass that pass through a cross-sectional area in the unit of time

$${\dot{M}}_{\text{wind}}=\mu n_H\pi R_w^2{v}_{\text{w}}$$(1)

where n_H is the total density and R_w the radius at the wind base.

For our rough estimate we assume that R_w is the emitting radius estimated as the Keplerian radius from the [O I] 630 nm line width. The total density is given by n_H = n_e /x_e, where n_e is the density at the wind base, that we derive from the [O I]557 nm/[O I] 630 nm ratio assuming T = 8000 K (see Sect. 7.2). We consider x_e = 0.1, which is the upper limit to the ionisation fraction derived in Sect. 7.3, while the wind velocity v_w is taken from the wind peak velocity corrected for the inclination angle. We measure the mass-loss rate only for those sources having a blue shifted velocity peak, outside the wavelength calibration error. In Fig. 16 we plot the so derived mass-loss efficiency, i.e. the Ṁ_wind /Ṁ_acc ratio, as a function of the source mass accretion rate. The mass-loss rate in most cases is between 0.1 and 1 times Ṁ_acc although a few sources show lower values of the order of 0.01, without a clear dependence on Ṁ_acc.

The uncertainties associated to the above derivation are however large. First of all the total density is poorly constrained since we have only an upper limit on the fractional ionisation. Giannini et al. (2019) estimated, from an analysis that takes into consideration a larger set of lines in DG Tau, HN Tau and DO Tau, that x_e in the LVC can be as low as 0.03, which would increase the final Ṁ_wind value by up to a factor of 3–4. On the other hand, the electron density giving rise to the [O I] 630 nm emission could be overestimated if the [O I] 557 nm line, with a critical density an order of magnitude larger than [O I] 630 nm, traces a much more compact region in wind with a large density gradient (see discussion in Sect. 9.1.1). Also, the estimated mass-loss rate would be lower if the emitting area is confined to a smaller section of the disk around the wind launching region. For example, if we assume that the region involved in the launching is ΔR/R ~ 0.1, then Ṁ_wind would decrease by about an order of magnitude. Overall the above uncertainties bring uncertainties on Ṁ_wind to more than one order of magnitudes on either sides.

We also consider that the line luminosity can be used to measure the volume of the oxygen gas giving rise to the observed emission (Natta et al. 2014, see e.g.):

$$V=\frac{L\left(\left[\text{O\hspace{0.17em}I}\right]630\right)}{n_Hϵ\left(\left[\text{O\hspace{0.17em}I}\right]630\right)\text{\hspace{0.17em}X}\left(\text{O}\right)}$$(2)

where ϵ([O I]630) is the emissivity of the [O I] 630 nm line in erg s⁻¹, and X(O) is the total abundance of neutral oxygen, for which we consider the solar value. With the approximation of a cylindrical geometry and a wind area at the base equal to $\pi R_w^2$, we have from here a rough estimate of the wind length above the disk. The estimated wind length is small, in many cases only a fraction of au, with a typical l_w/R_w ratio between 0.1 and 1.

Our estimate of Ṁ_wind for the NLVC is typically larger than the values derived by Fang et al. (2018), mainly due to the fact that we infer a l_w/R_w < 1 while Fang et al. (2018) assume that the wind height is larger than the corresponding Keplerian radius. We have however pointed out that given the large caveats in the determination, any estimate for Ṁ_wind in the NLVC can span about two order of magnitudes, from ~0.05 to more than ~1 the Ṁ_acc, value.

Fig. 16 The ratio Ṁloss /Ṁacc, plotted as a function of the source Ṁacc. Upper panel: Ṁloss for the narrow low velocity component estimated as in Sect. 9.1 and assuming a fractional ionisation xe = 0.1. Lower panel: One sided Ṁloss for the HVC. In cases where the HVC can be decomposed into two additional components at medium and high-velocity, the mass-loss has been separately estimated, with the medium velocity indicated with blue symbols.

Fig. 17 Map of the electron density of a photo-evaporative wind model from W20, with overlaid contours of the 80% emission regions of [O I] 630 nm (green) and [O I] 557 nm (magenta).

9.1.1 Empirical determination of Ṁ_wind from a PE model

In order to check how well our determination of Ṁ_wind from the [O I] 630 nm line can reproduce the actual mass-loss in the wind, we performed a test on the models presented in Weber et al. (2020, hereafter W20), who post-processed the output of the X-ray PE wind model by Picogna et al. (2019). To this aim, we apply the procedure described in Sect. 9.1 to the synthetic line profiles and luminosities given as the output of the model, and verify a posteriori whether the derived parameters are consistent with the theoretical physical and kinematical structure of the wind as predicted by the model.

As described in details in W20, the various forbidden lines trace different regions of the wind, depending on their excitation conditions. To illustrate this point, we show in Fig. 17 a density map for the model output that considers irradiation from a source with accretion luminosity L_acc = 2.6 × 10⁻² L_* (model PE1 in Table 1 of W20). Green and magenta lines delineate the wind regions where the [O I] 630 nm and [O I] 557 nm lines emit 80% of their luminosity, respectively. It can be seen that the [O I] 557 nm line is emitted from a region much more compact than [O I] 630 nm and thus it traces only the very base of the wind. Table 5 reports the physical parameters obtained by averaging the values within the emitting regions of the two lines, for two models with different L_acc, i.e. L_acc = 2.6 × 10⁻² L_* and L_acc = L_*, which also roughly correspond to the extreme values that we estimate for our sample.

The average densities significantly differ for the emitting regions of the two lines, and consequently the n_e value estimated from the [O I] 557 nm/[O I] 630 nm ratio is in between the two values, and thus overestimating the density in the [O I] 630 nm emitting region, typically by a factor of ten.

To test the method described in Sect. 9.1, we have derived R_w from the Gaussian fit of the synthetic line profile and the electron density from the [O I] 630 nm/[O I] 557 nm predicted luminosities, and applied Eqs. (1) and (2) to derive Ṁ_wind and the length of the emission region. For this test, the W20 model that assumes a total wind mass-loss rate Ṁ_wind = 10⁻⁸ M_⊙ yr⁻¹ and a disk inclination angle of 60 degree has been considered.

Table 6 summarises the output of this test for the two accretion luminosity models. For the model with low accretion luminosity, assuming the density from the [O I] 557 nm/[O I] 630 nm ratio and x_e = 0.1, the total density is overestimated by a factor about 10 with respect to the average in the model. On the other hand, the Keplerian radius estimated from the line width gives a slightly underestimated size of the emission region. Consequently, the final mass-loss rate derived from Eq. (1) is similar to the mass-loss in the [O I] 630 nm emitting region predicted by the model, i.e. ~2 × 10⁻¹⁰ M_⊙ yr⁻¹, but a factor of 100 lower than the total wind mass-loss rate of the model. This is because the wind is much more extended than the line emission (up to >100 au) and the [O I] 630 nm line traces only the inner part of it (12 au). If we then consider our estimate of the emitting volume and length for the [O I] 630 nm emitting region, we see that we underestimate both quantities, as we adopt a higher total density with respect to the average value in the model. In particular, the [O I] 630 nm emitting region extends up to about 3 au above the disk, while the measured length with our approximation is about 10 times smaller. This discrepancy is also due to the fact that the ‘empirical’ volume is compared to the 80% emission volume of the model under the assumption of uniform emission, when in reality the emission in the model is higher in the parts close to the star, which contribute only little to the volume, and lower in the outer parts of the emission region, which contribute the most to the volume.

Considering now the other extreme case of L_acc = L_*, we see that a higher accretion luminosity increases the extension of the [O I] 630 nm emission region (see also discussion in W20). In this case, the Ṁ_wind inferred from observations is about 15 times larger than the mass-loss rate in the [O I] 630 nm emitting region predicted by the model, but, noticeably, very close to the total mass-loss rate in the entire wind.

Therefore our conclusion is that the empirical estimates given in Sect. 9.1 can be considered as an upper limit on the mass-loss rate in the [O I] 630 nm emitting region, but a lower limit to total wind mass-loss rate. We remark that although this test has been performed in the assumption that the X-ray PE wind model by W20 gives a close to real representation of the physical structure of the observed low velocity wind, the results are more general. In fact, the main assumption here is that the physical parameters, and density in particular, have a large gradient in the wind, therefore the information inferred from a single, spatially unresolved line are intrinsically limited.

9.1.2 Limitations of [O I] 630 nm as a diagnostics for the LVC mass-loss rate

Given the above discussion, we summarise here the limits and uncertainty of the mass-loss rate determination based on the [O I] 630 nm line emission.

In general, the main limitation is that the wind geometry and extent are not known, which prevents to correctly estimate the length scale for the mass transfer. A second important limitation is the fact that large gradients in physical conditions are expected in the wind, and thus the lines whose ratio is used to derive the density intercept portions of the wind with different densities and temperatures. Finally, the choice of the tracer for the mass-loss is also important. In particular, we have shown that the [O I] 630 nm line might trace only the very inner wind region, which not necessarily gives the main contribution to the mass-loss.

Given the above considerations, our estimate for Ṁ_wind in the NLVC might easily span about two orders of magnitude, from ~0.05 to more than ~1 the Ṁ_acc, value. If one takes the X-ray PE model described in Weber et al. (2020) as a reference, these values should be considered only as lower limits to the global mass-loss in the wind. Therefore we cannot exclude the possibility that the mass-loss in winds ejected at radial distances larger than 0.5 au is of the same order of magnitude or even larger than the mass accretion rate.

Fang et al. (2018) suggested that the Ṁ_loss/Ṁ_acc in the BLVC is higher than that of the NLVC and thus that most of the wind mass-loss occurs close to the central star. This last conclusion is mainly driven by the higher luminosity and peak velocity of the BLVC with respect to the NLVC. While this is an objective observational difference, we think that the limitations discussed above do not allow one to draw solid conclusions on the difference in the mass-loss between the NLVC and the BLVC. In particular, the Fang et al. (2018) analysis is based on the assumption that the gas in both the NLVC and the BLVC is in LTE, which implies n_e > 10⁷ cm⁻³. However, if the density in the [O I] 630 nm emitting region is lower and different in the two components that trace different disk regions, then the Ṁ_wind estimated from the line luminosity will scale accordingly.

Also, the contribution to the mass-loss rate from different gas components, not traced by the optical lines, might be significant in the outer disk regions, for example in regions where the temperature drops below about 5000 K. For example, wide angle uncollimated molecular winds have been resolved in the H₂ near-IR lines (e.g. Agra-Amboage et al. 2011; Beck & Bary 2019) and show similar peak velocities as the [O I] 630 nm NLVC lines, but narrower FWHM, indicating they originate at larger disk radii (Gangi et al. 2020). Also, recent high resolution ALMA observations in HH30 and DG Tau B have shown the presence of cold CO molecular wide winds, whose estimated mass-loss rates are of the same order or even exceed the mass accretion rates for the sources (Louvet et al. 2018; de Valon et al. 2020). In conclusion, we think that our knowledge of the mass-loss rate in winds based on observations of optical forbidden lines is too uncertain to allow one to discriminate among different models and to evaluate the impact of the atomic winds on the overall disk dispersal.

9.2 Ṁ_loss in the HVC

We can use the luminosity and kinematical properties of the HVC to infer the mass-loss rate in the jet. This can be expressed as Ṁ_jet = M_gas × V_jet/l_jet, where V_jet and l_jet are the jet velocity and length, respectively. As for the LVC, here the jet velocity Vjet is estimated from the radial velocity of the line in exam, corrected for the jet inclination angle. We can however now remove the uncertainty on the geometry if we assume that the jet is more extended than the instrumental FoV, so to use the aperture size as the scale length in which M_gas is measured.

Unless estimates of jet width and total density is available for the mass measurement, M_gas is usually derived directly from the luminosity of the forbidden lines, assuming uniform temperature and density in the considered aperture. Hartigan et al. (1995) described in details the method for optical lines, which has been then applied in several subsequent studies.

Unless only a HVC red shifted component is detected, as in the case of BP Tau and GI Tau, we consider here only the blue shifted component and use half of the aperture width, which will give us the mass-loss rate only in the one-sided jet. In the assumption that the jet is symmetric but the red shifted is poorly visible due to a higher extinction or occlusion by the disk, the derived Ṁ_jet value needs to be multiplied by 2. For all our sources, the blue shifted HVC can be decomposed in one or two Gaussians. We separately derive Ṁ_jet for the observed components, that we generically name high- and medium- velocity components (as done in previous studies for resolved jets, as for example Agra-Amboage et al. 2011).

The multi-line analysis performed by the pilot study of Giannini et al. (2019), shows that the typical conditions pertaining to the HVC are n_e < 4 × 10⁴–10⁵ cm⁻³ and T_e ~ 10000–15000 K. We assume here n_e = 5 × 10⁴ cm⁻³ and T_e = 10000 K, on the basis that the line ratios in other objects do not differ within the error from those of this pilot study, and to be also consistent with the Ṁ_jet values measured in a similar way by Nisini et al. (2018) in a sample of Lupus and Chamaeleon sources. Considering however the emissivity variations in the range of density and temperature values, the uncertainty on the derived Ṁ_jet is almost a factor of ten.

Table 7 reports the derived Ṁ_jet for the sources showing a HVC, together with the derived values of the total jet velocity.

In Table 8 we collect Ṁ_jet derivations from other studies based on spatially resolved observations of the jet, to discuss whether they are comparable with our measurement of the mass-loss done based on individual spectra. Variations among different determinations for the same jet are typically a factor of a few and are mostly due to different methods, tracers and assumptions. For example the largest differences are observed on RW Aur A and CW Tau, where the Ṁ_jet in the resolved jet was estimated assuming the electron density from the ratio of the [S II] 673 nm doublet, which provides only a lower limit to the actual gas density (e.g. Maurri et al. 2014). We should however keep also in mind that differences in Ṁ_jet could be due to time variability of the mass-loss phenomenon (e.g. Takami et al. 2022), as seen in Sect. 5 for some of our sources. For example, the DG Tau [O I] 630 nm luminosity in the high-velocity component has more than doubled with vrespect to the value prior to 2012. For this specific case, however, the terminal velocity in 2012 was also higher than our measured value (-400 km s⁻¹ vs. 200 km s⁻¹). Therefore, the linear momentum M × V, and thus the Ṁ_jet, remained almost constant in time.

In Fig. 16, bottom panel, we plot the ratio Ṁ_jet/Ṁ_acc for the HVC with respect to the source mass accretion rate. We separate with red and blue symbols the mass-loss rates values relative to the high and medium velocity components, respectively. As also found by Nisini et al. (2018), the ratio spans more than an order of magnitude and it is mostly confined between 0.01 and 0.1, with a no clear dependence on the mass accretion rate.

10 Conclusions

We have analysed the profile and luminosities of optical forbidden lines in 36 CTTs of the Taurus-Auriga cloud whose disks span a range of different characteristics, and compared the properties of the inner atomic gas with the structure and properties of the outer disks. Following previous studies, we have decomposed the line profiles into HVC, NLVC, and BLVC components, focusing our analysis in particular on the NLVC. We used ratios of different forbidden lines to constrain the physical conditions of the atomic gas (temperature, density and ionisation fraction). For sources whose [O I] 630 nm profile peaks at the stellar velocity, we fitted the profile with a simple Keplerian disk model, giving constraints on the size of the emission line region. For those sources with blue shifted velocity centroids, indicative of an origin in slow winds, we constrained the wind mass-loss rate discussing the origin of all uncertainties in its determination. We finally derived the mass flux rate in the HVC and compared it to the values derived in the same sources from spatially resolved observations of the jets. Our main results can be summarised as follows:

• We find that in about 40% of our sources the peak velocity (V_p) of the NLVC is compatible with the stellar velocity. Such sources comprises all the transitional disks (TD) and are typically those having a single LVC, lower mass accretion rates and absence of a HVC. They therefore might represent later evolutionary stages where the emission from the disk is dominant with respect to the wind contribution. Full disks and disks with sub-structures do not show differences in the forbidden line properties;

• Assuming, for the lines with no-velocity shift, a Keple-rian disk model with a brightness scaling as a power-law index, we derive that the emission line regions extend from -^<0.01 au, and thus within the dust sublimation radius in some of the cases, up to several tens of au in a few sources. In these models, the characteristic radius Rkep, estimated from the HWHM of the line, represents the region within which on average 50% of the emission is confined;

• We find that the R_kep(OI) in TD does not increase with the size of the sub-mm dust cavity, therefore the clearing of inner warm atomic gas probing the upper disk surface layer does not follow the dispersal of cold molecular gas and dust during the late disk evolution. We have found, in addition, an anti-correlation between the [O I] 557/630 nm line ratio and R_kep, which suggests that the emitting region of [O I] is moving outwards in radius as the gas dominating the emission becomes cooler and less dense;

• We confirm previous findings that the line ratios observed in the LVC, if compared with a thermal single temperature and density model, implies n_e ~ 10⁶ – 10⁸ cm⁻³ and T_e ~ 5000–10 000 K, and we moreover constrain the ionisation fraction in the NLVC to be x_e < 0.1. We however discuss the limits of such a diagnostic if applied to the emission in spatially not resolved winds, where physical conditions have large spatial gradients, as the different lines probe different regimes of density and temperature;

• We constrain the mass-loss rate in the wind traced by the blue shifted LVC finding a_Ṁwind /Ṁ_acc, ratio that might span between ~0.01 and more than 1. When compared with synthetic [O I] 630 nm images from X-ray photoevaporation models, the estimated_Ṁwind represents a lower limit to the total rate of mass-loss of the model, indicating that the [O I] 630 nm line is probably not the best tracer to probe the mass-loss in winds;

• For the HVC, we found Ṁ_jet/Ṁ_acc ratios in the range 0.01–0.1, in line with previous findings in other star forming regions (for example Lupus and Cha, Nisini et al. 2018). The Ṁ_jet derived from the on-source spectra are consistent, within a factor of a few, with the values inferred in spatially resolved jet observations, providing that similar assumptions are used.

Two main conclusions can be drawn from this work. Firstly, we have shown that contribution from disk emission, in addition to disk winds, to the forbidden lines in T Tauri stars can be significant, especially for more evolved sources. We also conclude that without a better knowledge of the wind geometry and spatial extent, and given the limitation of the diagnostics, mass-loss rates in the wind traced by the blue shifted LVC cannot be constrained better than a factor of 100.

For a step forward in our understanding of this issue, one would need observations at higher spatial resolution, to significantly constrain the wind emitting region, but also of tracers at lower excitation likely probing winds at larger disk radii. In this respect, future JWST observations will allow us to investigate in detail the range of conditions, between 1000 and 5000 K, presently poorly traced by ground-based observations. In addition, the determination of the mass-loss rate should involve modelling of both line profiles and luminosities of the different tracers with appropriate wind models.

Acknowledgements

We thank Andrea Banzatti for providing the Keck [O I] 630 nm spectra and for useful discussions, and Davide Fedele for providing the model of line profiles from a Keplerian disk. We also thank Barbara Ercolano for helpful insights about photo-evaporative wind models. We finally thanks the anonymous referee for the useful comments and suggestions. This work has been supported by the projects PRIN-INAF 2019 “Spectroscopically Tracing the Disk Dispersal Evolution (STRADE)”, PRIN-INAF 2019 “Planetary systems at young ages (PLATEA)” and by the Large Grant INAF 2022 YODA (YSOs Outflows, Disks and Accretion: towards a global framework for the evolution of planet forming systems). C.F.M. is funded by the European Union (ERC, WANDA, 101039452). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. This work benefited from discussions with the ODYSSEUS team (HST AR-16129), https://sites.bu.edu/odysseus/.

Appendix A Profiles decomposition

Figures A.1–A.6 show, for each of the 36 sources of our sample, the continuum subtracted spectra of the five lines analysed in this paper, with overlaid the derived fit to the profiles obtained summing up one or more Gaussians.

Fig. A.1 Continuum-subtracted line profiles. Fluxes are in units of 10−15 erg s−1 cm−2 Å−1. The fit to the profile is shown in blue and it is obtained by adding the single or multiple Gaussians indicated in red. The two horizontal red dashed lines indicate ± 1 σ.

Appendix B Tables of line parameters

Tables B.1 and B.2 report the kinematical parameters of the line profile decomposition and the calculated luminosity for each component, corrected for the extinction values given in Table 1. The velocity V_p reported for each component has been corrected for the stellar radial velocity given in Table 1, while the Full Width at Half Maximum (FWHM) is corrected for the instrumental profile broadening.

References

All Tables

All Figures

	Fig. 1 Representative continuum-subtracted spectra of the forbidden lines observed in our survey. Flux units are erg s−1 cm−2 Å−1. Red lines shows the Gaussian profiles of the individual fitted components, while the blue line shows the combined fitted Gaussians.
In the text

	Fig. 2 Examples of [O I] 630 nm profile types for the LVC. Type 1 is a single low velocity profile. Type 2 is a double low velocity profile, that could be fitted either with a double component (one narrow and one broad) or with a Lorentzian profile. In Type 3 profiles, the LVC is composed by more than one component, not necessarily one narrow and one broad, peaking at different velocities.
In the text

Fig. 3 Comparison of Vp and deconvolved FWHM between the [O I] 630 nm line and the [O I] 557 nm (upper panels) and [S II] 406 nm (bottom panels) lines for the NLVC. Red dots indicate sources where the NLVC has been clearly identified in both lines, while blue dots refer to sources where the LVC in the [O I] 630 nm line has been decomposed into a NLVC and a BLVC, while a single component has been identified in the other lines.

In the text

Fig. 4 Deconvolved FWHM vs. line peak velocity Vp for the NLVC. The four panels highlight the distribution of sources with different properties. Upper left: disk type, upper right: presence of a HVC, bottom left: types 1, 2 and 3 profiles as described in the caption of Fig. 2, bottom right: value of Ṁaccbelow (filled circles) or above (empty circles) the median value for our sample (log Ṁacc = −8.2).

In the text

	Fig. 5 Comparison of continuum-subtracted and peak normalised [O I] line profiles observed in this work (black) and in Banzatti et al. (2019; red).
In the text

Fig. 6 Difference between the peak velocity measured in this work (Vp) for the NLVC and that measured by Banzatti et al. (2019; Vp[B19]), plotted against Vp. The Vp of Banzatti et al. (2019) have been corrected by the slightly different adopted RV with respect to our work. The two dashed lines at ± 1.5 km s−1 indicate the velocity uncertainty of our measurements. Symbols colours refer to the disk type (red: full disks, blue: disk with sub-structures, green: TD) while squares indicate sources with a HVC.

In the text

Fig. 7 Relation between the deconvolved HWHM normalised by the stellar mass and sin (idisk) for the NLVC. Black lines indicate line widths as a function of inclination at disk radii of 0.2, 0.5, 2 and 5 au, in the assumption of Keplerian motion. In the upper panel, blue and cyan symbols distinguish sources where the Vp of the line is consistent with the origin in the disk or consistent with a wind, respectively. In the bottom panel, the distribution with respect to the disk type is shown, with colours code as in Fig. 4.

In the text

	Fig. 8 Comparison between the [O I] 630 nm line profiles having Vp close to zero, and a model of emission from a Keplerian disk. In the observed profiles the NLVC (or the NLVC+BLVC in the cases of Type 2 profiles) has been isolated by subtracting the Gaussian fit derived for the other components at different velocities. The fitted parameters of the model are indicated in each panel.
In the text

	Fig. 9 Keplerian radius of the [O I] 630 nm NLVC emission region for the TD of our sample is plotted against the size of the dusty disk inner cavity. For sources where multiple determinations of the inner cavity are given in literature, they are both displayed and connected with a line.
In the text

Fig. 10 Keplerian radius of the [O I] 630 nm NLVC emission region is plotted against the estimated and/or measured radius of the dusty disk inner rim, following different approaches. Color circles refer to the dust sublimation radius (colour code as in left upper panel of Fig. 4), squares refer to the magnetospheric radius, while asterisks are the dust inner rim radius measured through near-IR interferometric observations (see text for more details).

In the text

	Fig. 11 Continuum-subtracted spectra of the [O I] and [S II] lines in DO Tau and UX Tau. The NLVC is well detected in the [O I] lines but not in the [S II] lines, with a [O I] 630 nm/[S II] 406 nm ratio >10.
In the text

Fig. 12 Slab model prediction for the [S II] 406 nm/[O I] 630 nm ratio for collisions with electrons. Solid lines refer to gas temperatures of 6000, 8000, and 10 000 K. The cyan area indicates the location of the range of values observed in the LVC (both NL and BL) for the sub-sample of sources where the [S II] 406 nm line has been detected. The grey area shows the location of the upper limits measured in all the other sources. Horizontal lines indicate the more stringent upper limits.

In the text

Fig. 13 Slab model prediction for the [O I] 557 nm/[O I] 630 nm ratio for collisions with electrons. Solid lines refer to gas temperatures of 5, 8 and 10 × 103 K. The cyan area indicates the location of the range of values observed in the NLVC for the sub-sample of sources where the [O I] 557 nm line has been detected in this component. The grey area indicates instead the location of the detected ratios in the HVC.

In the text

Fig. 14 The [O I] 557 nm/[O I] 630 nm ratio observed in the NLVC with no velocity-shifted peak is plotted against the estimated [O I] 630 nm emission radius assuming Keplerian rotation. Green symbols indicate TD sources. The dashed line represents the linear fit through the points, whose analytical expression is indicated in the figure. The Pearson correlation coefficient is also displayed in the bottom right of the figure.

In the text

Fig. 15 Ionization equilibrium model predictions for the [N II] 658 nm/[O I] 630 nm ratio for collisions with electrons, plotted as a function of the xe = ne/nH ratio. Different lines refer to different values of electron temperature and density as indicated in the legend. In the upper panel, horizontal lines indicate some of the more stringent upper limits of the [N II] 658 nm/[O I] 630 nm ratio in the NLVC, while in the lower panel the coloured region indicates the range of the [N II] 658 nm/[O I] 630 nm ratio observed in the HVC.

In the text

Fig. 16 The ratio Ṁloss /Ṁacc, plotted as a function of the source Ṁacc. Upper panel: Ṁloss for the narrow low velocity component estimated as in Sect. 9.1 and assuming a fractional ionisation xe = 0.1. Lower panel: One sided Ṁloss for the HVC. In cases where the HVC can be decomposed into two additional components at medium and high-velocity, the mass-loss has been separately estimated, with the medium velocity indicated with blue symbols.

In the text

	Fig. 17 Map of the electron density of a photo-evaporative wind model from W20, with overlaid contours of the 80% emission regions of [O I] 630 nm (green) and [O I] 557 nm (magenta).
In the text

	Fig. A.1 Continuum-subtracted line profiles. Fluxes are in units of 10−15 erg s−1 cm−2 Å−1. The fit to the profile is shown in blue and it is obtained by adding the single or multiple Gaussians indicated in red. The two horizontal red dashed lines indicate ± 1 σ.
In the text