© ESO, 2014

1. Introduction

The disk accretion phase assumes a pivotal role in the scenarios of early stellar evolution and planetary formation. Circumstellar disks are the ubiquitous result of the earliest stages of star formation, surrounding the vast majority of solar-type protostars at an age of ~1 Myr. Mass accretion from the inner disk to the central star regulates the star-disk interaction over the few subsequent million years; this process has a direct, long-lasting impact on the basic properties of the system, by providing at the same time mass and angular momentum. Angular momentum transport and mass infall have a most important part in dictating the dynamical evolution of protoplanetary disks, hence setting the environmental conditions which eventually lead to the formation of planetary systems (Armitage 2011).

The widely accepted paradigm for mass accretion in young stellar objects (YSOs) is that of magnetospheric accretion (Camenzind 1990; Königl 1991). A cavity of a few stellar radii extends from the stellar surface to the inner disk rim, as an effect of the inner disk truncation by the stellar magnetosphere (B ~ 2 kG). Accretion columns, channeled along the magnetic field lines, regulate the transport of material from the disk inner edge to the central star; as the accreting material impacts onto the stellar surface at near free-fall velocity, accretion shocks are produced near the magnetic poles. This picture is strongly supported by its predictive capability of the main observational features of accreting young solar-type stars (classical T Tauri stars, or CTTS): UV, optical and IR excess compared to the photospheric flux; spectral “veiling”; broad emission lines; inverse P Cygni profiles; pronounced spectroscopic and photometric variability (e.g. Bouvier et al. 2007b).

Several studies (e.g. Muzerolle et al. 2003; Rigliaco et al. 2011) have shown that the rates (Ṁ_acc) at which mass accretion occurs tend to scale with the mass of the central object (M_∗). Notwithstanding this general trend, a large spread in the Ṁ_acc values at a given M_∗ appears to be a constant feature (e.g. Hartmann et al. 2006); this is indicative of a complex and dynamic picture, shaped by a variety of parameters and concurrent processes. Similarly, observations suggest a progressive decrease of Ṁ_acc over the age of the systems (Hartmann et al. 1998; Sicilia-Aguilar et al. 2010), following the time evolution and dispersal of circumstellar disks (Haisch et al. 2001; Hillenbrand 2005). The large spread associated with the effective relationship attests that a multiplicity of elements concur to the disk dissipation; moreover, the initial conditions in which stars are formed (e.g. Dullemond et al. 2006) and local environmental properties (e.g. Guarcello et al. 2007) are likely to significantly affect the evolutionary pattern of individual star-disk systems.

Statistical surveys of the accretion properties of young stars provide essential insight into the dynamics of accretion. Large, homogeneously characterized samples of objects, spanning several ranges of parameters, are crucial to achieve a proper understanding of the physics governing the disk accretion history of forming stars and to explore different facets of the accretion mechanisms. To this respect, extensive mapping surveys of well-populated star-forming regions at UV wavelengths (e.g. Rigliaco et al. 2011; Manara et al. 2012) provide a most interesting standpoint, as they allow to directly probe accretion for large samples of stars from the excess emission produced in the shocked impact layer at the base of the accretion column.

Such single epoch surveys efficiently trace the instantaneous spread of accretion properties in a given star-forming region; however, they are unable to probe the nature of such a spread, i.e. to determine whether this is due to variability of individual objects or to an intrinsic spread amongst objects. And yet, a key aspect of the accretion dynamics is the intrinsic variability profile of the process (i.e., the amplitude and characteristic timescales of the variability in Ṁ_acc). The enhanced variability typical of young stars, cornerstone of the very definition of the T Tauri class (Joy 1945), is manifest over a broad range of wavelengths (X-rays, UV, optical, IR) and of time baselines (hours, days, years); a variety of processes concur to shape the variability pattern of individual objects, with a significant contribution coming from the geometry of the system. Understanding to which extent and on which timescales the accretion process is intrinsically variable, and characterizing its variability on a statistical basis compared to individual cases, is of utmost importance in order to accurately interpret the information provided by large accretion surveys and thus infer a detailed picture of the accretion process.

Multi-epoch surveys hence provide a more accurate description of the accretion process, as they allow to probe the time evolution of single snapshots of the accretion properties and, thus, to investigate how accretion occurs. This has been addressed in a few recent studies (e.g. Nguyen et al. 2009; Fang et al. 2013), which explored the variability of spectral accretion diagnostics for tens of objects with a few observing epochs over several months. Such efforts provided very valuable constraints on the average accretion variability shown by typical classical T Tauri stars (CTTS). However, a much tighter time sampling is needed in order to recover an estimate of the effective intrinsic accretion variability, isolated from geometric contributions which are simply due to a varying visible portion of the accretion spots during stellar rotation.

We have recently conducted an extensive observational program specifically aimed at addressing the issue of YSOs variability over the full spectrum. The coordinated synoptic investigation of NGC 2264 (CSI 2264; Cody et al. 2013) project was devised as a coordinated multi-wavelength exploration of the hours-to-weeks variability of the pre-main sequence (PMS) population of the star-forming region NGC 2264. The space telescopes CoRoT and Spitzer provided the backbone of the optical and IR investigations, with simultaneous continuous monitoring over 40 and 30 days respectively; the Flames multi-object spectrograph at the Very Large Telescope (VLT) provided a full set of spectra covering 20 different epochs for ~100 young stars, while multi-band optical+UV photometric monitoring obtained at the Canada-France-Hawaii Telescope (CFHT) provided contemporaneous, synoptic measures of the accretion rates and their variability from the direct diagnostics of the UV excess, with an overall time coverage of 14 consecutive days and several measurements per night. The targeted region, NGC 2264, has long been a benchmark for star formation studies, thanks to its youth (~3 Myr), its richness (~1500 known members, both CTTS and non-accreting, weak-lined T Tauri stars, or WTTS), its relative proximity (d ~ 800 pc), the low extinction on the line of sight and the association with a molecular cloud complex that significantly reduces the contamination from background stars (see Dahm 2008 for a recent review on the region).

As part of the CSI 2264 campaign, this paper specifically focuses on the issue of accretion properties and accretion variability in the PMS population of NGC 2264. Section 2 describes the photometric dataset obtained at CFHT, on which the study reported here is based. Section 3 reports on the characterization of photometric accretion signatures, which allowed us to perform a new, accretion-driven census of NGC 2264 members; an extensive, homogeneous investigation of individual stellar properties is also reported. In Sect. 4, we describe the procedure adopted to convert the UV excess measurement to a Ṁ_acc estimate; we analyze the Ṁ_acc over M_∗ distribution thus inferred and compare this picture with the amount of variability registered over a baseline of two weeks; we further probe the origin of this variability and infer a specific estimate of the intrinsic accretion variability as opposed to rotational modulation detected during the monitoring. Results are discussed in Sect. 5; our conclusions are synthetized in Sect. 6. This provides a complete picture of a whole star-forming region and its several hundred members at short wavelengths, from a specifically accretion-oriented perspective. A more detailed characterization of the variability properties monitored at CFHT and their relation to the underlying physics will be addressed in a forthcoming paper (Venuti et al., in prep.).

2. Observations and data reduction

2.1. CFHT MegaCam surveys of NGC 2264: data reduction and photometric calibration

We surveyed NGC 2264 in two multiwavelength observing campaigns at the CFHT, using the wide-field optical camera MegaCam; the camera has a field of view (FOV) of ~1 deg² (Boulade et al. 2003), thus fitting the whole region in a single telescope pointing. Five broad-band filters (u^∗,g′,r′,i′,z′) are adopted at the camera; their design, close to the Sloan Digital Sky Survey (SDSS) photometric system (Fukugita et al. 1996), ensures high efficiency for faint object detection and deep sky mapping.

2.1.1. Mapping survey

The first MegaCam survey of NGC 2264, held on Dec. 12, 2010, consisted of a deep u^∗g′r′ mapping of the region. In the r′-band, five short exposures (10 s) were obtained, using a dithering pattern in order to compensate for the presence of bad pixels and small gaps on the CCD mosaic. In the u^∗-band and g′-band, the same 5-step dithering procedure was repeated in two different exposure modes, short (10 s) and long (60 s), in order to obtain a good signal/noise (S/N) for the whole sample of objects without saturating the brightest sources. All images were obtained over a continuous temporal sequence within the night. On Feb. 28, 2012, concurrently with the monitoring survey reported in Sect. 2.1.2, we additionally obtained a deep i′-band mapping, in order to reconstruct a full 4-band picture of NGC 2264 together with the u^∗g′r′ survey of Dec. 2010. Similarly to previous observations, a dithering pattern was used, along with two different exposure times: 180 s, to detect all sources in the FOV, and 5 s, to recover the bright sources saturated in the long exposure.

Raw data have been preliminarily processed at the CFHT, in order to correct for instrumental signatures across the whole CCD mosaic (i.e., bad pixel clusters or columns, dark current, bias, non-uniformity in the response), and a first astrometric and photometric calibration has been performed. The subsequent photometric processing of the images obtained during the survey has been handled as described below.

For each spectral band and each exposure time, single dithered images have been combined, in order to discard spurious signal detections (e.g., cosmic rays) and correct for small CCD gaps. The photometric images thus assembled (one in the r′-band and two for each other filter, corresponding to the two different exposure times) have then been processed using the SExtractor (Bertin & Arnouts 1996) and PSFEx (Bertin 2011) tools, in order to identify true sources in the FOV and extract the relevant photometry. At first, we run a preliminary aperture photometry procedure with a detection threshold of 10 σ above the local background, in order to extract only brighter sources; this selection of objects has been used as input to run a pre-PSF fitting routine, aimed at testing the optimal parameters for the PSF function to be used for analyzing the images. A final, complete single-band catalog has eventually been produced from each image processed with SExtractor by performing a new PSF-fitting run with a lower detection threshold of 3σ, thus recovering the fainter component of the population.

A first filtering of each catalog population, based on the properties of the extracted sources and on the relevant extraction flags, allowed us to discard the stellar component potentially affected by saturation or other detection problems. Long and short exposure catalogs corresponding to the same spectral band, when present, have then been merged, in order to obtain a unique photometric catalog for each spectral band and avoid duplicates. The sample of objects common to both the long-exposure and short-exposure catalogs has been used to evaluate the photometric accuracy, to correct the short-exposure photometry for small magnitude offsets with photometry obtained from the long exposures and to define a confidence magnitude threshold below which photometry could begin to be affected by saturation. Objects have been preferentially taken from the long-exposure catalog; sources brighter than the long-exposure threshold have been retrieved from the short-exposure catalog, when available; sources brighter than the short-exposure threshold have been discarded.

