© The Authors 2023

1. Introduction

Star formation (SF) begins with the gravitational collapse of giant molecular clouds into molecular cores (Mac Low & Klessen 2004), which eventually evolve into protostars through mass accretion onto a central object. Once the main accretion phase is over, we are left with pre-main sequence (PMS) stars surrounded by circumstellar disks (Bouvier et al. 2007). These disks then dissipate within a few million years, as the stars continue their contraction toward the main sequence (MS). Understanding the mechanisms that drive the initial collapse and that regulate the SF process requires the combination of theoretical efforts with detailed observations of star-forming regions (SFRs) at different evolutionary stages and in different environments (Lada & Lada 2003). Indeed, the precise nature of the involved physical and dynamical processes have yet to be clarified. For example, it is still unclear whether SF is a slow, quasi-static process, regulated by the magnetic field and ambipolar diffusion (Tassis & Mouschovias 2004), or whether it occurs on shorter timescales, as expected if controlled by turbulence (Ballesteros-Paredes et al. 2007). Furthermore, recent studies suggest that, rather than a single SF episode, a molecular cloud may undergo several successive formation bursts, forming hierarchical spatial and dynamical structures, before becoming gravitationally relaxed. Also, the role of filaments is currently being actively investigated as the preferred SF locus (e.g., Bally et al. 1987; André et al. 2014; Montillaud et al. 2019a; Schisano et al. 2020; Naranjo-Romero et al. 2022). One approach to obtain a clearer observational picture of SF is to try to reconstruct the SF history of stars in active SFRs and to trace their stellar and cloud dynamics. However, for unbiased results, the full spatial structure of the region must be considered, including the outer regions often overlooked because of observational constraints. We may thus assess the presence of subclusters and trace them to temporally distinct and-or related SF episodes, and follow the dynamic evolution of the structures, such as the merging of filaments. Low-density haloes around SFR are also of great interest to obtain a full picture of SF. They might, for example, trace the evaporation of low-mass stars (leading to mass segregation), and thus be characterized by low stellar masses and characteristic space velocities, or they might trace and reveal older and now slowly expanding populations. Finding and characterizing these extended and low-density populations can be challenging due to the large contamination from field stars in the outskirts of SFRs. The problem is even more severe if most of these stars are old enough to have dissipated their circumstellar disks, making tracers of disks and accretion particularly biased and incomplete. The present investigation aims at improving our understanding of the spatial structure, SF history, and dynamical evolution of NGC 2264, which is one of the most prominent SFRs in the solar neighborhood.

NGC 2264 (Dahm et al. 2008) is a young stellar cluster populated by more than a thousand stellar-mass members; it is located about 760 pc from the Solar System (Park et al. 2000) and SF activity is still ongoing within its parental cloud. The median age of the cluster is approximately 3 Myr with an age dispersion of at least 5 Myr (Dahm et al. 2008; Venuti et al. 2019). Thanks to its relative proximity and low foreground extinction, NGC 2264 is one of the most accessible SFRs in our Galaxy. The MS and PMS populations of NGC 2264 are thus well-characterized, thanks to a large number of observations at all wavelengths. Infrared and X-ray observations have identified a complex spatial structure and two prominent sites of current SF activity in the central areas of the cluster (see, for example, Peretto et al. 2006). Sung et al. (2008) have shown that the SFR is preferentially elongated in the north-south direction and extends at least up to 3–5 pc from the field center or 2–3 pc from the S Mon O VII star and the Cone Nebula (see Fig. 1). Venuti et al. (2018) further investigated the structure and SF history of NGC 2264, using spectroscopic and photometric data from the Gaia-ESO and CSI-NGC 2264 surveys (Randich et al. 2013; Gilmore et al. 2012; Cody et al. 2014). They selected 655 likely members based on Li-, Hα-, IR- and UV- excesses, and on X-ray emission, confirming its hierarchical structure and an age spread of 4–5 Myr. Isochronal ages, gravity indexes, disk and accretion fractions, indicated multiple SF events: stars in the outer regions appear to be more evolved, having possibly been formed earlier and migrated outward. Moreover, SF was found to have started in the north and then propagated to the south, where it is still ongoing. An age difference between disk-bearing and accreting stars was also detected, indicating that disks outlast the end of the disk-accretion process. Disks around the O-type star S Mon were found to be shorter-lived than in the south possibly due to external photo-evaporation by the hot star.

Fig. 1. Optical DSS2 image of the 2.5×2.5 deg square field centered on NGC 2264 to which our study is limited. The locations of two prominent features, the O7 star S Mon and the Cone Nebula, are indicated by orange arrows. The areas shaded in blue and orange-brown are combined exposure maps of all available Chandra and XMM-Newton observations, respectively. Identifiers for the XMM-Newton field are also shown by red symbols in their respective lower-right corner. The two orange squares indicate the FoV of the CoRoT observations discussed in Sect. 2. The yellow circles and ellipse indicate the subregions identified by Sung et al. (2009). Contours of the four new spatial substructures identified in this study (Sect. 4.3), three inner ones and the outer polygon, are drawn in cyan.

Gaia DR2 data was employed by Kuhn et al. (2019) to investigate the internal dynamics of several young regions, including NGC 2264. Unlike most of the other regions, NGC 2264 showed no evidence of expansion or of subcluster coalescence. Buckner et al. (2020), focusing on the spatial clustering and dynamics of the central regions of NGC 2264, found evidence of a prolonged SF activity. Also the clustering of stars was found to be related to the evolutionary stage, the least evolved being more clustered. This implies that stars form in compact regions and then disperse. The compact subcluster of preferentially Class III stars around S Mon, was taken as evidence of a significant stellar feedback (photo-evaporation), confirming the finding of Venuti et al. (2018).

All of the studies described above were limited to the “classical” boundaries of the star-forming region, that is the “Halo” defined by Sung et al. (2008) and indicated in Fig. 1 by the outer yellow ellipse. Montillaud et al. (2019a,b) focused on G02.3+2.5, a north-south filamentary region, ∼0.2–1.0 degree north of the O7 star S Mon and studied its gas and dust emission and dynamics. Their multiscale analysis, from the compact sources, to cores and filaments, indicate that the recent and ongoing SF in this region is likely due to the collision of two filamentary structures. The collision itself might be due to the northward expansion of the H II region or gravitational collapse toward the main cluster to the south. The wider-field objective-prism H_α survey of Reipurth et al. (2004) uncovered a small number of dispersed accreting members, poorly correlated with the gas cloud and which they therefore suggested might have already drifted away from their formation locus. Venuti et al. (2019) also surveyed a broad region. They were able to identify the late-type population with limited field contamination, using exclusively Pan-STARRS1 and UKIDSS photometry. They retrieved the spatial structures of Sung et al. (2009), but also identified a more widespread population of ∼100 stars up to 10–15 pc away from the cluster center. As described below, clear indications of a larger extent of the region have been also derived independently by us, on the basis of the optical variability of ∼8000 stars in the NGC 2264 field, obtained with the CoRoT satellite (Baglin et al. 2006). This statistical study, briefly described in Sect. 2, was indeed our initial motivation for obtaining new X-ray observations of the outskirts of NGC 2264 so to identify and characterize its dispersed population, eventually aiming at a better understanding of how SF proceeded in the region. Since YSOs have X-ray luminosities 2–4 orders of magnitude higher than field stars (Walter et al. 1988), X-ray observations can easily distinguish high-probability YSO candidates from the field stars. Furthermore, X-ray studies are complementary to IR photometry: while this latter allows the identification of YSOs with excesses due to the presence of disks and envelopes, it does not identify PMS stars that have dissipated their disks.

In this work we also take advantage of a large body of literature data and, in particular, of the Gaia early 3rd data release (eDR3), both for defining membership and for characterizing young stars. As a byproduct of this effort we also test our approach for selecting candidate members of a SFR on the basis of optical variability. This approach might indeed prove useful in view of the present and future availability of suitable and extremely large databases of optical time-series for huge number of stars (such as Kepler, TESS, Rubin-LSST, Gaia, etc.).

The article is structured as follows: in Sect. 2 we give more details about CoRoT observations. In Sect. 3 we present literature data and the original X-ray observations we have used. For the X-ray observations, we describe the data reduction, source detection, and cross-identification procedures used to identify the optical and IR counterparts of X-ray sources. Section 4 introduces our primary results: the identification of samples of candidate cluster members, complete with estimations of their field-star contamination, and the definition of new spatial structures. In Sect. 5 we make use of our new member catalog to address a number of issues: first we discuss the efficiency of member selection via X-ray and variability criteria (Sect. 5.1); we then discuss the average physical properties of stars and their spatial variation. More specifically we asses density, distances, kinematics, masses and ages, and disk fractions. Finally, Sect. 6 summarizes our work and its main conclusions.

2. Motivation – Optical variability with CoRoT

CoRoT was a European satellite devoted to continuous optical observations of star fields for asteroseismological and planetary transit studies. We observed NGC 2264 with the exo-planetary camera of CoRoT in 2008 with an “additional program” (PI: F. Favata), obtaining 23.5 days of accurate and uninterrupted optical photometry for about 8000 preselected stars in a two adjacent 1.7 sq.degree fields of view (FoVs), each imaged by a different CCD (see Fig. 1). The northern region encompassed the whole SFR, while the southern one did not contain any members of NGC 2264 known at the time of the observation. The CoRoT targets were selected to include ∼270 known NGC 2264 members down to I ∼ 16, or M = 0.3 − 0.4 M_⊙ (see Alencar et al. 2010; Flaccomio et al. 2010), but the vast majority were randomly selected field stars in the same magnitude range. Their lightcurves are, as expected, statistically much less variable than those of known NGC 2264 members. However, some of these “field stars” appeared to be as variable as members. We thus embarked in a systematic search for variable stars that might be previously unrecognized members of NGC 2264. Variable field star contaminants were taken into account adopting, as a control field, the FOV of the CCD pointed to the south of NGC 2264, which is certainly dominated by field stars¹.

