© ESO 2021

1. Introduction

The final mass of low-mass stars (< 2 M_⊙) is mostly assembled during the protostellar phase. Therefore, it is during this phase that it is crucial to investigate how matter is accreted onto the forming star from the surrounding disk and envelope. The widely adopted model for mass accretion in young stellar objects (YSOs) is the magnetospheric accretion model (e.g., Hartmann et al. 2016). In this frame, the mass accretion rate not only regulates the final stellar mass build-up, but also the evolution of the protoplanetary disk, determining the initial conditions for planet formation (e.g., Manara et al. 2019). In particular, the mechanism driving the evolution of disks is responsible for the transport of material through the disk, and thus, in the end, for the mass accretion onto the central star.

According to the viscous evolution model (e.g., Lynden-Bell & Pringle 1974; Hartmann et al. 1998), protostars follow a stage of vigorous accretion during the initial 10⁵ yrs of life, that is during the Class 0 and I phases, when the protostar is still surrounded by its original envelope, while the mass accretion rate decreases significantly during the following pre-main sequence phase (Class II, ∼10⁶ − 10⁷ yrs) during which the disk is dispersed (e.g., Ercolano & Pascucci 2017) and planetary systems form, until the source approaches the zero age main sequence (ZAMS), that is to say the Class III phase. This broad classification has been empirically based on the spectral index of the spectral energy distribution (SED) of YSOs in the range between 2 and 24 μm (Evans et al. 2009). However, this simple classification might not reveal the complexity of such young stellar objects. For example, detailed 2D models (Whitney et al. 2003; Robitaille et al. 2006) show that, for a correct interpretation of the SED, the geometry of the system and its inclination angle need to be considered.

Observational verification of the above picture is still the subject of several studies. Most determinations of the mass accretion rate in YSOs have been obtained in classical T Tauri (CTT) stars, in other words Class II sources whose SEDs are mostly dominated by the stellar photosphere with little IR excess due to the disk contribution. These sources are accessible to UV and optical observations; therefore, the accretion luminosity (L_acc) and mass accretion rates (Ṁ_acc) can be directly measured from the UV-excess over the stellar photosphere caused by the accretion shock (e.g., Calvet et al. 2004; Herczeg & Hillenbrand 2008; Herczeg et al. 2009; Rigliaco et al. 2012; Ingleby et al. 2013; Alcalá et al. 2014; Fairlamb et al. 2015; Rugel et al. 2018). In addition, correlations between the L_acc and the luminosity of optical and IR permitted lines, such as the Balmer, Paschen, and Brackett lines, have been found (Muzerolle et al. 1998; Natta et al. 2004; Mohanty et al. 2005; Alcalá et al. 2017). Adopting these methods, several surveys of the accretion properties in YSOs of nearby star-forming regions with ages spanning between 1 and 10 Myr have been performed (e.g., Natta et al. 2006; Antoniucci et al. 2008, 2014; Manara et al. 2012, 2015, 2017a, 2020; Biazzo et al. 2012, 2014; Venuti et al. 2014, 2019; Alcalá et al. 2014, 2017; Kim et al. 2016), which have in particular addressed the connection between accretion rates and the stellar parameters.

Very little is known, however, on the accretion properties of younger and more embedded sources, such as the Class I YSOs, whose large extinction, caused by their residual envelope, prevents observations in the UV and optical. These sources are best observed in the IR, although a lack of a good determination of extinction, which can reach values of up to A_V = 40 − 50 mag, can hamper a correct interpretation of even the IR emission lines.

One of the most intriguing open problems on the accretion properties of younger and more embedded YSOs is the luminosity problem: The luminosity of these sources is in general too low, by as much as an order of magnitude, with respect to that expected if accretion proceeds at the high mass accretion rates predicted by the standard steady-state collapse model of Shu (1977). This problem was first noted by Kenyon et al. (1990) and later confirmed by the large core-to-disk (c2d) Spitzer photometric surveys (e.g., Evans et al. 2009). These IR observations suggest that most of the mass was accreted during the protostellar phase or in a non-steady frame (Evans et al. 2009), but direct accretion estimates are needed to confirm this. Assuming that the accretion decreases with time (e.g., Hartmann et al. 2006), studying the accretion rates in these young and poorly characterised objects can help us in understanding the luminosity problem itself and the role of the Class I stage in the overall picture of the accretion process.

Class I objects are also interesting to study as they still have enough mass in their protoplanetary disk to allow planet formation (e.g., Tychoniec et al. 2020), contrary to Class II objects (e.g., Manara et al. 2018). Stellar properties of relatively large samples of Class I sources were obtained by Doppmann et al. (2005) and Connelley & Greene (2010), who analysed the weak photospheric features of these objects which can be detected when they are not too ‘veiled’ by a large IR continuum excess.

Mass accretion rate estimates were obtained on small samples of Class I objects. In most studies, the luminosity of the Brγ and Paβ lines was used as a proxy for the accretion luminosity (Muzerolle et al. 1998; Nisini et al. 2005; Antoniucci et al. 2008, 2011). White & Hillenbrand (2004) derived both the accretion and stellar parameters of about 15 Class I stars through the analysis of their optical scattered light. These works highlight the large diversity of behaviour of the investigated Class I objects. In fact, few studied sources show high levels of accretion, which is compatible with a scenario in which their bolometric luminosity is dominated by the accretion luminosity (the so-called accretion dominated young objects or ADYOs, Antoniucci et al. 2008), while the majority of the samples have mass accretion rates comparable with those of Class II stars of a similar mass. These results, however, refer to small and inhomogeneous samples of sources located in different environments. Hence general conclusions on the accretion properties of young embedded sources with respect to more evolved objects in the same star-forming cloud can hardly be drawn from these studies. Nowadays, thanks to sensitive, high multiplex instrumentation in the IR, it is possible to perform more systematic surveys of embedded YSO populations in young star-forming clouds, thus contributing to the understanding of the accretion process in the early stages of star formation.

In this work we present the results of a spectroscopic survey of YSOs in the young (< 1 Myr) NGC 1333 cluster, located at 293 ± 22 pc (Ortiz-León et al. 2018, which takes into account recent Gaia results) in the Perseus cloud. This survey was performed with the K-band Multi Object Spectrograph (KMOS, Sharples et al. 2013) at the ESO Very Large Telescope (VLT).

The Perseus complex is a suitable laboratory for the study of low mass star formation because of YSOs at different evolutionary stages. The complex comprises several clusters, the largest being the IC 348 and NGC 1333 clusters (Sargent 1979; Kirk et al. 2006; Enoch et al. 2007).

The prestellar cores and the YSO population of Perseus have been extensively investigated at different wavelengths: in the NIR and MIR (Rebull et al. 2007; Jørgensen et al. 2007; Evans et al. 2009; Enoch & Evans 2009; Gutermuth et al. 2011; Young et al. 2015; Mercimek et al. 2017), in the FIR and sub-millimetre (Enoch et al. 2006; Hatchell et al. 2007a,b; Sadavoy et al. 2014; Tychoniec et al. 2020; Pezzuto et al. 2021), and at millimetre and radio wavelengths (Pech et al. 2016; Segura-Cox et al. 2018; Andersen et al. 2019).

The NGC 1333 cluster shows the highest fraction of Class I, in other words sources with α > 0.3. Even if the age of the main clusters, IC 348 and NCG 1333, is still debated, there is a general agreement on the fact that star formation in NGC 1333 has taken place considerably later than in IC 348. In particular, the original work from Sargent (1979) estimated the ages of IC 438 and NGC 1333 to be about 1.3 Myr and less than 1 Myr, respectively. Those ages have been confirmed by the latest studies in which the age of NGC 1333 has been measured to be < 1 Myr (Wilking et al. 2004; Young et al. 2015; Yao et al. 2018).

Optical and IR spectroscopy of YSOs in the NGC 1333 cluster has been performed by several authors and used to derive the stellar parameters, such as the spectral type, stellar luminosity, and mass, of the members of the clusters (the latest works are Winston et al. 2009; Itoh et al. 2010; Luhman et al. 2016; Esplin & Luhman 2017). In particular, Luhman et al. (2016) have reported the discovery of new members not previously observed with Spitzer, thus providing a brown dwarf catalogue of this cluster, later analysed by Esplin & Luhman (2017). Accretion properties of the NGC 1333 YSOs, however, have never been derived before.

The aim of our work is to fill in this gap, providing the first homogeneous study of the accretion properties of a sample of Class I and II sources in the NGC 1333 cluster through near-infrared (NIR) spectroscopy. The results of our survey are then compared with those derived in older regions to obtain information about the evolution of the accretion luminosity and the mass accretion rate with the age and the evolutionary stage of the YSOs.

The paper is organised as follows. In Sect. 2 we present the selection of our sample, the observations performed with KMOS, and the data reduction steps; in Sect. 3 we report how we measured the Paβ and Brγ lines, which are the lines we used to derive the accretion luminosity; in Sects. 4 and 5 we describe how we computed the stellar parameters and the accretion rates for Class II and I YSOs, respectively; in Sect. 6 we discuss our results; finally, in Sect. 7 we provide the summary and conclusions of this work.

2. Sample, observations, and data reduction

2.1. The sample

We selected our sample from the catalogue provided by Young et al. (2015), which contains 369 candidate YSOs in the Perseus cloud, classified using Spitzer photometry up to 70 μm and 2MASS photometry. In this catalogue 103 sources belong to the NGC 1333 cluster (40 Class 0+I/Flat, 58 Class II, and 5 Class III). To be able to observe those sources with KMOS, we selected only the targets which have a 2MASS counterpart (88 targets) and m_K < 14.5 mag, that is to say 84 targets (23 Class I/Flat, 56 Class II, and 5 Class III). We note that by considering sources that are bright in the NIR, we automatically removed any Class 0 from our sample. Therefore we refer to the sources of our sample as Class I objects.

A prioritised list of targets was provided as the input catalogue to the KMOS ARM Allocator (KARMA)¹ software. The allocation of targets to arms follows the priority of the targets (Priority One) is assigned to Class I/Flat, Priority Two to Class II, and Priority Three to Class III) and the positional constraints due to the mechanical properties of the arms. We identified three fields covering the largest number of sources (Fig. 1). Within these fields, KARMA selected 61 targets, that is to say 60% of the initial sample. These targets are listed in Table A.1.