A complete calibration of the CFHT instrumental magnitudes to the SDSS photometric system, including a correction for both offset and color effect, has then been performed. The calibration procedure has not been anchored to a sample of standard stars in the FOV; instead, we used an ensemble of around 3000 stars, present in the CFHT FOV, for which SDSS photometry was available from the seventh SDSS Data Release (Abazajian et al. 2009), to statistically calibrate the CFHT magnitudes to the SDSS system. The matching of sources between the SDSS and the CFHT catalogs has been performed using the TOPCAT (Taylor 2005) tool, fixing a matching radius of 1 arcsec and retaining only closest matches. The statistical ensemble of stars used for the calibration has been selected among the objects with complete CFHT u^∗g′r′i′ photometry and that matched the magnitude limits for 95% detection repeatability for point sources defined for the SDSS photometric survey (Abazajian et al. 2004). The calibration procedure adopted is based on the results of the photometric calibration study of the CFHT/MegaCam photometry to the SDSS system performed by Regnault et al. (2009). We compared the distribution of SDSS colors over the ensemble of stars with the SDSS dwarf color loci tabulated in literature (Covey et al. 2007; Lenz et al. 1998) and retained, for calibration purposes, only objects whose colors matched the tabulated ranges. After checking their consistency with the observed trends, Regnault’s linear calibration relations have then been used to fit the (${\mathit{m}}_{\mathit{\lambda }}^{\mathrm{CFHT}}\mathrm{-}{\mathit{m}}_{\mathit{\lambda }}^{\mathrm{SDSS}}$) vs. (${\mathit{m}}_{\mathit{\lambda }}^{\mathrm{SDSS}}\mathrm{-}{\mathit{m}}_{\mathit{\lambda }\mathrm{+}\mathrm{1}}^{\mathrm{SDSS}}$) distributions in u^∗, g′ and r′, by keeping the angular coefficient fixed and adapting the intercept value. In the i′-band, for which the filter in use at MegaPrime has been changed, due to accidental breakage, after the epoch of acquisition of data used in Regnault et al. (2009)’s study, no evident color dependence of the offset between the two photometric systems has been observed; the i′-band photometry has therefore been uniquely corrected for a constant offset. The system of 4 calibration equations thus obtained (one for each spectral band; Eqs. (1)–(4)) has then been solved to obtain the final conversion from CFHT instrumental magnitudes to SDSS photometry.

$\begin{array}{ccc}& & {\mathit{u}}_{\mathrm{CFHT}}\mathrm{-}{\mathit{u}}_{\mathrm{SDSS}}\mathrm{=}\mathrm{-}\mathrm{0.211}\hspace{0.17em}\left({\mathit{u}}_{\mathrm{SDSS}}\mathrm{-}{\mathit{g}}_{\mathrm{SDSS}}\right)\mathrm{-}\mathrm{0.63}\\ & & {\mathit{g}}_{\mathrm{CFHT}}\mathrm{-}{\mathit{g}}_{\mathrm{SDSS}}\mathrm{=}\mathrm{-}\mathrm{0.155}\hspace{0.17em}\left({\mathit{g}}_{\mathrm{SDSS}}\mathrm{-}{\mathit{r}}_{\mathrm{SDSS}}\right)\mathrm{-}\mathrm{0.51}\\ & & {\mathit{r}}_{\mathrm{CFHT}}\mathrm{-}{\mathit{r}}_{\mathrm{SDSS}}\mathrm{=}\mathrm{-}\mathrm{0.03}\hspace{0.17em}\mathrm{(}{\mathit{r}}_{\mathrm{SDSS}}\mathrm{-}{\mathit{i}}_{\mathrm{SDSS}}\mathrm{)}\\ & & {\mathit{i}}_{\mathrm{CFHT}}\mathrm{-}{\mathit{i}}_{\mathrm{SDSS}}\mathrm{=}\mathrm{-}\mathrm{0.02.}\end{array}$

A complete u^∗g′r′i′ catalog of NGC 2264, containing ~9000 sources in the FOV, has been built by cross-correlating the single-band catalogs on TOPCAT. A matching radius of 1 arcsec has been fixed for identifying common sources detected in the u^∗,g′,r′ concurrent observations, while a larger matching radius of 2 arcsecs has been introduced to match u^∗g′r′ sources with their i′-band counterpart, to account for a lower astrometric accuracy (lower angular resolution) characterizing the i′-band field, due to poorer seeing conditions at the time of acquisition of the i′-band images. Sources have been identified by retaining only best matches; we address potential problems in the derived colors associated with misidentifications and multiple matches in Sect. 2.2. Final SDSS u′g′r′i′ (hereafter ugri) photometry has then been associated to each source in the catalog using the derived calibration relations.

2.1.2. Monitoring survey

The second MegaCam survey of NGC 2264, performed in Feb. 2012, consisted in monitoring the u^∗-band and r′-band variability of the entire stellar population of the region over a mid-term timescale (Feb.14-28, i.e., 2 weeks), as a part of the CSI 2264 project. On each observing night during the run, the region was imaged repeatedly, with a temporal cadence varying from 20 min to 1.5 h. Each u^∗r′ observing block was performed using a 5-step dithering pattern, with single exposures of 3 s in the r′-band and 60 s in the u^∗-band.

The procedure adopted for the photometric processing of the images obtained during the u^∗r′ monitoring survey is similar to that applied to the mapping survey images, but each exposure has been individually processed in order to retrieve the luminosity variations over different timescales. The i-band image (Sect. 2.1.1) has been used as the astrometric reference exposure, as the nebulosity component was the least conspicuous in this image compared to the other filters; for each source, an aperture has been placed down on each u^∗ and r′ image at the position predicted from the i-band source detection and the amount of flux within this aperture has been extracted to build the light curve.

A preliminary quality check on the CFHT light curves consisted in closely examining a sample of NGC 2264 members with known periodicities, in order to attempt to recover their periods from the phase-folded CFHT light curves. This allowed us to ascertain the global accuracy of the u^∗r′ time series photometry, albeit occasionally affected by isolated discrepant points and non-negligible scatter due to poorer observing conditions. In order to locate the potentially problematic observing sequences, we examined the time series photometry for the whole sample of field stars in the CFHT FOV and derived, for each observing sequence, the distribution of median zero-point offsets referred to the master frame. The zero-point offset, for a given object, is defined as the mag difference between the frame of interest and the master frame; in case of large deviation from zero, the zero-point offset distribution associated to an observing sequence thus signals that the observation has been carried out through clouds. This translates to significantly lower than average S/N for a given source, which implies that, even for the brightest stars, no accurate measurements can be inferred from that specific observation; hence, epochs matching this description have not been used for any subsequent analysis.

As done earlier for the ugri snapshot survey, the instrumental light curve photometry has been calibrated to the SDSS system, referring the procedure to the calibrated photometry from the first CFHT survey. No offsets and color effects have been noticed in the r′-band (hereafter r-band) photometry, while a correction for both contributions has been applied to the u^∗-band photometry from the (${\mathit{u}}_{\mathrm{CFHT}\mathrm{\_}\mathrm{lc}}^{\mathrm{\ast }}\mathrm{-}\mathit{u}$) vs. (u − r) plane: ${\mathit{u}}_{\mathrm{CFHT}\mathrm{\_}\mathrm{lc}}^{\mathrm{\ast }}\mathrm{-}\mathit{u}\mathrm{=}\mathrm{-}\mathrm{0.1343}\hspace{0.17em}\mathrm{(}\mathit{u}\mathrm{-}{\mathit{r}}^{\mathrm{)}}\mathrm{-}\mathrm{0.1009.}$(5)

2.2. The CFHT catalog of NGC 2264

Figure 1 shows a MegaCam color picture of the field imaged at CFHT during the NGC 2264 campaign. The star-forming region extends over the central part of the field, with the most active sites of star formation located at the center of the image, northward of the Cone Nebula (cf. Dahm 2008).

We obtained a complete ugri dataset and photometric monitoring for ~9000 sources in the field, exploring the area projected onto the star-forming region as well as the periphery of the cluster and background/foreground sky. Table 1 provides some details on the photometric properties of the population of the CFHT catalog. The survey is complete down to u ~ 21.5 and r ~ 18.5 and the monitored objects span a range of ~7 mag. The photometry is globally accurate up to a few × 10^-2 mag and the relative accuracy rises up to order of 10^-3 mag for brighter objects.

In order to evaluate the probability of misidentifications as a result of the procedure adopted for matching sources in different catalogs, and thus evaluate the potential impact on the photometric properties inferred for individual sources, we reconsidered all multiple identifications within the allowed matching radii and closely examined their spatial distribution. The sources with potentially multiple identifications amount to ~2.9% of the population of the catalog. The difference in offset distance between the best match and the next best match is less than 0.5 arcsec for only 23% of these cases. Less than 20% of the multiply matched sources are located within the region occupied by most NGC 2264 members. For the statistical purposes of this study, this component is assumed to be negligible.

3. The PMS population of NGC 2264

A global picture of the photometric properties of the stellar population can be achieved from the analysis of color–color diagrams, as shown in Figs. 2–3.

Empirical dwarf SDSS color sequences (Covey et al. 2007; Kraus & Hillenbrand 2007) have been overplotted to the diagrams in order to follow the observed color distribution and investigate the relationship between color loci and stellar properties.

In Fig. 2, earlier-type field stars are located along the diagonal ellipse, with foreground or slightly reddened objects on the lowermost part of the distribution and background stars at the uppermost end; M-type stars are predominantly located on the horizontal branch, whose g − r position is in good agreement with the value of ~1.4 mentioned in previous studies (Covey et al. 2007; Finlator et al. 2000; Fukugita et al. 1996; Annis et al. 2011). A number of objects can be located to the left (i.e., blueward) of the standard locus in r − i; only for a part of them this displacement may be ascribed to photometric inaccuracies. The most likely explanation for these objects is that they are highly variable, and the inferred r − i color is spurious due to the large epoch difference between the r-band and i-band photometry (the over one year time lapse between the ugr and the i-band mappings of the region, as discussed in Sect. 2).

Similarly, a progression of increasing spectral types and reddening effects can be observed from left to right along the ellipse in Fig. 3; an interesting locus is defined by the point dispersion at 0.5 ≤ g − r ≤ 2.0, below the main body of the color distribution, as it shows u − g color excesses presumably linked with the YSO accretion activity, manifest at UV wavelengths with a flux excess compared to the photospheric emission. An offset of a few tenths of a mag is observed between the mean locus of field stars in our survey and the sequences of Covey et al. (2007) and Kraus & Hillenbrand (2007). This offset cannot be explained with imprecisions in our calibration. A possible explanation may lie in a lower than solar standards metallicity of the field population probed here; at any rate, this disagreement has no impact on our analysis, as we uniquely used our internal reference sequence to probe the photometric accretion signatures of members, as illustrated in the following sections.

We further pursued the color signatures of different YSO types in the SDSS filters by analyzing the color properties of known members of the region with respect to the loci traced on the diagrams by the large population of field stars. Known members were identified amidst the population of the region based on a wide collection of photometric/spectroscopic data available in the literature from various surveys of the cluster; membership and further information on the WTTS (no disk evidence) vs. CTTS (disk-accreting) nature of selected members have been inferred based on one or more of the following criteria: i) H_α emitter from narrow-band photometry and variability from the data of Lamm et al. (2004) and following their criteria; ii) X-ray detection (Ramírez et al. 2004; Flaccomio et al. 2006) and location on the cluster sequence in the (I, R − I) diagram when R,I photometry is available; iii) H_α emitter from spectroscopy (H_αEW > 10 Å, H_α width at 10% intensity >270 km s^-1); iv) radial velocity member according to Fűrész et al. (2006); v) member according to Sung et al. (2008; cf. criteria enumerated earlier) and Sung et al. (2009; Spitzer Class I/Class II). It is worth remarking that this preliminary analysis of membership and TTS classification is not based on direct accretion criteria.

Around 700 known members have been matched in the CFHT catalog within a matching radius of 1 arcsec; among these, ~60% are classified as weak-lined T Tauri stars (WTTS). Figures 4 and 5 show the color properties of these groups of objects respectively in the (r − i, g − r) and (g − r, u − g) diagrams.

In redder filters, as shown in Fig. 4, both CTTS and WTTS define a single color sequence following the one traced by field stars. A few CTTS appear well blueward of this sequence, reflecting multi-year variability that is evidently more common in the CTTS than the WTTS. Conversely, accreting and non-accreting stars define two clear distinct distributions on the (g − r, u − g) diagram in Fig. 5, with the latter displaying colors consistent with the photometric properties of field stars and the former located at smaller (i.e., bluer) u − g values compared to the main color locus defined by dwarfs.