For each CoRoT lightcurve we measured: i) the standard deviation of the flux, and ii) the height of the maximum peak in the Lomb-Scargle Normalized Periodogram (LNP, Scargle 1982). The standard deviation, divided by the median flux, is a measure of the variability amplitude; the maximum of the LNP is instead related to the degree of periodicity, irrespective of the period. These two quantities are meant to preferentially single out Classical and Weak-line T Tauri stars (CTTS and WTTS), respectively. Most CTTS, indeed, have irregular or quasi-periodic lightcurves with large amplitudes, while WTTSs show regular modulations and smaller amplitudes (Alencar et al. 2010; Cody et al. 2014).

We plot the two quantities derived from the CoRoT lightcurves², one versus the other, in Fig. 2. Empty blue circles indicate a sample of likely cluster members, selected using literature data and at least one of the following criteria: X-ray detection, mIR excess, strong emission in the H_α line (from low-resolution spectroscopy or narrow-band photometry). Members clearly cluster in the upper-right corner, having large-amplitude and periodic lightcurves. The locus of field stars, from ∼3600 stars in the control field, is indicated by the grayscale map, which is a smoothed density map (black stands for maximum density). The cyan contour encircles 90% of field stars (in the smoothed map, 92% of actual data-points). Seventy-four per cent of the 93 known members lie outside of this contour and in the upper-right quadrant. We thus considered as candidate new members all the stars in the northern CoRoT field that fall in this region of Fig. 2 and that were not previously known as members (filled red circles). This selection results in a total of 170 candidate members (including previously known ones). Notably, the new candidate members have a broader spatial distribution with respect to the known members, and a noticeable higher concentration toward the west of the cluster. On the basis of our control field, however, 3.2% of the field stars are estimated to fall in the member-selection region of Fig. 2. If we make the simplifying assumption that a negligible fraction of new members is present among the 2089 stars in the north field that are not known members, we estimate ∼67 field star contaminants (3.2% of 2089), and a contamination fraction of 39% (67/170)³. This fraction might actually be lower since the CCD covering NGC 2264 also includes the densest part of the molecular cloud and is thus expected to contain fewer background stars than the control field, down to any given magnitude. In Sect. 5.1 we reexamine the variability-selected sample to assess the effectiveness of the chosen variability-based membership criterion.

Fig. 2. Height of the most significant LNP peak vs. standard deviation of light curves normalized to the median flux value. The grayscale map represents the density of stars in the control field, with the cyan contour encircling 90% of these stars. Blue open points are known members of NGC 2264, while red points indicate the new candidate members.

3. Observations and data analysis

Stimulated by the CoRoT observation discussed in the previous section, we obtained new X-ray data of the outskirt of the region, mainly to the west of the known cluster area. X-ray detection will indeed serve as one of our primary membership criteria in the outskirts of the region. We will also extensively exploit, for the same purpose, the Gaia astrometric data, from the early data release 3 (Gaia-eDR3). Additional optical and IR data from the literature will also be used, both as help in establishing membership, and to subsequently characterize the selected members. In the following two subsections we will discuss the literature data and the new X-ray data, respectively.

3.1. Literature data

We collected a large body of available data for the region, in order to define membership and characterize NGC 2264 stars. The largest field that we consider is a 2.5×2.5 deg square centered on RA = 06:40:48, Dec = 9:42:00, which likely contains the whole NGC 2264 population. An optical DSS2 image of the region is shown in Fig. 1. The yellow circles and ellipse toward the center of the field delimit the subregions defined by Sung et al. (2009). The outer ellipse, in particular, is indicated by Sung et al. (2009) as the boundary of the cluster “Halo”. The four subregions drawn in cyan are newly introduced in this paper (Sect. 4.3). We began with the master catalog of cross-identified objects also used by Cody et al. (2014). This contains several photometric and spectroscopic catalogs, optical and IR, mostly limited to the inner part of the region. Specifically, we include: SDSS system photometry from Venuti et al. (2014), optical Johnson-Cousins photometry (Dahm & Simon 2005; Lamm et al. 2004; Rebull et al. 2002; Sung et al. 2008), UKIRT nIR photometry (King et al. 2013), Spitzer mIR photometry (Sung et al. 2008; Cody et al. 2014), rotational periods (Lamm et al. 2004; Rebull et al. 2002; Venuti et al. 2017), and spectral types (Walker 1956; Makidon et al. 2004; Dahm & Simon 2005). To this catalog we added data from the Gaia-ESO survey (DR5) and the following wide area surveys covering the full 2.5×2.5 deg field in Fig. 1: Gaia (early 3rd data release, Gaia Collaboration 2021), Pan-STARRS1 (1st data release, Chambers et al. 2016; Flewelling et al. 2020, PS1 hereafter⁴), IPHAS (2nd data release, Barentsen et al. 2014), 2MASS (point source catalog, Skrutskie et al. 2006), and allWISE (Wright et al. 2010).

Each of these additional catalogs was added in sequence. Objects in our master catalog were matched to those in each of the catalogs using identification radii based on the provided positional uncertainties (limited to a minimum radius, chosen for each catalog) and the spatial resolution of the survey. A thorough visual examination was then performed at each step to resolve uncertain and-or ambiguous identifications. The allWISE catalog received an ad-hoc treatment. Because of the poor spatial resolution of the WISE mIR images, the appropriate identification radii were rather large and visual inspection of the ∼1.3 × 10⁵ catalog entries in our FOV turned out impractical. We therefore decided to limit identifications to isolated cases, that is allWISE sources with a single counterpart in our master catalog and for which the second closest object was farther away than both 3 times the identification radius and five arcseconds. A total of 7.6×10⁴ allWISE sources (∼57% of the full catalog) where thus matched with our catalog. Following a preliminary membership analysis (Sect. 4.1) we also identified the remaining allWISE sources with all candidate members (member sample “a” in Sect. 4.1). Following visual inspection, 2376 additional allWISE sources were matched with 3151 stars in our master catalog. This latter step was helpful both to assign mIR photometry to our candidate members and to confirm the membership status of stars on the basis of their mIR excesses.

NGC 2264 has been observed in X-rays multiple times with both Chandra and XMM-Newton. Figure 1 shows the spatial coverage of these X-ray observations. Four different sets of Chandra observations performed with the ACIS-I detector, have been presented by Ramírez et al. (2004), Flaccomio et al. (2006, 2010), and Cody et al. (2014; see also Flaccomio et al. 2018), covering the central S Mon and Cone regions with total exposure time between ∼50 ks (toward S Mon) and ∼460 ks (north of the Cone Nebula).

Most of the Chandra fields cover the southern cores, that is the Cone Nebula, IRS1, and IRS2 (or Spoke-cluster) regions (PIs: Sciortino, Micela), while a single shallower exposure (PI: Stauffer) covers the S Mon region. They all fall inside the Halo region as defined by Sung et al. (2008, 2009), the outermost yellow ellipse in Fig. 1. An analysis of all the above Chandra observations was presented by Townsley et al. (2019), who found 3373 point sources. However, we here prefer to use our own more conservative analysis of the same data, to be fully described in Flaccomio et al. (in prep.). The resulting catalog of 1043 X-ray sources was included in our master catalog. Basic information on these sources is listed in Table A.1.

Finally, results for the two central XMM-Newton fields were published by Dahm et al. (2007), but we did not include the respective source lists in our master catalog since we reanalyze these observations consistently with the new ones here.

3.2. XMM-Newton X-ray data

As mentioned above, two XMM-Newton exposures (PI: Simon) were taken, respectively, on 2001 March 20 and 2002 March 17, pointed toward the inner regions of the cluster (Dahm et al. 2007). They lasted ∼40 ks and were performed using a thick filter to minimize the contamination of the ultraviolet light from the massive cluster stars. Hereafter, they are named field N (for north) and field S (south). We obtained four additional observations between 2011 and 2014 each lasting about 50 ks (P.I. E. Flaccomio). These observations were pointed to the outer regions of the cluster and were performed with a medium filter. Hereafter these fields are named field 1, 2, 3, and 4. Figure 1 shows the footprints of the observed regions in the sky. The first five columns in Table 1 list, for each field, observation IDs (in the XMM-Newton archive), pointing coordinates, and observation dates.

3.2.1. Data reduction and source detection

The data of each EPIC detector, MOS1, MOS2, and PN, were reduced with the XMM-Newton SCIENTIFIC ANALYSIS SYSTEM (SAS; Gabriel et al. 2004) version 14.0.0. We reprocessed the PN and the MOS data to obtain tables of photons with astrometry, energies and arrival times. We then filtered out artifacts due to bad pixels, bad columns and cosmic rays. In order to maximize the signal-to-noise ratio of weak sources, background levels during each observation were monitored and data taken at times of high background were removed. To this effect we followed the standard SAS data analysis threads⁵. Figure 3 shows, as an example, the lightcurve of background photons detected by the EPIC-pn for field 3, the one with the worst background flares and, thus, the largest fraction of rejected data. The lightcurve refers to events with energy between 10 and 12 keV, pattern = 0, and standard PN flag filtering. In this case, time intervals with background count rates greater than 0.5 cnt/s were discarded. This rejection threshold was chosen individually for each observation and EPIC camera (PN, MOS1, and MOS2), by inspecting the relative background lightcurves. Columns 6–8 in Table 1 list the total and background-filtered exposure times (on-times) for each observation and detector.

Fig. 3. Background light curve of EPIC PN-detected events in field 3 (Obs.ID 0721400201); the horizontal red line is the count-rate threshold that we chose to define and filter-out high-background time intervals, thus maximizing the S/N of faint sources.

