© ESO, 2017

1. Introduction

The circumstellar disk around young stars is the birthplace of planets, and the architecture of the forming planetary system depends on the evolution with time of the surface density of gas and dust in such disks (e.g., Thommes et al. 2008; Mordasini et al. 2012). Different processes are at play during the disk evolutionary phases, such as accretion of matter through the disk and on the central star (e.g., Hartmann et al. 2016), removal of material through winds driven by the high-energy radiation from the central star, or the result of magnetic torques on disk material (e.g., Alexander et al. 2014; Gorti et al. 2016; Armitage et al. 2013; Bai 2016), or even loss of material due to external disturbances, such as binarity and encounters (e.g., Clarke & Pringle 1993; Pfalzner et al. 2005) or external photoevaporation (e.g., Clarke 2007; Anderson et al. 2013; Facchini et al. 2016). Each of these processes modifies the distribution of material in the disk and is thus relevant for the planet formation process. Models aiming at explaining the observed properties of exoplanets or our own solar system need constraints on the contribution of each of these disk evolution processes (e.g., Adams 2010; Bitsch et al. 2015; Pfalzner et al. 2015). These constraints are obtained by studying properties of young stars and their disks at different ages, evolutionary stages, and with a wide span of stellar properties. Here we focus on measuring the rate at which material is accreted onto the central star, which is crucial information for all the aforementioned processes, as a function of stellar mass and luminosity for a complete sample of young stars that are surrounded by disks in the nearby (d = 160 pc¹) ~2–3 Myr old Chamaeleon I star-forming region (Luhman 2008).

Measurements of the mass accretion rate on the central star (Ṁ_acc) are taken from the excess emission in spectra of young stars that is due to the shock of gas infalling from the disk onto the central star along the stellar magnetic field lines (e.g., Calvet & Gullbring 1998). This excess emission is especially strong in the ultraviolet (UV) spectral range, in particular in the Balmer continuum, as well as in the optical (e.g., Fischer et al. 2011). Modern instruments mounted on 8 m class telescopes, such as the X-shooter spectrograph (Vernet et al. 2011) at the ESO Very Large Telescope (VLT), allow us to access the Balmer continuum region at λ ≲ 346 nm with high sensitivity while simultaneously obtaining medium-resolution flux-calibrated spectra that cover up to λ ~ 2.5μm. This has allowed researchers to derive Ṁ_acc for large samples of objects in different star-forming regions by modeling the Balmer continuum excess in the spectra of young stars (e.g., Taurus, Herczeg & Hillenbrand 2008; Lupus, Alcalá et al. 2014; σ-Orionis, Rigliaco et al. 2012). None of these studies has been performed on a statistically complete sample of disk-bearing young stars. A complete sample is instead crucial to firmly constrain the effects of different mechanisms on the disk evolution.

In a previous work including ~40% of the young stars that are surrounded by disks in the Chamaeleon I region (Manara et al. 2016a, hereafter MFH16) we have derived stellar and accretion properties using VLT/X-shooter spectra. In this sample Ṁ_acc scales with stellar mass (M_⋆) as a power law with exponent ~2, in agreement with previous results in other star-forming regions (e.g., Muzerolle et al. 2003; Mohanty et al. 2005; Natta et al. 2006; Herczeg & Hillenbrand 2008). Moreover, the young stars in the sample of MFH16 whose disk present a dust-free region in their inner part, the so-called transition disks (TDs), have similar Ṁ_acc as full disks with the same M_⋆. This is consistent with previous results (e.g., Manara et al. 2014; Keane et al. 2014). The spread of values of Ṁ_acc at any given M_⋆ was instead found to be smaller than those found in ρ-Ophiuchus (Natta et al. 2006) or in Taurus (Herczeg & Hillenbrand 2008), but larger than the sample comprising ~50% of the stars with disks in Lupus (Alcalá et al. 2014). However, complete samples, as analyzed here, are needed to constrain the real extent of the spread of Ṁ_acc and the dependence of Ṁ_acc on M_⋆, to be then compared with expectations from different models of disk evolution.

In the following, we present our sample and the data acquisition and reduction in Sect. 2, then we describe the stellar and accretion parameters derived for the sample in Sect. 3. In Sect. 4 we study the dependence between the accretion parameters and the stellar mass and luminosity, and we discuss the implications of our findings in Sect. 5. Finally, we present our conclusions in Sect. 6.

2. Data collection

2.1. Sample selection

The sample of objects discussed here is based on the selection made for the accompanying ALMA survey of the disk population of the Chamaeleon I region (Pascucci et al. 2016). The ALMA sample includes all the objects displaying excess emission with respect to the photosphere in more than one infrared wavelength ranging from the near-infrared, that is, 2MASS data, to mid-infrared, which means Spitzer or WISE data, and when available, far-infrared, that is, data obtained with Herschel (Luhman et al. 2008; Szűcs et al. 2010; Olofsson et al. 2013). This selection excludes all the objects classified as Class III, or disk-less, by Luhman et al. (2008). Objects still surrounded by an optically thick envelope, the Class 0 and Class I targets, were also excluded from the sample. These selection criteria result in a sample of 93 disk-bearing objects. According to previous spectral type classifications (Luhman et al. 2008; Luhman & Muench 2008; Luhman 2007), this sample is complete down to M6, which roughly corresponds to M_⋆~ 0.1 M_⊙, and it includes three targets with later spectral type, as discussed by Pascucci et al. (2016). This sample includes some binary stars, as reported in Pascucci et al. (2016).

The spectroscopic survey presented here similarly targets all the known young stellar objects in the region harboring a disk. Of the ALMA sample, only 5 targets were never observed with X-shooter²: 4³ because they are too faint to obtain a high enough signal-to-noise ratio in the optical spectrum, and one (CHXR30B) as it was not in the slit when observing CHXR30A. Two of these targets have previously derived spectral types M6 and M8, which means that they are in the range where the sample is incomplete. All other targets were observed with X-shooter. In particular, 43 ALMA targets observed in Pr.Id. 095.C-0378 (PI Testi) and 8 in Pr.Id. 090.C-0253 (PI Antoniucci) are presented here for the first time. Of these, 1 target is a newly discovered binary, as discussed in Sect. A.1. The remaining targets in the ALMA sample were observed with X-shooter in the past. In particular, 35 stars were observed in Pr. Id. 084.C-1095 (PI Herczeg) and analyzed by Manara et al. (2014, 2016a). Three of the ALMA targets studied by MFH16 are in binary systems that were resolved and studied separately. The latter sample contained mainly solar-mass objects, while the more recent spectra analyzed here are mostly focused on the lower-mass objects and selected to complete the initial sample. Two other targets were observed in Pr.Id. 085.C-0238 and 089.C-0143 (PI Alcalá). One of them is a well-known edge-on target, ESO-Hα 574, which was analyzed by Bacciotti et al. (2011) and is not discussed here, since the stellar and accretion parameters of edge-on targets are too uncertain to be included in the analysis. The other, ISO-ChaI 217, is a a well-known brown dwarf with a jet (Whelan et al. 2014) and is included in our sample. With respect to the ALMA sample, our sample includes 2 additional targets. One is Sz18, which was not included in the ALMA sample because it was originally classified as a Class III target by Luhman et al. (2008), but was reclassified as a transition disk by Kim et al. (2009) and thus studied by Manara et al. (2014). The other addition, Cha-Hα1, is a brown dwarf and in the spectral type range where the ALMA survey was incomplete, and it was studied by MFH16.

With respect to the ALMA sample, which is complete down to M6, we therefore lack only 3 targets with spectral type earlier than M6. We accordingly consider our total sample of 94 objects, which includes resolved components of binary systems, to comprise 97% of the disk-bearing young stellar objects with spectral type earlier than M6 in the Chamaeleon I region.

2.2. New observations