These distinctive features are well observed in Fig. 6, which shows how members distribute on the (u − r, r) color–magnitude diagram. There, the cluster sequence is clearly traced by WTTS, while CTTS appear broadly scattered to the left as a result of their u-band excess, indicative of active accretion activity, compared to WTTS. A number of presumed non-accreting members show a non-negligible displacement blueward/redward of the main locus; these two groups of objects are discussed further in Sect. 3.2.2.

Two main inferences can be drawn from these diagrams:

• 1.

  CTTS and WTTS members do differentiate in the SDSS colors on the basis of photometric accretion signatures at short wavelengths (UV flux excess);

• 2.

  the UV excess linked with accretion, detected in the u-band, does not affect significantly the observations in filters at longer wavelengths, as emphasized by the substantial agreement of the CTTS and WTTS color sequences in Fig. 4 and the color saturation at g − r ~1.4, a behavior common to field stars.

These photometric properties, observed for different ensembles of stars, enable us to perform a new membership and population study on an individual basis, from an accretion-driven perspective. This analysis has the twofold purposes of investigating the presence, in the CFHT FOV, of additional objects that are new candidate members, and of re-examining the classification of known members, looking in both cases for classification outliers in the diagram loci dominated by accretion.

3.1. UV excess vs. different accretion diagnostics

In order to ascertain the coherence of the classification based on the UV excess detection on these diagrams and define some confidence limits for a reliable identification of accreting members, we compared the u-band excess information drawn from CFHT photometry with different diagnostics commonly used to investigate accretion, such as the equivalent width (EW) and the width at 10% intensity in the Hα emission line, probing the accretion funnel, or the mid-infrared emission, that signals the presence of material in the inner disk when a flux excess compared to the photospheric level is detected.

Figure 7 shows the same color–color diagram as in Fig. 5 with additional information on the Hα EW values measured for the subsample of objects having VLT/Flames spectra from the CSI 2264 campaign (A. Sousa, UFMG). A threshold value commonly adopted to identify accreting YSOs from Hα emission is HαEW ≥ 10 Å (Herbig & Bell 1988), while a more robust scheme relating the threshold value to the stellar spectral type has been introduced by White & Basri (2003). As can be observed in Fig. 7, the UV excess and the Hα EW diagnostics are fully consistent, with the objects having the least Hα EWs located on the color locus traced by WTTS and field stars, while a region of the diagram dominated by accretion can clearly be delimited below the field stars distribution and identified from the detection of a UV excess.

A useful tool for characterizing presence and properties of a circumstellar disk is the examination of the infrared spectral energy distribution (SED) of the system. A particularly interesting indicator of the SED morphology is the measure of its slope in the spectral range of interest (i.e., the spectral index α first introduced by Lada 1987): smaller amounts of material in the inner disk result in smaller contributions, i.e. smaller flux excesses, at mid-IR wavelengths and hence more negative values of the α-index, which approaches the slope of a reddened blackbody, while positive values of the mid-IR SED slope are specific to deeply embedded, Class I sources. In Lada et al. (2006), the α-index diagnostics for disks is applied to the mid-IR emission sampling between 3.6 and 8 μm performed with the Spitzer instrument IRAC for the young cluster IC 348; the same tool has been adopted by Teixeira et al. (2012) to characterize the evolutionary state of NGC 2264 members based on the inner disk properties. We compare their results to the colors of NGC 2264 YSOs in the (g − r, u − g) diagram in Fig. 8. Following the classification proposed in Table 2 of Teixeira et al. (2012), values of α_IRAC ≲ −2 are indicative of objects with anaemic disks or naked photospheres, while objects with thick disks (hence likely to be actively accreting) are characterized by α_IRAC > −2. Indeed, as can be observed in Fig. 8, objects in the first α_IRAC group are mainly located on the WTTS/dwarf locus in u − g vs. g − r, while objects showing a UV excess (and thus actively accreting stars) have the largest values of α_IRAC, as globally expected.

3.2. UV census of the NGC 2264 population

3.2.1. New CTTS candidates in NGC 2264

A most interesting locus to identify accreting members is the data points dispersion below the main body of the distribution in Fig. 5. In fact, this allows easy detection and measurement of the u-band excess with respect to WTTS. The inferred color excess estimate is independent of distance and has no significant dependence on visual extinction A_V, as WTTS define a color sequence nearly parallel to the reddening vector down to ~M0 spectral type. For later-type stars, the color sequence traced by WTTS turns down, saturating at g − r ~ 1.4 (cf. Fig. 4); the implications of this for the UV excess measurement are discussed in Sect. 4. In order to minimize the contamination, we traced a straight line, parallel to the reddening vector, bordering the accretion-dominated area below the lower envelope of the dwarf color locus (Fig. 9); we then selected as new CTTS candidates all the objects lying in this region and marked as field stars. For all these objects, a preliminary check of the photometry on the CFHT field images has been performed, in order to discard the sources whose flux measurements could be affected by detection issues (e.g., too faint source, presence of two close sources, partial projection of the source on CCD gaps or bad pixel rows).

A second group of new CTTS candidates has been extracted from the color distribution in Fig. 4, namely among the color outliers located to the left and, to a smaller extent, below the main color distribution. Similarly to the procedure adopted for the selection of objects on the diagram in Fig. 5, a check on the photometry quality for all the objects of interest has preceded the identification of new candidates.

A third color–color diagram of interest for the analysis of membership is the one depicting u − g vs. r − i (Fig. 10); as a “transitional” diagram between u − g vs. g − r (Fig. 5; direct UV excess detection) and g − r vs. r − i (Fig. 4; nearly horizontal photospheric color sequence for late-type stars, with no color saturation on r − i), this diagram provides a greater sensitivity to the UV excess for spectral types ≥M0. Additional, later-type CTTS candidates have then been identified on this diagram, and their photometric properties inspected similarly to what done for the candidate members extracted earlier.

A second aspect of our selection of new candidate, actively accreting YSOs concerns the group of already known members which had previously been considered WTTS (that is, YSOs without evidence of significant, ongoing accretion at the time of the previous surveys). We examined the locations of these putative WTTS in our several diagrams (Figs. 4 through 9) and identified a significant number whose colors fell generally within the UV-bright regions – for example, the small group of blue points located below the normal dwarf locus in Fig. 5. This procedure allowed us to re-classify 19 WTTS members as CTTS, which add to the 50 newly accretion-identified candidates; this varied status does not necessarily reflect an erroneous previous classification, but may as well be the result of a long-term (~years) variability in the accretion activity of the objects.

The results of our UV-based membership investigation are provided in Table A.1 (new CTTS candidates) and Table A.2 (re-classified CTTS). This revised census of accreting members we infer for NGC 2264 likely encompasses all objects showing a strong UV excess, hence a significant amount of accretion. However, our sample might not be exhaustive regarding the weakly-accreting CTTS component. Indeed, a non-negligible overlap between the CTTS and WTTS populations can be seen in Figs. 5, 6, which cannot be efficiently probed with our diagnostics. Hence, we do identify some new CTTS, but we might be missing a number of lower accretors, showing a smaller UV excess that cannot be easily separated from the color properties of WTTS/field stars.

The wide FOV of CFHT/MegaCam allowed us to probe a quite extended region around the cluster, encompassing most of the FOVs adopted in previous surveys. It is thus of interest to examine the spatial distribution of the new CTTS candidates, selected as described in Sect. 3.2.1, compared to the distribution of known NGC 2264 members on the CFHT FOV. This comparison is shown in Fig. 11. A fraction of the object positions are projected onto the main cluster region, while a significant part are located on the periphery or beyond. This provides evidence for a wider spatial spread of the PMS population associated with NGC 2264, and attests to the presence of actively accreting members far from the original sites of star formation, currently known to be located towards the north and the center of the field imaged at CFHT (Dahm 2008).

3.2.2. Field stars contamination in the NGC 2264 sample

Based on the discussion in Sect. 3.2.1, another potentially interesting location for our accreting members census is the data point dispersion to the left of the cluster sequence (WTTS distribution) in Fig. 6. A direct search for candidate members in this color–magnitude locus would hardly be straightforward, due to the distance-dependence and to the vast preponderance of field stars; on the other hand, this diagram provides an additional test bench for investigating the nature of individual groups of objects.

Similarly to Fig. 5, we can identify on Fig. 6 a group of putative non-accreting members with inferred u − r colors bluer than expected. This feature, suggestive of ongoing accretion activity for a young cluster member, is consistently observed on different color-magnitude diagrams (g vs. u − g, i vs. u − i) for several of these objects; interestingly, these objects populate the thin lower branch of the WTTS distribution, below the main WTTS locus, on the lower part of the dwarf locus in Fig. 5 (u − g vs. g − r diagram). Little spectroscopic information is available for this group of stars; contrary to what is expected for objects clearly showing accretion activity, the level of variability exhibited by many of these objects during our monitoring is not significant over the noise level defined, as a function of magnitude, by field stars in our survey (see Venuti et al., in prep.). Additionally, several of these putative members would be located close to or onto the MS-turnoff (or below) in the HR diagram of the region (Fig. 12). We therefore suspect at least a part of these objects to possibly be misclassified dwarfs instead. A list of objects with questioned membership, of the relevant CFHT photometry and of the original selection criteria is provided in Table A.3.

Another interesting group of objects in Fig. 6 is the sample of WTTS displaced redward of the main cluster sequence. Most of these objects consistently appear significantly reddened through all CFHT photometric diagrams; in several cases, they display little variability. While contamination by field stars might be partly reflected in the observed properties, most of the objects in this stellar group would actually look like very young, or equivalently overluminous, objects on the HR diagram of the region (i.e., they are located on the upper right envelope of the data point distribution in Fig. 12). This would rather suggest binarity and/or earlier (i.e., more embedded) evolutionary status as possible explications of the distinctive properties of these objects.

3.3. NGC 2264 members in the CFHT survey: individual extinction and stellar parameters

The sample of PMS stars retrieved or newly identified in the CFHT survey amounts to about 750 objects. Based on the multiband photometry obtained during our campaign, as well as on the collection of previous literature data in different spectral bands, we performed an extensive investigation of individual stellar parameters (like extinction, mass and radius); this is aimed at characterizing the properties of members in a uniform and self-consistent manner.

Extinction towards NGC 2264 is known to be low (cf. Dahm 2008, for a review of literature results), with a typical value of ≲0.4 mag for the visual extinction A_V. Indeed, the good agreement of the cluster sequence with the color locus traced by foreground stars on Fig. 4 and the color saturation at ~1.4 in g − r (cf. discussion at the beginning of Sect. 3) is consistent with a low value of the average A_V. We use the photometric properties displayed in this diagram to derive an estimate of the individual A_V for late-type stars, i.e. the group of members located on the horizontal branch in the color locus. For this group of objects, A_V is determined by displacing each point along the reddening direction until the lowermost envelope of the branch distribution defined by field stars is reached and computing the corresponding amount of reddening. This procedure has the advantage of relying solely on observed quantities, without involving the comparison with an external reference sequence that would intrinsically introduce an additional source of uncertainty on the result.

For earlier-type stars (≲M0), the cluster sequence is nearly parallel to the reddening direction in Fig. 4, thus preventing a direct A_V determination from the object location on the diagram. For this group of objects, we derive a first estimate of the individual A_V from the analysis of the infrared photometry in the JHK_S filters, as measured in the 2MASS survey (Skrutskie et al. 2006). This consists in comparing their infrared colors with the JHK_S color sequence for dwarfs (Covey et al. 2007) and with the CTTS infrared locus defined in Meyer et al. (1997) and converted to the 2MASS photometric system as in Covey et al. (2010) (K. Covey, private communication). Additional A_V estimates for a part of these sources are available from previous, similar studies (e.g., Rebull et al. 2002, where A_V is investigated from the color excess on R − I, or Cauley et al. 2012, where A_V is computed by fitting stellar spectra with spectral templates) and/or derived from our optical photometry collection (from the measurement of the color excess on R − I, V − I or B − V in order of preference). In the absence of a direct A_V estimate from CFHT photometry, all available A_V derivations have been examined and their consistency checked on the color properties displayed in the CFHT diagrams, in order to discard discrepant estimates; an average value has then been adopted as a final estimate.