Source detection was performed for events in the 0.3–7.9 keV energy band, appropriate for coronal sources, using the PWXDetect wavelet-based code developed at the INAF – Osservatorio Astronomico di Palermo by Damiani et al. (1997). The code combines the MOS and PN data so to maximize the detection efficiency for faint sources. Significance thresholds for source detection were determined so to yield, on average, one spurious detection per field due to background fluctuations. These thresholds were determined individually for each of our six datasets, by simulating 100 blank (background-only) fields and running the detection code on these fields. The blank fields were simulated adopting the actual exposure map and number of background events, this latter determined subtracting the number of source events, from a preliminary detection, from the total number of detected events. The final significance thresholds, 4.9 or 5.0σ, depending on field, are listed in Col. 9 of Table 1 followed, in Col. 10, by the number of detected X-ray sources (1081 in total). Detections were inspected individually so to spot obviously spurious sources produced by residual instrumental artifacts that were not filtered out by the data processing described above. Thirty-five sources were then removed from the detection lists: twenty-nine in field 2, four in field N, and two in field S.

We finally created a merged XMM-Newton source catalog, by cross-identifying the six individual lists from each observation. Because of the (small) overlap between the fields, some sources are indeed detected in more than one observation. We employed an iterative merging process: we started by cross-identifying sources detected in field N and field S. The resulting source list was then cross-identified with that of field 1, and so on until all lists were merged. The cross-identification radii of individual detections were arbitrarily chosen as times the formal positional uncertainty provided by the PWXDetect code, which is based on the observed count statistic and photon spatial distribution. After merging each catalog, new positions and uncertainties were computed as the uncertainty-weighed averages of those of the matched detections. After each step, all identifications were visually checked and compared with the spatial distribution of photons. In eight cases we manually cross-identified source pairs that had not fulfilled our matching criterion. This can be quite reasonable, since our identification radii, based on positional uncertainties, can be significantly smaller than the spatial extent of X-ray sources (for example for extended sources or, more often in our case, groups of unresolved point sources). Our final catalog of unique XMM-Newton X-ray sources, 944 entries, is presented in Table 2. We list: sky positions, uncertainty-based identification radii, detection significances, background-subtracted source counts, and photon fluxes, as estimated by the PWXDetect code⁶.

3.2.2. Cross-identification with the master catalog

We matched the X-ray sources with our master catalog, described in Sect. 3.1, using a positional criterion. Identification radii were chosen as the quadrature sum of the identification radii for the XMM-Newton sources and those for the master catalog. For XMM-Newton sources the radii are listed in Table 2 and range between 1.0″ and 10.7″ (median: 3.6″). For the master catalog, radii were taken, for each object, as the smallest positional uncertainty among those adopted for the original catalogs in which the object is present. They range between 0.3″ and 8.1″ (median: 0.5″).

A significant fractions of X-ray sources end up with multiple counterparts in our master catalog. This is due, on one hand, to the sizable identification radii for the X-ray sources, and on the other, to the large density of objects in our master catalog, the large dynamic range covered, and, quite importantly in parts of the FoV, the presence of many obviously spurious entries in some of the constituent optical-IR catalogs. The IPHAS and UKIRT catalogs we have adopted, in particular, clearly include a large number of spurious sources, many of which can be ascribed to the nebulosity in the region or to the wings of bright stars. For many XMM-Newton sources, therefore, some of the potential optical-IR counterparts can be discarded with confidence.

We thus resorted to inspecting each identification by eye and evaluating the multiple optical-IR counterparts of any given X-ray source based on: differences in magnitudes, in one or more bands, and in positional offsets from the X-ray source, their position relation to the actual spatial distribution of X-ray photons, deep optical images (such as those from Pan-STARRS), which help identify likely spurious sources. In particular, we most often discarded counterparts that were > 5 magnitudes fainter than another counterpart. We also often discarded faint counterparts only listed in the UKIRT or IPHAS catalogs, when other brighter and more credible counterparts, often detected in multiple optical-IR catalogs, were also present. At the end we discarded 1866 counterparts to 570 XMM-Newton sources, out of the initial 3050 counterparts to the 944 XMM-Newton sources.

In several cases, we also added counterparts that formally did not match the positional identification criterion. In some cases, a prominent optical-IR source fell close to the edge of the identification circle, whose size is based on the positional uncertainty with a given confidence level. In these cases both statistical and unaccounted-for systematic uncertainties (for example due to the complicated PSF of XMM-Newton) can explain the observed offsets. In a number of cases, two (or more), often similarly bright, optical-IR sources may lie at opposite sides of the identification circle. Looking at the X-ray photons spatial distributions, some of these X-ray sources are more extended than expected for a point source. In some cases the X-ray source can even be seen to be marginally resolved into two separate components. In these cases, in which the uncertainty-based identification radii are not very informative, we added all the reasonable counterparts outside of the identification circles. A total of 129 counterparts were added for these reasons to 103 XMM-Newton sources.

The result of the above changes to the initial identifications is that the number of X-ray sources with unique identifications is increased from 200 to 680 and the number of X-ray sources with no optical-IR identification is reduced from 108 to 75.

Finally, we inspected once more the 264 X-ray sources that, after the above steps, were left with more than one possible counterpart. In doing so we searched for additional optical-IR counterparts that were somewhat unlikely to be responsible for the observed X-ray emission, for example because of their larger offset and significantly smaller optical-IR luminosity with respect to other counterparts. Contrary to the discarded counterparts described above, these identifications were kept but were flagged as unlikely. To this effect we examined again positional offsets from the X-ray position and magnitude differences, and considered as unlikely counterparts those that had both significantly larger offsets and were > 2 magnitudes fainter than the brighter and closer counterpart(s). We also examined the diagrams that provide membership information and which are fully discussed in Sect. 4.1. In cases in which at least one of the counterparts appeared to be an NGC 2264 member, we considered as unlikely counterparts those that had significant indication of being nonmembers (field stars or background AGNs, statistically expected to be fainter than members in X-rays). In a couple of cases we considered as unlikely a counterpart which appear to be a spurious source in the original literature catalog, based on deep archival optical images and other plausibility criteria (for example extremely close pairs of sources that could not have been resolved in the literature data). For XMM-Newton sources also detected at higher spatial resolution by Chandra, we also considered as unlikely counterparts optical-IR objects that were > 2 mag fainter than the unique counterpart of the Chandra source and that were clearly not detected in the Chandra data. Disregarding all the above unlikely counterparts, the number of XMM-Newton sources with a single likely counterpart increases from 680 to 735.

4. Results

The master catalog we have assembled lists ∼6.5×10⁵ objects throughout our 2.5×2.5 deg FoV. We now proceed to use the information gathered from the X-ray and optical-IR catalogs to identify members of NGC 2264, estimate field-object contamination, and discuss the spatial extent and structure of the cluster.

4.1. Candidate member selection

The first row of panels in Fig. 4 shows, for all objects in our catalog, the IPHAS r-H_α versus r − i color-color diagram, and three Gaia-based diagrams utilizing: Gaia magnitude and colors (G versus BP − RP)⁷, the two components of the proper motion vector (PM_RA versus PM_Dec), and parallaxes (PLX versus G). Given the large number of plotted points (first number in the bottom-right of each panel) we plot smoothed grayscale maps of the point density. The maps are binned into 128×128 pixels and convolved with a Gaussian kernel (FWHM = 3 px). Dashed blue lines in the plots indicate member loci, defined below, from which candidate members are selected, with varying degree of contamination from field objects (dotted blue lines, when shown, indicate more conservative, less-contaminated, definitions). Solid red lines indicate the boundaries of field loci, from which likely field objects are identified (in two of the panels plotted here the boundaries of the field loci coincide with those of member loci but this is generally not the case). The second row of panels in Fig. 4 shows the same four plots for the subsample of objects falling within the XMM-Newton field of view. Field object contamination within all defined member loci is clearly severe for both the full sample and for this latter subsample. The cluster can, however, be clearly identified in the proper motions and parallaxes plots referring to the smaller field. The third row refers to X-ray sources, detected by Chandra and-or XMM-Newton, and with an unique (or preferred) optical-IR counterpart. Comparison with the preceding row clearly demonstrates the power of X-ray observations in selecting candidate young members of a star-forming region. Indeed, using these plots for X-ray sources, we are able to define member and field-object loci for the full sample. Finally, the fourth row of panels in Fig. 4 refers to our working candidate member sample (sample “c”), which is introduced below.

Fig. 4. Object density in four different diagrams (one per column) and four different samples (one per row). The diagrams are, from left to right: the r-Hα vs. r − i color-color diagram using IPHAS data, G vs. BP − RP color magnitude diagrams (with three reddening vectors, at different colors), PMDec vs. PMRA, and parallax vs. G magnitudes (the last three diagrams using Gaia data). The four samples, from top to bottom, are: all objects in our catalog, objects within the XMM-Newton field of view, objects unambiguously detected in X-rays, candidate members in our sample “c” (Sect. 4.1). Density map are smoothed using a Gaussian Kernel. The grayscales range from black at maximum density to white at zero density, scaling with the square root of the density, rather than linearly, to reveal more details at low-densities. Isolated gray dots in low-density areas indicate individual stars. Dashed and solid lines delineate member and field loci, respectively, and, in the first row, “m” and “f” labels help clarify the distinction between the two. Dotted lines in the last two diagrams refer to more conservative member loci. The first number in the bottom-right corners indicate the total number of plotted objects. The following ones refer to the number of stars in the member- and field-loci, with 3- and 2-σ significance, respectively.

In addition to the plots described above, using IPHAS and Gaia measurements, several other diagrams providing membership information can be constructed with our data. We used 19 additional plots, presented in Appendix B. In 16 of these we defined field-star loci, while in 9 we defined member-loci. As mentioned, for each of the panels in Fig. 4 and B.1, we selected as candidate members all the objects that fall in the member locus (delimited by dashed blue lines), with a 3σ confidence, where the confidence was estimated from the uncertainties in the plotted quantities and Monte Carlo methods⁸. The same procedure was followed to identify likely field objects, with the difference that the confidence level for this selection was lowered to 2σ (eventually resulting in a more aggressive vetting of members).