The new data presented in this work were obtained in two different observing runs using the VLT/X-shooter spectrograph (Vernet et al. 2011). The targets of Pr.Id. 095.C-0378 (PI Testi) were observed using both narrow (see Table D.1) and large (5.0′′ × 11′′) slits. The narrow-slit observations were carried out by nodding the target along the slits on two slit positions A and B, with one exposure per position. The nodding cycle was ABBA, with a total of four exposures. The exposure times and slit widths were chosen to obtain enough signal based on the brightness of the targets, and are reported in Table D.1. These exposures with narrow slits were obtained to achieve the highest spectral resolution, typically R ~4000–10 000 in the UVB arm (λλ ~300–550 nm), and R ~ 6700–18 000 in the VIS arm (λλ ~ 550–1000 nm). Simultaneous observations in the near-infrared (NIR) arm (λλ ~ 1000–2500 nm) will be described in a following paper. The large-slit observations were carried out in stare mode on the target immediately after the narrow-slit exposures, and they lasted for ~10% of the total exposure time of the observations taken with the narrow slit. These large-slit observations led to spectra with a much lower resolution, but without slit losses, and are thus crucial for a correct flux calibration of the spectra obtained with the narrow slit.

The eight targets observed during Pr.Id. 090.C-0253 (PI Antoniucci) were similarly observed by nodding the telescope and using narrow slits, with widths reported in Table D.1. All these spectra have a resolution of R ~ 10 000 in the UVB arm and R ~ 18 000 in the VIS arm. However, no exposures with the large slit were obtained in this observing run. We therefore used non-simultaneous photometric data to calibrate these spectra, as discussed in the next section.

2.3. Data reduction

Data reduction was performed with the ESO X-shooter pipeline (Modigliani et al. 2010) version v.2.5.2 run through the Reflex workflow (Freudling et al. 2013). The pipeline performs the usual reduction scheme, including bias subtraction, flat-fielding, wavelength calibration, flexure, and atmospheric dispersion correction, background removal (in stare mode) or combination of spectra obtained in a nodding cycle, and spectrum extraction. The latter is performed by the pipeline on the 1D background-subtracted spectra using a large extraction window. In order to maximize the signal in the UVB spectra, we manually extracted the spectrum using the apall task in IRAF⁴. The flux calibration of the spectra was also performed in the pipeline by deriving a spectral response function using a flux standard star observed during the same night. Typically, this procedure leads to a flux calibration accuracy of ~2% (e.g., Alcalá et al. 2014; Manara et al. 2016a). However, the spectra obtained using the narrow slit suffer slit losses. To correct for this effect, we scaled the flux-calibrated spectra obtained with the narrow slit to the flux-calibrated spectra obtained with the large slit using the median of the ratio between the two spectra as a correction factor in each arm when the narrow slit was 0.9′′ wide, or more, and a wavelength-dependent linear fit of the ratio in each arm when the narrow slit was 0.5′′ or 0.4′′ wide. The latter is needed to account for wavelength-dependent differential slit losses that are due to the high airmass of the targets at the time of the observations. Typically, the correction factors are ~1.5–2, and they are ~3–5 in a few cases when observations were performed with the narrowest slits and with large seeing. The agreement between the flux calibrated spectra in the overlapping region of the VIS and UVB arm is excellent. This same procedure to correct for slit losses was not possible for the eight targets observed in Pr.Id. 090.C-0253, as these were observed only with the narrow slits. We thus used available photometry to correct these spectra (see Table A.1). We first matched the spectra from the three arms and then scaled all of them to minimize the scatter to the photometric points. The correction factors found with this method are ~3 − 5, compatible with those found when the large slits were also adopted. Finally, telluric correction was performed in the VIS arm using telluric standard stars observed close in time and airmass to the target. The spectra of these standard stars were continuum normalized and the photospheric absorption features were removed before correcting the telluric lines using the telluric task of IRAF.

The sample includes some unresolved close binaries, T5, CHXR30A, CHXR71, Cha-Hα2, T45, T46, CHXR 47, Hn13, CHXR79, and T43. Their spectra thus include both components. On the other hand, the two components of 2MASS J11175211−7629392 were separated, and we manually extracted the two spectra. Their analysis is discussed in Sect. A.1.

Three objects needed a particular procedure to reduce their spectra. CHXR71 was observed during a cloudy night, and the spectrum obtained with the large slit is fainter than the spectrum with the small slit. We thus corrected for slit losses using a correction factor of 3, which leads to a spectrum matching the available photometry. The ratio of the UVB spectra with the large and small slits for both T48 and T27 could not be fitted with a linear dependence to wavelength, and we thus used a second-order polynomial to fit this ratio and correct the small-slit spectra.

3. Data analysis

The spectra analyzed here contain multiple features typical of young stellar objects. We checked in particular whether they present the Hα line in emission and the Li absorption line at λ670.8 nm. These features are usually considered a confirmation of the nature of young stars for the target. All objects except for ISO-ChaI 79 present a clear Hα line in emission. This line is visible in the spectrum of this target only with some smoothing, as the spectrum is very noisy, probably because the disk is viewed edge-on, obscuring the star (see Sect. 3.2) The lithium absorption line is detected in 36 targets in the sample. These are the objects whose spectrum has signal-to-noise S/N ≳ 10 in the continuum at λ ~ 700 nm (see Table D.1). We thus confirm that all the targets are young stellar objects, since the spectra of 36 targets present the lithium line in absorption and the Hα line in emission, while the other targets show the Hα line in emission and the lithium line is not detected due to the low signal-to-noise ratio of the spectra.

In the following, we describe the method used to derive stellar and accretion parameters for the targets.

3.1. Determining the stellar and accretion parameters

The spectral type of the target, its stellar luminosity (L_⋆), extinction (A_V), and accretion luminosity (L_acc) are obtained by modeling the spectrum of the targets from λ ~ 330 nm to λ~715 nm following the automatic procedure described by Manara et al. (2013b). This same procedure was used by MFH16, Manara et al. (2014), and Alcalá et al. (2014, 2017), among others. Briefly, the observed spectrum is modeled as a sum of a photospheric template and of a slab model to reproduce the excess emission due to accretion (as in, e.g., Valenti et al. 1993; Herczeg & Hillenbrand 2008), which are reddened assuming a typical extinction law with R_V = 3.1 (Cardelli et al. 1989). The best match is found by minimizing a distribution over the spectral type of the template, the parameters of the slab model, and the values of A_V. With respect to MFH16, here we include additional photospheric templates of non-accreting young stars in addition to those by Manara et al. (2013a, 2014). They have spectral types between G5 and K6, and between M6.5 and M8 (Manara et al. 2017). The spectral type of the best-matching photospheric template is adopted for the target, and the flux ratio between the template and the target spectrum corrected for extinction using the best-fit A_V is used to derive L_⋆ (Manara et al. 2013b). These values are reported in Table 1. Typical uncertainties on A_V are <0.5 mag, on the spectral type the uncertainties are ±0.5 subclass for M-type young stars and ±1 subclass for earlier-type stars, and on logL_⋆ the uncertainties are 0.2 dex. The final T_eff for the targets is obtained from the spectral type using the relation by Luhman et al. (2003) for M-type stars and the relation by Kenyon & Hartmann (1995) for earlier-type stars. Finally, L_acc is derived by directly integrating the flux of the best-fit slab model, with typical uncertainties on L_acc of ~0.25 dex.

These results from the automatic procedure are also validated against independent methods. The spectral type is compared with the type derived using different spectral indices by Riddick et al. (2007), Jeffries et al. (2007), and Herczeg & Hillenbrand (2014). Of the indices by Riddick et al. (2007), we selected the same as have been used by Manara et al. (2013a), while we used the indices G-band, R5150, TiO-7140, TiO-7700, and TiO-8465 from Herczeg & Hillenbrand (2014) depending on the spectral type of the target. The spectral types derived with these indices are reported in Table 3, where we list the mean value when multiple indices from the same author are adopted. In general, these values agree within ±1 subclass with our estimates, which is within their validity ranges. The estimated L_⋆ agree within ~0.2 dex with those obtained using the bolometric correction by Herczeg & Hillenbrand (2014). Finally, L_acc from the method just described, which is directly obtained from the continuum UV-excess emission, is always in agreement with the one derived by converting the luminosity of different emission lines into L_acc using the relation by Alcalá et al. (2014) within the uncertainties.