To evaluate the uncertainty on our individual A_V estimates, we collected all available A_V determinations for any object in our sample, and measured the amount of scatter observed among values inferred from different authors/methods. This comparison allowed us to conclude that different A_V estimates are typically consistent within a radius of a few 0.1 mag, which then provides an order-of-magnitude uncertainty on the A_V values adopted in this study.

Of the sample of CFHT members, spectral type is known for around 50% and has been retrieved, in order of preference, from the studies of Dahm & Simon (2005), Rebull et al. (2002) or Walker (1956). For the remaining objects, the spectral type has been derived from the comparison of the dereddened colors with the empirical optical color sequence of Covey et al. (2007). We used the dereddened photometry of members with known spectral type, earlier than M1, to recalibrate the empirical spectral type–color sequence onto CFHT photometry. Spectral types were converted to effective temperatures T_eff following the scale of Cohen & Kuhi (1979), checked for consistency against the scale of Luhman et al. (2003).

Bolometric luminosities have been derived from the dereddened J-band photometry retrieved from the 2MASS catalog, adopting T_eff–dependent J-band bolometric corrections (BC_J) obtained as a fit to the T_eff–BC_J scales of Pecaut & Mamajek (2013) and Bessell et al. (1998). A cluster distance value of 760 pc (Sung et al. 1997) has been adopted for the conversion to absolute magnitudes; this commonly adopted estimate of distance has been recently strengthened by the results of Gillen et al. (2014), who performed a detailed characterization of a newly identified, well sampled low-mass PMS eclipsing binary in NGC 2264, with strong evidence of membership and an inferred distance of 756 ± 96 pc. For a few tens of objects, for which J-band photometry from the 2MASS survey was not available, L_bol has been derived using an empirical calibration relationship between the dereddened r-band photometry and L_bol, which has been inferred from the distribution of dereddened r-band magnitudes vs. L_bol values derived from the J-band photometry over the whole sample.

Stellar masses have been determined by placing each object on the Hertzsprung-Russell diagram and comparing their position with the mass tracks of the PMS model grid of Siess et al. (2000; Fig. 12). A track-fitting tool provided by L. Siess on his webpage has been used to interpolate the mass value between two mass tracks. CFHT members span a wide range of masses, varying from 0.1 to 2 M_⊙. Stellar radii have been determined from the estimates of L_bol and T_eff as $\begin{array}{ccc}{\mathit{R}}_{\mathrm{\ast }}\mathrm{=}\sqrt{\frac{{\mathit{L}}_{\mathrm{bol}}}{\mathrm{4}\mathit{\pi \sigma }{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}}}& & \end{array}$where σ is the Stefan-Boltzmann constant.

Photometric data and stellar parameters for the whole population of NGC 2264 members monitored at CFHT are reported in Tables 2 and 3, respectively.

4. UV excess and accretion in NGC 2264

The investigation of YSOs at short wavelengths offers a unique window on the accretion mechanisms: it provides one of the most direct probes to the hot emission from the impact layer where the accreting material hits the stellar surface at near free-fall velocities, producing a shocked area up to several ×10³ K hotter than the stellar photosphere. This enhanced luminosity at UV wavelengths is reflected in the characteristic color properties of accretion-dominated objects compared to non-accreting young stars (cf. Figs. 5 and 6) and hence provides a direct proxy to the accretion luminosity, which in turn enables the investigation of the rates of mass accretion in individual systems.

The u-band flux excess is trivially given by ${\mathit{F}}_{\mathit{u}}^{\mathrm{exc}}\mathrm{=}{\mathit{F}}_{\mathit{u}}^{\mathrm{obs}}\mathrm{-}{\mathit{F}}_{\mathit{u}}^{\mathrm{phot}}\mathit{,}$(6)where ${\mathit{F}}_{\mathit{u}}^{\mathrm{obs}}$ is the measured flux in the u-band and ${\mathit{F}}_{\mathit{u}}^{\mathrm{phot}}$ is the amount of flux that would be expected from a purely photospheric emission. It is important to remark that, in order to isolate the component of additional flux produced in the accretion shock and thus probe accretion onto the star, the flux excess (compared to dwarfs) due to the enhanced chromospheric activity of T Tauri stars has to be included in the photospheric emission. This issue is addressed by taking the photometric properties observed for the WTTS population as the reference for the definition of the UV excess. The color excess will thus correspond to the difference between the observed stellar colors and the expected photospheric colors based on the properties of the non-accreting counterparts of the objects of interest: $\mathit{E}\mathrm{(}\mathit{u}\mathrm{-}\mathit{m}\mathrm{)}\mathrm{=}\mathrm{(}\mathit{u}\mathrm{-}\mathit{m}{\mathrm{)}}_{\mathrm{obs}}\mathrm{-}\mathrm{(}\mathit{u}\mathrm{-}\mathit{m}{\mathrm{)}}_{\mathrm{phot}\mathrm{+}\mathrm{chrom}}\mathit{.}$(7)As discussed in Sect. 3, no significant impact of the accretion luminosity is observed on the stellar flux detected at filters redder than the u-band; hence, the color excess of Eq. (7) basically corresponds to the u-band excess revealing the presence of accretion, E(u) ≃ E(u − m).

We measured the UV color excess of accreting stars in two different ways, from the properties displayed on the u − g vs. g − r and r vs. u − r diagrams (Figs. 5 and 6, respectively), as detailed below.

• In the first case (Fig. 5), we referred the color excess measurement to a straight line nearly parallel to the reddening vector and following the upper part of the WTTS distribution, as shown in Fig. 9. E(u) is thus measured as $\mathit{E}\mathrm{\left(}\mathit{u}\mathrm{\right)}\mathrm{=}\mathrm{(}\mathit{u}\mathrm{-}\mathit{g}{\mathrm{)}}_{\mathrm{obs}}\mathrm{-}\mathrm{(}\mathit{u}\mathrm{-}\mathit{g}{\mathrm{)}}_{\mathrm{ref}}\mathit{,}$(8)where (u − g)_ref is the reference (i.e., expected) color at the observed g − r level, (u − g)_ref = 1.312 (g − r)_obs + 1.03. The WTTS color sequence is in good agreement with the reference line until spectral type ~M0, and hence the adoption of this reference sequence allows us to derive a UV excess estimate essentially unaffected by reddening. For later-type stars, however, the real WTTS color sequence turns down on the diagram; this implies that the extrapolation of the WTTS trend at spectral types earlier than M0 to the whole sample yields an excess overestimation (i.e., more negative values of E(u)) for M-type stars. As the bulk of M-type WTTS lies about 0.2 mag below the reference line, we corrected for this bias by adding a constant offset of 0.2 mag to the E(u) estimates derived for this group of objects.

• In the second case, as can be observed on Fig. 6, WTTS nicely trace the cluster sequence over the whole range of magnitude. We thus dereddened the u,r photometry and defined the sequence of reference colors as the least-squares fit polynomial to the resulting WTTS distribution. The UV excess is thus provided by $\mathit{E}\mathrm{\left(}\mathit{u}\mathrm{\right)}\mathrm{=}\mathrm{(}\mathit{u}\mathrm{-}\mathit{r}{\mathrm{)}}_{\mathrm{obs}}\mathrm{-}\mathrm{(}\mathit{u}\mathrm{-}\mathit{r}{\mathrm{)}}_{\mathrm{ref}}\mathit{,}$(9)where (u − r)_ref is the reference color at brightness r_obs. Due to the small number of points and the large scatter affecting the definition of the sequence at the brighter end, we decided to restrict our UV excess analysis to objects fainter than 14.5 in r. This second method, compared to the first, has the disadvantage of being intrinsically more uncertain, as it is subject to errors in A_V determination (while the first method is essentially A_V-independent); on the other hand, the application of this method to u- and r-band observations taken during the monitoring survey allows us to derive both a global picture of accretion in NGC 2264 from the average photometry, and a variability range for the UV excess on a timescale of a few weeks (cf. Sect. 4.2).

An average rms error of ~0.16 mag has been associated to the UV excess determinations, in order to account for the scatter of the WTTS distribution around the reference sequence.

From Eq. (6) and $\mathit{E}\mathrm{\left(}\mathit{u}\mathrm{\right)}\mathrm{=}{\mathit{u}}_{\mathrm{obs}}\mathrm{-}{\mathit{u}}_{\mathrm{phot}}$(10)we obtain ${\mathit{F}}_{\mathit{u}}^{\mathrm{exc}}\mathrm{=}{\mathit{F}}_{\mathit{u}}^{\mathrm{obs}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{10}}^{\mathrm{+}\mathrm{0.4}\hspace{0.17em}\mathit{E}\mathrm{\left(}\mathit{u}\mathrm{\right)}}\mathrm{)}\mathrm{=}{\mathit{F}}_{\mathit{u}}^{\mathrm{0}}\hspace{0.17em}{\mathrm{10}}^{\mathrm{-}\mathrm{0.4}\hspace{0.17em}\mathit{u}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{10}}^{\mathrm{+}\mathrm{0.4}\hspace{0.17em}\mathit{E}\mathrm{\left(}\mathit{u}\mathrm{\right)}}\mathrm{)}\mathit{,}$(11)where u is the (dereddened) apparent u-band magnitude of the object, E(u) is the UV excess computed as in Eqs. (8)–(9) and ${\mathit{F}}_{\mathit{u}}^{\mathrm{0}}$ is the SDSS u-band zero-point flux (${\mathit{F}}_{\mathit{u}}^{\mathrm{0}}$ = 3767.2 Jy).

We used Eq. (11) (cf. also Rigliaco et al. 2011) to measure the u-band flux excess over the whole sample and then converted ${\mathit{F}}_{\mathit{u}}^{\mathrm{exc}}$ to ${\mathit{L}}_{\mathit{u}}^{\mathrm{exc}}$ using a distance of 760 pc.

The results of a similar, single epoch photometric survey to detect and measure mass accretion in NGC 2264 from the UV excess diagnostics have been reported in Rebull et al. (2002). The authors measured the UV excess displayed by accreting stars from U-band and V-band photometry, by comparing the observed colors with colors expected based on the spectral type of individual objects. In Fig. 13, we compare their measured UV excess luminosities (L. Rebull, priv. comm.) with both the median ${\mathit{L}}_{\mathit{u}}^{\mathrm{exc}}$ values and the relevant variability ranges we detect here, on a timescale of a couple of weeks, for CTTS targeted in both surveys (about 100). As can be observed, data points result to be well distributed around the equality line; the rms scatter measured about the line, amounting to 0.5 dex, is quite consistent with the average amount of variability on ${\mathit{L}}_{\mathit{u}}^{\mathrm{exc}}$ (about 0.6 dex) detected across the sample from CFHT monitoring.

4.1. Accretion luminosity and mass accretion rates

The u-band excess luminosity represents a single component of the total excess emission over the whole spectral range. In order to characterize the properties of the accretion process, it is therefore necessary to also take the flux excess at unobserved wavelengths into account, i.e., to apply a bolometric correction.