We end-up with twelve membership criteria from these plots, plus X-ray detection, and twenty-one exclusion (field object) criteria. Because of the widely variable spatial coverage and depth of the diverse datasets we have collected, only a fraction of the ∼6.5×10⁵ objects in our full catalog have the data to test each of these criteria. The numbers of objects placed in each plot are given in the lower-right corner (first figure), and range between 550 for the spectroscopic Hα data and 3.4×10⁵ for the J vs. J − H plot. A significant fraction of objects can be placed in one or more of the plots providing membership criteria (68%) or field-object criteria (86%). For 23% (67%) of the objects three or more membership (field-object) criteria can be defined.

In both Fig. 4 and B.1, the numbers of candidate members and field-objects selected by each criterion are shown by the second and third figures in the lower-right corner of the relative plot. The number of candidate members range between 173, for the IPHAS Hα data, and 6418 for the PM data. Most, if not all, of the selected member subsamples are contaminated by field-objects to varying degrees. For example, it is clear from the PM plots that, over the whole field, candidate members selected from PM alone are actually dominated by field objects belonging to the broad distributions of points roughly centered at PM_RA = 0, PM_Dec. = 0. Indeed, 79% of these candidate members also fall in the field-object locus (with 2σ confidence) of at least one other plot. The same considerations apply to the 2167 candidate members selected in the parallax versus G plane, 78% of which also fall in at least one field-object locus⁹. Lower but significant “contamination fractions” (as approximated above), are also estimated for other membership criteria, with the lowest estimates, 25%, 28%, 31%, and 33%, found for criteria based on H_α EW, H_α FWZI, X-ray detection, and the Spitzer [3.6]–[5.8] versus [3.6]–[4.5] color-color diagram, respectively.

Conversely, actual NGC 2264 members are likely to fall into (and contaminate) several field-object loci, for example because of the very large extinctions of stars highly embedded in the molecular cloud, or large systematic errors in the photometry in one or more bands (for example because of contamination from neighboring objects or nebulosity), or even because of more physically interesting causes. One example are the so-called below-main-sequence stars, which can easily be mistaken for field stars. These are often found in young clusters and their peculiar magnitudes or colors may be due to the edge-on viewing geometry of the star-disk systems, and the associated attenuation and scattering of their photospheric emission by the disk material (see, for example, Guarcello et al. 2010; Bonito et al. 2013).

Given the significant contamination of individual membership and field-object criteria, we try to combine them to define samples of likely members that may serve different purposes. Our sample “a” tries to include most members, regardless of contamination; sample “b” tries to include only bona-fide members, minimizing contamination; sample “c”, finally, provides a reasonable compromise between completeness and contamination.

Sample “a” includes the 10267 stars for which we have any indication of membership, that is those that fall in at least one of the member loci or are X-ray detected (with unique optical-IR identification). All of these stars are listed in Table 3, available in its entirety on-line. For each of these objects we present sky coordinates along with Gaia, PS1, 2MASS, CSI, Chandra, and XMM-Newton identifiers. Two sequences of flags then show the criteria indicating member and field-object status. In the following, we do not further discuss sample “a”, focusing instead on subsamples with smaller contamination from field objects.

After extensive experimentation, and guided by the estimates of field star contamination discussed in the next section (Sect. 4.2), we define samples “b” and “c” as comprising stars that: (i) fall in more member loci than field-object loci, or, (ii) fall in at least three member loci, or (iii) are indicated as members by at least two of the “stronger” criteria, that is loci indicative of accretion or presence of disks, and X-ray detection. This latter provision has the advantage of being able to retrieve several of the below-main-sequence members discussed above. Sample “b”, which yielded the lowest field star contamination fraction among the investigated samples, was defined as the 1971 stars fulfilling the above criteria when considering the most stringent member loci for PMs, parallaxes, and RVs, (the two rightmost panels in Fig. 4, and fifth panel in Fig. B.1, respectively). Sample “c”, finally, comprises the 2257 stars that fulfill the same criteria when considering the larger, less conservative, member loci. The last column in Table 3 lists the samples each star in sample ‘a’ belongs to. Stars in sample ‘c’ are plotted in the fourth row of Fig. 4 and in Fig. B.2.

Since most membership criteria can only be applied to a (small) fraction of the considered field, in order to avoid spatial biases and, crucially, to discuss the overall spatial distribution of stars in the region, we also introduce a further sample, defined similarly to sample “c”, but only using information and catalogs available for the full field: 2MASS, Gaia, PS1, IPHAS, and WISE. This latter sample includes 1828 star and we refer to it as “c-wide” (indicated as ‘w’ in the last column of Table 3). In spite of the better homogeneity across the plane of the sky, it is also biased in several ways, for example because it includes a smaller fraction of the most embedded (and youngest) members which are best selected through deep X-ray and Spitzer mIR data.

Figure 5 shows sky density maps for samples “b”, “c-wide”, and “c”.

Fig. 5. Spatial distributions for candidate member samples “b”, “c-wide”, and “c” (from left to right). The density maps are smoothed with a Gaussian kernel as in Fig. 4, and dotted lines show the boundaries of the regions discussed in the text (Sect. 4.3) and also shown in Fig. 1. The figures at the bottom right of each panel indicate the number of objects in the sample.

4.2. Field-object contamination

Because of the complex procedure with which we built our candidate member samples, estimating the degree of field-object contamination is not straightforward. Assuming that the SFR is fully contained within our 2.5×2.5 deg field, we decided to estimate the sky density of field objects from the density of member candidates far from the cluster center. More specifically, we start by constructing density profiles of the candidate-member samples. Considering the overall north-south elongation of the association and its asymmetry with respect to the north-south direction, as apparent from Fig. 5, we construct our density profiles considering elliptical annuli (rather than circular) and four separate quadrants (north, south, east, and west). The elliptical annuli, in 10′ steps, are centered and have the same ellipticity (a/b = 0.7) as the Halo region defined by Sung et al. (2008). Figure 6 shows, in different panels, these profiles for samples “b”, “c-wide”, and “c”. The insets in the upper-right corners show the spatial distribution of stars in each sample, color-coded according to the quadrants, and superimposed on the elliptical annuli used to compute densities. The main plots show the stellar density as a function of major axis radius (R), separately for the four quadrants.

Fig. 6. “Radial” distributions for samples “b”, “c-wide”, and “c”. Different colors refer to different “quadrants”. Spatial distributions are shown at the top-right of each panel, with stars color-coded according to their quadrant. The elliptical annuli where the densities for the main plot were computed are also shown. Points with error bars, based on Poisson uncertainties, are slightly shifted along the x-axis to avoid overlap among the four independent samples. The dotted black line refer to the combination of the four quadrants. The thick horizontal magenta segments toward the bottom right indicate the average densities of stars in the left quadrant for r > 50 arcmin.

The density falls-off more rapidly in the eastern direction (magenta curve), as also clearly observed in Fig. 5, and appears to level off for R > 50′. We thus estimate the density of contaminants as the mean density of objects in the eastern quadrant for R > 50′, indicated in Fig. 6 by a horizontal magenta segment. We report, at the bottom of each panel, the number or stars in the sample, the estimated contaminants, both in number (density × area of the full field) and as a percentage of the total, and the estimated number of real members within the sample. We must keep in mind, however, that these estimations are subject to several biases. First of all, the actual density of field objects, at a given sensitivity, is significantly greater in the outskirts rather than at the field center, due to the obscuring cloud. Indeed, using the same procedure as above (plots not shown), we find that: (a) the density of all Gaia-PS1-2MASS-WISE objects is ∼2–3 times larger in the outer region than in the center¹⁰ (b) the density of field objects, here taken as all those that meet any of the criteria available for the full field, is ∼5 times larger in the outer region wrt. the center of the field. This is particularly relevant for sample c-wide, which, being based on full-field data, is our simplest case. Because of the above, the contamination estimated from the outer regions for sample c-wide is most likely an upper limit.

Samples “b” and “c”, however, are selected using also datasets which are limited to the field center, and with varying coverage. Indeed, the total density of objects in our catalog is ∼7 times higher in the center than in the outskirts¹¹. Consequently, the number of both membership and field-object criteria used is greater in the center of the field than in the outer regions, which is where we estimate the density of contaminants. It is not easy to predict what the outcome in terms of contaminating fraction is: while each additional membership criteria may introduce some spurious member candidate, additional field-object criteria help reduce contamination. Overall, however, the overwhelming majority of the contaminants are expected to originate from the proper motion and parallax criteria, based on Gaia data, which are available for the full field. Fewer contaminants are actually expected from objects that are prominently missing from the Gaia catalog, such as embedded objects (toward the field center) selected as members through X-ray detection and-or mIR excesses. We thus argue that our estimations of contaminants, ranging from 5.7% for sample “b” to 10.2% for sample “c-wide”, are reasonable approximations and, possibly, upper limits.

4.3. Spatial extension and structure

We are now in the position to discuss the extension of NGC 2264 in the plane of the sky, starting from the density maps for samples “b”, “c-wide”, and “c” shown in Fig. 5. In spite of the differences in completeness and-or contamination levels of the three samples, in all the maps we tentatively identify a new extended structure, marked by the outer polygon (dotted). Crucially, the region is visible even for the “c-wide” sample which, if anything, is likely to have a bias against including embedded stars toward the field center. We name this structure the “Extended Halo” to distinguish it from the Halo, which is the smaller elliptical region within the polygon. We refer to the region outside the extended halo, as the “Field”.