This method is only applicable when the spectrum has signal down to λ ~ 330 nm, and it cannot be applied to 19 late-M brown dwarfs with non-negligible extinction. For these young stars, marked in Table 1 with “m”, the stellar parameters are derived using the spectral indices by Riddick et al. (2007), Jeffries et al. (2007), and Herczeg & Hillenbrand (2014) and by comparing the observed spectrum with the spectra of photospheric templates reddened by increasing values of A_V in steps of 0.1 mag, until a best match is found visually. The emission line fluxes are then used to derive L_acc using the relations between the line luminosity and L_acc by Alcalá et al. (2014).

In the following sections, we discuss the results for the whole Chamaeleon I sample, including the targets studied by MFH16, Manara et al. (2014), and Whelan et al. (2014). For these objects we also report the spectral type, T_eff, A_V, L_⋆, and L_acc in Table 2. These values have been derived with the automatic method described in this section, or with the other method we described in the case of ISO-ChaI 217 (Whelan et al. 2014).

3.2. Determining stellar mass and mass accretion rate

In order to derive M_⋆ for the observed targets, pre-main sequence evolutionary models are necessary. The choice of model has an effect on the estimate of M_⋆ (see, e.g., Stassun et al. 2014; Rizzuto et al. 2016). We show in Fig. 1 the Hertzsprung-Russel diagram (HRD) for the whole Chamaeleon I sample and two different sets of isochrones from the models by Siess et al. (2000) and Baraffe et al. (2015). These models lead to similar values of M_⋆, while the isochronal ages are different. We decide to use the more recent models by Baraffe et al. (2015) to derive M_⋆ by interpolating the models to the position of the targets on the HRD. These models are an improvement from past models in their ability to reproduce the observed distribution of low-mass targets between 3200–5000 K on the HRD for different young associations. However, these models are limited to M_⋆ < 1.4 M_⊙ and can therefore not be used for five targets whose position on the HRD indicate higher M_⋆. For these targets we use M_⋆ derived using the evolutionary models by Siess et al. (2000), since they lead to similar M_⋆ as in Baraffe et al. (2015) in the overlapping mass range, and they extend to higher M_⋆. We do not consider the individual isochronal ages of the targets in this work, therefore the large differences in the age estimates using one or the other model do not affect our analysis. The values of M_⋆ derived with these models are reported in Tables 1 and 2 and have typical uncertainties of ~0.1 dex. The method used to derive M_⋆ as well as the evolutionary model adopted for M_⋆ > 1.4 M_⊙ differ from the choice of the companion paper that analyzes the ALMA data for Chamaeleon I objects (Pascucci et al. 2016). There we used a Bayesian approach to derive the stellar mass and its uncertainty for each object using the Baraffe et al. (2015) models complemented with those from Feiden (2016). However, the differences in M_⋆ are very small, typically Δ(log M_⋆) < 0.05 dex, with only seven objects with differences up to 0.07 dex.

The targets in the Chamaeleon I region studied here are distributed in the HRD along the 3 Myr isochrone of the Baraffe et al. (2015) models, but with a wide spread of up to ~1 dex in logL_⋆ at all T_eff. This spread is wider than the typical uncertainties on L_⋆ and would in turn imply a wide spread in ages. However, ages derived from evolutionary models are uncertain, therefore it is not straightforward to derive an age spread from the HRD (e.g., Soderblom et al. 2014).

The positions on the HRD of three objects are clearly below the main locus of the Chamaeleon I sample and below the 30 Myr isochrone of the Baraffe et al. (2015) model. They are marked in the plots with open squares and in Table 1 as “UL”, which stands for underluminous. Two of these targets are known to present an edge-on disk (ESO Hα 569, Robberto et al. 2012; 2MASS J10533978−7712338, Luhman 2007). We also consider the third object, ISO-ChaI 79, to be seen behind an edge-on disk, although future studies are needed to confirm this finding. These objects will not be considered in the following analysis, as their stellar and accretion parameters are more uncertain because of the high gray extinction on the line of sight caused by the edge-on disk. One additional object, ESO-Hα 574, also has an edge-on disk, but is not included in our sample or discussed further. These four objects, together with the target 2MASS J11082570-7716396, which was not observed with X-shooter, are the same five underluminous objects in the ALMA sample as were reported by Pascucci et al. (2016). The fact that stars seen through an edge-on disk have L_⋆ corresponding to a position on the HRD below the 30 Myr isochrone is consistent with what is found in the Lupus star-forming region by Alcalá et al. (2014, 2017).

One target, Hn13, has L_⋆ significantly higher than the luminosity corresponding to the youngest isochrone for its temperature according to the models by Baraffe et al. (2015), and we adopt the value of M_⋆ obtained using the models by Siess et al. (2000) for this target. This choice is possible since M_⋆ derived with the two models is very similar in this stellar mass range.

The values of L_acc measured by fitting the observed spectrum are finally combined with M_⋆ and the stellar radius (R_⋆), derived from L_⋆ and T_eff, to derive Ṁ_acc with the usual relation Ṁ_acc = L_accR_⋆/ (0.8 · GM_⋆) (e.g., Hartmann et al. 1998). These values of Ṁ_acc are reported both for the sample analyzed here and for the remaining targets in Chamaeleon I in Tables 1 and 2, and their typical uncertainties are ~0.35 dex.

3.3. Non-accreting targets

Young stars are known to have a very active chromosphere. The emission from this chromosphere contributes in the spectra of young stars to both the continuum emission (e.g., Ingleby et al. 2011) and the line emission (e.g., Manara et al. 2013a). The choice of using spectra of non-accreting young stars as photospheric templates for the analysis is also driven by the fact that this allows us to consider the approximate chromospheric contribution to the observed continuum emission of the target star. When objects are strongly accreting, the chromospheric emission represents a negligible contribution to the total excess emission. On the other hand, there are objects for which the excess emission is small and where most, if not all, of this emission may come from a pure chromosphere. Manara et al. (2013a) have shown that the chromospheric emission measured from the luminosity of emission lines in the spectra of young stars can be converted into a typical bias, or noise, on the accretion luminosity (L_acc,noise). The typical intensity of the chromospheric emission, when converted into L_acc,noise and measured as the ratio with L_⋆, decreases for stars with later spectral type (Manara et al. 2013a), while it is constant for late-G and K-type young stars (Manara et al. 2017).

In order to constrain the importance of chromospheric emission with respect to the total measured excess emission, we compare the values of log(L_acc/L_⋆) measured for our targets with the expected chrosmopheric emission at the same T_eff, log(L_acc,noise/L_⋆) in Fig. 2. The values of log(L_acc,noise/L_⋆) expected at each T_eff are shown with a dashed line for M-type stars (Manara et al. 2013a) and a dotted line for earlier-type stars (Manara et al. 2017). The vast majority of the targets have L_acc/L_⋆ much higher than the typical L_acc,noise/L_⋆. However, 13 targets have L_acc comparable or below the expected emission of an active chromosphere. This number includes the two components of the binary 2MASS J11175211-7629392, which are then excluded in the following analysis. It is possible that the excess emission, or the emission line luminosity, measured for these targets is mainly due to chromospheric emission. However, it is also possible that the chromosphere of these targets is less active than typical young stars, and thus this emission really originates from the accretion of matter onto the central star. Since it is not possible to distinguish between these possibilities, we refer to these targets as “dubious accretors” or low-accretors in the text, and we mark them in Tables 1 and 2, and with open symbols in the following plots. The measured values of L_acc for the dubious accretors are not upper limits on their accretion rate, but a measurement that is most likely strongly contaminated by another process. We also show in Appendix B.2 that considering their measured L_acc as detection or upper limit does not affect our results. We stress that because of our selection criteria, all these objects show evidence of infrared excess due to the presence of warm dust around them, and that millimeter continuum emission from a disk has been detected for five of them⁵ by Pascucci et al. (2016). Interestingly, objects with an ALMA counterpart are mostly located on the L_acc/L_⋆ vs. T_eff plane very close to the threshold for being considered accretors, while most of the objects not detected with ALMA are found even at much lower values of L_acc/L_⋆ than this threshold. We note that objects of this type, that is, those showing infrared excess and/or millimeter continuum emission from a disk while not showing a signature of accretion, are well known in the literature (e.g., Mohanty et al. 2005; McCabe et al. 2006; Fedele et al. 2010; Wahhaj et al. 2010; Furlan et al. 2011; Hernández et al. 2014). It is as yet unclear, however, whether they represent a typical evolutionary stage of disk evolution or are peculiar. However, this discussion is beyond the scope of this paper.