Gullbring et al. (1998) showed that a direct correlation exists between the total accretion luminosity L_acc and the excess luminosity observed in the U-band. This is indeed expected, since the excess emission likely dominates the observed flux at wavelengths <4500 Å (i.e., below the g-band filter). In order to achieve this result, they measured the excess emission from spectra covering the range 3200–5400 Å, for a sample of 26 CTTS. They modeled the excess spectrum from the accretion spot assuming a region of constant temperature and density (slab model) and used this model fit to estimate the amount of excess emission at unobserved wavelengths. The authors eventually resolved to adopt a uniform correction factor, across the whole sample, for the ratio of total-to-observed excess flux, corresponding to a T ~ 10 000 K. For each object, they subsequently synthetized a U-band magnitude from the observed spectrum and computed the corresponding flux excess. The comparison of the two quantities (L_acc and ${\mathit{L}}_{\mathit{U}}^{\mathrm{exc}}$) enabled the authors to yield a well-determined calibration equation, $\log \mathrm{(}{\mathit{L}}_{\mathrm{acc}}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}\mathrm{)}\mathrm{=}\mathrm{1.0}{\mathrm{9}}_{-0.18}^{\mathrm{+}\mathrm{0.04}}\hspace{0.17em}\log \mathrm{(}{\mathit{L}}_{\mathit{U}}^{\mathrm{exc}}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}\mathrm{)}\mathrm{+}\mathrm{0.9}{\mathrm{8}}_{-0.07}^{\mathrm{+}\mathrm{0.02}}$, that was ultimately consistent with the model predictions of Calvet & Gullbring (1998), who introduced a detailed description of the accretion shock region.

Here we adopt a simple phenomenological approach to derive a direct relationship between the u-band excess luminosity and the total accretion luminosity L_acc. In a first approximation, we describe the stellar emission as a blackbody emission at the photospheric temperature T_phot. Similarly, the emission from the accretion spot is approximated by a blackbody at T = T_spot. The excess emission linked with accretion is thus given by the difference between the flux emitted by the distribution of accretion spots and the flux emitted by a corresponding area, of the same extension, on the stellar photosphere. In terms of blackbodies, this becomes the difference between a blackbody flux density at T = T_spot and a blackbody flux density at T = T_phot, both integrated over the same area. In this picture, the ratio r of total-to-observed flux excess in the u-band, i.e. the corrective factor to be introduced to derive the total L_acc from ${\mathit{L}}_{\mathit{u}}^{\mathrm{exc}}$, can be estimated as $\mathit{r}\mathrm{=}\frac{\mathit{\sigma }{\mathit{T}}_{\mathrm{spot}}^{\mathrm{4}}\mathrm{-}\mathit{\sigma }{\mathit{T}}_{\mathrm{phot}}^{\mathrm{4}}}{\mathit{\pi }{\mathrm{\int }}_{\mathit{u}}\mathrm{\left(}{\mathit{B}}_{{\mathit{T}}_{\mathrm{spot}}}\mathrm{\right(}\mathit{\lambda }\mathrm{)}\mathrm{-}{\mathit{B}}_{{\mathit{T}}_{\mathrm{phot}}}\mathrm{(}\mathit{\lambda }\mathrm{\left)}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\lambda }}\mathit{,}$(12)where the integration is performed over the u-band window. No assumptions on the filling factor are needed in this description.

In order to investigate the relationship between ${\mathit{L}}_{\mathit{u}}^{\mathrm{exc}}$ and L_acc, we selected a subsample of 44 CTTS, covering the spectral type range M3.5:K2 (or equivalently the mass range ~0.2–1.8 M_⊙) and whose variability, monitored in the u-band (λ_eff = 3557 Å) and r-band (λ_eff = 6261 Å) at CFHT, is likely accretion-dominated. An extensive spot modeling (cf. Bouvier et al. 1993) of the simultaneous u + r variability amplitudes has been performed in order to assess the dominating features of the variability displayed by individual objects; accretion-dominated sources have been identified as well-represented variables in terms of a surface spot distribution ~10³ to several × 10³ K hotter than the stellar photosphere. A full description of the variability analysis for the CFHT NGC 2264 members will be provided in a forthcoming paper (Venuti et al., in prep.).

For each of the selected objects, we adopted the individual T_spot estimate inferred from spot models, that corresponds, for a given object, to the color temperature of the spot distribution that best reproduces the simultaneous flux variations observed at different wavelengths; T_phot is derived from the spectral type as described in Sect. 3.3. We then derived, for each object, the ratio r of Eq. (12), through a numerical integration of B(λ) in the spectral range 3257–3857 Å (cf. Fukugita et al. 1996); an uncertainty on the value of r has been inferred correspondingly to the uncertainty on the spot temperature T_spot. Based on these estimates, the percentage of the total excess flux detectable in the u-band is on average ~10%. The accretion luminosity L_acc is then computed as ${\mathit{L}}_{\mathrm{acc}}\mathrm{=}\mathit{r}\hspace{0.17em}{\mathit{L}}_{\mathit{u}}^{\mathrm{exc}}$.

A least-squares fit to the $\log \left({\mathit{L}}_{\mathit{u}}^{\mathrm{exc}}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}\right)$ vs. log (L_acc/L_⊙) distribution inferred as previously described provided us with the following calibration relationship: $\log \left(\frac{{\mathit{L}}_{\mathrm{acc}}}{{\mathit{L}}_{\mathrm{\odot }}}\right)\mathrm{=}\mathrm{(}\mathrm{0.97}\mathrm{\pm }\mathrm{0.03}\mathrm{)}\hspace{0.17em}\log \left(\frac{{\mathit{L}}_{\mathit{u}}^{\mathrm{exc}}}{{\mathit{L}}_{\mathrm{\odot }}}\right)\mathrm{+}\mathrm{(}\mathrm{1.09}\mathrm{\pm }\mathrm{0.07}\mathrm{)}\mathit{.}$(13)Figure 14 compares our result with the calibration inferred from the study of Gullbring et al. (1998); the parameters of the two calibrations are fairly consistent within the error bars.

The locus traced by our point distribution is globally consistent with the prediction of Gullbring et al.’s study, albeit showing a systematic shift above their calibration line. Our estimates result to be larger by a factor of ≲2 to 4, corresponding to a difference from ~0.2 dex at the highest Ṁ_acc regimes (~10^-7 M_⊙/yr) to ~0.6 dex at the lowest Ṁ_acc regimes (~10^-10 M_⊙/yr).

This difference between the two conversions is significantly reduced if we adopt, in our approach, a constant spot temperature of 10 000 K, an assumption closer to the uniform correction, at T ~ 10 000 K, that was introduced in Gullbring et al.’s method. Indeed, individual T_spot estimates adopted in our study cover a wide range of values, with a mean around ~6500 K, a spread of ~2000–3000 K at a given stellar temperature and nominal T_spot that can reach up to ~9000–11 000 K and down to ~5000 K. As described earlier, spot temperatures, in our approach, are inferred case by case by fitting simultaneous optical (r-band) and UV (u-band) variability amplitudes with spot models. Hence, these correspond to the color temperatures that best reproduce the observed, λ-dependent contrast of the spot distribution with the stellar reference emission level (which is the effect we aim at accounting for when introducing a bolometric correction). As the peak of the emission of a blackbody in the u-band occurs around 10 000 K, lower T_spot estimates naturally imply that a larger correction factor is needed in order to recover the total flux emitted by the spot. While, for a typical photosphere (T ~ 4000 K), 12% of the excess emission of a 10 000 K spot is observed in the u-band, this percentage reduces to ~10% for a 7000 K spot and to less than 8% for a 6000 K spot. In order to evaluate how much of the discrepancy could actually be due to this effect, we repeated our procedure assuming a constant T_spot = 10 000 K over the whole subsample of objects; as shown in Fig. 14, this indeed would allow us to mostly cancel the offset between the two calibrations.

We are aware of the simplified nature of the blackbody approximation for stellar and spot fluxes; at the same time, we stress that the results inferred from this description are not inconsistent with those obtained from more sophisticated pictures. The discrepancy, albeit systematic, between L_acc estimates from Eq. (13) and those from Gullbring et al. (1998)’s calibration is smaller than, or comparable to, the usual estimates of the uncertainty affecting the accretion rate measurements (~0.5 dex), and a non-negligible contribution to this offset may actually derive from the different assumption on the correction factors adopted across the samples of objects considered in each study. Figure 10 of the theoretical study performed by Calvet & Gullbring (1998), though it provides full endorsement of Gullbring et al.’s relationship, yet depicts the complexity of achieving a detailed physical representation of the accretion shock. The range of predictions from model configurations they explored would lie, on Fig. 14 of our work, across the region delimited by our (Eq. (13)) and Gullbring et al. (1998)’s relationship.

Considering this, and the fact that each model has its own limitations, we believe it beneficial to explore whether, and to which extent, information deduced from data might actually be sensitive to different model assumptions. We therefore adopt the calibration in Eq. (13) to estimate the total L_acc from the measured u-band luminosity excess, and will refer at the same time to Gullbring et al.’s, where relevant, for comparison purposes.

If we assume that the accretion energy is reprocessed entirely into the accretion continuum (cf. Herczeg & Hillenbrand 2008), the mass accretion rate Ṁ_acc can be derived from L_acc as $\mathit{Ṁ}\mathrm{acc}\mathrm{=}{\left(\mathrm{1}\mathrm{-}\frac{{\mathit{R}}_{\mathrm{\ast }}}{{\mathit{R}}_{\mathrm{in}}}\right)}^{-1}\frac{{\mathit{L}}_{\mathrm{acc}}{\mathit{R}}_{\mathrm{\ast }}}{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}\mathrm{~}\mathrm{1.25}\hspace{0.17em}\frac{{\mathit{L}}_{\mathrm{acc}}{\mathit{R}}_{\mathrm{\ast }}}{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}\mathit{,}$(14)where M_∗ and R_∗ are the parameters of the central star, R_in is the inner disk radius and we are assuming that the accretion funnel originates at the truncation radius of the disk, at R_in ~ 5 R_∗ (Gullbring et al. 1998).

Figure 15 compares the Ṁ_acc distributions inferred from each of the two UV derivations described at the beginning of Sect. 4. The diagram reports the Ṁ_acc determinations for every accreting object, in the sample, that displays a UV excess according to the identification procedure previously described. The two distributions are on average consistent, within a small mass-dependent offset likely reflecting a residual overestimation of the (u − g) color excess from Fig. 9 (see discussion around Eq. (8)). Variability is likely to affect, at least partly, the scatter observed within ±0.25 dex around the fitting line (see Sect. 4.2).

Figure 16 shows the mass accretion rate distribution inferred for a population of 236 high-confidence NGC 2264 members, as a function of stellar mass, whose range extends down to ~0.1 M_⊙ and up to ~1.5 M_⊙. These objects have been classified as CTTS, i.e., accreting objects, as opposed to WTTS among the PMS population analyzed so far, based on spectroscopic criteria (H_α EW, H_α width at 10% intensity) and/or photometric criteria (significant UV excess from the CFHT sample, α_IRAC value consistent with Class II sources from the data of Teixeira et al. 2012, large variability from the CFHT photometry; see Venuti et al., in prep., for this last point). The accretion parameters for these objects are reported in Table 4. The nominal values of the accretion rates are computed from the photometry corresponding to the median UV excess detected during the CFHT/MegaCam 2-week long monitoring and are estimated to be accurate within a factor of ~3 typically.

A confidence detection threshold for Ṁ_acc has been computed by evaluating the average scatter of the WTTS around the reference color sequence and testing the corresponding level of “noise” on the UV excess – and consequently accretion rate – determination. Namely, the magnitude-dependent rms scatter of WTTS around the reference sequence has been measured to compute a (magnitude-dependent) minimum UV excess level to be detected in order to be confident that accretion is indeed the dominant source of excess. This allowed us to trace a conservative boundary separating the objects with an unambiguous Ṁ_acc detection from the objects for which it is not possible, a priori, to evaluate the relevance of spurious contributions (i.e., chromospheric activity) to the apparent UV excess. For the accreting objects for which the actually measured UV excess is less than the corresponding boundary value, only an upper limit at the detection threshold has been reported for Ṁ_acc.