The map in Fig. 7 shows the center of the cluster, roughly within the halo, using sample “c”. In this case the binned density map was smoothed using the adaptive kernel smoothing algorithm of Ebeling et al. (2006)¹². In addition to the structures already identified by Sung et al. (2008) we tentatively identify three new substructures: two within the “S Mon” region and one within the Cone(C) region. The first two, named “S Mon(C)” and “S Mon(ref)”, correspond to two moderate but observable density enhancements; subregion S Mon(C), where C stands for core, is centered on the bright S Mon star, while subregion S Mon(ref) is named after the reflection nebula which is prominent at its location in optical images. The third new subregion, Cone(C-IR), lies within Cone(C) and corresponds to a compact member overdensity, which is very prominent in nIR and mIR images including a large fraction of young and embedded objects (also called IRS1).

Fig. 7. Spatial distribtion for sample “c” (center), as in the rightmost panel of Fig. 5, zoomed in by a factor of four (a 37.5 arcmin square). We here employed an adaptive kernel smoothing (see text).

In addition to these regions, close inspection of the density maps in Fig. 5 may also reveal hints of two filamentary structures: a vertical one close to the northern edge of our field, at RA ∼ 100.3 deg, which corresponds to the G02.3+2.5 region, and another one departing at a ∼45 degree angle from the western edge of the “halo” (Dec ∼ 9.5 deg) and reaching the edge on the extended halo. We do not further discuss these possible filaments.

In the following, in order to investigate how SF proceeded in the region, we characterize the properties of the stellar populations that lie, in projection, within all of these regions. Field-object contamination, which varies from region to region, can adversely affect our results. Adopting the average sky densities from Sect. 4.2, it is however small for all but the outermost subregions. Table 4 lists, separately for samples “b”, “c-wide”, and “c”, and for each subregion, the number of candidate members within the region, the expected net number of members after subtracting the estimated contamination, and the fraction of contaminants. For the subregions that contain nested subregions, we provide, separately, figures referring to the areas that exclude the inner subregions, and to the whole area within their respective perimeters (indicated by an asterisk following the region name). The indentations of the region names in the first column of the table indicate their spatial hierarchy: regions with a given indentation level contain those listed below it with a higher indentation level.

Contamination is clearly an issue mainly for the “Field” subregion (excluding the Extended Halo) where estimates range between ∼30% and ∼35%, depending on sample. With these levels of contaminations, moreover subject to uncertainties, characterizing this population is difficult and more data, such as future Gaia releases, might be needed to determine to what degree the cluster extends in these outer regions. On the other hand, the population of the other regions appear to be relatively unpolluted. Even for the Extended Halo^*, our estimates for the contamination are only 1–2%, rising to 3–5% when excluding the inner Halo^* region. The small “core” regions (Cone(C)^*, Spokes^*, S Mon(C)^*, and S Mon(ref)^*) all have negligible contaminations. This analysis also supports the use of sample “c” for the investigation of the mean (or median) properties of each subregion.

5. Discussion

We now proceed to exploit our membership determinations, described above, for our science goals. Our member samples are significantly richer-more complete with respect to previously available ones. We may, for example, compare our sample “c”, comprising 2257 stars, with that assembled by the CSI project (Cody et al. 2014). That sample, based on a early version of our current catalog, listed 1444 “very likely” members. Of these, 322 stars are not included in our sample “c”¹³, which, however, includes 1135 new candidates. Overall, we thus almost doubled the known population of NGC 2264.

Before focusing on our study of the SF process in the region and on the early stellar evolution, we make use of our selected member sample to address the usefulness of the X-ray detection and optical variability as proxies of stellar youth and, thus, as membership indicators for SFRs.

5.1. X-ray detection and optical variability

We can exploit our final member selection, sample “c”, to investigate the effectiveness of individual membership criteria. We here address X-ray detection, which is well-established, and “strong” optical variability, based on CoRoT lightcurves (see Sect. 2 and Fig. 2), for each of our subregions separately¹⁴. Selection efficiency and contamination are both estimated for each of the two criteria, the former being defined as the fraction of sample “c” members selected by the criterion, and the latter as the fraction of selected candidate members that are not included in sample “c” (that is of those that were rejected as nonmember).

For these statistics, in order to take into account the limited sensitivity of our X-ray and CoRoT observations, we restrict all the involved samples to stars brighter than a given limiting optical magnitude. From the G versus BP − RP diagrams of counterparts to X-ray sources (Fig. 4) we set this limit for X-ray detected stars to G < 19. We note, however, that the depth of the X-ray observations varies considerably across the field. For the variability-based membership we adopt G < 17 as the completeness limit, from a similar plot for stars with available CoRoT lightcurves (not shown). We caution that the X-ray detections were actually used to define our reference member list (sample “c”); assuming that our reference sample is reasonably complete, down to the adopted limiting magnitude, this should not overly bias the statistics. On the other hand, the CoRoT data were not used to define sample “c”. However, the targets observed by CoRoT were selected in advance and, in the central part of the field, roughly within the Halo region, are biased toward known NGC 2264 members. This may well bias the estimated contamination fractions in these regions, but should not affect the selection efficiency.

Table 5 lists, for each of our subregions (including the whole considered area, that is the Field^*): the number of candidate members in the samples (X_N and CoRoT_N), the selection efficiencies (X_sel.f and CoRoT_sel.f) and the contamination fractions (X_rej.f and CoRoT_rej.f), both for X-ray detection and for the variability analysis based on CoRoT data. For X-ray detection the selection efficiency ranges from 81–92% in the dense Cone, Cone^*, Cone(C), and Spokes regions, which are deeply imaged in X-rays, to ∼50% in the Halo and Extended Halo regions, which are instead mostly or exclusively covered by the relatively shallow XMM-Newton observations (cf. Fig. 1). Indeed, Table 5 also shows the result of limiting the sample to G < 17, instead of G < 19, yielding higher and more consistent selection efficiencies. Contamination from objects rejected as nonmembers by our selection procedure is typically < 10% in the densest regions, but reaches as high as ∼32% in the low-density Extended Halo.

From the selection efficiencies listed in Table 5 we infer that member selection based on optical variability (Sect. 2) is as effective as X-ray detection. As for contamination fractions, because of the mentioned target selection bias, the only reliable figure is probably that for the Extended Halo region, in which no members were known at the time of the CoRoT target selection. Although relatively high, 42%, it is comparable with that for X-ray detection in the same region (32%). Even considering that these results are based on a specific analysis of CoROT lightcurves, having specific and favorable characteristics in terms of cadence and duration, these estimates suggest that optical variability may prove a powerful method for selecting young stars in SFRs, and may be particularly useful when applied to forthcoming large datasets, such as those from Gaia and Rubin-LSST (see, for example, Bonito et al. 2018).

5.2. Physical properties

We now investigate the properties of the NGC 2264 population, both globally and for each of the subregions we have defined. As described in Sect. 4.3, we have tentatively identified four new substructures in addition to those of Sung et al. (2009): the Extended Halo, significantly enlarging the known extension of the cluster in the plane of the sky (×2.7 the area of the Halo), two substructures within the S Mon region (the S Mon “core”, S Mon(C)^*, around the bright O-type star, and S Mon (ref)^*, an apparently looser subclustering, southwest of S Mon, corresponding to a reflection nebula in the optical images), and the Core(C-IR)^* region within Cone(C), corresponding to a clear compact density enhancement, also prominent in nIR and mIR images. In the following we mostly ignore our “Field” subregion which, as discussed in Sect. 4.3, may or may not contain a significant population of NGC 2264 members. We begin assessing the spatial and kinematic properties of these structures and then discuss stellar masses, ages, and circumstellar disk frequencies. We will often refer to Table 5, which lists several average or median properties of stars in each subregion.

5.2.1. Spatial distribution and densities

Our substructures span a wide range of average stellar densities (cf. Table 5), a factor of 28(55) between the Extended Halo and the Spokes^* cluster (Cone(C-IR)^* region). Figure 8 shows the density profile within the Extended Halo^* region, computed in elliptical annuli concentric with the Halo perimeter (cf. Sect. 4.3 and Fig. 6). Two profiles are shown, both including and excluding the denser S Mon^* and Cone^* regions. A break in the density profiles, at a radius similar to the boundary between the Halo and the Extended Halo, possibly suggests that the two structures might indeed be distinct, in spite of the fact that most of the other average stellar properties, discussed below, appear to be compatible. Note that the density of the plateau at radii larger than 45′, at ∼0.08 stars/arcmin², is ∼10 times the density of the field-object contaminants, as estimated from Fig. 6.

Fig. 8. Density profile of sample “c” members within the Extended Halo* subregion. Densities plotted as a solid line are computed in the intersection between the Extended Halo* region and elliptical annuli, concentric with those used in Fig. 6. The dotted line is the result of excluding the dense S Mon* and Cone* region. The vertical red line indicates the boundary of the Halo region.

The external parts of the two main subclusters, the S Mon and the Cone regions, have similar densities, ∼1 star/arcmin², but the two cores within the Cone region are slightly denser than those in the S Mon region (4.0–4.8 versus 2.8–3.8 stars/arcmin). The densest structures are indeed the Cone(C-IR)^* and Spokes^* subcluster(s) dominated by an embedded and presumably very young population (9.3 and 4.8 stars/arcmin², respectively on average, with a central peak of ∼30 stars/arcmin² in both cases).

5.2.2. Distances

We estimate the distance to the NGC 2264 region from the median of the individual Gaia parallaxes of 1337 members (in our sample “c”) within the Ex.Halo^* region. We choose not to use the other candidate members, in the “Field” region, because of the potentially significant contamination from nonmembers. Our best median distance estimate is thus of 724±8 pc. We also estimate a distance of 722±2 pc, excluding 211 “outliers”, with parallaxes outside of the 1.1–1.7 mas range¹⁵.