MFH16 used a different criterion to asses the status of accretor for the targets, namely the width at 10% of the peak and the equivalent width of the Hα line, which is a criterion usually adopted in the literature (e.g., White & Basri 2003; Muzerolle et al. 2003; Mohanty et al. 2005). This leads to some discrepancies with the results just discussed. In particular, only T54-A and ISO-ChaI 52 are considered dubious accretors here and also by MFH16. Three other objects that seem to be non-accreting according to the width of their Hα line, that is, T45a, T4, and T33-A are located very close to the L_acc,noise value for their T_eff. Finally, Hn17 and Hn18, whose Hα line width indicates that their are not accreting, are instead found to be well above the typical L_acc,noise for their T_eff, and they are thus accreting. The method of comparing the measured L_acc with L_acc,noise is less affected by uncertainties on the real peak value of the Hα line, which is needed to derive the width at 10% of the peak, and by the effects of a strong photospheric absorption line, which would modify the measured equivalent width of the emission line. Moreover, the width of the Hα emission line is subject to the rotation of the star itself, and non-accreting objects could present very wide Hα emission lines if they rotate very fast. This is the case, for example, of Sz121 (Manara et al. 2013a). In the following we consider as dubious accretors only those classified comparing L_acc with the typical L_acc,noise.

4. Results

This section is focused on the relationships between the accretion and stellar properties for the whole sample of Chamaeleon I young stars. Different statistical tests are run to determine the shape of the L_acc–L_⋆ and Ṁ_acc–M_⋆ relations. Additional tests are reported in Appendix B. In the following subsections, objects whose position on the HRD is well below the 30 Myr isochrone (see Sect. 3.2) are excluded from the analysis and are not shown on the plots.

4.1. Accretion luminosity and stellar luminosity dependence

The logarithmic dependence of L_acc on L_⋆ is shown in Fig. 3. In general, L_acc increases with L_⋆. However, objects with L_⋆≳ 0.1 L_⊙ reach higher ratios of the accretion to stellar luminosity (L_acc/L_⋆), up to ~1, than lower luminosity stars. The lines of constant L_acc/L_⋆ ratio are overplotted in Fig. 3 for increasing values from 0.01 to 1. No targets with L_⋆ ≲ 0.1 L_⊙ have ratios L_acc/L_⋆ ≳ 0.1. In contrast, ~50% of the objects with L_⋆ ≳ 0.1 L_⊙ have an accretion luminosity higher than 10% of the stellar luminosity. This suggests a change in the accretion properties for objects with a stellar luminosity above and below L_⋆~ 0.1 L_⊙.

We first explore our data using non-parametric tools, which have the advantage that they do not rely on prior assumptions on the real model describing the data. Both a Nadaraya-Watson (spatial averaging) and a local-polynomial fit carried out using the PyQt-fit module on the L_acc vs. L_⋆ data show two different regimes in the L_acc–L_⋆ relation: a steep increase for log(L_⋆/L_⊙) ≲ −0.6, and a flatter relation at higher L_⋆, with a slope ~1. A similar result is derived when computing the median values of log(L_acc/L_⊙) as a function of log(L_⋆/L_⊙) by dividing the sample into bins containing the same number of objects. The results from these tests, discussed in Appendix B, suggest that the dependence of L_acc to L_⋆ is possibly more complex than a single power-law distribution, and that two power laws, a segmented line on the logarithmic plane, are a possible representation of the observations. Therefore we consider in the following both hypotheses, and we try to quantify which model is a better fit using different statistical tests.

The fit is performed using two different Python modules, linmix by Kelly (2007) and scipy.optimize.curve_fit. The former allows one to consider the uncertainties on the measurements on both axes, fits the data with a fully Bayesian analysis, but can only be used to fit linear relations. The latter does not include a treatment of the measurement errors, but is able to fit more complex relations between the quantities. It uses the nonlinear least-squares Levenberg-Marquardt algorithm to fit a user-defined function to the data. We show in Appendix B that the two tools lead to results that are compatible with each other since the scatter in the values dominates the outcome of the fit more than the measurement errors. Moreover, we also show in Appendix B that treating the accretion rate measured for the dubious accretors as detections is equivalent to considering these values as upper limits in the fit.

The two models considered here have a different number of free parameters. The single power law is described by two parameters, a slope and intercept in the logarithmic plane. The double power law is a segmented line in the logarithmic plane, which is described by the following function: (1)where θ_i (i = 0,1,2) and x_c are free parameters of the fit, y is log (L_acc/L_⊙) and x is log (L_⋆/L_⊙). Hence, this model has four free parameters. To properly quantify whether the single or double power-law is to be preferred, we use three different information criteria to estimate the goodness of fit, namely R², the Akaike information criterion (AIC), and the Bayesian information criterion (BIC), as discussed by Feigelson & Babu (2012), for instance. While the former does not consider the number of free parameters in the model, both AIC and BIC take this aspect into account. In general, the model that better reproduces the data should maximize R², and it should minimize both AIC and BIC⁶. As discussed for example by Raftery (1995) and Riviere-Marichalar et al. (2016), a difference in the values of BIC between 2 and 6 shows that the model with lower BIC is more plausible, while a difference of 10 or more excludes the model with the higher value of BIC with a high probability. Similarly, Murtaugh (2014) discussed that differences in AIC values of 14 or more generally firmly exclude the model with the higher AIC value, while differences of ~4–8 correspond to a p-value of 0.05, thus to a distinction between the models at a lower significance.

We first fit the data with scipy.optimize.curve_fit. The slope of the best fit with a single line is 1.75, and the intercept is −0.96. The fit with a segmented line leads to a lower AIC (208 vs. 212), a higher R² (0.78 vs. 0.70), and the same value of BIC (217). According to the AIC criterion, the segmented line is preferred, and R² also points toward a better fit with the segmented line. However, the latter does not consider the number of degrees of freedom in the model, and the BIC criterion does not prefer any of the two models. Therefore, the statistical tests on this relation are not conclusive. The best fit with the segmented line is shown in Fig. 3, and more discussion of the results with this method are reported in Appendix B. The best-fit parameters for the segmented line model are θ₀ = − 0.55, θ₁ = 2.08, θ₂ = − 1.45, and x_c = − 0.34, which imply (2)The linmix tool is also used to test the two models. The single-line fit is shown in Fig. 3 and is (3)with R² = 0.69, BIC = 219, and AIC = 214, and a 1σ dispersion of 0.67 ± 0.08 around the best fit. We then use the value of log(L_⋆/L_⊙) at which the scipy.optimize.curve_fit module finds a change in the slope of the segmented line fit as initial parameter to divide the sample into two subsamples. These two subsamples are then fitted separately with the linmix module (see Appendix B). The two slopes for high and low L_⋆ are compatible with those obtained with scipy.optimize.curve_fit, namely 2.5 ± 0.2 below log(L_⋆/L_⊙) = − 0.34, and 1.1 ± 1.3 above. However, the value of the slope for higher L_⋆ is not constrained by the fit, since the uncertainty is too large. This is probably due to the large scatter of ~0.9 dex in this region of the plot, which then results in a low correlation (r = 0.3 ± 0.3), in contrast to the very high correlation found below the break value (r = 0.91 ± 0.04). This fit with a segmented line leads to a slightly higher R² = 0.71, to a higher value of BIC = 222, and to a lower AIC = 212. Therefore, in this case the statistical tests are not conclusive either in the choice of best model to describe the data, since the better difference in AIC values is not significant, and the BIC statistics would prefer the linear model with a low confidence.

We can then conclude that the fit with a segmented line, with an exact value of the slope in the range of value log(L_acc/L_⋆) ≥ −0.34, which is very uncertain, is slightly preferred by some statistical tests, but not by others, and never with results that clearly exclude the single power-law model. Therefore, both models are plausible.