Little extensive investigation of mass accretion rates in NGC 2264 has been reported in the literature prior to the survey we are currently performing. Here we refer again to the study of Rebull et al. (2002). Those authors reported Ṁ_acc estimates for about 75 CTTS in the region; their accretion rates are derived from the measured UV excess luminosities by adopting Gullbring et al. (1998)’s relationship to obtain the total L_acc values. Results of the comparison of our Ṁ_acc estimates with theirs, for the few tens of sources in common (that cover nonetheless about 2 orders of magnitude in Ṁ_acc), are quite consistent with what deduced from the comparison of the measured ${\mathit{L}}_{\mathrm{UV}}^{\mathrm{exc}}$ (Fig. 13). When using the same ${\mathit{L}}_{\mathit{u}}^{\mathrm{exc}}$–L_acc calibration, Ṁ_acc estimates from the two studies distribute around the equality line, with an rms scatter of about 0.5 dex, pretty consistent with the average amount of Ṁ_acc variability we detect on week timescales during our monitoring (see Sect. 4.2).

Some important features can be observed on Fig. 16. The Ṁ_acc distribution globally traces a mass-dependent trend, showing an average increase of the mass accretion rates with stellar mass. The existence of such a trend has been investigated from observations over the last decade, albeit with a large uncertainty on the quantitative estimate of the actual relationship, which is likely to be dependent, among other factors, upon the probed stellar mass range and the mean age and age spread across the cluster (cf., e.g., Hartmann et al. 2006 and references therein, Rigliaco et al. 2011, Manara et al. 2012). Scattering of objects at a given stellar mass (see discussion later in the Section) also contributes to increasing the uncertainty on the actual relationship. Another important point to this respect is to understand the actual role played by detection biases in the observed mass-Ṁ_acc trend. As previously discussed, the statistical sample of accretion rates is limited at the lower Ṁ_acc end (Fig. 16), due to the impossibility of an accurate detection below a certain threshold; for this group of objects, only an upper limit to the actual Ṁ_acc value can be inferred. Contrast effects favor the detection of smaller accretion rates on cooler, i.e. lower-mass stars over higher-mass stars (as can be seen in Fig. 16), thus unavoidably introducing a mass dependence. It is thus of interest to assess the degree of disentanglement between the observed mass-Ṁ_acc dependence and the effects of censoring in the data.

A statistical approach that allows us to infer a robust measure of correlation in a given sample, in presence of multiple censoring (both ties, here points at the same x-value but with different y-values, and non-detections, i.e. upper limits), is a generalization of Kendall’s τ test for correlation (see Feigelson & Babu 2012; Helsel 2012). The test consists in pairing the data points (measurements and upper limits) in all possible configurations and, for each pair, measuring the slope (positive or negative) of the line connecting the two points. A pair of data points having the same x-value will be counted as an indeterminate relationship, likewise a pair of upper limits; upper limits will uniquely contribute a definite relationship when paired with actual detections at higher y-value. The number of times positive slopes occur is then compared to the number of negative slopes, and their difference weighted against the total number of available pairs. This defines the statistical correlation coefficient, τ. ln this picture, upper limits have the effect of lowering the likelihood of a correlation, if this is present in the data, as they introduce a number of indeterminate relationships among the tested pairs. Conversely, no correlation will be found in a pure sample of upper limits, even though these trace a well-defined trend, since the number of pairs with definite relationship will be zero. In the null hypothesis of no correlation, τ is expected to follow a normal distribution centered around zero; by measuring the deviation from zero of the actually found value, in terms of σ, it is thus possible to assess a confidence level for the presence of correlation in a given sample.

The application of this statistical tool to our data allowed us to establish the presence of a correlation (τ ~ 0.28) that appears to hold across the whole mass range studied here with a confidence of more than 6σ. The robustness of this result was tested against a similar search for correlation in randomly generated distributions of Ṁ_acc in the log range [− 10.5 , − 6.5] for the same stellar population with the same individual detection thresholds. The τ distribution inferred from 100 iterations is peaked around τ = 0.00 (no correlation) with a σ of 0.04; this suggests that indeed an intrinsic correlation is observed in our data to a significance of >6σ. We then inferred a statistical estimate for the slope following the approach of Akritas-Theil-Sen nonparametric regression (see Feigelson & Babu 2012):

• 1.

  we derived a first guess for the slope from a least-squares fit to the data distribution and explored a range of ±0.5 around this first estimate with a step of 0.005;

• 2.

  for each value x, we subtracted the xM_∗ trend from the initial Ṁ_acc distribution and repeated the Kendall’s τ computation procedure on the distribution of residues thus inferred.

The best fitting trend is the one that, subtracted to the point distribution, produces τ = 0, while an error bar can be inferred corresponding to the interval of values that produce τ< |σ|. We applied this procedure to both the mass-Ṁ_acc distribution inferred from the u − r excess distribution (see discussion around Eq. (9)) and the one inferred from the u − g excess distribution (see discussion around Eq. (8)), obtaining respectively: $\begin{array}{ccc}\mathit{Ṁ}\mathrm{acc}& \mathrm{\propto }& {\mathit{M}}_{\mathrm{\ast }}^{\mathrm{1.5}\hspace{0.17em}\mathrm{\pm }\hspace{0.17em}\mathrm{0.2}}\\ \mathit{Ṁ}\mathrm{acc}& \mathrm{\propto }& {\mathit{M}}_{\mathrm{\ast }}^{\mathrm{1.3}\hspace{0.17em}\mathrm{\pm }\hspace{0.17em}\mathrm{0.2}}\mathit{.}\end{array}$The two trends are consistent within 1σ; the smaller exponent in Eq. (16) compared to Eq. (15) may reflect the small, mass-dependent offset between the Ṁ_acc distributions inferred from the two methods (cf. Fig. 15). A slope of 1.7 ± 0.2, instead of the estimate in Eq. (15), would be inferred if we converted u − r excesses to L_acc values following Gullbring et al. (1998). The two estimates are still consistent within 1σ; this comparison illustrates that a (small) range in values around the true slope may be induced by external factors such as the ${\mathit{L}}_{\mathit{u}}^{\mathrm{exc}}$–L_acc conversion.

The analysis of correlation described above refers to the homogeneous sample of Fig. 16 as a whole. While the relationship we infer seems to be robust across the whole sample, this likely provides a general, average information on the global trend of Ṁ_acc with mass, not a punctual description of the trend. We notice, for instance, that little mass dependence is observed in the upper envelope of the Ṁ_acc distribution above M_∗ ~ 0.25 M_⊙. We attempted a more detailed investigation of the observed Ṁ_acc vs. M_∗ relationship by probing separately two different, similarly populated mass bins, M_∗ ≤ 0.3 M_⊙ and M_∗ > 0.3 M_⊙. In both subgroups, we found some evidence of correlation, albeit with different values and significance levels: a correlation with a slope of 3.3 ± 0.7 was inferred, to the 4σ level, in the lower-mass group, while a correlation with a slope of 0.9 ± 0.4 was obtained for the higher-mass group with a significance of 2.4σ. These results might suggest that the dependence of Ṁ_acc on M_∗ is not uniform across the whole M_∗–Ṁ_acc range, but shallower at higher masses and higher accretion rates; on the other hand, a closer investigation of the M_∗–Ṁ_acc relationship from separate, small subsamples is necessarily more uncertain than the assessment of a general correlation trend observed across the whole sample.

Another remarkable feature on Fig. 16 is the large dispersion in the Ṁ_acc values observed at each stellar mass, covering a range of ~2 to 3 dex. Similar results have been inferred from previous statistical surveys of accretion in a given cluster (see, e.g., the recent studies of Rigliaco et al. 2011 on σ Ori and of Manara et al. 2012 on the ONC). This implies that only the average behavior of the actual trend of Ṁ_acc with M_∗ can be deduced, and suggests that different accretion regimes may co-exist within the same young stellar population. To this respect, we highlight in Fig. 16 the accretion properties of two, morphologically well-distinct, classes of YSOs in NGC 2264, both actively accreting. The first corresponds to the newly identified (Stauffer et al. 2014) class of objects whose light curves are dominated by accretion bursts. The second class includes objects whose light curves show periodic or quasi-periodic flux dips likely resulting from a rotating inner disk warp partly occulting the stellar emission (Bouvier et al. 2007a; Alencar et al. 2010). As can be seen in Fig. 16, both classes of objects show a level of accretion on average larger than the mean level of accretion activity detected across the population of NGC 2264; however, the two classes are well distinguished, with the first showing a mean accretion rate about three times larger than the second.

4.2. Mid-term Ṁ_acc variability

The scatter in Fig. 16 could in principle partly arise from the huge photometric variability exhibited by T Tauri stars, both on short-term (hours) and mid-term (days) timescales (e.g., Ménard & Bertout 1999). Short-term variability can affect the Ṁ_acc estimates inferred for single objects, but its contribution is likely negligible for our purposes (i.e., for a statistical mapping of accretion throughout the population of the cluster). On the other hand, mid-term variability, comprising the geometric effects linked with stellar rotation and the intrinsic accretion variability (timescales relevant for hot spots), might effectively blur the Ṁ_acc distribution inferred from the mapping of the region, since each object is observed at an arbitrary phase.

In order to assess the contribution of mid-term variability to the Ṁ_acc spread, we probed the variability of the detected Ṁ_acc, over timescales relevant to stellar rotation, by measuring the UV excess and correspondingly computing the mass accretion rate from each observing epoch during the CFHT u-band and r-band monitoring. On average, 16–17 points distributed over the 2-week long survey have been considered to probe variability. This procedure allowed us to robustly associate a range of variability to the individual average Ṁ_acc, inferred for each object as described in Sect. 4.1.

Figure 17 shows the results inferred from this analysis; each variability bar encompasses all values of Ṁ_acc detected for the corresponding object during the monitoring (i.e., it traces the difference between the highest and the lowest detected Ṁ_acc). Objects for which only an upper limit to the average Ṁ_acc had been derived have not been considered for this step of the analysis; for objects whose nominal Ṁ_acc value fall below the detection threshold (Sect. 4.1) at certain epochs, an upper limit to the lowest apparent Ṁ_acc has been drawn.

As can be observed in Fig. 17, variability properties are not uniform across the sample. No clear trends emerge with stellar mass or mass accretion rates; the geometry of individual systems might have a more important contribution in determining the apparent Ṁ_acc variability on week timescales, as suggested by the very few points on the diagram, maybe observed in a pole-on configuration, that show a negligible variation in the measured mass accretion rate. The (a)symmetry of individual bars, albeit potentially affected by the specific temporal sampling, may reflect the way accretion proceeds onto the central object, as well as the dominant features of the system: a strongly asymmetric bar with the median Ṁ_acc point close to the lower end can identify an object showing sudden, short-lived accretion bursts, while the opposite case well describes a system whose typical luminosity state is interspersed by short duration flux dips. These properties are indeed observed for some of the objects marked in Fig. 16.

The detected amplitudes of variability in Ṁ_acc range from small fractions to ~1 order of magnitude. In order to better visualize the average effects, we subdivided the sample of objects in Fig. 17 in three broad mass intervals (M_∗ < 0.4 M_⊙; 0.4 M_⊙ ≤ M_∗ ≤ 1 M_⊙;M_∗> 1 M_⊙) and for each of them computed the mean Ṁ_acc and amplitude of variability. These results are reported in Table 5, while a typical variability bar is represented in red in the lower right part of Fig. 17. On average, a variability of ~0.5 dex in Ṁ_acc is observed over a couple of weeks (which correspond to ~2–3 rotational cycles for most of the objects; see Affer et al. 2013), against a typical Ṁ_acc spread of ~2 dex. This comparison clearly shows that the observed Ṁ_acc spread at each mass is significantly larger than what can be accounted for by variability on week timescales, be it intrinsic (non-steady accretion) or geometric (rotational modulation).