In order to investigate possible differences in distance across the region, Fig. 9 shows a binned color-coded map of the region, with each pixel indicating the average parallax of the stars within that pixel (most pixels in the outer sparse regions contain only one star). Because of the relatively large uncertainties on parallaxes, we limit the sample to stars for which these uncertainties are smaller than 0.05 mas. No outstanding pattern is visible across the whole region. However, the bulk of the stars in the Halo^* region share a similar parallax (∼1.4 mas, in green), while the stars in the Field region seem to be either farther away (blue) or closer (red). Within the Halo we also note that the Cone^* and the S Mon^* regions contain several clustered blue and black areas, with the main concentration corresponding to the embedded Spokes^* cluster.

Fig. 9. Spatial distributions for sample “c”, color-coded according to four different variables. The top row refers to our full FoV, while the bottom row shows the central area in better detail. Within each row, the panels show, from left to right: parallax, PMRA, PMDec., and radial velocity. Color are assigned using the average value of the variable for all stars within each pixel and following the scale on the right-hand side of each panel. Most pixels on the outskirts contain a single star. The numbers on the bottom left indicate the number of stars used for the panel.

Figure 10 shows the median distance for all the defined subregions, here computed excluding outliers (parallaxes outside of the 1.1–1.7 mas range) as a function of average stellar density. The distances of most of the subregions are compatible with the adopted NGC 2264 distance, that is that of the Ex.Halo^* region, indicated in the figure by the filled gray circle. The most notable exception is the Spokes^* regions, which seems to be ∼20 pc farther away than the rest of the cluster. When considering the statistical uncertainties on the medians, the discrepancy between the Spokes^* and Ex.Halo^* distances is significant at the 2.1σ level, rising to 2.7σ when comparing with the Ex.Halo region (the outer population) which might be at a slightly closer distance.

Fig. 10. Distance vs. average stellar density for all the defined subregions. Distances are computed as the median of all stars with Gaia parallaxes in each subregion, excluding outliers (see text). Error bars refer to uncertainties on the median values. Filled symbols refer to the whole area within a region perimeter, while empty symbols to “haloes”, i.e. excluding nested regions, if present. The color coding is as in the legend.

Similar conclusions can be drawn comparing the distributions of the distances of individual stars within each subregion, again excluding distance outliers, using two-sided Kolmogorov-Smirnov (KS) tests. The probability that the distribution of distances for stars in the Spokes^* and Ex.Halo^* regions are drawn from the same parent population is 0.52%, which is reduced to 0.02% when comparing the Spokes^* and Ex.Halo (outer population only) regions. Similar but slightly less significant null probabilities, 1.4% and 0.03%, are obtained when also including parallax outliers, but excluding values with uncertainties > 0.3 mas.

We also note that these results are obtained in spite of the small fraction of available parallaxes among the members in the Spokes^* region: 33% and 36% for the no-outlier and uncertainty-restricted samples, respectively, significantly lower than those for the corresponding samples in the Ex.Halo^* region, which are 63% and 64%, respectively. This is easily understandable since we expect that the most embedded stars, which are optically faint, are less likely to have measured and good-quality Gaia parallaxes. The true median distance of the Spokes^* region might thus be even larger than estimated with the available parallaxes.

We conclude that the Spokes^* cluster is very likely behind the main NGC 2264 population. Given the uncertainties, the exact depth separation is only loosely constrained. It is, however, of the same order as the extent of the star-forming molecular cloud in the plane of the sky, as traced by the population we have selected: ∼1 deg or ∼12.5 pc.

5.2.3. Kinematic properties

The second, third, and fourth panels of Fig. 9 show binned spatial maps of members, color coded according to the mean RA and Dec proper motions, and radial velocity, respectively. The proper motion scale was converted to km s⁻¹, assuming a distance of 722 pc. Patterns and trends are clearly observable, but obtaining a clear picture of the structures and the nature of their motion is not straightforward. Velocities do not appear to be strictly associated with the subregions we have defined on the basis of the spatial distributions. We tentatively recognize signs of expansion in all three directions, more markedly in the north-south direction, consistent with the overall elongation of the cluster. It also seems that the Halo region is expanding toward the viewer with respect to the denser inner subregions (Cone^* and S Mon^*) at 2–3 km s⁻¹.

Rotation is also a possibility, mostly in the plane of the sky. On average, the southern region (Cone^*) moves toward the west at ∼0.5 mas yr⁻¹ (1.7 km s⁻¹) relative to the northern S Mon^* region. This effect can be significant, as it translates to 25 arcmin (5.5 pc) in 3 Myr (the approximate age of the cluster). The region might thus have rotated significantly since collapse started.

Figure 11 shows the evolution of stellar density in the last million years as derived by simply back-tracing the coordinates of members with measured proper motions. This extrapolation assumes rectilinear motions and no intervening SF. Only stars with PM uncertainties < 0.4 mas yr⁻¹ are used and, in the four panels representing epochs from 1 Myr ago to the present time, the red cross in the upper left corner indicate the maximum 1σ uncertainty in individual stellar positions. We observe the overall counter-clockwise rotation discussed above. Probably more interestingly, we see clear signs of collapse: two structures, at least one with a filament-like shape, appear to fall toward each other in the southern regions of the cluster, forming the present day substructures we named Cone(C) and Stokes. Both of these regions contain a rich embedded and extremely young population, very likely much younger than 1 Myr. We argue that the compression, or density increase, due to the collapse of these two structures is the cause of this recent formation outburst.

Fig. 11. Spatial densities predicted at three past epochs: 1, 0.67, and 0.33 Myr (from left to right) in the past, plus the present one (on the right). As indicated in the panels, only stars with PM uncertainties < 0.4 mas yr−1 are used and the maximum 1σ uncertainties for the reprojected position are indicated by the red cross in the upper-left corner. The density maps are smoothed with an adaptive kernel as in Fig. 7.

This scenario is somewhat similar to that inferred by Montillaud et al. (2019a) for G02.3+2.5, ∼0.2–1.0 deg north of S Mon, from the radial motions of its gas and dust filaments. Although the region is for the most part included in our field, we do not have evidence of this collision from stellar motions on the plane of the sky. It might be more evident in stellar radial velocities, which are not available in this region or, alternatively, we might simply have too few stars with proper motions in this region to notice the effect, especially since many young stars may be embedded and lack good Gaia astrometry.

5.2.4. Masses and ages

We estimated masses and ages for our candidate members using the MIST (v. 1.2) evolutionary tracks (Choi et al. 2016; Dotter 2016). Four different sets of masses and ages were derived by interpolating the evolutionary tracks for rotating (v/v_crit = 0.4), solar abundance stars, in different planes: photometric CMDs from the IPHAS, Pan-STARRS, and Gaia surveys (i versus r − i, for the first two surveys and G versus BP − RP for Gaia), and the fundamental plane (L_bol versus T_eff). For this latter, spectral types (available for 433 members in sample “c”) were used to estimate effective temperatures, bolometric corrections, and, together with photometry, extinction values for each star¹⁶. These individual extinctions values, were also used to derive mean values, which we adopted to de-redden the stars in the photometric planes before interpolating the tracks.

The four sets of masses and ages were then compared among each other, obtaining reasonable agreements: the 1-σ logarithmic dispersion between all the different pairs of estimates range from 0.06 to 0.14 dex for masses and between 0.18 and 0.40 dex for ages. We finally merged the four estimates, assigning the ‘best’ available values to each star, with the following precedence with respect to the origin of the estimates: the L_bol versus T_eff spectro-photometric plane (with individual extinction corrections), Pan-STARRS, IPHAS, and, finally, Gaia photometry. While we expect these estimates to be reasonable for the majority of NGC 2264 members, we warn that, for unresolved binaries and for heavily absorbed or embedded stars without spectral types, and therefore with no individual extinction estimates, estimated masses and ages might be rather biased. An additional small bias might also be introduced by variation in distances, for example making the embedded and more distant stars in the Spokes region look fainter and, therefore, older than they really are.

Figure 12 shows the mass and age distributions for the whole sample “c” (within the Extended Halo) and for some of the subregions. No striking differences are seen for masses. However, KS tests suggest, with a ∼99% confidence, that the S Mon(C) subregion might harbour more massive stars when compared to the Halo, the S Mon, and the Cone^* subregions. This is in spite of the fact that the most massive star in the region, the O7 star S Mon, is not included in the mass distribution.

Fig. 12. Logarithmic mass and age distributions (left and right panels, respectively) for selected subregions, color coded as in the caption in the upper right corner. For regions with nested subregions, dashed lines refer to the whole area within their perimeters, and solid lines to the ‘haloes’. Thick solid lines in the left panel show the Salpeter’s Initial Mass Function (IMFs, Salpeter 1955, slope −1.35) in orange, and the Chabrier (2003) system and single stars IMF, for M < M⊙, in cyan and magenta, respectively. The short vertical lines near the x-axes, at log M ∼ −0.5 and log Age ∼ 6.5, indicate the median values of each samples.

The age distributions of most subregions are statistically compatible, when comparing them using KS tests. The exception seems to be Cone(C), which does appear to be somewhat younger than other subregions, with null probabilities < 1% when compared to the Ex. Halo^*, the Ex.Halo, the Halo, and S Mon(C)^* (n.p. = 0.9%, 0.8%, 0.2%, and 0.3%, respectively, when comparing the distributions of ages plotted in Fig. 12)¹⁷. For some of the subregions, notably for the Spokes^* and the Cone (C-IR)^* regions (not shown in Fig. 12), the fraction of stars with age estimates is small, 41% and 33%, respectively, and their median ages are likely overestimated. Indeed, their embedded and likely young populations may well be under-represented, and, as noted above, the ages that are available might be adversely affected by the unaccounted-for larger extinctions and distances. Moreover, all of our age estimates may be affected by the presence of accretion disks, through the luminosity excesses and color shifts produced by accretion spots (in the blue part of the spectrum or photometric bands), thermal emission by inner disks (in the red bands), and scattering of stellar emission by the disks. Since, as shown in the next section, the fraction of stars with accretion disks varies significantly in different subregions, even the relative age scale might be skewed. We thus consider our median age estimates with extreme caution.