4.2. Mass accretion rate and stellar mass dependence

While the relation between L_acc and L_⋆ is derived independently of evolutionary models, the physical quantity that characterize the accretion process is instead the amount of material accreted onto the central star per unit time, Ṁ_acc. It is thus important to study how this varies with M_⋆ to constrain models of disk evolution. We show the logarithmic relation between these two quantities in Fig. 4. Stars with higher M_⋆ have generally higher Ṁ_acc. Similarly to the findings by Manara et al. (2014, 2016a), as well as others (e.g., Sicilia-Aguilar et al. 2009; Fang et al. 2013; Keane et al. 2014; Alcalá et al. 2017), transition disks are well mixed with objects with a full disk. Non-parametric analyses (see Appendix B) suggest that a break in the distribution of Ṁ_acc for increasing M_⋆ is present, with a steeper slope for log(M_⋆/M_⊙) ≲ −0.5 than above this value, where Ṁ_acc is almost constant. Here we proceed similarly as in Sect. 4.1 to quantify whether a single power-law or a double power-law model better describes the data.

The fit with scipy.optimize.curve_fit is performed first, using both a single linear relation in the logarithmic plane and a segmented line equivalent to the line used for the L_acc–L_⋆ relation (Eq. (1)). We find that the segmented line leads to a better fit to the data with respect to a single line, with a higher value of R² (0.6 vs. 0.5), and both a lower AIC (223 vs. 236) and BIC (233 vs. 241). The single-line fit has a slope of 2.2 and an intercept of −8.2. The best-fit values for the segmented line are θ₀ = − 6.45, θ₁ = 4.31, θ₂ = − 3.75, and x_c = − 0.53, and this best fit is shown in Fig. 4. This implies two lines with equation: (4)We then fit the data using the linmix tool, starting with a single power-law fit. The best fit in this case is shown in Fig. 4 and is (5)and the 1σ dispersion around the best fit is 0.8 ± 0.1. This dispersion is consistent with the one reported by MFH16 for a subsample of Chamaeleon I targets. Moreover, the slope of this relation is slightly steeper than those reported in the past (e.g., Mohanty et al. 2005; Natta et al. 2006; Herczeg & Hillenbrand 2008; Alcalá et al. 2014; Manara et al. 2016a), but still compatible within the uncertainty with the typical values of 1.8–2 reported in the past. This small difference in the slope is only partially due to the evolutionary model chosen here. For comparison, using the models by Baraffe et al. (1998) on this same dataset would result in a best-fit slope of 2.2 ± 0.2. This fit results in values of R² = 0.5, BIC = 241, and AIC = 236. Finally, we test the segmented line fit also using the linmix tool on the two subsamples, considering either the objects with log(Ṁ_acc/M_⋆) ≥ −0.53, or lower than this value, where this is selected based on the results of the fit with scipy.optimize.curve_fit. The result is as follows. The correlation between the two quantities is strong in the sample of objects with M_⋆ lower than the value of the break (r = 0.9 ± 0.1), with a slope of 5 ± 1. Conversely, the spread is wide at higher values of M_⋆, and the fit results in a slope of 0.7 ± 0.8, with r = 0.2 ± 0.2. This fit with a segmented line is a better representation of the data according to all the three information criteria we use, as it leads to a higher R² = 0.6, and both a lower BIC = 238 and AIC = 228 than the single-line fit. However, the differences in the values of these information criteria do not rule out any of the two models.

Even if the statistical improvement of the segmented line is small, the statistical information criteria seem to prefer this model, which leads to a steeper slope at lower M_⋆ and a lack of correlation at higher M_⋆. However, none of the two description of the data can be ruled out.

There is a discrepancy in the results of the statistical tests between the best model to describe the Ṁ_acc–M_⋆ and the L_acc–L_⋆ relations. The model with a double power-law is slightly preferred by all the information criteria in the case of the Ṁ_acc–M_⋆ relation, while this is not necessarily the case for the L_acc–L_⋆ relation. This discrepancy could be attributed to the fact that both M_⋆ and Ṁ_acc are derived from the measured quantities L_⋆ and L_acc using evolutionary models, in this case, from Baraffe et al. (2015) and Siess et al. (2000). In particular, it is possible that the conversion from L_⋆ and T_eff into M_⋆ enhances the break in the Ṁ_acc–M_⋆ relation. Several studies have shown the limits of the various evolutionary models in estimating M_⋆ (e.g., Stassun et al. 2014; Herczeg & Hillenbrand 2015; Rizzuto et al. 2016), and this issue may add a word of caution on the double power-law model, but also on the values of the single power-law fit. Furthermore, we checked that there is a small difference in the typical values of A_V in the lower M_⋆ and higher M_⋆ subsamples, with the latter having a higher median A_V by ~0.9 mag. This could explain the wider spread in value of Ṁ_acc at higher M_⋆. However, the distributions of A_V at low and higher M_⋆ are similar, and thus this difference in median values is probably not the origin of the possible bimodal distribution in the Ṁ_acc–M_⋆ plane. Finally, we do not expect strong effects on our results because of the incompleteness in the sample. We discussed in Sect. 2.1 that our sample comprises nearly all the stars with disks in the Chamaeleon I region down to M_⋆~ 0.1 M_⊙. Unless we systematically missed the strongly accreting brown-dwarfs, for example, when they are still embedded in the parental cloud and are thus classified as Class I or are undetected at infrared wavelengths, the sample incompleteness would not affect the result. This said, at the present time we cannot rule out either of the two models, that is, whether a single or a double power-law are the best description of the data, and the implications of both possibilities are discussed in the next section.

5. Discussion

We have shown in Sect. 4 that the sample of stars with a disk in Chamaeleon I shows an increase in L_acc with L_⋆, and likewise in Ṁ_acc with M_⋆. We have discussed that these relations can be modeled either with a single or double power-law. The latter leads to a steeper slope at lower L_⋆ or M_⋆ than for higher L_⋆ or M_⋆ objects. Our statistical tests cannot firmly rule out any of the two hypotheses, although the double power-law description is preferred in particular for the Ṁ_acc–M_⋆ relation. In the following, we thus discuss the implications of our results in light of both descriptions of the data. We then discuss the overall distribution of data on the Ṁ_acc–M_⋆ plane, and how this affects our understanding of disk evolution.

5.1. Single power-law describing the relation of accretion to stellar parameters

When fitting the logarithmic relations between accretion and stellar parameters with a single linear fit, we derive L_acc ∝ L_⋆^1.9±0.1 and Ṁ_acc ∝ M_⋆^2.3±0.3, which is consistent with previous studies (e.g., Muzerolle et al. 2003; Natta et al. 2006, 2014; Mohanty et al. 2005; Herczeg & Hillenbrand 2008; Rigliaco et al. 2011a; Alcalá et al. 2014, 2017; Venuti et al. 2014; Kalari et al. 2015; Manara et al. 2012, 2016a). These values are in broad agreement with expectations from different theoretical models.

The L_acc–L_⋆ relation was studied by Tilling et al. (2008). The authors showed that in the context of viscously evolving disks the slope of this relation can be related to the index of the exponential decay of Ṁ_acc with time in the self-similar late phase of disk evolution. This index depends on the scaling of the viscosity (ν) with the disk radius (R). A scaling in the form ν ∝ R (Hartmann et al. 1998) implies L_acc ∝ L_⋆^1.7. A steeper slope of the L_acc–L_⋆ relation would be reproduced by a steeper dependence of the viscosity on the disk radius. However, the authors argue that different assumptions on the scaling of the viscosity with radius similarly fill the whole distribution of the observed L_acc–L_⋆, thus it is not possible to derive constraints on the viscosity law from this distribution, unless there are regions of the allowed parameter space that are empty of observed points. We return to this point below.