4.3. Accretion variability: intrinsic vs. modulated

As discussed in Sect. 4.2 (see also Sect. 5.2), an average variability of ~0.5 dex is derived when monitoring the observed Ṁ_acc for individual objects over days-to-weeks timescales. This amount of variability is expected to incorporate both the intrinsic variability associated with the accretion process and the geometric effect of varying accretion shock projection along the line of sight as the star rotates.

An illustration of this is shown in Fig. 18.

The underlying sinusoidal profile describes the geometric modulation of the accretion spot linked with stellar rotation, which is reflected in the measured UV excess and, consequently, in the value of Ṁ_acc inferred at a given phase. Two additional effects contribute to the variability picture: the intrinsic variability of accretion on days timescales, reflected in the variation of the UV excess measured at the same phase in two different cycles, and the intrinsic variability of accretion on hours timescales, reflected in the rms variation of the UV excess measured within the same rotational cycle around a given phase (i.e., measurements obtained a few hours apart during the same night).

In order to probe the relative significance of these different contributions to the monitored Ṁ_acc variability and characterize the intrinsic accretion variability over hours-to-days timescales, we selected a subsample of 35 accreting members of NGC 2264, subject to the following conditions:

• 1.

  being an actual detection in Ṁ_acc (i.e., not an upper limit);

• 2.

  displaying a detectable periodicity, known from the study of Affer et al. (2013);

• 3.

  having rotational period small enough to be able to detect at least a partial superposition of different cycles in the phased diagram.

Affer et al. (2013)’s study of rotation in NGC 2264 is based on a 23 day-long optical monitoring of the star-forming region performed with the CoRoT satellite in March 2008. Accreting stars with detected rotational modulation amount to about 16% of the CTTS in our sample. These are distributed over the whole Ṁ_acc range at all masses in Fig. 16, except at the lowest Ṁ_acc regimes (≲5 × 10^-10 M_⊙/yr), typically populated by objects fainter than 17 in r and hence not attainable with CoRoT.

For each of these objects selected as described above, we phased the light curve of the UV excess and measured independently the average rms variation of detections obtained within the same observing night and the average rms variation of the average UV excess measured in a given phase bin of 0.05 at different rotational cycles.

The average peak-to-peak variation in the color excess throughout the sample amounts to ~0.47 mag, which corresponds to a variation of ~0.52 dex in Ṁ_acc, consistent with what we found in Sect. 4.2. 50% of the objects show an average intrinsic rms variation of <0.1 mag at a given phase from cycle to cycle (80% at ≤0.15 mag and 90% at ≤0.25 mag). Hence, on average a ratio of ~0.25 between the intrinsic rms variation on day timescales and the peak-to-peak variation can be inferred. The average rms variation on hour timescales is smaller, amounting to ≤0.05 mag. On the assumption that the two contributions (accretion variability on hours vs. days timescales) are mutually independent, we conclude that these two sources of variability together contribute on average just ~28% of the observed peak-to-peak variation in the UV excess. For a typical ΔṀ_acc of 0.52 dex, intrinsic accretion variability is thus expected to contribute for about 0.15 dex, while a major contribution arises from the geometric modulation of the accretion shock.

Interestingly, in the picture of modulation-dominated variability and assuming a homogeneous phase sampling, we would expect to statistically observe symmetric variability bars around the average value of Ṁ_acc; this is indeed consistent with the results reported in Sect. 4.2 and synthetized in Table 5.

5. Discussion

Accretion is a key ingredient in early stellar evolution: it regulates the star-disk interaction over the first few million years, with a major impact on the dynamical evolution of both the central object and the protoplanetary disk. Inferring a complete picture of accretion within large, well-characterized samples of young stars allows us to explore the accretion properties over a variety of parameters (e.g., stellar mass, age, disk properties), and hence represents a crucial issue for understanding how accretion evolves on an individual vs. statistical basis.

5.1. The M_∗ − Ṁ_acc relationship

Understanding what relationship links disk accretion rates and the mass of central objects has been a major focus for studies of accretion in YSOs over the last years. A prime motivation for this lies in that the M_∗–Ṁ_acc relationship likely reflects the actual mechanisms governing the accretion process.

Hartmann et al. (2006) reviewed a number of different mechanisms which may be relevant to disk accretion and its variability and to the M_∗–Ṁ_acc relationship. As noted by those authors, several unrelated mechanisms may contribute to establishing this relationship, and those mechanisms may act differently in different mass regimes. A number of previous studies have attempted to better understand the mass dependence of the M_∗–Ṁ_acc relationship (Vorobyov & Basu 2008; Rigliaco et al. 2011; Fang et al. 2013); these analyses suggest a diminution of the steepness of the Ṁ_acc–M_∗ dependence as M_∗ increases from a few tenths or hundredths of solar mass to a few solar masses. This is consistent with what we observe here for NGC 2264; on the other hand, the overall trend observed in the subsolar-to-solar mass regimes appears to be significantly shallower than that observed in the intermediate (1–10 M_⊙) mass regime (Ercolano et al. 2014).

In this study, we examined the distribution of mass accretion rates over mass for a homogeneous population of 236 confirmed, accreting members in NGC 2264, covering the mass range ~0.1–1.5 M_⊙. We inferred the presence of a correlation between Ṁ_acc and M_∗ and further showed that this result is not significantly affected by a mass-dependence in the Ṁ_acc detection limits. The same fraction of accreting stars over the total young population (~40%) is roughly observed at lower (<0.4 M_⊙) and higher (>1 M_⊙) masses; unless the true fraction of accreting objects is considerably larger at higher masses than at lower masses, this attests that no significant selection biases towards lower-mass objects affect our accretion survey. Hence, we obtain $\mathit{Ṁ}{}_{\mathrm{acc}}\mathrm{\propto }{\mathit{M}}_{\mathrm{\ast }}^{\mathit{\alpha }}$, with α ~ 1.4 ± 0.3 (Eqs. (15)–(16)). This exponent is smaller than the α = 2 dependence inferred by Muzerolle et al. (2003) over the mass range 0.04–1 M_⊙; the trend they obtained, however, is strongly influenced by the substellar component of the PMS population, which is not probed in our study. Indeed, results similar to Muzerolle et al.’s findings have been subsequently inferred from studies probing a comparable mass range in the substellar down to planetary mass regimes (e.g., Mohanty et al. 2005 for a composite sample of members from different star-forming regions, Natta et al. 2006 in ρ Ohiuchi, Herczeg & Hillenbrand 2008 in Taurus, Alcalá et al. 2014 in Lupus). On the other hand, studies addressing more specifically the T Tauri mass regime have recently proposed somewhat lower estimates for the quantitative dependence of Ṁ_acc on M_∗. Rigliaco et al. (2011) probed a mass range, in σ Ori, similar to the one investigated in Muzerolle et al.’s study (0.06–1 M_⊙) and inferred a value of α = 1.6 ± 0.4. Barentsen et al. (2011) studied a population of 158 accreting stars, covering the mass range 0.2–2 M_⊙, in IC 1396, and inferred a M_∗–Ṁ_acc relationship with α = 1.1 ± 0.2. All these results are quite consistent with the estimate of slope here inferred. The range of α from these studies may in part reflect the inhomogeneous nature of the methodologies and stellar samples; however, taken at face value, the results suggest a possible dependence of α on mass, with the higher mass regime showing a shallower slope.

5.2. Accretion variability: what contribution to the Ṁ_acc spread at a given mass?

A ubiquitous result of the studies of Ṁ_acc distributions for large samples of young stars is the significant spread observed in Ṁ_acc at a given mass, enveloping the average trend over 2 to 3 orders of magnitude. CTTS are known to display a strong variability over different timescales; this takes contribution both from geometric effects linked with the rotational modulation of surface starspots and disk features (Herbst et al. 1994) and from the intrinsic variability of the accretion process. Variability could thus be partly responsible for the observed scatter.

In order to probe this effect, we monitored the UV excess and Ṁ_acc variability for about 150 NGC 2264 accreting members over two full weeks, with typically 16–17 epochs for each object; this allowed us to show that variability on week timescales determines on average a variation of the detected Ṁ_acc over ~0.5 dex, significantly smaller than the extent of the average Ṁ_acc scatter.

This study complements the results inferred from recent studies that explored the variability in mass accretion rates on longer timescales. Nguyen et al. (2009) analyzed the Ṁ_acc variability of a sample of 40 members in Taurus-Auriga and Chamaeleon I, from multi-epoch optical spectra distributed over ~10 months, with typically 4 epochs per object on a baseline varying from hours to months. They probed the accretion properties from Hα 10% width and Ca II flux and inferred a median variability estimate amounting respectively to 0.65 dex and 0.35 dex. Fang et al. (2013) probed the Ṁ_acc variability of several tens of members in L1641, monitored on timescales of 1, 10 and 22 months; measuring the amount of accretion from the full width of the H_α line at 10% intensity, the authors inferred a typical Ṁ_acc variability of ~0.6 dex.

A recent spectroscopic survey of the Ṁ_acc variability in T Tauri and Herbig Ae stars, addressing shorter (≲hour to days) timescales, has been performed by Costigan et al. (2014). The authors probed accretion properties and their variability from H_α EW measurements for a sample of 14 objects; targets have been monitored over ≲1 h blocks, in some cases iterated within a single night, and repeated over a few nights. Combining the results of this short-term variability monitoring with multi-year information for a few objects, as well as with results from different surveys, the authors found that the dominant contribution to the observed variability in Ṁ_acc arises from the days timescale, with a spread in accretion rates ranging from 0.04 to 0.4 dex.

All these results robustly show that variability on day-to-month timescales does not impact significantly on the large spread observed in Ṁ_acc around the average Ṁ_acc–M_∗ trend; hence, the spread is real and possibly linked to the different evolutionary stages of individual objects (see Sect. 5.3). Additionally, a strong overlap in the age distributions for the populations of accreting (CTTS) and non-accreting (WTTS) members seems to hold; modulo the issue of an accurate age determination for young stars, this suggests that an intrinsically large range in Ṁ_acc exists at a given mass and age, which may reflect various other parameters (like, e.g., the properties of the stellar magnetic field or the initial mass of the disk).

Remarkably, the order-of-magnitude estimates for the Ṁ_acc variability inferred on day-to-week timescales (this study; Costigan et al. 2014) and on month-to-year timescales (Nguyen et al. 2009; Fang et al. 2013) are quite consistent; this suggests that the same mechanisms dominate variability over these baselines, with a major contribution arising from the mid-term variability component. This picture is supported by the results of the comparison, shown in Fig. 13, between the ${\mathit{L}}_{\mathrm{UV}}^{\mathrm{exc}}$ measurements obtained from ours and Rebull et al. (2002)’s survey of NGC 2264, respectively; what we observe in this case is that, indeed, UV excess luminosities measured at the two epochs are globally consistent, with an rms scatter, about the equality line, of the same order of magnitude as the average amount of variability detected during our 2 week-long CFHT monitoring.

5.3. What is the origin of the intra-cluster spread in Ṁ_acc?

A star-to-star variety in the evolutionary stage and evolutionary scenario may provide an important contribution to the interpretation of the large range of accretion properties observed within the region. The cluster shows a hierarchical structure (see Dahm 2008 and references therein), with regions of active, current star formation embedded in a more widespread population of somewhat older stars (Sung et al. 2009). Several episodes of star formation may have occurred, as proposed by Adams et al. (1983). Despite the uncertainty on a quantitative estimate, there is a general consensus about an age spread of up to several Myr within the population of NGC 2264. It is thus of interest to attempt a characterization of the accretion properties relative to age indicators and disk evolutionary stages.