5.2.5. Disk and accretion fractions

We investigate differences in disk and accretion fractions across the SFR, using both binned spatial maps and focusing on the subregions we have defined. We define the subsample of members with disks as the stars that fall in one of the two mIR-based member loci ([3.6]–[5.8] versus [4.5]–[8.0] and W1 − W2 versus W2 − W3, in Fig. B.2), with 2σ confidence¹⁸. Likewise, we define as accreting members those that lie, with 2σ confidence, in at least one of the five members loci based on Hα emission, as observed either photometrically or spectroscopically (Figs. 4 and B.2). Note that having independent accretion assessment is useful, not only to enlarge the sample of stars with appropriate data, but also to account for the variable nature of accretion.

Figure 13 shows, in the left-hand panel, the binned map of candidate members with sufficient mIR data to determine disk presence; the central panel refers to the subset of stars with disks, and the one on the right shows the fraction of stars with disks in each spatial bin. Although the former two maps cover the whole region, in order to reduce Poisson uncertainties, we only show the disk fractions for the spatial bins containing a minimum of 5 members per bin (in the left-most panel). The map, which is therefore only defined within the dense cores, shows that the Cone region has larger disk fractions than the S Mon region. Figure 14 similarly shows maps for accreting stars. Note that the heavily absorbed or embedded regions, such as the Spokes cluster, are quite clearly under-sampled by optical observations and that, consequently, their accretion fraction are likely underestimated.

Fig. 13. Binned spatial maps of (from left to right): i) stars placed in either of the two mIR based diagrams from which we assess the presence of disks, ii) stars with disks, iii) the fraction of stars with disks. In the rightmost plot the gray scale follows the color-bar on the right-hand side, and the spatial bins in which the denominator of the fraction contained less than 5 stars are plotted in cyan.

Fig. 14. Same as Fig. 13, but for the accretion indicators. The leftmost panel shows stars placed in any of the five Hα based diagrams, the central one those lying in any of the five accretion loci, and the rightmost one the ratio between the former two maps.

Highly significant differences are seen when comparing the disk fractions of our subregions: from 80–88% for the Spokes^* and Cone(C-IR)^* regions to ∼27% for the S Mon(C) region. This is shown in Fig. 15, where we plot the disk fractions versus the mean stellar densities of the subregions. A correlation is seen, with the possible outlier of the S Mon region, which appears to have a factor of ∼2 lower disk fraction with respect to the general trend. An even clearer correlation is seen when plotting accretion fractions (instead of disk fraction) versus density (not shown, data in Table 5), again with a noticeable ∼30% reduction in accretors fraction for the S Mon region.

Fig. 15. Mid-IR excess vs. stellar density for independent regions. Symbols and colors as in Fig. 10. The numbers to the right of the region names are the relative fraction of stars for which the disk indicators are available.

Indeed, the accretion fractions estimated from Hα emission correlate well with the disk fractions from the mIR data (Fig. 16). Spearman’s ρ and Kendalls’s τ tests give 0.5–1.2% probabilities that the two fractions are uncorrelated. The Spokes^* region seems to have somewhat fewer accretors than stars with disks, which might be explained as a selection effect, that is by a lower fraction of available Hα data among the embedded (young) stars. Similarly, the deficit of accreting stars in the Ex. Halo might be attributed to the fact that the outer regions are covered by only one Hα survey (IPHAS) while the central region is covered by up to 5 datasets. When using the same samples, that is stars with both disk and accretion information, the correlation between the two fractions becomes even more significant, in spite of the resulting larger uncertainties (null probabilities: 0.16% and 0.67%). Moreover, all accretion and disk fractions become compatible at the ∼1.5σ level. We can conclude that our data is insensitive to possible differences between the timescales for the evolution of accretion and disks.

Fig. 16. Disk fractions vs. accretion fractions for independent subregions. Symbols and colors as in Fig. 10. The two numbers beside each region name are the relative fraction of stars for which indicators of accretion (left) and disks (right) are available. The identity line is shown in gray and the results of Spearman’s ρ and Kendalls’s τ correlation tests are shown in the bottom right.

We tried to correlate disk and accretion fractions with stellar ages but, unfortunately, as discussed (Sect. 5.2.4) median ages for the subregions are affected by large uncertainties, both statistical and systematic. Disregarding the two heavily embedded regions, for which unaccounted-for extinctions and selection effects probably dominate, the subregion with the highest disk and accretion fractions, Cone(C), might indeed be the youngest, while the S Mon(C)^* region, the one with the lowest disk and accretion fractions, might be the oldest (cf. Fig. 12 and KS tests discussed in Sect. 5.2.4), as expected but rarely, if ever, observed within a single star-forming region. However, we do not report this as a result, as we are concerned that age determinations might be affected by the presence of disk- and accretion-related excess, indeed spuriously producing such a correlation.

6. Summary and conclusions

Following the recognition that a significant number of stars far from the known NGC 2264 members show strong optical variability, similarly to the known members, we embarked on a novel investigation of the structure of this well studied star-forming region. First, we have obtained new X-ray imaging data, collected with XMM-Newton, significantly enlarging the area in which X-ray sources, which are likely to be young stars, have been detected. We have reduced and analyzed these original data, along with previously published X-ray observations, in a consistent manner. Our list of XMM-Newton sources appreciably enriches the known population of X-ray emitting stars in the region. Then, we have collected a large catalog of photometric, spectroscopic, and astrometric data for an even wider field, a 2.5×2.5 degree area centered on the known cluster. Most notably, the Gaia eDR3 astrometric data have been extremely useful, albeit not per-se sufficient, to derive a reasonably accurate census of NGC 2264 members. Optical and IR photometry have also been needed to clean-up the astrometrically selected candidates from field stars. Spitzer and WISE photometry have allowed us to identify stars with a mIR excess and therefore likely surrounded by circumstellar disks. We have derived indication of ongoing mass accretion from narrow band H_α photometry, as well as from observation of the same line with low and medium resolution spectroscopy, these latter data being limited to the central regions.

We have assembled a catalog of more than ten-thousand candidate members, according to at least one of several criteria, some subject to significant field-star contamination. Out of this, we have extracted several less contaminated samples, the most useful of which, our sample “c”, comprises 2257 candidate members, less than 8% of which are expected to be contaminants, mostly located in the outskirt of our search area. With this catalog, almost twice as large than previous determinations, we reach a number of conclusions.

– We demonstrated that the optical variability of young stars, a manifestation of magnetic activity (like the X-ray emission) and of circumstellar accretion, is as powerful a membership criterion as X-ray emission. The effectiveness of X-ray member selection was indeed confirmed, even in the sparse outskirts of the region where, however, further membership confirmation is useful to reduce field-star contamination.

– We have defined four new substructures with respect to Sung et al. (2008), based on the surface density map of candidate members. The compact S Mon(C) region may be a physically recognizable structure, being possibly older than the surrounding region and-or having a distinctly lower fraction of stars with disks and undergoing accretion. The embedded Cone (C-IR) is, instead, a recently formed embedded subcluster. The Extended Halo and the S Mon(ref) region, on the other hand, do not constitute obviously coherent structures, either kinematically or physically.

– The identification of the population in the Extended Halo significantly enlarges the known extent of NGC 2264. The cluster is likely even broader, but a less contaminated membership is needed in order to characterize the outer population.

– NGC 2264 is certainly not a dynamically relaxed cluster. We see signs of ordered bulk motions such as expansion and rotation. No mass segregation, neither spatial or in velocity-space, is clearly observed, although stars very close to the O-type star S Mon might be, on average, slightly more massive than elsewhere. Notably, we observe the collapse of filamentary structures, which may have recently triggered new SF activity in the southern regions.

– Interpreting disk and accretion fractions as proxies of stellar ages, stars in the southern regions are younger than those close to S Mon in the north. We thus speculate that SF started ∼4 Myr ago in the S Mon region and extended progressively toward the south, where it is presently continuing owing to the aforementioned collapse.

– The observed trend of the increasing disk fraction with stellar density confirms previous suggestions that stars are preferentially born in high density regions that later disperse as they age and dissipate their disks. The lower disk and accretion fractions of the S Mon(C)^* region might be due to disk evolution, implying a slightly older age with respect to other regions, or to the strong UV radiation of the O7 star evaporating the disks of close-by stars. If the region were older, its high density (for the age) could be explained by a higher-than-average density in the original cloud, or by the gravitational pull of S Mon. A more precise estimation of stellar ages is required to choose between the two possible explanations for the lower disk and accretion fractions.

Our study confirms that NGC 2264 is one the best templates to investigate the formation mechanisms of stars and clusters. More precise astrometric data, as those foreseen for the future Gaia data releases, along with future spectroscopic characterization of the new candidate members found here, will allow us to obtain a more complete picture of the present status of the region and of its recent evolution, to be compared with theoretical models of cloud collapse and of star-disk evolution.

Acknowledgments

E.F. would like to thank Francesco Damiani for suggesting to trace the position of stars back in time. R.B. acknowledges financial support from the project PRIN-INAF 2019 “Spectroscopically Tracing the Disk Dispersal Evolution” and the “Preparing for Astrophysics with LSST” Program, funded by the Heising Simons Foundation through grant 2021-2975, and administered by Las Cumbres Observatory.

References

Appendix A: Chandra ACIS sources