The dependence of Ṁ_acc to M_⋆ with a power law with exponent ~2 has been discussed by several authors. Alexander & Armitage (2006) and Dullemond et al. (2006) argued that this dependence is a result of the imprint of initial conditions on the subsequent viscous evolution of disks. The former authors were able to reproduce the observed Ṁ_acc ∝ M_⋆² relation under the assumption that the viscous timescale (t_ν) scales with the inverse of M_⋆. This implies a longer viscous timescale for lower mass objects, thus, in a viscously evolving system, a longer disk lifetime. However, this could be in contrast with our results for a fit with a double power-law, as we discuss in the following. Dullemond et al. (2006) instead assumed that cores with significantly different masses rotate similarly to breakup rate ratios, and this resulted in a power-law dependence with exponent 1.8 ± 0.2 for the Ṁ_acc–M_⋆ relation, which is compatible with our result. As a consequence of viscous evolution, the authors suggest that Ṁ_acc and the disk mass (M_disk) must also be correlated with their assumptions on initial conditions. This correlation was recently observed in both the Lupus star-forming region (Manara et al. 2016b) and in the Chamaeleon I region (Mulders et al., in prep.).

Another theoretical argument proposed by Padoan et al. (2005) to explain the Ṁ_acc∝M_⋆² relation is that accretion onto the star is a consequence of Bondi-Hoyle accretion of the gas in the surrounding star-forming region onto the disk-star system. However, the actual rate of accretion onto the central star may be lower than the Bondi-Hoyle rate (e.g., Mohanty et al. 2005), and a dependence of Ṁ_acc on the properties of the gas surrounding the young stars has never been observed (e.g., Mohanty et al. 2005; Hartmann et al. 2006).

Photoevaporation by high-energy radiation from the central star can also explain the dependence of Ṁ_acc on M_⋆. In this context, our results are more in agreement with expectations from X-ray-driven photoevaporation (Ercolano et al. 2014, Ṁ_acc ∝ M_⋆^1.6−1.9) than UV-driven photoevaporation (Clarke & Pringle 2006, Ṁ_acc ∝ M_⋆^1.35).

Finally, the dependence of Ṁ_acc on M_⋆ can also be explained if the ionization of the disk, and thus the magnetorotational instability that generates the viscosity driving the accretion in the disk, is strongly dependent on M_⋆ (Mohanty et al. 2005). As discussed by Hartmann et al. (2006), in a disk model where only the surface layer is ionized, and hence accretion is driven in this layer, Ṁ_acc is independent on M_⋆. By including additional heating by irradiation from the central star in a disk model with ongoing layered accretion, the authors were able to predict a dependence Ṁ_acc ∝ M_⋆. They suggest that this could explain the properties of disks around solar-mass stars. They then suggest that disks around very low-mass stars may be magnetically active either because of a very small initial disk radius that is due to a strong dependence of the disk radius on stellar mass (), or because the viscosity parameter α is high. This scenario leads to a steep dependence of Ṁ_acc on M_⋆, which is compatible with observations. In this scenario, very low-mass stars and brown dwarfs will have a small and magnetically active disk that quickly evolves viscously and settles at a lower rate of Ṁ_acc. This predicts a faster evolution for low-mass than for solar mass objects. We note that recent observations of disks around brown dwarfs with ALMA have found that these disks are small in the young ρ-Ophiuchus region (e.g., Testi et al. 2016). As we discuss in the following, this could be in line with our results.

5.2. Implications of a bimodal distribution in the relation of accretion to stellar parameters

The possibility that the L_acc–L_⋆ and Ṁ_acc–M_⋆ relations are described by a double power-law with a steeper slope at M_⋆ ≲ 0.3 M_⊙ is considered here. This possibility has been little explored in the past and could have strong implications on our theoretical understanding of disk accretion and evolution. We have presented some possible caveats to this finding in Sect. 4.2, but based on our current knowledge, we cannot exclude that this bimodal distribution is indeed the correct representation of the data. We thus qualitatively discuss some possible explanation for this bimodality. We also note that a similar behavior is observed in the complete sample of young stars with disks in the Lupus star-forming region that has been studied with the same method by Alcalá et al. (2017), and is possibly present in the survey of NGC 2264 by Venuti et al. (2014). Both regions have a very similar age as Chamaeleon I (Alcalá et al. 2014; Venuti et al. 2014, and references therein). Moreover, Fang et al. (2013) discuss that the same behavior of a steeper dependence at lower M_⋆ is compatible with their data in the similarly young L1641 region, although they did not perform a fit of this relation. Vorobyov & Basu (2009) also suggested that the observed Ṁ_acc–M_⋆ relation, when very many objects are considered, is better reproduced by a segmented line with a break around M_⋆ ~ 0.2 M_⊙. While this qualitatively matches our observed values, they derived a shallower slope than we obtained, although the two values are compatible. Finally, we note that the linear slope derived for the higher L_⋆ and M_⋆ subsample of objects is consistent with the almost linear relation between L_acc and L_⋆ observed in samples of Herbig Ae-Be stars (e.g., Mendigutía et al. 2015), which have stellar masses higher than the objects studied here, but are also in general older.

A possible interpretation for this bimodal distribution on the Ṁ_acc–M_⋆ plane is a different evolutionary timescale for disk accretion around stars with different masses. In this view, disks around stars with M_⋆ ≲ 0.3 M_⊙ will faster evolve to lower values of Ṁ_acc. The complete U-band photometric surveys in σ-Orionis (Rigliaco et al. 2011a) and in the Orion Nebula Cluster (Manara et al. 2012) and the spectroscopic survey of L1641 by Fang et al. (2013) came to the similar conclusion of a faster evolutionary timescale of accretion for lower-mass stars than for solar-mass stars. This is a possible explanation for the difference with the results obtained by Manara et al. (2015) with a similar procedure on a sample of very low-mass stars and brown dwarfs in the younger ρ-Ophiuchus star-forming region. Manara et al. (2015) found that the observed Ṁ_acc–M_⋆ relation for the objects in ρ-Ophiuchus follows the same Ṁ_acc ∝ M_⋆^~1.8 relation as was found for objects of higher stellar masses derived with an incomplete sample in the Lupus star-forming region by Alcalá et al. (2014). Thus, Manara et al. (2015) argued that the Ṁ_acc–M_⋆ relation is the same from brown dwarfs up to solar-mass stars. However, these objects in ρ-Ophiuchus have higher Ṁ_acc than those derived here for objects in the Chamaeleon I region with similar M_⋆, and a simple explanation is then that targets in ρ-Ophiuchus are younger than those in the Chamaeleon I or in the Lupus region. However, the sample studied by Manara et al. (2015) is highly incomplete and biased toward stronger accretors, therefore additional data are needed to confirm this hypothesis.

This hypothesis slightly contradicts the well-established result that the dusty inner disks evolve on longer timescales for very low-mass stars than for solar- and higher-mass stars (e.g., Carpenter et al. 2006; Bayo et al. 2012; Ribas et al. 2015). It should be considered that here we included only objects with M_⋆ ≲ 2 M_⊙, which means that they are in a different mass range than was considered by Ribas et al. (2015). In this context, the result of Ribas and collaborators of a shorter evolutionary timescale for disks around stars with M_⋆ > 2 M_⊙ cannot be compared with our findings. However, different studies in various star-forming regions have shown that several objects still show evidence of the presence of a dusty inner disk from infrared excess, but show no signatures of ongoing accretion (e.g., Mohanty et al. 2005; McCabe et al. 2006; Fedele et al. 2010; Wahhaj et al. 2010; Furlan et al. 2011; Hernández et al. 2014). Similarly, the dubious accretors discussed here show no evidence of ongoing accretion, but they all have non-negligible infrared excess, and their disks are detected with ALMA in five cases (see Sect. 3.3). Therefore, this discrepancy can be ascribed to a different evolutionary timescale for the dusty inner disk with respect to the timescale of the process of accretion of material onto the central star.

A faster evolution for lower-mass stars is opposite to predictions by Alexander & Armitage (2006), who postulated that the viscous timescale increases with stellar mass to explain the Ṁ_acc ∝ M_⋆² relation. As discussed, a linear relation between Ṁ_acc and M_⋆ similar to the relation we observe for M_⋆ ≥ 0.3 M_⊙ is predicted in the context of centrally irradiated accretion disks around solar-mass stars with an active accretion layer (e.g., Mohanty et al. 2005; Hartmann et al. 2006). Intriguingly, the models by Hartmann et al. (2006) imply a faster evolution of disks around low-mass stars and brown dwarfs, which is in line with our results.