Looking at how points distribute on the HR diagram of the cluster (Fig. 12) and comparing this to model isochrones provides a debated indication on nominal ages. In order to probe the evolution of accretion properties with stellar age, we inferred an age estimate for each source from isochrone fitting, adopting the models of Siess et al. (2000). The results for accreting sources are displayed vs. the Ṁ_acc estimates in Fig. 19. Only objects which are actual detections in Ṁ_acc and whose age estimate from models is between ~10⁶ yr and 10⁷ yr have been retained for the comparison.

An overall decreasing trend of the average Ṁ_acc with increasing stellar ages is traced in Fig. 19, consistently with previous results regarding the average evolution of accretion with stellar age (e.g. Hartmann et al. 1998; Sicilia-Aguilar et al. 2010; Manara et al. 2012). The anticorrelation trend is well described in the M_∗ < 1 M_⊙ regime, reaching a significance of over 5σ according to a Spearman’s rank order test (Press et al. 1992). A linear regression to this group of objects yields $\mathit{Ṁ}\mathrm{acc}\mathrm{\propto }{\mathit{t}}^{\mathrm{-}\mathit{\eta }}\mathit{,}\hspace{0.17em}\mathit{\eta }\mathrm{\simeq }\mathrm{1.51}\mathit{,}$(17)quite consistent with the η = 1.5 estimate inferred by Hartmann et al. (1998). Little information can be achieved on the accretion properties of the higher-mass group (M_∗ > 1 M_⊙) from our UV excess analysis.

The actual Ṁ_acc-age relationship may vary depending on the mass range probed, as suggested by the looser correlation that can be inferred from the least massive group of objects (M_∗ < 0.4 M_⊙) compared to the 0.4–1 M_⊙ group. This might reflect truly different evolutionary timescales, but also be affected by uncertainties on the age estimates inferred from models. Indeed, a clear-cut decreasing trend is observed on the lower-left envelope of the data point distribution in Fig. 19, with the youngest group being among the most luminous sources in the sample. While limits on the detection of small accretion rates in luminous stars might somewhat affect the observed trend, a more significant bias here is likely introduced by the model age estimates. We notice that lower-mass stars appear to be systematically older than higher-mass stars, based on model ages. If this is not a physical effect, but rather an artifact driven by systematics in the model, the actual anticorrelation trend that can be derived from Fig. 19 is, at best, doubtful.

The reliability of such time evolution analyses of accretion is subject to the assumption that the interpretation of the luminosity spread for a given PMS population on the HR diagram in terms of a true age spread is overall correct. However, this interpretation is subject to a severe controversy. Indeed, recent studies have questioned the accuracy of individual age estimates from the stellar luminosity. Baraffe et al. (2009) claimed that the luminosity spread observed on the HR diagram of a given star-forming region does not correspond to an actual age spread, but merely reflects the individual episodic accretion history, whose dynamics may mimic an age spread of up to a few Myr for coeval stars. While no consensus has obviously been achieved on this issue, an intrinsic age spread on the order of Myr, if present, is likely to be reflected in different properties of the region, such as the spatial distribution of members. This point is discussed in Sect. 5.4.

As mentioned in Sect. 3.1, information on the circumstellar disk morphology, hence on the evolutionary status, is conveyed by the mid-IR slope of the SED of the system. In Fig. 20, the accretion properties for three different mass groups of members are compared to the α_IRAC classification from Teixeira et al. (2012). Smaller values of α_IRAC correspond to more evolved disks; the average accretion rates decrease as disks become more evolved, albeit within a significant spread for each mass group.

These comparisons show that several parameters, beyond stellar mass, regulate the accretion evolution of individual systems and eventually contribute to shape the complex picture of a star-forming region such as we may infer from an extensive accretion survey.

5.4. A spatial mapping of accretion properties in NGC 2264

Characterizing the spatial distribution of accreting members according to their accretion properties is a useful tool to investigate the evolution dynamics of the region.

Figure 21 shows a spatial mapping of accretion in NGC 2264. Most accretors are projected onto the cluster; a part of the population of accreting or mildly accreting members are instead located on the peripheral regions. Some interesting features can be observed: the strongest accretors tend to be concentrated around the two known active sites of star formation within the cloud, a few parsecs apart, while the population of more moderate accretors tends to be more evenly distributed throughout the whole cluster. The halo of members is spread over a few parsecs around the star-forming sites.

This differential distribution of mass accretion rates may be the product of a sequential star formation process: the youngest stars, likely the most actively accreting, are indeed expected to be found close to their birth sites, while dynamical evolution over a few Myr (the average age of the cluster) can explain displacements up to a few parsecs relative to the original location where stars were formed. If this is the case, the distribution of Ṁ_acc also reflects the presence of an intrinsic age spread within the region, hence corroborating the interpretation in terms of a time evolution over different ages.

We attempted to pursue further this scenario by investigating the spatial distributions of distinct, equal-sized groups of objects at the two age extremities in Fig. 19. If our interpretation is correct, we would expect the presumed younger member group to be spatially more compact than the older one. A statistical method to probe spatial compactness is the minimum spanning tree (MST) method (Prim 1957). This (see also Allison et al. 2009; Parker et al. 2011) consists in computing the minimum path-length required to connect all points of a given sample in a network of direct point-to-point links with no closed loops. Shorter MST lengths will thus identify, to a certain degree of confidence, more compact populations. The significance of the result can be tested by randomly selecting two groups of objects of the same size and belonging to the same population as the two groups of objects under study, and computing the MST of the first random group, the MST of the second random group and the ratio between the two. The iteration of this procedure allows us to reconstruct the normal distribution of MST₁/MST₂ values when no significant difference exists between the spatial distributions of two given groups of objects selected among the population of interest.

For our population of accreting objects, sampled in random groups of 20 to 30 elements, we derived a normal distribution of MST ratios centered on 1.0 (homogeneous spatial distribution) with a σ of 0.2. Little evidence for an actual difference in the spatial distributions of “younger” vs. “older” members in Fig. 19 can be inferred when comparing the result of their MST test to this distribution. Somewhat smaller values than the mean of 1.0 are obtained from the MST_younger/MST_older ratio, oscillating around a significance of 1.5σ that this difference is truly meaningful. This might suggest that the younger population is to some extent more compact than the older population; however, no conclusive evidence can be drawn from this analysis, and the small statistics at lower and higher ages in Fig. 19 yields some dependence of the specific result on the specific selection of objects for the two groups. On the other hand, the complexity of the cluster, with its clumpy structure, the presence of at least two separate poles of star formation, the association with several massive OB stars, intrinsically translates to a non-linearity in the cluster evolution and, consequently, to a certain degree of ambiguity in the information provided by the spatial distribution of more evolved members.

6. Conclusions

We have performed an extensive photometric investigation of accretion in the young open cluster NGC 2264. Our sample, homogeneous in the dataset and in the characterization of individual properties, includes ~750 members, spanning in mass the range ~0.1–2 M_⊙; of these, about 40% show ongoing accretion. Accretion activity has been traced from the detection and monitoring of the u-band excess emission relative to non-accreting young stars; we showed that this excess is a reliable indicator of accretion, matching the properties observed at different wavelengths (H_α emission, IR excess), with the advantage of providing alone a direct proxy to the accretion luminosity.

We re-examined the membership of the cluster from an accretion-related perspective, exploring a wider area around the spatial distribution of currently known members within the star-forming region. This allowed us to expand the census of accreting members, with the identification of 50 new candidates based on their accretion signatures at short wavelengths. We detected UV color excess ranging from a few tenths of a mag to ~3 mag, with a sensitivity of ~0.2 mag deriving from the definition of the reference color locus from the properties of non-accreting young stars. Following a new empirical prescription for the conversion of the UV excess into accretion luminosity, we inferred a complete picture of accretion in the cluster and monitored its variability over a few rotational cycles. A broad interval of Ṁ_acc values has been inferred within the region, ranging from ~10^-7 to a few × 10^-10 M_⊙/yr; at any given stellar mass M_∗, a large spread in Ṁ_acc values, typically extending over 2 orders of magnitude, is observed.

The average Ṁ_acc shows a statistically significant correlation with M_∗; more uncertain is the functional form of this dependence, which in turn appears to be mass-dependent and likely time-evolving, as highlighted by the results inferred from previous studies for various young stellar populations. A statistical slope of 1.4 ± 0.3 for the log Ṁ_acc over log M_∗ relationship is here derived for NGC 2264 in the mass sample probed; while the absence of the lower-mass (<0.1 M_⊙) stellar component may contribute to explain this less steep dependence than the value of ~2 proposed in several studies, the large spread of Ṁ_acc values observed at a given mass attests the mere indicative value of such a relationship, which weighs the contributions of several trends as well as of undetected regimes.

Little contribution to the spread in the Ṁ_acc vs. M_∗ distribution arises from Ṁ_acc variability, which is found to amount to ~0.5 dex on the average, at least on days to months timescales; hence, this spread is physically relevant and provides valuable insight on the complexity of the accretion and evolutionary picture for a given stellar population. A differential distribution of evolutionary stages among members, such as we may deduce from the HR diagram of the cluster and/or from disk evolution tracers (α_IRAC), would naturally impact the global accretion picture we may infer for the cluster as a whole. Indeed, Ṁ_acc is shown to anticorrelate with such intrinsic/evolutionary age indicators, although these trends are in turn not exempt from a large amount of scatter. The detected spread in Ṁ_acc may thus reflect the dynamics of cluster evolution and the presence of a non-negligible age spread among its members, be it the result of sequential star formation episodes and/or of a multiplicity of evolutionary patterns. Moreover, a diversity in the accretion regimes at play, from episodic, short-lived accretion bursts to steadier, funnel-flow accretion, may concur to shape the features of different stellar groups and result in a range of observed accretion rates, from the upper envelope of the Ṁ_acc vs. M_∗ distribution to lower accretion rate regimes.

An important contribution toward the achievement of a complete physical picture of the accretion process is provided by a detailed analysis of the variability properties. Geometric effects constitute an important filter to the effective variability profile for a given system, contributing up to 75%, on the average, of the detected Ṁ_acc variability; however, much information on the underlying physics is conveyed by the specific photometric signatures of this variability, and a comparative analysis over intrinsically different stellar groups may be of remarkable interest to discern different scenarios. This aspect will be fully addressed for the young stellar population of NGC 2264 in a forthcoming paper (Venuti et al., in prep.); this is intended to provide a complete characterization of the accretion process and its impact onto the dynamics of the cluster.

Acknowledgments

We would like to thank Christian Veillet, former director of CFHT, for granting discretionary time to perform the mapping survey in December 2010, and the Terapix center at Institut d’Astrophysique de Paris, and in particular Yannick Mellier, for the prompt processing of the MegaCam images obtained during this run. We also thank Nadine Manset and Jim Thomas at CFHT for efficient run scheduling and data retrieval procedures. We warmly thank Luisa Rebull for retrieving for us the L_acc values she derived for her 2002 paper on disk-bearing objects in NGC 2264. We thank Kevin Covey for discussions on SDSS dwarf sequences and Lynne Hillenbrand for discussions on bolometric corrections scales; we also acknowledge useful discussion on the Ṁ_acc–M_∗ relationship with Beate Stelzer. This publication makes use of data products from the Sloan Digital Sky Survey and the Two Micron All Sky Survey. This project was in part supported by the grant ANR 2011 Blanc SIMI5-6 020 01. S.H.P.A. and A.P.S. acknowledge support from CNPq, CAPES and Fapemig.

References

Appendix A: Accretion-oriented revision of the census of NGC 2264 members: new CTTS candidatemembers and possible field contaminants

All Tables

All Figures