Fig. 1. 2MASS mosaic map of the NGC 1333 cluster achieved through the Image Mosaic Service (http://montage.ipac.caltech.edu/). Big magenta circles correspond to the KMOS pointings of 7.2′ diameter. Circles refer to the sample of YSOs with mK > 14.5 mag. The colour coding is the following: red (blue) and green (yellow) symbols indicate observed (not-observed) Class II/III and I/Flat sources, respectively.

Fig. 2. Histogram of the K-band magnitude of the Spitzer identified YSOs in the NGC 1333 cluster (black, Young et al. 2015) and the sample observed by KMOS (red). The black-dashed vertical line represents the cut-off at 14.5 mag due to KMOS sensitivity.

According to the classification given by Young et al. (2015), 20 of them are Class I/Flat, 39 are Class II, and two are Class III objects. Therefore, in total we successfully observed 50% of the Class I/Flat and 67% of the Class II sources. Our sample is 87% and 70% complete within m_K < 14.5 mag for Class I/Flat and Class II, respectively (for more details see Figs. A.1 and A.2).

2.2. KMOS observations

The observations were carried out between October 2 and 31, 2019 in service mode, with KMOS. KMOS is able to perform integral field spectroscopy in the NIR of 24 objects simultaneously within a 7.2′ patrol field. The targets were selected by 24 pick-off arms which place small mirrors in the focal plane to send a field of view of 2.8″ × 2.8″ to an image-slicing integral field unit (IFU).

For each of the three fields, shown in Fig. 1, we performed two observations, one optimised for bright targets (m_K < 11 mag) and the other for the faintest ones (11 mag < m_K < 14.5 mag).

Each observation was performed with a dithering sequence (Freedither mode) with one cycle of a six-point dither pattern on target. The dither step size is 0.72″ on right ascension (RA) and declination (Dec). This pattern is interleaved with sky observations with an ABAA configuration for bright targets and ABA for dim targets, where A is the pointing on the target and B is the pointing on the sky. At each position of the pattern, we performed a long (∼280 − 300 s) or a short (∼60 − 80 s) exposure in order to have S/N > 50 for all the faint and the bright targets, respectively. Observations were taken, using this strategy for each J, H, and K band grating. The resolving power at the band centre is about ∼3600, ∼4000, and ∼4200, respectively.

Three reference stars are acquired in each field. A set of calibration images was taken in addition to each of the observations. This set consists of dark frames, flat fields, and an arc-lamp exposure. A G-type telluric standard star from the HIPPARCOS catalogue is observed in each field.

Further observational details are reported in Appendix A.

2.3. Data reduction

Data reduction was performed with the kmos-2.1.0 pipeline (Davies et al. 2013) in combination with the ESOReflex pipeline, version 2.9.1 (Freudling et al. 2013). The pipeline routines were used to calibrate and reconstruct the 3D data cubes from the raw data of the scientific targets and standard star. The procedures used by the pipeline of KMOS are mostly standard for NIR spectroscopy: flat field correction, wavelength calibration, sky sybtraction, extraction of the spectra, telluric correction, and instrumental effects correction. We corrected the atmospheric transmission using a model of telluric absorption. This model was obtained using the molecfit code (Smette et al. 2015) applied to the standard star spectrum, normalised for a black body function and the instrumental response. This telluric model was then applied to all the spectra of the targets to remove the telluric absorption effect. The reconstruction of the 3D data cubes was performed by interpolating from tables which provide the (x, y, z) location in the final cube for each pixel on the detector.

From the calibrated 3D data cubes, we then extracted the 1D spectra (see Appendix B). In order to maximise the signal-to-noise ratio (S/N) and be sure to avoid the contribution of spurious features in the spectrum, we extracted the spectra in a region of 0.6″ × 0.6″ (3 × 3 pixels), centred on the peak of the source. This aperture is typically smaller than the airmass-corrected seeing, which is usually ∼0.8 − 1″ at the wavelength range and time of the observations. This assures us to be dominated by stellar emission, minimising the background contribution. For those sources which are not well centred in the small IFU field of view, the spectrum has been extracted from a narrower aperture centred on the source peak. Since the KMOS spectro-photometry was not reliable, we decided to calibrate the flux scale using the 2MASS photometry (see Table 1²) for all targets. The median scale factors (2MASS/KMOS photometry) are 2.4, 2.7, and 2.6 in J, H, and K band, respectively, and they range between 0.94 and 5.10.

On the basis of the observed spectra, eight sources have not been considered as YSOs. As a result, we removed them from the analysis in this work (see Appendix C).

3. Measurement of Paβ and Brγ lines

The brightest lines tracing the accretion process in the NIR are the Paβ at 1.28 μm and the Brγ at 2.16 μm lines. We measured the flux of these two emission lines in each target using the Image Reduction and Analysis Facility (IRAF, Tody 1993). The line flux was determined by simply summing the flux of the pixels contained between the continuum and the line. We repeated the same measurement three times per line, considering the highest, lowest, and the middle position of the local continuum, depending on the local noise in the spectrum. We used the mean and the standard deviation of the measurements as the best estimate for the line flux and its error. We considered a line detected when its S/N > 3.

For the non detected lines, we measured the root mean square (RMS) of the observed flux between 1.275 μm and 1.285 μm for Paβ and between 2.160 μm and 2.170 μm for Brγ. Then we computed the upper limit of the line flux ($F_{\mathrm{line}}^{\mathrm{upp}}$) as in Nisini et al. (2018):

$$\begin{array}{c}\hfill F_{\mathrm{line}}^{\mathrm{upp}}=3\times \mathrm{rms}\times ({\lambda }_{\mathrm{line}}/R)\end{array}$$(1)

where λ_line is the central wavelength of Paβ or Brγ and R is the resolution of the instrument in the J or K band, respectively.

As a result, 12 Class II YSOs show both Paβ and Brγ, seven show only Paβ, and one shows only Brγ; in three Class I YSOs, we detect both Paβ and Brγ, and in seven of them we detect only Brγ. Therefore ∼55% of Class II (20 out of 36) and ∼59% of Class I (ten out of 17) have at least one accretion line tracer above the threshold noise.

Detected line fluxes and upper limits are reported in Table 2. Taking into account the uncertainty on the absolute flux due to the variability of CTT stars (∼0.5 mag, Lorenzetti et al. 2013), the fluxes we measured could vary by about 50%.

4. Analysis of the Class II sources

In this section we present the analysis of the spectra of the Class II sources. Source #178 is a Class III YSO, therefore is not included in the discussion. However, its Paβ and Brγ fluxes and stellar parameters are reported in Tables 2–4, marked with an asterisk (*).

The analysis performed consists in measuring the stellar parameters of the sources (i.e., spectral type, extinction, star luminosity, and mass) and their accretion parameters (accretion luminosity and mass accretion rate) from the luminosity of the Paβ and Brγ emission.

4.1. Stellar parameter determination

To determine the spectral type (SpT) and the extinction (A_V) of each target, we compared our spectra with a grid of photospheric templates consisting of observed spectra of non-accreting Class III YSOs (Manara et al. 2013, 2017b), which cover spectral types from G4 to M9.5, with a typical step of one spectral sub-class for G and K types, and 0.5 spectral sub-classes for the M type. For each source, we estimated the SpT and A_V by varying the template and the extinction and by matching the shape and the molecular features (mostly H₂O in H band) of our data with the ones of the Class III template (see Appendix B of Manara et al. 2013, for a list of NIR molecular bands). This is a procedure typically used in the literature (e.g., Alcalá et al. 2014; Manara et al. 2015). To compare the templates to KMOS spectra, we reddened and normalised them within a window of 0.02 μm around λ = 1.60 μm, where no molecular features are present. The comparison was made by eye, taking into account the match with the molecular bands and the overall shape of the spectrum in the J and H bands. We did not consider the K band for matching due to the IR excess expected in this band from to the disk emission. An example of a spectral typing result is shown in Fig. 5.

The J band of the KMOS spectra is in general noisier than the other bands, and for several sources the S/N is too low to compute the scale factor in this band. We verified that the spectral typing changes at most within one sub-Class and one mag in A_V, scaling the template flux to the flux of our data in the J or H band for the sources whose J band is not too noisy. Therefore, we decided to use H band for the overall sample. This method can provide some small degeneracies in SpT and A_V, because the presence of certain molecular bands is related to the spectral class.

We adopted ±1 sub-class and ±1 magnitude for the uncertainties of spectral type and extinction, respectively.

We compared our SpT and A_V with the ones estimated by Luhman et al. (2016) for the sources in common, finding a general agreement and verifying that different methods give consistent results. We discuss this comparison more in detail in Appendix E.

We further estimated A_V with other methods and compared the results. The first method is based on the use of the [J − H] vs. [H − K] colour–colour diagram. This is shown in Fig. 3 together with the observed sample and the extinction-vector, which indicates the direction of the colour difference of a star extincted by A_V = 5 mag according to the Cardelli et al. (1989) extinction law, with the standard value R_V = 3.1. We estimated the A_V from the colour difference between the position of the observed targets and the CTTS locus (Meyer et al. 1997), following the extinction-vector.

Fig. 3. Colour–colour diagram of our sample. Class I, II, and III YSOs are plotted as green, red, and blue crosses, respectively. Triangles are lower limits. Class I YSOs which show HI in emission are surrounded by a black circle. The main sequence (MS, Luhman et al. 2010) and CTT locus (Meyer et al. 1997) are shown in black and purple, respectively. The green arrow in the left-top corner defines the reddening direction, following the Cardelli et al. (1989) extinction law. Its length corresponds to AV = 5 mag.

Fig. 4. HR diagram for the Class II YSOs of NGC 1333 (green stars), Lupus (orange crosses), and Cha I (purple asterisks). The green, orange, and purple lines correspond to the median of NGC 1333, Lupus, and Cha I samples (sub-luminous objects are not included), respectively. The birthline (Palla & Stahler 1993) is plotted as a dot-dashed line. Black solid and dashed lines show the isochrones and evolutionary tracks by Siess et al. (2000).

Fig. 5. Example of the spectral typing result for Class II source #156. The observed spectrum is shown in red. The dereddened spectrum (black) matches the spectrum of an M6.5 Class III (green) from Manara et al. (2013). The difference in the K band is due to the fact that source #156 is accreting, while the Class III is not.

Another method is to estimate A_V from the accretion luminosity obtained with tracers at different wavelengths. Indeed, as we show in the next section, for those sources in which both Paβ and Brγ are detected, we can estimate the value of extinction that makes $L_{\mathrm{acc}}^{{\mathrm{Pa}}_{\beta }}=L_{\mathrm{acc}}^{{\mathrm{Br}}_{\gamma }}$.

Table 3 summarises A_V estimates derived from all methods. The A_V derived from the comparison with the reddened spectral templates are generally consistent with those derived from the colour–colour diagram within 2 mag for the majority of sources. Larger differences occur in a few cases, including sources # 165, #180, #197, and #213. The extinction derived from the colour–colour diagram is on average smaller than that estimated through the other methods, which can be explained by the presence of scattered light in the source environment.

In the following, we refer to A_V as the extinction derived by comparing the spectra with the reddened templates ($A_V^{\mathrm{SpT}}$). The SpT and A_V for each Class II YSO is reported in Table 4.

Knowing the SpT, we can then estimate the star luminosity (L_⋆) from the observed magnitudes corrected for the extinction and assuming a bolometric correction (BC). We used bolometric correction values in the J band (BC_J) of 5 − 30 Myr stars from Table 6 of Pecaut & Mamajek (2013) for spectral types down to M5, and from Herczeg & Hillenbrand (2015) for SpT from M5.5 to M7.5. From the same works, we also took the effective temperature (T_eff) corresponding to the estimated SpT. In several cases, the SpT class or sub-class of our targets is not present in these works. In such cases, we obtained BC and T_eff by interpolating with the available SpTs.

Then, L_⋆ is given by:

$$\begin{array}{c}\hfill log\left(\frac{L_{\star }}{L_{\odot }}\right)=0.4[M_{\mathrm{bol},\odot }-M_{\mathrm{bol}}],\end{array}$$(2)

where M_bol, ⊙ is the bolometric magnitude of the Sun (Mamajek et al. 2015), and the bolometric magnitude of a source is M_bol = m_J − 5 log(d/10[pc]) + BC_J where m_J is the extinction corrected magnitude. The typical uncertainty in computing the luminosity in this way is of 0.2 dex in log L_⋆/L_⊙, taking into account the uncertainty of the observed magnitude and of the SpT.

After deriving L_⋆ and T_eff for each source, we used the evolutionary tracks by Siess et al. (2000) to determine the corresponding M_⋆. Since Siess et al. (2000) models are available down to 0.1 M_⊙, for objects with lower masses, we used the Baraffe et al. (2015) evolutionary tracks instead. As investigated in Alcalá et al. (2017) and Manara et al. (2017b), the results of both stellar parameters and accretion rates are compatible within each other by using different evolutionary models. The range of stellar masses derived for our sample is 0.08 − 0.72 M_⊙, with 0.1 dex in log(M_⋆/M_⊙) of uncertainty (see Table 4).

We plotted the locus of the NGC 1333 stars in the HR diagram in Fig. 4. For comparison, in this diagram, we also plotted the values for Lupus and Chamaeleon I (Cha I) samples taken from Alcalá et al. (2019) and Manara et al. (2019), respectively. We note that there is one target (#224) below the 30 Myr isochrone. This is highlighted in Fig. 4 with a red box. The spectrum of this target shows HI emission lines and features typical of a star as young as the others in the sample. Therefore, we suggest that this target is a candidate sub-luminous object, where the low apparent luminosity can be due to the presence of an edge-on disk which occults the central star (e.g., Alcalá et al. 2014).

By looking at the median stellar luminosity of the sources, we observe that the Class II targets in NGC 1333 have a systematically higher stellar luminosity than the targets in the Lupus and Cha I regions. This is an indication that, as expected and in agreement with previous estimates of this cluster age, NGC 1333 is younger than Lupus and Cha I.

Assuming that the central object emits as a black-body, we computed the stellar radius:

$$\begin{array}{c}\hfill R_{\star }=\frac{1}{2T_{\mathrm{eff}}^2}\sqrt{\frac{L_{\star }}{\pi \sigma }}\end{array}$$(3)

where σ is the Stefan-Boltzmann constant. The mean uncertainty on the radius is about ±0.6 R_⊙. All the stellar parameters are listed in Table 4.

4.2. Accretion rate estimate

As discussed in the introduction, the accretion luminosity (L_acc) of the sources can be measured by adopting the relationships that connect L_acc with the luminosity of the Paβ and Brγ lines (L_line).

We computed the line luminosity as L_line = 4πd²F_line, where F_line is the extinction corrected line flux and d is the distance to the NGC 1333 region. Then, we estimated the accretion luminosity by using the relations from Alcalá et al. (2017)

$$\begin{array}{c}\hfill log\left(\frac{L_{\mathrm{acc}}}{L_{\odot }}\right)=alog\left(\frac{L_{\mathrm{line}}}{L_{\odot }}\right)+b,\end{array}$$(4)

where a = 1.06 ± 0.07 and b = 2.76 ± 0.34 for Paβ, and a = 1.19 ± 0.10 and b = 4.02 ± 0.51 for Brγ. The value of the accretion luminosity for each Class II source is reported in Table 4. For the targets observed twice, we considered the mean value of the two L_acc estimates. We computed the uncertainty on the accretion luminosity by propagating the error of the coefficients and L_line, obtaining a typical error of 0.4 dex in log(L_acc/L_⊙). When there is no detection of Paβ nor Brγ, we report in Table 4 the larger 3σ upper limit. The accretion luminosity in our sample is −3.57 ≤ log(L_acc/L_⊙) ≤ −0.55.

From the derived L_acc, M_⋆, and R_⋆, we computed the mass accretion rate Ṁ_acc using the relation

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

where R_in is the inner-disk radius which we assume to be R_in ∼ 5R_⋆ (Hartmann et al. 1998). By propagating the uncertainties, we estimated the error on Ṁ_acc, finding an average value of Δ Ṁ_acc/Ṁ_acc = 0.08. The mass accretion rate of our sample ranges between ∼10⁻¹¹ and 10⁻⁷ M_⊙ yr⁻¹. The derived Ṁ_acc are listed in Table 4.

5. Analysis of Class I YSOs

Class I objects have a large residual envelope that significantly contributes to the NIR emission and cannot be separated from the photospheric and disk contribution. Consequently, for those sources, it is not possible to adopt the same classification method we applied to the Class II sample. Indeed, estimating the extinction from the colour–colour diagram (Fig. 3) for this sample is not possible because half of the sources have upper limits in J magnitudes and, in principle, we cannot assume that the CTTS locus is also valid for Class I YSOs. Therefore, for the observed Class I objects, we adopted a different method, already applied by Antoniucci et al. (2008), based on a procedure which we describe in the next section.

5.1. Method

The self-consistent method we used is based on the computation of the stellar luminosity in two equations which both involve extinction. Firstly, the stellar luminosity can be computed by assuming that the bolometric luminosity of YSOs is due to the sum of the accretion and the stellar contribution

$$\begin{array}{c}\hfill L_{\star }=L_{\mathrm{bol}}-L_{\mathrm{acc}}.\end{array}$$(6)

As shown in Sect. 4.1, the stellar luminosity depends on the extinction and the bolometric correction. In Class I sources, we also need to take the contribution of the emission into account due to the envelope and the disk, which we have considered negligible in the H-band of Class II sources. This contribution is generally described with the veiling parameter r = F_ex/F_⋆, representing the excess emission above the photosphere. Therefore for Class I sources for which we used the Brγ line flux, in Eq. (2), the bolometric absolute magnitude can be given as follows:

$$\begin{array}{c}\hfill M_{\mathrm{bol}}={\mathrm{BC}}_K+m_K+2.5log(1+r_K)-A_K-5log(d/10\mathrm{pc}).\end{array}$$(7)

Secondly, as discussed in Sect. 4.2, it is possible to compute the accretion luminosity from the flux of HI emission lines (Eq. (4)), which depends on the extinction as well.

The observed parameters that enter into Eqs. (6) and (7) are the Brγ flux, the observed K-band magnitude (m_K), the bolometric luminosity, and the distance, while assumptions on the K-band veiling (r_K) and on the age of the object, which set the spectral type and thus the bolometric correction (BC_K) in Eq. (7), need to be made. If all these parameters are determined, then Eqs. (6) and (7) only depend on A_K, so that we can find the value of the extinction which consistently provides the same L_⋆ from the two computations.

For our analysis, we assumed that our Class I sources are located between the birthline, as defined by Palla & Stahler (1993), and the 1 Myr isochrone of the Siess et al. (2000) models. Appendix D displays the corresponding evolutionary tracks and how the sources are displayed on them.

Once the previous quantities were determined, the other stellar and accretion parameters, namely M_⋆, R_⋆, and Ṁ_acc, were derived as described in Sect. 4.1.

5.2. Results

Our sample comprises 17 confirmed Class I/Flat sources (see Table 2). In only ten of them, the Brγ line was detected at minimum. We notice that there is only once Flat source (namely source #202) in our sample that shows accretion features. We analysed only sources with at least one HI line detection.

Our adopted method relies on a good determination of the source bolometric luminosity. For this reason, we computed L_bol by integrating the observed spectral energy distributions (SEDs) of the sources including, in addition to NIR and Spitzer photometry, also more recent Herschel photometric measurements and millimetre measurements, when available. Details on the construction of SEDs are given in Appendix F, while the derived L_bol for the Class I accretors are listed in Table 5.

As discussed in the previous section, the method also assumes a value for the K-band veiling of our sources. The effect of the veiling is to fill the photospheric absorption lines, decreasing their equivalent widths and making them undetectable for extreme veiling values. From the inspection of the K-band spectra of our ten accreting sources, we see that photosperic lines typical of late-type stars, such as the Na doublet at 2.209 μm or the Al I at 2.117 μm (see e.g., Nisini et al. 2005; Doppmann et al. 2005), are detected in four sources (namely #202, #215, #231, and #232), while in the other ones such lines are either very weak or not detected. For Class I sources where photospheric features have been detected, we assume r_K values ranging between 1 and 3, as estimated in the samples of Class I objects in the literature (e.g., Doppmann et al. 2005; Connelley & Greene 2010). The same works estimate r_K ≳ 5 in sources where photospheric lines are not or are barely detected. In the bright source #191, where the continuum was measured at the highest S/N, weak photospheric lines of Al and Mg are detected at a 5σ level (Fig. 6). Comparing the equivalent width of these lines with templates of stars with SpT between K7 and G1, we estimate that the veiling should be between 3 and 6. To be conservative, we explored a range of veiling between 3 and 8 for the source #191 and the other six sources where photospheric lines could not be detected.

Fig. 6. K-band spectrum of SVS13 A, normalised to the continuum. Red-dashed lines correspond to photospheric lines as labelled.

The results of our analysis are shown in Table 5, where we report the range of values for the stellar and accretion parameters that we derived from our assumptions.

All Class I YSOs but three (namely sources #186, #191, and #200) are late-type stars (K-M), with stellar parameters (mass and luminosities) in the range of those estimated for the Class II sources of the region. Three sources are, however, more massive and luminous. In particular, source #191 (SVS13 A) and #186 have masses in the range of 2.5 − 3.5 M_⊙, and stellar luminosities between 12 L_⊙ and 40 L_⊙. These results are due to the high bolometric luminosities of the two sources and also to the derived accretion luminosity, which can account for a maximum of 60% of L_bol.

The derived extinction is high in all the sources, and it ranges between 12 and 51 mag. Finally, the accretion luminosity and mass accretion rate range are 10⁻² − 10² L_⊙ and 10⁻⁹ − 10⁻⁶ M_⊙ yr⁻¹, respectively. We discuss this last issue in Sect. 6.2 and the peculiarity of source #191 in Sect. 5.4.

5.3. Upper limits on L_acc in Class I with no line detection

Upper limits on the L_acc for seven Class I sources show that no HI emission lines have been derived from the 3σ upper limits of the Brγ line (see Table 2), and that applying a reasonable estimate of the extinction, based on the spectral features such as IR excess and photospheric lines.

Sources #147, #190, #196, and #237 show small or no IR excess. Therefore, it is possible to estimate A_V and SpT by comparing their spectra with Class III templates, as we did for the Class II sample (see Sect. 4).

Sources #161 and #183 show some IR excess, but also photospheric lines in the K-band, as MgI, AlI, and SiI. For these sources, we assumed the mean of the values we found as extinction for the low-veiling sources (A_V = 20 mag).

Source #218 shows a large IR excess and no photospheric absorption lines. For this source we assume the typical extinction value we find for the highly veiled sources (A_V = 40 mag).

The derived upper limits for the above sources under these assumptions are listed in Table 6.

5.4. The case of SVS 13A

The most luminous source of our sample, source #191, is the well known variable Class I object SVS 13A (or V512 Per), driving the chain of the Herbig Haro objects HH7-11. The source is part of a multi-object system composed of SVS 13A, B, and C (Looney et al. 2000). SVS 13A itself is a close radio binary separated by 0.3″ (VLA 4A and 4B Anglada et al. 2000), while another radio object has been identified at about 5″ (SVS13 A2 Tobin et al. 2016). The NIR source responsible for the HH7-11 outflow has been associated with VLA 4B (Hodapp & Chini 2014; Lefèvre et al. 2017).

SVS 13A had a burst of about 1.5 mag in the K-band during 1990 (up to a K-mag of about 8 mag), followed by a slow decline in which it never reached the pre-burst value of about 9.5 mag (Hodapp & Chini 2014). In this paper, we use photometric values that were taken in 2015 (Giannini, priv. comm.) where the K-band magnitude was 8.7 mag, for example, about 0.5 mag weaker than the 2MASS value taken in 2000. The quick rise and slow decline of brightness resemble the characteristics of FUor-like objects; however, the small amplitude of the burst and the spectrum rich in emission lines exhibited by this source is more typical of EXor-type outbursts, as noted by Hodapp & Chini (2014). High resolution spectral images of this source (Hodapp & Chini 2014) have identified an atomic micro-jet surrounded by a series of expanding bubbles, whose dynamical time suggests that the most recent of them can be associated with the photometric outburst from 1990.

Our analysis suggests that SVS 13A is an intermediate mass star of ∼3 M_⊙, accreting at a rate of a few 10⁻⁶ M_⊙ yr⁻¹. This mass accretion rate is in line with the upper values of Ṁ_acc derived on Herbig stars of a similar mass (e.g., Wichittanakom et al. 2020). We exclude that this result could be affected by the binarity of the source. First of all, the companion to VLA 4B is not detected at NIR wavelengths, as testified by the high spatial resolution observations of Hodapp & Chini (2014). And also, observations in the I band performed with AstraLux at Calar Alto did not detect the companion, indicating that the difference in brightness with respect to the primary should be at least 3 mag (Comeron, priv. comm.). Secondly, although at 3.6 cm, VLA 4A and B have comparable brightness, VLA 4B has a much steeper SED at shorter wavelengths and it is about five times brighter than VLA A at 1.4 mm (Anglada et al. 2004; Lefèvre et al. 2017). This suggests that the contribution from VLA 4A to the source bolometric luminosity should be negligible.

From spatially resolved observations of the [FeII] micro-jet of this source, Davis et al. (2011) estimated a mass ejection rate of Ṁ_loss = 1.4 × 10⁻⁶ M_⊙ yr⁻¹, which would imply a Ṁ_acc/Ṁ_loss ∼ 1, that is to say about a factor of 5 − 10 smaller than expected. A factor of a few larger mass accretion rate would be derived in our analysis only by assuming extreme values of the veiling in the K-band, in other words larger than 20. However, the detection of weak photospheric features in the K-band spectrum of SVS13 A clearly indicates that the IR excess of this source is not extreme, and very high veiling values are ruled out. More likely, the Ṁ_loss derived from the micro-jet gives a value of the mass ejected averaged in time; therefore, it takes into account periods of higher accretion activity of the source.

6. Discussion

6.1. Accretion properties of Class II sources

Previous studies on CTTs in different star-forming clouds have shown that a correlation exists between accretion properties (L_acc and Ṁ_acc) and properties of the central star (L_⋆ and M_⋆), see for example, Hillenbrand et al. (1992), Mohanty et al. (2005), Natta et al. (2006), Venuti et al. (2014), Hartmann et al. (2016), Alcalá et al. (2017), Manara et al. (2017b), Rugel et al. (2018), Arun et al. (2019) and Biazzo et al. (2019). These studies have been mainly concentrated on relatively old clouds (i.e., 2 − 10 Myr). Here, we want to investigate these correlations for our sample of Class II stars in NGC 1333 and how the results in this cluster are related to those of older regions.

6.1.1. Accretion luminosity

Figure 7 shows the distribution of L_acc − L_⋆ for our sample. For comparison, we report on the same plot and also the values relative to the samples of the Lupus and Chamaeleon I clouds (Alcalá et al. 2017; Manara et al. 2017a, rescaled for Gaia distances by Alcalá et al. 2019; Manara et al. 2019). The figure highlights that, for the stars of NGC 1333, L_acc correlates with L_⋆ in a similar way as the sources of the other clouds, roughly following a linear slope in a logarithmic scale. Using the IDL routine linmix_err, which takes into account errors on both axes and the presence of upper limits (Kelly 2007), a linear fit through the data points like

$$\begin{array}{c}\hfill log{L}_{\mathrm{acc}}=\gamma log{L}_{\star }+\beta \end{array}$$(8)

Fig. 7. Accretion luminosity as a function of stellar luminosity for the NGC 1333 Class II sample (green filled circles). Orange and purple stars represent Lupus and Cha I sources, respectively. Upper limits are displayed as triangles. Solid lines show the best fit for each star-forming region as specified in the legend. The standard deviation from the best fit is shown with dashed lines. The subluminous target is highlighted with a red box.

yields γ = 1.5 ± 0.2 and β = −1.0 ± 0.1, with a standard deviation of σ = 0.3 ± 0.1. The same values for the fit were obtained excluding subluminous object #224. In order to compare this relationship with those of the Lupus and Chamaeleon I data, we applied the same fitting procedure to these data, considering upper limits, finding the following fit parameters: γ_Lup = 1.6 ± 0.1, β_Lup = −1.3 ± 0.1 with a standard deviation of σ_Lup = 0.6 ± 0.1, and γ_ChaI = 1.9 ± 0.1, β_ChaI = −1.0 ± 0.1, with a standard deviation of σ_ChaI = 0.8 ± 0.1. These values are consistent with the parameter fits given by Alcalá et al. (2017) and Manara et al. (2017b), before correcting the luminosities for the Gaia distance. The best fit relationships, with their uncertainties, are displayed in Fig. 7 for comparison.

We note that the fit through the data of NGC 1333 has a larger uncertainty, by up to a factor of two in the slope, with respect to the other regions, which is mainly due to the lower number of data points and higher fraction of upper limits. This figure also shows that the NGC 1333 and the Lupus samples have a similar slope of the L_acc vs. L_⋆ correlation, while the fit through the Cha I data is steeper, although compatible within the uncertainties. To compare the distribution of L_acc among the three regions, we plotted the histograms of the luminosity ratio L_acc/L_⋆ in Fig. 8. We note that, similarly to the Lupus and Cha I samples, in NGC 1333 there are no sources where the accretion luminosity is larger than the stellar luminosity. This result is in agreement with previous studies and is theoretically justified by the reasonable assumption that at this evolutionary stage the disk and star lifetimes coincide, as explained in Clarke & Pringle (2006) and Tilling et al. (2008). In addition, the NGC 1333 histogram is similar in shape to the Lupus one and peaks at a roughly higher L_acc/L_⋆ value, that is to say log L_acc/L_⋆ between −2.3 and −1.8.

Fig. 8. Histograms of the luminosity ratio Lacc/L⋆ of NGC 1333, Lupus, and Cha I samples in green, orange, and purple, respectively. The green-dashed line corresponds to the histogram of NGC 1333, including the upper limits.

We note that the NGC 1333 sample has a lower spread in L_acc/L_⋆ with respect to Lupus and Cha I (Fig. 8). We think this is mostly due to the smaller number of sources in the NGC 1333 sample and to the lower sensitivity of our observations in detecting low L_acc values. To better discuss this latter point, we also consider the histogram of the L_acc/L_⋆ upper limits in Fig. 8 (green dashed line). This shows that the derived upper limits are not stringent enough to confirm or disregard the possibility that the stars of the NGC 1333 cluster do indeed have a smaller spread. Therefore, we cannot exclude the possibility that this small spread is related to the younger age of the cluster with respect to Lupus and Cha I. A similar effect of an increasing L_acc spread with age is indeed predicted by the protostellar evolution model of Padoan et al. (2014), in which the infall rates from turbulent clouds drive the accretion rates from the collapse to the pre-main sequence (PMS) phase.

In conclusion, our analysis suggests that the relationship between L_acc and L_⋆ in NGC 1333 is similar to that of older regions, although with a smaller spread in the data. Our smaller statistical sample, however, does not allow us to infer whether this small spread is true and due to evolutionary effects or to a limited sensitivity.

6.1.2. Mass accretion rate

The relation between the mass accretion rates vs. the stellar mass for the NGC 1333 sources is shown in Fig. 9 together with the objects in the Lupus and Cha I samples, for comparison. Similarly to what we have done for the L_acc − L_⋆ relationship, we performed a linear fit considering the upper limits on Ṁ_acc in the following form:

$$\begin{array}{c}\hfill log{\dot{M}}_{\mathrm{acc}}=\eta log{M}_{\star }+\zeta .\end{array}$$(9)

Fig. 9. Mass accretion rate as a function of stellar mass for the NGC 1333 Class II sample (green). Detections are displayed as filled circles, while upper limits are shown as triangles. Orange and purple stars represent Lupus and Cha I sources, respectively. Solid lines show the best fit for each star-forming region as specified in the legend, and the standard deviation from the best fit is shown with dashed lines. The subluminous target is highlighted with a red box.

The linear fit of NGC 1333 results in η = 2.6 ± 0.9 and ζ = −7.3 ± 0.6 with a standard deviation of σ = 1.1 ± 0.1. By repeating the fit excluding the subluminous source #224, we obtain the same values for the fit. The only difference is in the fact that the standard deviation decreases to a value of σ = 0.9 ± 0.1. We used this value of the standard deviation in the following. The slope obtained for the sample in the NGC 1333 region is compatible within the error with the slope obtained here for both Lupus and Cha I (η_Lup = 2.1 ± 0.2 and η_ChaI = 2.3 ± 0.3). The fit obtained for the sample in NGC 1333 appears slightly shifted towards higher Ṁ_acc values with respect to the other clouds (whose intercepts are ζ_Lup = −8.2 ± 0.1 and ζ_ChaI = −8.0 ± 0.1). The standard deviation of the best fit for Lupus and Cha I are σ_Lup = 0.6 ± 0.1 and σ_ChaI = 1.0 ± 0.1, respectively. We can thus state here that in this plot the Cha I linear fit is compatible, within the error, with the fit of the other two regions.

We also note that the slope of the linear fit for NGC 1333 is driven by the large number of upper limits for log M_⋆ < −0.6. In fact most of the detections are located above the linear relationship. This is better shown in Fig. 10, which displays the distribution of the Ṁ_acc/M_⋆ values in the three clouds. In NGC 1333, these values show a spread in the range ∼10⁻¹¹ − 10⁻⁷, without a defined peak. This range is in agreement with the accretion values of the CTT sample investigated by Natta et al. (2006) in ρ-Ophiucus (∼1 Myr). However, more recent observations set a better and nearer estimate of the distance of this region, and the corresponding decrease in the luminosity implies an older age of this cloud, as recently confirmed by Esplin & Luhman (2020), which set the ρ-Ophiucus age at about ∼6 Myr. Indeed, more recent estimates of ρ-Ophiucus accretion rates provided by Rigliaco et al. (2011), Manara et al. (2015) are smaller by an order of magnitude than the NGC 1333 CTT sample.

Fig. 10. Histograms of the ratio Ṁacc/M⋆ for NGC 1333, Lupus, and Cha I samples in green, orange, and purple, respectively. The green-dashed line corresponds to the histogram, which also includes upper limits.

When comparing the NGC 1333 Class II sample with Lupus and Cha I regions, the Ṁ_acc/M_⋆ ratio ranges in the same interval, peaking at ∼10⁻⁹. Although the lack of low Ṁ_acc/M_⋆ values in NGC 1333 can be caused by the limited sensitivity of our IR observations with respect to the optical surveys of Lupus and Cha I, this histogram suggests that in NGC 1333 there is a larger fraction of sources with high Ṁ_acc/M_⋆ with respect to the older clouds.

Alcalá et al. (2017) and Manara et al. (2017b) suggest that the log Ṁ_acc vs. log M_⋆ relationship in Lupus and Cha I can be better reproduced by a double-slope power law, as the higher mass (> 0.2 M_⊙) stars show a flatter distribution of Ṁ_acc as a function of M_⋆ with respect to the sources at a lower mass. Moreover, theoretical models by Vorobyov & Basu (2009) predicted a bimodal distribution to be due to the action of large gravitational instabilities in the more massive disks, limiting the disk masses and the corresponding mass accretion rates.

Here, we can investigate if a similar bimodal distribution is present in our NGC 1333 sample by addressing the role of possible evolutionary effects at the origin of such a distribution. Unfortunately, due to the large number of upper limits for M_⋆ < 0.2 M_⊙, the statistics that we have is not sufficient to draw firm conclusions. We note, however, that all the Ṁ_acc values measured in sources with masses of M_⋆ < 0.2 M_⊙ (7 values) are above the fitted single linear relationship and are higher than the corresponding Ṁ_acc of the Cha I and Lupus in the same mass range.

This result therefore suggests that a change in the slope at the lower masses is not present in the sample of Class II of NGC 1333. A similar result is found by Manara et al. (2015) for the very low mass stars in ρ-Ophiucus, which is a star-forming region younger than Lupus and Cha I. A possible interpretation of this absence of the bimodal distribution can be that disks evolve with different timescales around stars with different masses, as already noticed by Manara et al. (2015). In particular, the mass accretion rate of very low mass stars (M_⋆ < 0.2 M_⊙) evolve faster than solar-type stars. This is in agreement with the models by Hartmann et al. (2006) and is opposite to what was predicted by Alexander & Armitage (2006).

In conclusion, our analysis of the mass accretion rate as a function of the stellar mass gives indications of larger Ṁ_acc/M_⋆ values, on average, for the sample of YSOs in NGC 1333 with respect to the older regions Lupus and Cha I, despite the poor statistics and the upper limits on the low mass regime that prevent us from drawing firm conclusions.

6.2. Accretion properties of the Class I sample

According to the viscous disk accretion model, the mass accretion rate is high during the first 3 × 10⁵ yr, after which it declines with time (Hartmann et al. 1998; Alexander & Armitage 2009; Gorti & Hollenbach 2009). Therefore we would expect that objects in the Class I stage, which are assumed to be younger than Class II sources, would show higher accretion rates with respect to the latter. In the following we discuss whether our results are in agreement with this prediction or not through the comparison of the accretion properties of Class I and II objects of the same cloud.

Figures 11 and 12 show the distributions of L_acc − L_⋆ and Ṁ_acc − M_⋆, respectively, for both the Class I and the Class II samples of NGC 1333. As discussed in Sect. 5, the inferred spread of stellar and accretion parameters for Class I objects is large, and it is caused by the assumptions made on the possible range of veiling and age for these sources. Consequently, the Class I results in these plots are not displayed as individual data points with errors, but rather as areas covering the possible range of values. The object number is also indicated to facilitate the identification of the different sources. The colour of each source is related to the veiling. Sources with 1 < r < 3 are displayed in magenta, while highly veiled (r > 3) sources are in black. Less veiled sources show a lower accretion luminosity than the Class I YSOs with r > 3. In also considering the Class II sources, which have a smaller veiling in general, our results suggest that the more embedded and veiled the objects, the higher the accretion luminosity is, as already noted in previous works (e.g., Calvet et al. 2004; Fischer et al. 2011).

Fig. 11. Accretion luminosity as a function of stellar luminosity for NGC 1333. The Class II sample is shown as in Fig. 7; Class I sources with veiling between 1 and 3 are plotted in magenta, while those with veiling between 3 and 8 are in black.

Fig. 12. Mass accretion rate vs. stellar mass of the NGC 1333 cluster. Class II sample as in Fig. 9. Class I sources with veiling between 1 and 3 are plotted in magenta, and the ones with veiling between 3 and 8 are in black. The black dashed-dotted line represents a linear relation with a unitary slope.

Figure 11 shows a diversity of behaviour among the various sources. The stars with a mass > 0.6 M_⊙ and the sources #231 and #232 are consistent with the relationship found for the Class II objects, within the associated uncertainty. On the contrary, the other sources have L_acc appreciably higher than Class II sources of the same stellar luminosity, up to a factor of 2.

The spread in age that we assume for the Class I analysis is at the origin of the even higher uncertainty in the range of possible Ṁ_acc and M_⋆ values, which is displayed in Fig. 12. A wider spread in the stellar mass is present in the Class I sample, with respect the CTTs.

Given the small number of sources and the large uncertainty on the parameters, no linear fit was performed on this sample. Nevertheless, we can infer that low-mass sources have, on average, mass accretion rates higher than Class II YSOs of the same mass, up to more than two orders of magnitude. In fact, the distribution of values is compatible with a relationship Ṁ_acc ∝ M_⋆, which is indicated in the plot with a dashed line. Hartmann et al. (2006) suggested that this form of relationship is valid in the limit of pure irradiation with pure viscous heating, for stars with a mass up to 2 M_⊙.

The degeneracy of the parameters for source #200 does not allow us to infer if this is a low mass source with a mass accretion rate of the order of 10⁻⁶ M_⊙ yr⁻¹ or a 2 M_⊙ source with a factor almost 100 lower value of Ṁ_acc. On the contrary, the mass accretion rates of the two more massive stars #186 and #191 are quite constrained between 10⁻⁷ and a few times 10⁻⁶ M_⊙ yr⁻¹, respectively, that is to say in agreement with the relationship found for the Class II sources (assuming that this relation is also valid at higher masses).

The mass accretion rates of these two sources are consistent with the upper values found by Wichittanakom et al. (2020) in Herbig stars of similar masses, as already noted for source #191 in Sect. 5.4.

We therefore conclude that the accretion luminosities and mass accretion rates derived for the Class I YSOs are, on average, higher than those of Class II stars of the same mass. However, this finding does not necessarily imply that the sources we have investigated are in their main accretion phase, that is the phase in which a protostar is accreting most of its mass and its luminosity is dominated by the accretion luminosity. Previous studies on the accretion properties of Class I sources have shown that indeed only a small fraction of them have L_acc/L_bol > 0.5 (e.g., White & Hillenbrand 2004; Nisini et al. 2005; Antoniucci et al. 2008). Figure 13 shows L_acc/L_bol as a function of the bolometric luminosity for our sample. We can see that, although the large uncertainty on the L_acc value, only up to four out of the ten Class I sources of our sample (i.e., ∼22%) might have their bolometric luminosity dominated by accretion.

Fig. 13. Lacc/Lbol ratio plotted as a function of the bolometric luminosity of our Class I sample. Low-veiled sources are plotted in magenta, and high-veiled ones are in black. Blue triangles are the upper limits for the Class I/Flat YSOs which show no HI detection. The horizontal dashed line represents the accretion luminosity fraction above which a source is in the main accretion phase.

This evidence suggests that, assuming a steady-state condition, the objects considered should have already accreted most of their mass. This can be tested from our estimated values for the mass accretion rates in these sources. Assuming a lifetime of τ_I = 0.54 Myr (Enoch & Evans 2009; Dunham et al. 2014) for the Class I phase and that the mass accretion rate remains constant during this time, we can easily find that the sources in our sample should have accreted only 4% − 35% of their mass during this phase. Consequently, either they have acquired most of their mass already during the Class 0 phase, or the accretion process proceeded in a non-steady framework.

This result is in agreement with the so-called luminosity problem. In an attempt to solve the luminosity problem, different models have been developed during the last years to try to overcome the limitations of the standard Shu collapse model. The main feature of such models is the assumption that the accretion is not a steady process, but it is characterised by variability of the mass accretion rate. In particular, Dunham & Vorobyov (2012) propose a collapse model in which the accretion declines with time but has also short-term variability and episodic bursts caused by disk gravitational instability and fragmentation. In their simulations, Dunham & Vorobyov (2012) foresee that the duration of the embedded phase depends on the final stellar mass produced, while the mass accretion rate, averaged over the full lifetime, is independent of the stellar mass. The range of mass accretion rates that we estimate for our sources, that is from 10⁻⁸ to 10⁻⁶ M_⊙ yr⁻¹, would be in agreement with those predicted by these models in quiescent periods; while during bursts, Ṁ_acc increases up to three orders of magnitude. Our derived values are also not in contradiction with the possibility that Ṁ_acc does not correlate with M_⋆, given the associated uncertainty on the estimated values. The fact that none of our sources are observed in bursts is expected, as their occurrence is estimated to be once every 10⁴ yrs during the embedded phase (Scholz et al. 2013).

A different model proposed to solve the luminosity problem is the one by Padoan et al. (2014). This model assumes that the accretion rate, during both the embedded and the Class II phases, is controlled mainly by the mass infall from the large-scale envelope and the turbulent cloud, and not by the disk, whose role is only to mediate accretion. Following the evolution of the infall rates during the formation of a protostar, Padoan et al. (2014) find that they are comparable with accretion rates inferred from protostellar luminosities or measured in pre-main sequence stars. Interestingly, in these simulations, large variations in the infall rates during the protostellar evolution are found due to the stochastic nature of turbulent flow and not to accretion bursts. Consequently, large spreads of accretion rates are predicted at any age and for any mass, with infall rate values again being consistent with the mass accretion rate measured in our sample.

6.3. Class I with no line detection

We have so far discussed the results found on the ten Class I sources of our sample for which we were able to give an estimate of the accretion properties.

Seven other sources do not show any detected line emission in their spectra, which suggests that they should be low accretors, as shown by the upper limits on the accretion luminosity derived in Sect. 5.3. In this section, we discuss the characteristics of these ‘candidate weak accretors’.

Table 3 shows that four of them are classified as Flat, with −0.25 ≤ α_IR ≤ 0.06, and three are classified as Class I, with 0.60 ≤ α_IR ≤ 1.56. Thus, the IR spectral index of this sub-sample does not have a peculiar trend with respect the Class I strong accretors.

Most of these sources have a SED (see Fig. F.2) similar to the other accreting Class I of the sample. However, the SEDs of sources #161 and #237 decrease at wavelengths longer than 2 μm and increase again at λ > 5 μm. This might suggest that the NIR source is not the counterpart of the Spitzer source.

Also, by inspection at the colour–colour diagram in Fig. 3, we note that all these sources lie in the same region of the diagram in which the Class II sample is located. In addition, because of the lack of IR excess in the K band (see Fig. B.2), in four of these sources, we are able to determine the spectral type and the stellar parameters in the same way we have done it for the Class II sample (see Sect. 4). Indeed, these sources lie in the same region of the Class II sample.

Moreover, because Class I protostars are defined by the condition that L_acc dominates their luminosity, in Fig. 13, we plotted the corresponding L_acc/L_bol upper limits for a comparison with the other Class I objects. The derived upper limits suggest that these sources have L_acc/L_bol more than two orders of magnitude lower than those of the Class I targets which show Brγ lines. This suggests that the observational classification based on SEDs’ spectral index value (Class 0, I, II, III) does not match with the corresponding evolutionary ‘Stage’ (e.g., Stage 0, I, II, III), which instead refers to a source true physical nature (Jayawardhana et al. 2001; Enoch & Evans 2009; Dunham et al. 2014). We conclude that this sample is composed of Class I YSOs in the evolutionary stage, Stage II.

Spectra with typical features of more evolved objects than Class I, such as photospheric absorption lines, were also found by Carney et al. (2016) for a sample of sources classified by Spitzer as YSOs, meaning that the spectroscopic follow-up is probably crucial to determine the physical nature of the evolutionary stage of YSOs.

Impact on the lifetime estimates of Class I/Flat sources

As we have seen in previous sections, the Class I sample contains strong accretors which show strong HI emission lines, and more evolved YSOs whose accretion, if any, is weak, hence not detected over the continuum. It is however instructive to highlight the results found on the global sample of embedded objects.

From the original sample of YSOs in NGC 1333, we observed 20 Class I/Flat sources, that is ∼55% of all the sources of these classes identified in NGC 1333 in Young et al. (2015). Of these 20 sources, ten are strong accretors, seven are YSOs with an evolutionary stage compatible with Class II, and three have been disregarded as YSOs on the basis of their spectra (see Appendix C). In conclusion, we derive that only 85%, 17 out of 20 of our investigated Class I sample, are accreting YSOs, and only 50% are strong accretors. Assuming that the same fraction applies to the sources of the Young et al. (2015) sample that we have not observed, we might conclude that the number of embedded accreting stars is smaller than estimated from the photometric information alone. This evidence can have implications for the estimated protostellar embedded phase lifetime. In fact, Evans et al. (2009) derive the lifetime of the Class I and Flat phases from the number of Class I/Flat relative to the Class II sources, assuming that the Class II lifetime is of the order of 2 ± 1 Myr. Our results show, however, that the contamination of non-accreting sources, likely more evolved objects, in the sample of embedded YSOs identified through Spitzer can be high, and consequently the lifetime of the Class I/Flat phase might be overestimated by up to a factor of two.

7. Summary and conclusions

In this paper we present the analysis of IR VLT/KMOS spectra of 52 young stellar objects in the < 1 Myr old NGC 1333 cluster in the Perseus star-forming region. Our goal is to investigate the evolution of accretion properties of YSOs from Class I to Class II.

The stellar (A_V, L_⋆, M_⋆, R_⋆) and accretion (L_acc, Ṁ_acc) parameters of the sample have been derived from the NIR spectra. For the Class II sample, we derived A_V and SpT through a comparison with a grid of non-accreting Class III stellar templates. We then computed the accretion luminosity by using the known correlation between L_acc and the luminosity of HI lines (Paβ and Brγ). For the Class I sample, where the presence of a large IR excess in the K band makes use of the same spectral typing method impossible, we used a procedure that allowed us to measure the stellar and accretion luminosity in a self-consistent way, based on a range of veiling values. Mass accretion rates Ṁ_acc were then measured once stellar masses and radii were estimated, adopting suitable evolutionary tracks. This work represents the first unbiased study of the accretion properties of 17 Class I YSOs belonging to the same cloud. In particular, we were able to measure the mass accretion rates of ten sources of our sample, where HI lines have been detected, while we set upper limits for the remaining seven sources.

We compared our results of the Class II NGC 1333 sub-sample with the samples of CCTs in the older Lupus and Chamaeleon I regions (Alcalá et al. 2017; Manara et al. 2017a).

The main results of this analysis are the following:

• The accretion luminosity of the Class II sample correlates with the stellar luminosity, as found for the other regions. The linear fit in the log plane is log L_acc = (1.5 ± 0.2)log L_⋆ + ( − 1.0 ± 0.1), which is in agreement within the errors with the relationship found on the Lupus sources and shallower with respect to the Cha I relation. The mass accretion rates of the Class II sample correlate with the stellar masses, following the relationship: log Ṁ_acc = (2.6 ± 0.9) log M_⋆ + (−7.3 ± 0.7), in agreement within the error with both Lupus and Cha I relationships. However, at variance with what is found for the log Ṁ_acc − log M_⋆ distribution in Lupus and Cha I, we do not find evidence for a double power-law trend in NGC 1333. If confirmed, this might suggest that the double-slope found in older regions is due to a faster evolution of the Ṁ_acc in the low mass regime.

• The accretion luminosity of the Class I sample (∼10⁻² − 10² L_⊙) is higher by up to two orders of magnitude with respect to values found on the Class II sample (∼10⁻⁴ − 10⁰ L_⊙). Similarly, the mass accretion rate of the Class I sample (∼10⁻⁹ − 10⁻⁶ M_⊙ yr⁻¹) is higher than the one found on the Class II stars up to two orders of magnitude (∼10⁻¹¹ − 10⁻⁷ M_⊙ yr⁻¹). We find that the higher accretion rates are associated with the sources with the larger IR veiling. Given the large uncertainty on the derived parameters, due to the assumptions needed by our method, and the low number of Class I targets, we are not able to define whether a correlation exists between the stellar parameters and L_acc or Ṁ_acc.

• We computed the L_acc/L_bol ratio for the Class I sample, finding that for the majority of the sources (13 out of 17), this ratio is smaller than 0.5. This indicates that for most of the sources, the bolometric luminosity is not dominated by the accretion luminosity and thus the sources are not in their main accretion phase. This result is in agreement with the well known luminosity problem and suggests that those objects have already accreted most of their mass or the accretion occurs in a non-steady fashion.

• We observed 55% of the Class I/Flat YSOs listed in the Young et al. (2015) catalogue of the NGC 1333 cluster. The spectroscopic follow-up of them shows that three are not YSOs and seven show no HI emission lines. These results highlight the importance of more extended spectroscopic surveys of YSOs in different star-forming clouds in order to confirm their classification and infer a more solid estimate of their lifetimes and of the star formation efficiency.

In conclusion, our analysis suggests that early stage YSOs have accretion rates on average higher than Class II stars of the same cloud. However, such values are still not sufficiently high to be responsible for the assembly of most of the stellar mass during the source lifetime, giving support to theories predicting that mass accretion does not occur in a steady fashion. To provide more solid conclusions, observations of larger samples of YSOs in young star-forming regions with sensitive IR instrumentation, as will be made possible with James Webb Space Telescope (JWST) in the near future, are highly demanded.

Acknowledgments

This work has been supported by the project PRIN-INAF-MAIN-STREAM 2017 “Protoplanetary disks seen through the eyes of new-generation instruments” and by the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 823823 (DUSTBUSTERS). We also acknowledge the support by INAF/Frontiera through the “Progetti Premiali” funding scheme of the Italian Ministry of Education, University, and Research and by the Deutsche Forschungs-Gemeinschaft (DFG, German Research Foundation) - Ref no. FOR 26341/1 TE 1024/1-1. E. Fiorellino was supported by an ESO studentship during the preparation of this work.

References

Appendix A: Technical details about observations

In this section, we report some technical details of each OB. In Table A.1, the name of each OB (Col. 1) and the date and time when the OB was observed (Col. 2) are reported. The exposure time on target is the product of the DIT time in seconds (Col. 3), times the number of DIT (NDIT, Col. 4), times the number of exposures in the A position (i.e., ‘on sky’, see Sect. 2, Col. 5). The airmass (Col. 6) and the seeing corrected for the airmass (Col. 7) are reported as an average value. Column 8 shows the number of KMOS arms allocated to science targets for each OB.

In the OB named 1bright, the following sources were observed: #184, #188, #192, #197, #200, #201, #203, #204, #205, #208, #209, #218, #219, #222, #226, #228, #231, #232, #233, #235, #237, and #241.

In the OB named 1dim, the following sources were observed: #180, #181, #186, #188, #190, #196, #200, #202, #203, #204, #205, #215, #218, #224, #227, #231, #232, and #237.

In the OB named 2bright, the following sources were observed: #152, #153, #156, #159, #160, #161, #162, #164, #166, #173, #174, #177, #178, #180, #186, #191, #194, and #196.

In the OB named 2dim, the following sources were observed: #147, #149, #153, #161, #165, #168, #174, #175, #183, #197, and #198.

In the OB named 3bright, the following sources were observed: #167, #168, #175, #176, #183, #195, #198, #199, #213, #221, #223, and #234.

In the OB named 3dim, the following sources were observed: #164, #194, #195, #199, #213, and #221.

We also added information about completeness for each class, showing the histograms of Class I (Fig. A.1) and Class II (Fig. A.2) sub-samples.

Fig. A.1. Histogram of the K-band magnitude of the Spitzer identified Class I in the NGC 1333 cluster (black, Young et al. 2015) and the sample observed by KMOS (red). The black-dashed vertical line represents the cut-off at 14.5 mag due to KMOS sensitivity.

Fig. A.2. Histogram of the K-band magnitude of the Spitzer identified Class II in the NGC 1333 cluster (black, Young et al. 2015) and the sample observed by KMOS (red). The black-dashed vertical line represents the cut-off at 14.5 mag due to KMOS sensitivity.

Appendix B: Observed spectra

We report here the complete JHK spectra and the individual spectra of Paβ and Brγ lines for our sample.

Figure B.1 shows the spectra of the Class I YSOs for which we detect at least one HI line in emission. Figure B.2 shows the spectra of the Class I YSOs for which no HI line is detected. Figures B.3–B.8 show the Class II YSOs spectra from early K to late M spectral type. We shifted the spectra in flux by a suitable amount to optimise the visualisation of the YSOs spectra in the figure.

Fig. B.1. J, H, and K band spectra of Class I YSOs which show at least one HI line. All the spectra are suitably shifted in flux for clarity.

Fig. B.3. J, H, and K band spectra of Class II YSOs with SpT from K7 to M2. All the spectra are suitably shifted in flux for clarity.

Fig. B.4. J, H, and K band spectra of Class II YSOs with SpT M3 and M3.5. All the spectra are suitably shifted in flux for clarity.

Fig. B.5. J, H, and K band spectra of Class II YSOs with SpT M4 and M4.5. All the spectra are suitably shifted in flux for clarity.

Fig. B.6. J, H, and K band spectra of Class II YSOs with SpT M5 and M5.5. All the spectra are suitably shifted in flux for clarity.

Fig. B.7. J, H, and K band spectra of Class II YSOs with SpT M6 and M6.5. All the spectra are suitably shifted in flux for clarity.

Fig. B.8. J, H, and K band spectra of Class II YSOs with SpT from M7 to M9.5. All the spectra are suitably shifted in flux for clarity.

Fig. B.2. J, H, and K band spectra of Class I YSOs which show no HI line. All the spectra are suitably shifted in flux for clarity.

Figure B.9 shows Paβ lines, while Fig. B.10 shows Brγ lines for Class II YSOs. Figures B.11 and B.12 show Paβ and Brγ lines, respectively, for Class I YSOs. Some sources were observed two times (dim and bright OBs, see Table A.1), but since there is no variability between them, we report only the spectra observed with a longer integration time.

Fig. B.9. Paβ lines of Class II objects. The flux is in erg s−1 cm−2 μm−1.

Fig. B.9. continued.

Fig. B.10. Brγ Class II lines. The flux is in erg s−1 cm−2 μm−1.

Fig. B.10. continued.

Fig. B.11. Paβ lines of Class I objects. The flux is in erg s−1 cm−2 μm−1.

Fig. B.12. Brγ Class I lines. The flux is in erg s−1 cm−2 μm−1.

Appendix C: Other sources

Here, we report some information about the eight sources for which an estimate of the accretion parameters was not possible on the bases of their spectra (Table C.1). We removed these objects from our sample for two reasons: the spectrum shows at least one HI line (Paβ and/or Brγ) in absorption, and/or it was not compatible with the one of YSOs.

In more detail, the spectra of sources #176 and #208 are not compatible with the YSOs ones. Source #176 shows HI lines in absorption and other absorption lines that are not compatible with the spectra of YSOs. Computing the extinction with the colour–colour diagram method (see Sect. 4), we find this source is not extinct (A_V = 0). Looking on the Simbad database for further information on #176, we found that its name is BD+30547, a G4V MS source (Frasca et al. 2018). Source #208 shows absorption lines typical of evolved stars, as is clearly visible in Fig. C.1.

Fig. C.1. J, H, and K band spectra of removed sources. All the spectra are suitably shifted in flux for clarity.

For what concerns sources #152, #160, #198, #223, and #228, they all show HI lines in absorption, so they cannot be analysed by using the empirical relationships of Alcalá et al. (2017). Source #223 is a non-accreting Class III YSO. For Class II objects (#160 and #198), we conclude they are non-accreting YSOs as well. For objects classified as Class I/Flat (see Table C.1), we checked the bolometric luminosity in the literature, finding L_bol = 1.9 L_⊙ and 0.07 L_⊙ for #152 and #228, respectively (Enoch & Evans 2009). Given the presence of HI in absorption and the low bolometric luminosity, we think these objects are non-accreting YSOs, suggesting they are not in the evolutionary stage expected for Class I YSOs. As a matter of fact, source #228 is spatially associated with ASR 7, a time-variable object exhibiting strong circular polarisation which suggests gyrosynchrotron emission (Rodríguez et al. 1999), despite the α_IR classifcation as a Class I YSO.

Appendix D: HR diagram of Class I

Figure D.1 shows the position in the HR-diagram of the ten Class I YSOs which show HI lines, in our assumption that the age of these sources is between the birthline and 1 Myr (see Sect. 5).

Fig. D.1. HR diagram of the Class I YSOs. Sources with veiling between 1 and 3 are plotted in magenta, while those with veiling between 3 and 8 are in green. The birthline (Palla & Stahler 1993) is plotted as a dot-dashed line. Black solid and dashed lines show the isochrones and evolutionary tracks by Siess et al. (2000).

From this diagram, we can appreciate how the stellar mass, given the luminosity derived from our method, changes by varying the age from the birthline to the 1 Myr. In particular, we note that sources #191 and #186 both remain in the intermediate-mass range if we put them on the birthline or in the 1 Myr isochrone, while source #200 has a larger uncertainty on the mass due to the kink of the birthline curve in the luminosity range of this source.

Appendix E: Comparison with Luhman et al. (2016)

Luhman et al. (2016) provide a census of YSOs and brown dwarfts (BDs) in the IC 348 and NGC 1333 clusters. In their work, the spectral characterisation of the sample is made through a comparison between the averages of dwarf and giant standards (Luhman 1999) spectra with the low-resolution spectra (R  =  100 − 800) of the sources taken with different spectrographs (SpeX@IRTF, IMACS@Magellan, GNIRS@Gemini, NIRC@Keck I). The extinction was estimated by dereddening the position of each source in the [J − H] vs. [I − K_s] diagram until it intersected the sequence of colours for dwarfs (Luhman et al. 2003).

We matched the sample analysed by Luhman et al. (2016) with ours, finding 40 sources in common. Of these, only 35 have both the extinction and the spectral type estimated by Luhman et al. (2016) (see Tables 4 and 6).

Figure E.1 shows a comparison between the SpT estimated by Luhman et al. (2016) and that derived in this paper. To be able to quantify the difference in the spectral type, we assigned a value in the following way to each spectral type: we set the M0 class as 0, then earlier types were associated to negative numbers with 0.1 steps (i.e., K7 corresponds to −0.1, K6 to −0.2, and so on), and later types to positive numbers (i.e., M1 corresponds to 0.1, M5 to 0.5, etc.). A general agreement is found between the two samples. As discussed in Sect. 4, the SpT and the extinction depend on each other. Therefore, in Fig. E.2, we plotted the variation in SpT and A_V between our results and those of Luhman et al. (2016). For only two sources, the discrepancy on A_V is > 5 and all the sources have ΔSpT contained within a spectral class. We also note that the most populated quadrant is the fourth (eight sources). This is a tentative trend, according to which, on average, Luhman et al. (2016) classification prefers higher extinction and early-type spectral types. This can be due to the low resolution of their observations which were less sensitive to the molecular features.

Fig. E.1. SpT of this work vs. the spectral type from Luhman et al. (2016) for the common sources.

Fig. E.2. Variation between this work and the Luhman et al. (2016) values of the extinction and the spectral type of our sample. We attributed an increasing number from early SpT to late SpT.

Appendix F: Class I SED’s and L_bol

Here, we present the spectral energy distributions of the Class I YSOs and the derivation of their bolometric luminosity.

To construct the SEDs, we used the IR photometry 2MASS and Spitzer of the Young et al. (2015) catalogue and, where possible, we added the measurements in FIR-submillimetre, from Herschel (Pezzuto et al. 2021), and millimetre from Bolocam (Enoch & Evans 2009). The derived SEDs are displayed in Fig. F.2.

Fig. F.1. Comparison between the bolometric luminosity computed by Evans et al. (2009) and the one computed in this work.

Fig. F.2. Spectral energy distributions (λFλ [erg s−1 cm−2] vs. λ [μm]) of the Class I YSOs. We use the 2MASS, Spitzer and Herschel, and Bolocam photometry (Young et al. 2015; Pezzuto et al. 2021).

In order to minimise the contamination from other stars in the relatively large Hershel FoV with respect to the Spitzer one, we inspected the Spitzer and Hershel 70 μm maps and disregarded the Herschel fluxes when a clear contamination is found. This is, for example, the case of #200.

We computed the bolometric luminosity (Table 5) by integrating the SEDs with an ad hoc IDL³ procedure, already used in Antoniucci et al. (2008). The integration was performed starting from the J-band, and considering straight lines in the Log(λ)−Log(λF_λ) plan relatives to available SED points; a final correction at the longest wavelengths was applied assuming that the emission decreases as 1/λ² after the last available observation. The derived values are shown in Tables 5 and 6, and they range between 0.06 L_⊙ and 58.8 L_⊙.

Figure F.1 shows the comparison between our estimates of L_bol and the ones computed by Evans et al. (2009) for the Class I sample. We find that for all the sources but #200 the bolometric luminosity we estimated is compatible or larger with the one estimated by Evans et al. (2009), as expected by adding Herschel fluxes in the FIR.

The results we show in this work are the first estimates of L_bol which include the information of Herschel photometry.

All Tables

All Figures

Fig. 1. 2MASS mosaic map of the NGC 1333 cluster achieved through the Image Mosaic Service (http://montage.ipac.caltech.edu/). Big magenta circles correspond to the KMOS pointings of 7.2′ diameter. Circles refer to the sample of YSOs with mK > 14.5 mag. The colour coding is the following: red (blue) and green (yellow) symbols indicate observed (not-observed) Class II/III and I/Flat sources, respectively.

In the text

	Fig. 2. Histogram of the K-band magnitude of the Spitzer identified YSOs in the NGC 1333 cluster (black, Young et al. 2015) and the sample observed by KMOS (red). The black-dashed vertical line represents the cut-off at 14.5 mag due to KMOS sensitivity.
In the text

Fig. 3. Colour–colour diagram of our sample. Class I, II, and III YSOs are plotted as green, red, and blue crosses, respectively. Triangles are lower limits. Class I YSOs which show HI in emission are surrounded by a black circle. The main sequence (MS, Luhman et al. 2010) and CTT locus (Meyer et al. 1997) are shown in black and purple, respectively. The green arrow in the left-top corner defines the reddening direction, following the Cardelli et al. (1989) extinction law. Its length corresponds to AV = 5 mag.

In the text

Fig. 4. HR diagram for the Class II YSOs of NGC 1333 (green stars), Lupus (orange crosses), and Cha I (purple asterisks). The green, orange, and purple lines correspond to the median of NGC 1333, Lupus, and Cha I samples (sub-luminous objects are not included), respectively. The birthline (Palla & Stahler 1993) is plotted as a dot-dashed line. Black solid and dashed lines show the isochrones and evolutionary tracks by Siess et al. (2000).

In the text

	Fig. 5. Example of the spectral typing result for Class II source #156. The observed spectrum is shown in red. The dereddened spectrum (black) matches the spectrum of an M6.5 Class III (green) from Manara et al. (2013). The difference in the K band is due to the fact that source #156 is accreting, while the Class III is not.
In the text

	Fig. 6. K-band spectrum of SVS13 A, normalised to the continuum. Red-dashed lines correspond to photospheric lines as labelled.
In the text

Fig. 7. Accretion luminosity as a function of stellar luminosity for the NGC 1333 Class II sample (green filled circles). Orange and purple stars represent Lupus and Cha I sources, respectively. Upper limits are displayed as triangles. Solid lines show the best fit for each star-forming region as specified in the legend. The standard deviation from the best fit is shown with dashed lines. The subluminous target is highlighted with a red box.

In the text

	Fig. 8. Histograms of the luminosity ratio Lacc/L⋆ of NGC 1333, Lupus, and Cha I samples in green, orange, and purple, respectively. The green-dashed line corresponds to the histogram of NGC 1333, including the upper limits.
In the text

Fig. 9. Mass accretion rate as a function of stellar mass for the NGC 1333 Class II sample (green). Detections are displayed as filled circles, while upper limits are shown as triangles. Orange and purple stars represent Lupus and Cha I sources, respectively. Solid lines show the best fit for each star-forming region as specified in the legend, and the standard deviation from the best fit is shown with dashed lines. The subluminous target is highlighted with a red box.

In the text

	Fig. 10. Histograms of the ratio Ṁacc/M⋆ for NGC 1333, Lupus, and Cha I samples in green, orange, and purple, respectively. The green-dashed line corresponds to the histogram, which also includes upper limits.
In the text

	Fig. 11. Accretion luminosity as a function of stellar luminosity for NGC 1333. The Class II sample is shown as in Fig. 7; Class I sources with veiling between 1 and 3 are plotted in magenta, while those with veiling between 3 and 8 are in black.
In the text

	Fig. 12. Mass accretion rate vs. stellar mass of the NGC 1333 cluster. Class II sample as in Fig. 9. Class I sources with veiling between 1 and 3 are plotted in magenta, and the ones with veiling between 3 and 8 are in black. The black dashed-dotted line represents a linear relation with a unitary slope.
In the text

	Fig. 13. Lacc/Lbol ratio plotted as a function of the bolometric luminosity of our Class I sample. Low-veiled sources are plotted in magenta, and high-veiled ones are in black. Blue triangles are the upper limits for the Class I/Flat YSOs which show no HI detection. The horizontal dashed line represents the accretion luminosity fraction above which a source is in the main accretion phase.
In the text

	Fig. A.1. Histogram of the K-band magnitude of the Spitzer identified Class I in the NGC 1333 cluster (black, Young et al. 2015) and the sample observed by KMOS (red). The black-dashed vertical line represents the cut-off at 14.5 mag due to KMOS sensitivity.
In the text

	Fig. A.2. Histogram of the K-band magnitude of the Spitzer identified Class II in the NGC 1333 cluster (black, Young et al. 2015) and the sample observed by KMOS (red). The black-dashed vertical line represents the cut-off at 14.5 mag due to KMOS sensitivity.
In the text

	Fig. B.1. J, H, and K band spectra of Class I YSOs which show at least one HI line. All the spectra are suitably shifted in flux for clarity.
In the text

	Fig. B.3. J, H, and K band spectra of Class II YSOs with SpT from K7 to M2. All the spectra are suitably shifted in flux for clarity.
In the text

	Fig. B.4. J, H, and K band spectra of Class II YSOs with SpT M3 and M3.5. All the spectra are suitably shifted in flux for clarity.
In the text

	Fig. B.5. J, H, and K band spectra of Class II YSOs with SpT M4 and M4.5. All the spectra are suitably shifted in flux for clarity.
In the text

	Fig. B.6. J, H, and K band spectra of Class II YSOs with SpT M5 and M5.5. All the spectra are suitably shifted in flux for clarity.
In the text

	Fig. B.7. J, H, and K band spectra of Class II YSOs with SpT M6 and M6.5. All the spectra are suitably shifted in flux for clarity.
In the text

	Fig. B.8. J, H, and K band spectra of Class II YSOs with SpT from M7 to M9.5. All the spectra are suitably shifted in flux for clarity.
In the text

	Fig. B.2. J, H, and K band spectra of Class I YSOs which show no HI line. All the spectra are suitably shifted in flux for clarity.
In the text

	Fig. B.9. Paβ lines of Class II objects. The flux is in erg s−1 cm−2 μm−1.
In the text

	Fig. B.9. continued.
In the text

	Fig. B.10. Brγ Class II lines. The flux is in erg s−1 cm−2 μm−1.
In the text

	Fig. B.10. continued.
In the text

	Fig. B.11. Paβ lines of Class I objects. The flux is in erg s−1 cm−2 μm−1.
In the text

	Fig. B.12. Brγ Class I lines. The flux is in erg s−1 cm−2 μm−1.
In the text

	Fig. C.1. J, H, and K band spectra of removed sources. All the spectra are suitably shifted in flux for clarity.
In the text

	Fig. D.1. HR diagram of the Class I YSOs. Sources with veiling between 1 and 3 are plotted in magenta, while those with veiling between 3 and 8 are in green. The birthline (Palla & Stahler 1993) is plotted as a dot-dashed line. Black solid and dashed lines show the isochrones and evolutionary tracks by Siess et al. (2000).
In the text

	Fig. E.1. SpT of this work vs. the spectral type from Luhman et al. (2016) for the common sources.
In the text

	Fig. E.2. Variation between this work and the Luhman et al. (2016) values of the extinction and the spectral type of our sample. We attributed an increasing number from early SpT to late SpT.
In the text

	Fig. F.1. Comparison between the bolometric luminosity computed by Evans et al. (2009) and the one computed in this work.
In the text

	Fig. F.2. Spectral energy distributions (λFλ [erg s−1 cm−2] vs. λ [μm]) of the Class I YSOs. We use the 2MASS, Spitzer and Herschel, and Bolocam photometry (Young et al. 2015; Pezzuto et al. 2021).
In the text