An additional possibility is that the bimodal distribution is a result of two different accretion regimes at different stellar masses, as initially suggested by Vorobyov & Basu (2009) to explain the bimodality of Ṁ_acc vs. M_⋆. The authors modeled the evolution of disks considering the self-consistently generated gravitational torques, which efficiently drive accretion onto the central star in solar-mass stars, as well as an effective turbulence implemented in the models, needed to model accretion in the very low-mass stars where the effects of disk self-gravity are weak. Their models result in slightly higher Ṁ_acc than those observed here, but they are compatible within the spread, in particular in the solar-mass range. The measured disk masses for this same sample of targets in the Chamaeleon I region shows that these disks are currently not gravitationally unstable (Pascucci et al. 2016). It is nevertheless possible that these disks were gravitationally unstable at earlier ages.

Further investigations of the effect of different magnetic field topologies in different ranges of M_⋆, or other differences between solar-mass and very low-mass stars, should be pursued to explain our findings.

5.3. Spread of the Ṁ_acc–M_⋆ relation

It was argued by several authors, including Clarke & Pringle (2006) and Tilling et al. (2008), that the observed Ṁ_acc fill the whole observable range of values. If this is the case, then it is not possible to derive constraints on the viscosity law or the dependence of disk on stellar mass from the L_acc–L_⋆ or Ṁ_acc–M_⋆ relations. It is thus relevant to test on this complete sample whether the whole range of observable Ṁ_acc is filled. This observable range extends from the upper boundary L_acc = L_⋆ (Tilling et al. 2008) to the lower boundary imposed by chromospheric emission. As discussed in Sect. 3.3 and by Manara et al. (2013a), the emission by active chromospheres in young stars prevents detecting accretion below the typical values of chromospheric emission. This is then an observational limitation. We therefore expect to be able to detect any accretion rate between these two boundaries.

We show in Fig. 5 the Ṁ_acc–M_⋆ relation together with the upper and lower boundaries. We include in this plot all the objects with disks in the Chamaeleon I region that have been studied with X-shooter and the similarly complete sample of stars with disks in the Lupus region that was studied by Alcalá et al. (2014, 2017). The Lupus region has a similar age to Chamaeleon I and a similar distribution of stellar masses. The data have been analyzed in the same way, therefore the inclusion of this sample allows us to examine a more statistically robust sample. The boundary L_acc = L_⋆ is calculated using the 2 Myr isochrone of the Baraffe et al. (2015) models for M_⋆ ≤ 1.4 M_⊙, and the 3 Myr isochrone by Siess et al. (2000) for higher M_⋆. The points slightly above this line correspond to the same targets above the L_acc = L_⋆ line in Fig. 3. In general, this line is a good upper boundary of the distribution of points for M_⋆ ≳ 0.2–0.3 M_⊙. However, fewer targets are present in the lower-mass range, where stars seem to have much lower values of L_acc/L_⋆, as already observed in Fig. 3. Again, this could be an evolutionary effect, in the sense that these objects already have a much lower accretion rate than the initial one, which could have been L_acc = L_⋆.

The lower boundary is taken from the typical value of the chromospheric emission (L_acc,noise) derived by Manara et al. (2013a) for ~3 Myr old objects, completed with the similar result for higher-mass stars by Manara et al. (2017). As expected, the dubious accretors nicely follow this lower boundary. Moreover, stars with M_⋆≳ 0.5 M_⊙ have values of Ṁ_acc that are spread throughout the observable range. Objects with M_⋆ ≲ 0.2 M_⊙ also show values of Ṁ_acc as low as this observational limit. This is not the case for stars with 0.2 M_⊙≲M_⋆ ≲ 0.5 M_⊙: an empty region just above the chromospheric noise limit is present in this stellar mass range, at Ṁ_acc ~ 10^-10M_⊙/yr. This region is highlighted in Fig. 5 with a blue dot-dashed box. This empty region is also present in the Ṁ_acc–M_⋆ relation derived in the Lupus star-forming region (Alcalá et al. 2017). This region of the plot that is empty of data is of particular interest because the lack of objects with Ṁ_acc and M_⋆ values within the highlighted box in Fig. 5 could imply that disks around stars with this mass evolve very fast once their Ṁ_acc decreases below a certain threshold. This is what photoevaporation predicts: a fast dispersal of the inner disk, thus a fast drop of Ṁ_acc, once Ṁ_acc ≲ 10^-9M_⊙/yr (e.g., Alexander et al. 2014; Gorti et al. 2016). Therefore, this empty region of the plot could be explainable by these models. Further work must be carried out to confirm this hypothesis.

6. Conclusions

We have presented a study of a sample of 94 young stars with disks in the Chamaeleon I star-forming regions. This sample represents 97% of the stars with disks and with M_⋆ ≳ 0.1 M_⊙ in this region.

We have analyzed the spectra of these objects obtained with the ESO VLT/X-shooter to derive the stellar and accretion parameters self-consistently. The objects are distributed on the HRD with a wide spread between the 1 Myr and 10 Myr isochrones of the Baraffe et al. (2015) models, with only three objects located well below the 30 Myr isochrone. Two of these three underluminous targets are known to have an edge-on disk, and we suggest that the third also has a similar disk-viewing geometry. By comparing the measured L_acc with typical chromospheric noise for young stars with similar spectral type as our targets, we found that the measured excess emission for 13 targets is compatible with being largely due to chromospheric emission.

We analyzed the logarithmic dependence of L_acc on L_⋆ and of Ṁ_acc on M_⋆ and found a positive correlation between these quantities. Moreover, stars surrounded by transition disks are well mixed with those harboring full disks in these plots. We have further investigated these logarithmic relations with different statistical tests. Non-parametric analyses suggest that the relation could be described with a model more complex than a single line, in particular, that a break might be present and lead to a steeper relation at lower L_⋆ or M_⋆, and a linear relation at higher L_⋆ or M_⋆. We then fit both relations with either a single line or a segmented line. In the former case, we obtain a slope of 1.9 ± 0.1 and 2.3 ± 0.3 for the L_acc–L_⋆ and the Ṁ_acc–M_⋆ relations, respectively. These values are consistent with previous results and with theoretical expectations from various studies. However, the segmented line fit with a steeper relation for L_⋆ ≲ 0.45 L_⊙, or M_⋆ ≲ 0.3 M_⊙, is a statistically slightly preferred description of the Ṁ_acc–M_⋆ relation according to different statistical estimators, in line with the findings in a similar survey with VLT/X-shooter in the Lupus star-forming region (Alcalá et al. 2017), although the single-line fit is not statistically excluded. We suggest that the steeper relation for lower-mass stars is due to a faster evolution of the accretion process around these objects, as was found with previous photometric studies of complete samples of young stars with disks in different star-forming regions (Rigliaco et al. 2011a; Manara et al. 2012). This result is in agreement with theoretical predictions of faster evolution of disks around very low-mass stars because these disks are smaller, highly ionized, and MRI active (Mohanty et al. 2005; Hartmann et al. 2006). Another possibility is that two different accretion regimes are present, with gravitational instability governing accretion in disks around solar-mass stars, and viscosity at lower masses, as predicted by Vorobyov & Basu (2009).

Finally, by exploring the distribution of measured accretion rates in comparison with the values that might be observed as a result of physical and observational boundaries, we find two main features. First, a lack of very low-mass stars with high accretion rates, possibly due to the same evolutionary effect that causes the steepening of the Ṁ_acc–M_⋆ relation. Second, a lack of targets with M_⋆ ~ 0.3–0.4 M_⊙ just above the observational limits imposed by chromospheric emission, as if the disks around these stars were rapidly dissipated once Ṁ_acc is below a certain threshold. This is what photoevaporation theory predicts, and the distribution of data in our sample may be a sign of ongoing photoevaporation.

Future theoretical work is needed to constrain the hypotheses on this steeper relation of accretion with stellar parameters, and to use this information to better describe the evolution of protoplanetary disks.

---

Acknowledgments