Table A.1 presents basic information on X-ray sources detected on the available Chandra ACIS data for NGC 2264 (Obs.ID: 2540, 2550, 9768, 9769, 13610, 13611, 14368, and 14369). We list: source number, coordinates, positional uncertainties, and background-subtracted source photons. Sources were detected with PWdetect (Damiani et al. 1997), and analyzed using the ACIS-EXTRACT package (Broos et al. 2010). The table is available in its entirety on-line. A full description of the data and analysis will be provided by Flaccomio et al. (in preparation).

Appendix B: Additional diagrams used for membership

Figure B.1 shows, for all the objects in our catalog, the 19 plots we used for selecting likely members and field objects, in addition to the four already shown in Fig. 4. The following figure, Fig. B.2, shows the same plots for our final and most-used selection of likely members (sample “c”).

Fig. B.1. Nineteen additional diagrams used for membership determination in addition to those in Fig. 4. All objects in our catalog are plotted. See Fig. 4 for a description of lines and symbols.

Fig. B.2. Same as Fig. B.1 for sample “c”

Within each of these figures, the first four panels in the upper row show indicators of the strength of the Hα line versus stellar color: the first two use narrow-band photometric indexes similar to the IPHAS-based one in Fig. 4 (Sung et al. 2008; Lamm et al. 2004), the third uses the spectroscopic Hα equivalent width (EW, from a variety of sources), and the fourth the Full Width at Zero Intensity (FWZI) values derived by Bonito et al. (2020) from the spectra obtained by the GAIA-ESO survey (see their Table 1; if more thans one spectra were available for a single source, we adopted the maximum measured value of FWZI). For this latter plot and for those using the photometric H_α indexes, we define three loci to distinguish low magnetic activity, characteristic of old stars, strong magnetic activity, expected from young PMS and MS stars, and strong accretion-related line emission, taken as indication of membership. This latter locus is also defined for the spectroscopic Hα EW versus V-I plot. The first panel in the second row shows spectroscopic radial velocities versus GAIA G magnitudes. We then plot ten different color-magnitude diagrams (CMDs). In addition to the GAIA G versus BP-RP diagram in Fig. 4, we also plot the absolute G (G - distance modulus, DM) versus BP-RP, restricted to objects whose parallaxes have S/N> 5. We then show: i versus r-i¹⁹, I versus R-I, I versus V-I, J versus J-H, K versus H-K, [3.6] versus [3.6]-[4.5], [8.0] versus [4.5]-[8.0], [24.0] versus [4.5]-[8.0], and [3.6] versus [3.6]-[24.0]. In all of these CMDs we define loci from which we select field object. Finally, we plot four color-color diagrams, J-H versus H-K, [3.6]-[5.8] versus [4.5]-[8.0] (Spitzer), W1-W2 versus W2-W3 (WISE), and i-z versus [3.6]-[4.5] (PS1 and Spitzer). We define member loci, based on the IR excesses produced by circumstellar disks, in all four diagrams, and field-object loci in the first three.

All Tables

All Figures

Fig. 1. Optical DSS2 image of the 2.5×2.5 deg square field centered on NGC 2264 to which our study is limited. The locations of two prominent features, the O7 star S Mon and the Cone Nebula, are indicated by orange arrows. The areas shaded in blue and orange-brown are combined exposure maps of all available Chandra and XMM-Newton observations, respectively. Identifiers for the XMM-Newton field are also shown by red symbols in their respective lower-right corner. The two orange squares indicate the FoV of the CoRoT observations discussed in Sect. 2. The yellow circles and ellipse indicate the subregions identified by Sung et al. (2009). Contours of the four new spatial substructures identified in this study (Sect. 4.3), three inner ones and the outer polygon, are drawn in cyan.

In the text

	Fig. 2. Height of the most significant LNP peak vs. standard deviation of light curves normalized to the median flux value. The grayscale map represents the density of stars in the control field, with the cyan contour encircling 90% of these stars. Blue open points are known members of NGC 2264, while red points indicate the new candidate members.
In the text

	Fig. 3. Background light curve of EPIC PN-detected events in field 3 (Obs.ID 0721400201); the horizontal red line is the count-rate threshold that we chose to define and filter-out high-background time intervals, thus maximizing the S/N of faint sources.
In the text

Fig. 4. Object density in four different diagrams (one per column) and four different samples (one per row). The diagrams are, from left to right: the r-Hα vs. r − i color-color diagram using IPHAS data, G vs. BP − RP color magnitude diagrams (with three reddening vectors, at different colors), PMDec vs. PMRA, and parallax vs. G magnitudes (the last three diagrams using Gaia data). The four samples, from top to bottom, are: all objects in our catalog, objects within the XMM-Newton field of view, objects unambiguously detected in X-rays, candidate members in our sample “c” (Sect. 4.1). Density map are smoothed using a Gaussian Kernel. The grayscales range from black at maximum density to white at zero density, scaling with the square root of the density, rather than linearly, to reveal more details at low-densities. Isolated gray dots in low-density areas indicate individual stars. Dashed and solid lines delineate member and field loci, respectively, and, in the first row, “m” and “f” labels help clarify the distinction between the two. Dotted lines in the last two diagrams refer to more conservative member loci. The first number in the bottom-right corners indicate the total number of plotted objects. The following ones refer to the number of stars in the member- and field-loci, with 3- and 2-σ significance, respectively.

In the text

	Fig. 5. Spatial distributions for candidate member samples “b”, “c-wide”, and “c” (from left to right). The density maps are smoothed with a Gaussian kernel as in Fig. 4, and dotted lines show the boundaries of the regions discussed in the text (Sect. 4.3) and also shown in Fig. 1. The figures at the bottom right of each panel indicate the number of objects in the sample.
In the text

Fig. 6. “Radial” distributions for samples “b”, “c-wide”, and “c”. Different colors refer to different “quadrants”. Spatial distributions are shown at the top-right of each panel, with stars color-coded according to their quadrant. The elliptical annuli where the densities for the main plot were computed are also shown. Points with error bars, based on Poisson uncertainties, are slightly shifted along the x-axis to avoid overlap among the four independent samples. The dotted black line refer to the combination of the four quadrants. The thick horizontal magenta segments toward the bottom right indicate the average densities of stars in the left quadrant for r > 50 arcmin.

In the text

	Fig. 7. Spatial distribtion for sample “c” (center), as in the rightmost panel of Fig. 5, zoomed in by a factor of four (a 37.5 arcmin square). We here employed an adaptive kernel smoothing (see text).
In the text

	Fig. 8. Density profile of sample “c” members within the Extended Halo* subregion. Densities plotted as a solid line are computed in the intersection between the Extended Halo* region and elliptical annuli, concentric with those used in Fig. 6. The dotted line is the result of excluding the dense S Mon* and Cone* region. The vertical red line indicates the boundary of the Halo region.
In the text

Fig. 9. Spatial distributions for sample “c”, color-coded according to four different variables. The top row refers to our full FoV, while the bottom row shows the central area in better detail. Within each row, the panels show, from left to right: parallax, PMRA, PMDec., and radial velocity. Color are assigned using the average value of the variable for all stars within each pixel and following the scale on the right-hand side of each panel. Most pixels on the outskirts contain a single star. The numbers on the bottom left indicate the number of stars used for the panel.

In the text

Fig. 10. Distance vs. average stellar density for all the defined subregions. Distances are computed as the median of all stars with Gaia parallaxes in each subregion, excluding outliers (see text). Error bars refer to uncertainties on the median values. Filled symbols refer to the whole area within a region perimeter, while empty symbols to “haloes”, i.e. excluding nested regions, if present. The color coding is as in the legend.

In the text

Fig. 11. Spatial densities predicted at three past epochs: 1, 0.67, and 0.33 Myr (from left to right) in the past, plus the present one (on the right). As indicated in the panels, only stars with PM uncertainties < 0.4 mas yr−1 are used and the maximum 1σ uncertainties for the reprojected position are indicated by the red cross in the upper-left corner. The density maps are smoothed with an adaptive kernel as in Fig. 7.

In the text

Fig. 12. Logarithmic mass and age distributions (left and right panels, respectively) for selected subregions, color coded as in the caption in the upper right corner. For regions with nested subregions, dashed lines refer to the whole area within their perimeters, and solid lines to the ‘haloes’. Thick solid lines in the left panel show the Salpeter’s Initial Mass Function (IMFs, Salpeter 1955, slope −1.35) in orange, and the Chabrier (2003) system and single stars IMF, for M < M⊙, in cyan and magenta, respectively. The short vertical lines near the x-axes, at log M ∼ −0.5 and log Age ∼ 6.5, indicate the median values of each samples.

In the text

Fig. 13. Binned spatial maps of (from left to right): i) stars placed in either of the two mIR based diagrams from which we assess the presence of disks, ii) stars with disks, iii) the fraction of stars with disks. In the rightmost plot the gray scale follows the color-bar on the right-hand side, and the spatial bins in which the denominator of the fraction contained less than 5 stars are plotted in cyan.

In the text

	Fig. 14. Same as Fig. 13, but for the accretion indicators. The leftmost panel shows stars placed in any of the five Hα based diagrams, the central one those lying in any of the five accretion loci, and the rightmost one the ratio between the former two maps.
In the text

	Fig. 15. Mid-IR excess vs. stellar density for independent regions. Symbols and colors as in Fig. 10. The numbers to the right of the region names are the relative fraction of stars for which the disk indicators are available.
In the text

	Fig. 16. Disk fractions vs. accretion fractions for independent subregions. Symbols and colors as in Fig. 10. The two numbers beside each region name are the relative fraction of stars for which indicators of accretion (left) and disks (right) are available. The identity line is shown in gray and the results of Spearman’s ρ and Kendalls’s τ correlation tests are shown in the bottom right.
In the text

	Fig. B.1. Nineteen additional diagrams used for membership determination in addition to those in Fig. 4. All objects in our catalog are plotted. See Fig. 4 for a description of lines and symbols.
In the text

	Fig. B.2. Same as Fig. B.1 for sample “c”
In the text