We thank the anonymous referee for the useful comments that helped improved the presentation of the results. C.F.M. gratefully acknowledges an ESA Research Fellowship and support from the ESO Scientific Visitor Programme. This work was partly supported by the Gothenburg Centre for Advanced Studies in Science and Technology as part of the GoCAS program Origins of Habitable Planets and by the Italian Ministero dell’Istruzione, Università e Ricerca through the grant Progetti Premiali 2012-iALMA (CUP C52I13000140001). I.P. acknowledges support from an NSF Astronomy & Astrophysics Research Grant (ID: 1515392). AN acknowledges funding from Science Foundation Ireland (Grant 13/ERC/I2907). D.F. acknowledges support from the Italian Ministry of Education, Universities and Research project SIR (RBSI14ZRHR). We thank E. Whelan for sharing additional information on ISO-ChaI 127. We acknowledge particularly insightful discussions with C. Clarke. This research made use of Astropy, a community-developed core Python package for Astronomy at http://www.astropy.org, and of the SIMBAD database, operated at CDS, Strasbourg, France.

References

Appendix A: Additional data

Additional information from the literature on the new targets discussed here is reported in Table A.1.

Appendix A.1: 2MASS J11175211−7629392: a newly identified binary

The target 2MASS J11175211−7629392 was identified as a member of the Chamaeleon I region by Luhman (2007) and the membership was also confirmed by Lopez Martí et al. (2013) based on proper motion analysis. However, the latter report the proper motion to be dubious.

In the acquisition image of the VLT and in the raw data we clearly see two components for this system with a separation of ~2′′, equivalent to ~320 au at the distance of Chamaeleon I. Since the target was observed at a seeing of ~0.9′′, the two components are resolved and we are able to extract the spectra of the two components separately. The fact that the system

is a binary can explain the dubious proper motions obtained by Lopez Martí et al. (2013).

The analysis of the spectra of the two components shows that both objects have the same spectral type, which is consistent with the type reported by Luhman (2007), that is, M4.5. None of the objects is accreting, and there is no sign of excess emission with respect to the photosphere in the near-infrared part of the spectrum. The only sign of excess emission is from Spitzer 12 μm photometry, but at the resolution of this telescope it is not possible to separate the emission from one or the other component. Moreover, the millimeter ALMA observation of this system did not detect any continuum emission (Pascucci et al. 2016). Therefore, the presence of a disk around these stars is not confirmed, and we thus do not include these objects in the analysis of the dependence of accretion properties with stellar properties.

Appendix B: Additional statistical tests

Here we present some additional tests that we carried out to study the dependence of the accretion on stellar parameters, as discussed in Sect. 4.

Appendix B.1: Non-parametric statistics

Non-parametric statistics allows one to explore the data without assuming an underlying model that describes the data. Therefore, we use this technique to understand what types of models we should use to fit the distribution of data in the L_acc–L_⋆ and Ṁ_acc–M_⋆ planes. In the following, we describe two different types of tests we carried out.

Appendix B.1.1: Medians

We divided the sample into bins comprising an equal number of objects. We carried out the analysis with bins comprising aboutfve or seven objects. We then computed the median value of both logL_acc and logL_⋆ in one case, and logṀ_acc and logM_⋆ in the other. As always, we excluded the edge-on targets from the analysis, while we included the dubious accretors in the sample. The results are shown in Fig. B.1 for logL_acc-logL_⋆, and Fig. B.2 for logṀ_acc-logM_⋆. In both cases the median values show a steeper slope at low L_⋆, or M_⋆, than at higher values, where the slope is almost linear. The break is located at log(L_⋆/L_⊙) ~−0.7 and log(M_⋆/M_⊙) ~−0.6, respectively, as derived from a fit of these median values with a segmented line with the same functional shape as Eq. (1). The best fit of these median values for a single or broken power-law is also shown in Fig. B.1 for L_acc–L_⋆, and in Fig. B.2 for Ṁ_acc–M_⋆, using different colors for the results using five or seen objects per bin. The fit with a broken power-law matches the median values better.

Appendix B.1.2: Non-parametric fit

We used two different non-parametric fit methods available with the Python package PyQt-fit, namely the Nadaraya-Watson, a spatial averaging technique, and the local-polynomial fit on the L_acc–L_⋆ and Ṁ_acc–M_⋆ planes. The local-polynomial fit is carried out using a quadratic and a cubic polynomial. The results are as follows.

For L_acc–L_⋆ there is a ~1 slope at logL_⋆ ≥ −0.6, and a steeper slope for lower L_⋆. In the plot (Fig. B.3) the best fit is shown with a solid line, where different colors correspond to different methods, as reported in the legend. The green filled region represents the 95% confidence level interval on the local-polynomial fit using a quadratic polynomial, derived with a bootstrap technique.

The fits of the Ṁ_acc–M_⋆ are more uncertain (see Fig. B.4), probably due to the larger scatter of points. However, a flatter slope is present for log(M_⋆/M_⊙) ≳ −0.5, and steeper below this value.

Appendix B.2: Treating dubious-accretors in the fit

As discussed in Sect. 3.3, 13 objects in our sample have an accretion luminosity below the typical chromospheric contribution to the excess emission with respect to the photospheric one, thus they are considered as dubious accretors. Here we show that considering their measured value of L_acc as an upper limit or as a detection does not affect the fit of the L_acc–L_⋆ and Ṁ_acc–M_⋆ relations.

First, we fit the data with a single power-law relation using the maximum-likelihood Bayesian tool linmix by Kelly (2007), which derives the linear dependence between logL_acc and logL_⋆ considering measurement uncertainties on both axes. We ran the fit with three different assumptions on the dubious accretors. First, we only included clearly accreting targets, second, we included the dubious accretors considering the L_acc measured for these objects as an upper limit, and finally, we included the dubious accretors considering their measured L_acc as detection. The three cases consistently led to the same slope of 1.9 ± 0.1 of the linear relation, while the intercept varies from − 0.6 ± 0.1 when only accreting objects are included, to − 0.8 ± 0.2 when non-accreting targets are included either as detection or upper limits, as shown in Fig. B.5. Therefore, the choice of including dubious accretors in the analysis is relevant only when determining the intercept of the relation, but the choice of considering them as upper limit or detection does not change the results.

The same result is found when analyzing the Ṁ_acc–M_⋆ relation. The slope of the best fit is compatible when considering dubious accretors as detection, or upper limits, or even neglecting them. These results are shown in Fig. B.6.

Similarly, the choice of considering the accretion rate for dubious accretors as detections is solid even when the sample is divided into two to fit the two power-laws (see Sect. 4), and the slope of the two parts of the segmented line are independent of how the dubious accretors are treated, as we show in Fig. B.7 for the L_acc–L_⋆ relation, and in Fig. B.8 for the Ṁ_acc–M_⋆ relation.

We then test how strongly the results obtained with scipy.optimize.curve_fit are affected by including or excluding the dubious accretors. This tool does not allow including upper limits in the test. The result is shown in Fig. B.9 for the L_acc–L_⋆ relation, and in Figs. B.6 and B.10 for the Ṁ_acc–M_⋆ relation. The slope of the single power-law fit is the same regardless of whether the dubious accretors are included, while small differences are present in the intercept for both the single power-law and the segmented line fit.

Appendix B.3: Comparison between methods considering or excluding measurement errors

We tested whether the uncertainties on the individual measurements contribute to the estimate of the best fit parameters more than the scatter of the data itself. We find that the scatter is the most important effect in determining the fit parameters following a similar approach to Pascucci et al. (2016). We fit the observations using linmix, where the uncertainties are considered, and two methods that derive the best fit without considering the uncertainty: scipy.optimize.curve_fit and cenken in the NADA R package. We discussed in Sect. 4 that the former leads to best-fit values that are compatible with linmix. The same is true for the latter, both when considering dubious accretors as detection or upper limits and when fitting the whole sample or the two subsamples at high and low L_⋆ and M_⋆ separately. Therefore, we conclude that the scatter dominates the fit determination.

Appendix C: Best fit

Here we show the best fit of the Balmer continuum region obtained with our method as described in Sect. 3.1.

Appendix D: Observation log

Appendix E: Stellar masses and mass accretion rates using different evolutionary models

All Tables

All Figures