© ESO 2020

1. Introduction

With the second instalment of the Gaia survey (DR2; Gaia Collaboration 2018), high-precision position and kinematic data became available for a large number of young clusters which previously lacked reliable parallax and proper motion measurements. Now, it is possible to probe the dynamical evolution and star formation history of these clusters through substructure, mass segregation, and relative dynamics studies of their young populations.

One of the fundamental questions in such analyses is to what degree stars “cluster” together and how does this change as the cluster evolves. The anwser to this question requires a combined study of the spatial intensity, correlation, and distribution of stars/clumps with the kinematic data. For this type of characterisation the use of local indicators (Anselin 1995) is suggested. Unlike global indicators (e.g. the two-point correlation function) that derive a single parameter for a group of stars as a whole, local indicators derive a parameter for each unique source such that variations and trends as a function of fundamental parameters (stellar mass, evolutionary stage, position, and individual dynamical histories) can be distinguished. Unfortunately local indicators have remained largely ignored in cluster analysis due to a distinct lack of appropriate astro-statistics tools, and the best understood methods from other fields cannot be easily applied to (or are simply invalid for) astronomical datasets.

Our aim with this paper series is the development and application of local statistic tools, optimised for stellar cluster analysis. In Paper I (Buckner et al. 2019) we introduced the tool INDICATE (INdex to Define Inherent Clustering And TEndencies) to assess and quantify the degree of spatial clustering of each object in a dataset, demonstrating its effectiveness as a tracer of morphological stellar features in the Carina Nebula (NGC 3372) using positional data alone. In this paper we demonstrate that when combined with kinematic data from Gaia DR2, INDICATE is a powerful tool to analyse the star formation history of a cluster in a robust manner.

Embedded in the Mon OB1 cloud complex, NGC 2264 is located at a Gaia DR2 determined distance of ${723}_{-49}^{+56}$ pc (Cantat-Gaudin et al. 2018). Structurally, the cluster is elongated along a NW-SE orientation with two sub-clusters C and D in the southern region (Fig. 1). There is an age spread of ∼3–4 Myr between the older star formation inactive northern region which contains the bright O-type binary star, S Mon, and the younger ongoing star formation southern region within the C and D sub-clusters (Mayne & Naylor 2008; Venuti et al. 2017). Recent studies have found NGC 2264 to be rich in young stellar objects (YSOs) of all evolutionary stages (Teixeira et al. 2012; Povich et al. 2013; Rapson et al. 2014; Venuti et al. 2018). Moreover, owing to its close proximity, this cluster is one of the best researched in the literature with numerous studies into its recent and ongoing star formation (for example Sung et al. 2009; Sung & Bessell 2010; Teixeira et al. 2012; Venuti et al. 2017, 2018; González & Alfaro 2017). As such we have chosen to focus our efforts on NGC 2264 as (i) the validity of our results can be checked against what is already known about the cluster, and (ii) its large YSO population makes this cluster an ideal candidate to show that INDICATE can successfully provide the rigorous spatial analysis necessary to validate and correctly interpret dynamical behaviours found with DR2 data for young clusters.

Fig. 1. Distribution of our sample overlaid on the Herschel 70 μm map of the region to clearly show the locations of S Mon and sub-clusters C, D. Colours/shapes: Class 0/I (red triangles), Class II (green plus signs), Class TD (purple stars), Class III/F (dark blue dots), Class II/III (grey crosses), Class Ambiguous (pink squares), and non-YSO (light blue pentagons).

The paper is structured as follows. We introduce our sample of YSOs in Sect. 2 and analysis methods in Sect. 3. In Sect. 4 we present the results of our spatial and kinematic analyses, which we discuss in Sect. 5. A conclusion is given in Sect. 6.

2. Sample of young stellar objects

In this section we describe our sample of young NGC 2264 members. The following terminology is employed to discriminate between the different evolutionary stages of YSOs:

• Class 0: protostars without dust emission.

• Class I: protostars with envelope and disc dust emission.

• Class II: Pre-main sequence (PMS) stars with circumstellar accretion discs.

• Class TD: transition discs; an intermediate stage between Class II and III where the disc has a radial gap.

• Class III: PMS stars without discs.

2.1. Catalogue selections

We draw our sample from two independent catalogues of the region. The first was constructed by Kuhn et al. (2014, hereafter K14) as part of the MYStIX project (Feigelson et al. 2013) which surveyed 20 OB-dominated young clusters using a combination of Spitzer IRAC (Fazio et al. 2004) infrared and Chandra (Weisskopf et al. 2000) X-ray photometry. The bounds of the ∼0.19° NGC 2264 region surveyed (Fig. 1) correspond to the limits of the Chandra mosaic coverage, as this is smaller than that of the Spitzer IRAC photometry. Identification and classification of YSOs were carried out by Povich et al. (2013) who used Spitzer IRAC, 2MASS (Skrutskie et al. 2006), and UKIRT (Lawrence et al. 2007) imaging photometry with spectral energy distribution fitting to flag sources as “0/I”, “II/III”, “non-YSO (stellar)” or “Ambiguous (YSO)”. The catalogue consists of 969 sources of which 139 are Class 0/I, 298 Class II/III, 413 non-YSO (stellar), and 119 Ambiguous (YSO).

The second catalogue is by Rapson et al. (2014, hereafter R14) who analysed 2MASS and Spitzer photometry of Mon OB1 East to identify YSO members. A three-phase classification method by Gutermuth et al. (2009) which utilised photometry in eight infrared bands (J, H, K, 3.6 μm, 4.5 μm, 5.8 μm, 8.0 μm, and 24 μm) was employed, resulting in 10454 potential candidates. We selected only sources in the smaller K14 region of which there are 1645 comprised of 70 Class 0/I, 307 Class II, 26 Class TD, 1189 Class III/F (where “F” denotes source is potentially a line-of-sight field star: see next section), and 53 contaminants (AGN, Shock, PAH).

We merged the two samples and performed a cross-match to remove duplicates, finding 1848 unique sources for the region. There is significant overlap between the two catalogues with a R14 counterpart for 766/969 K14 sources and a discrepancy in the assigned classifications of 24.4% (Table B.1). Interestingly, 101/119 of the duplicate sources flagged as Ambiguous (YSO) in K14 have a definitive classification from R14 (i.e. 0/I, II, TD, III/F, AGN, Shock, PAH). Sources classified by R14 form the majority of the merged sample so we adopted their classifications for all sources that appear in both catalogues. This is statistically justified as we are interested in the spatial behaviour of the YSO population as a whole, not individual sources. Thus (i) a single classification system should be utilised, where possible, to ensure the spatial analyses are systematic; and (ii) if an incorrect classification is assumed for an individual source, it would effectively be an outlier so does not have a statistically significant impact on our results as our methodology is robust against outliers (Sect. 3.1).

After removing the contaminants, our final sample contains 1795 sources. Table 1 details its composition and Fig. 1 its distribution. We create four sub-samples: S₁ – all sources (n₁ = 1795); S₂ – Class 0/I only (n₂ = 111); S₃ – Class II only (n₃ = 307); S₄ – Class III/F only (n₄ = 1189).

2.2. Field star contamination

The classification method employed by R14 distinguishes Class III sources from earlier type YSOs by their (lack of) 3.6 μm and 4.5 μm excess emission. However, as field stars in the line of sight also lack this excess they cannot be readily distinguished from true Class III cluster members. To estimate the number of field star contaminants, the authors calculated the expected number of field stars from comparison to two control regions neighbouring Mon OB1 East, determining a contamination of ∼29% in the region of NGC 2264. Thus 345 of the 1189 Class III/F sources in our sample are expected to be field stars and 844 Class III members.

For an independent measure, we cross-matched our sample with the Coordinated Synoptic Investigation of NGC 2264 (CSI 2264; Cody et al. 2014). Uniquely, CSI 2264 continuously observed the region for 30 days using Spitzer IRAC and the Convection, Rotation and Planetary Transits satellite (CoRoT; Baglin et al. 2006) simultaneously, with additional observations from 13 other telescopes. Subsequently the CSI 2264 photometric database is one of the most comprehensive of the region to date containing an impressive 146 855 sources. Membership of each source was assessed against photometric, spectroscopic, spatial, and kinematic criteria (see Appendix A.1 of Cody et al. 2014 for details) and flagged as “very likely member”, “possible member”, “likely field object” or “no membership information”. For our 1189 Class III/F objects, 775 are very likely or possible members, 91 are likely field objects, 320 have no membership information and 3 do not appear in the catalogue. Of the 320 with no membership information, 11 are identified by K14 as members of the cluster indicating a field contamination rate of 8 − 34%. This is in good agreement with the estimate of R14 of 29% which suggests the true contamination is towards the upper limit of this range. We therefore conclude the R14 contamination calculation is reasonable and assume this value for our analysis (but see Sect. 3.2).

2.3. Completeness of sample

We anticipate two sources of incompleteness in our sample. The first relates to the heavy extinction present in the cluster owing to its embedded nature, which can be seen in the Herschel 250 μm image of NGC 2264 shown in the left panel of Fig. 2. Despite compiling our sample from two catalogues of the region, in the absence of longer wavelength photometry it is reasonable to assume that they suffer from incompleteness due to dust extinction. In particular, we anticipate the majority of “missing” sources are located in regions of highest extinction (sub-clusters C and D) and for these to primarily be the most deeply embedded Class 0/I objects which have not been detected. Indeed, Class 0/I objects constitute only 16% of sources in the right panel of Fig. 2. To gauge how many Class 0/I sources are “missing” we consult submillimetre surveys of the sub-clusters (Peretto et al. 2006; Teixeira et al. 2007), and find there are at least 16 Class 0/I sources which are not included in our sample. We refrain from appending our sample to ensure it remains homogeneous (and thus results reliable) as these surveys only cover relatively small areas of NGC 2264’s southern region. However it should be noted that the inclusion of these highly concentrated sources would strengthen, rather than diminish, the trends found in Sect. 4 and thus our conclusions remain unchanged irrespective of our decision to exclude these objects.

Fig. 2. Left panel: sample overlaid on a Herschel 250 μm image of NGC 2264. White lines and crosses denote the borders and sources of the sample respectively. Right panel: zoomed-in plot of Fig. 1 in the NGC 2264-C and -D sub-cluster region, overlaid on a corresponding Spitzer 24 μm image. Colours and symbols as defined in Fig. 1.

The second relates to the point spread function (PSF) wings of a bright star at the centre of NGC 2264-C, as seen in the Spitzer MIPS 24 μm image (Fig. 2 right panel). There is significant angular dispersion of the PSF which likely occludes a number of fainter sources in 2MASS and Spitzer IRAC bands (from which both catalogues were derived).

2.4. Gaia DR2 kinematic data

We cross-match all 1795 sources with the Gaia DR2 database by colour and position, finding proper motions for 1268 sources and radial velocities for 20 sources (Table 1). As expected there is a distinct lack of Class 0/I objects with kinematic data in DR2 because the magnitudes of these deeply embedded objects are typically below the Gaia detection limit. Before proceeding we must consider the impact of systematic errors in the proper motion measurements caused by the Gaia scanning law as our sample occupies an area ≪1° and spatial scales < 1° are most affected with a root mean square amplitude of 0.066 mas yr⁻¹ (Lindegren et al. 2018). To ensure these errors do not dominate the measurements it is necessary to exclude any sources from our kinematic analysis in Sect. 4.3 with a proper motion (within error bounds) of < 0.066 mas yr⁻¹. A search of the sample reveals 29 sources that meet the exclusion criteria. We further exclude 672 objects with r_hi <  674 pc or r_lo >  779 pc, where r_hi is the upper error bound and r_lo the lower error bound on the distance estimate determined by Bailer-Jones et al. (2018) from DR2 parallaxes. The cut-off distance values correspond to the upper/lower DR2 distance estimate of ${723}_{-49}^{+56}$ pc for NGC 2264 (Cantat-Gaudin et al. 2018).

Two perspective corrections on the proper motions are needed prior to analysis: first, the radial motions of members cause NGC 2264 to appear to contract and, second, members appear to move towards a point of common convergence as they are part of the same stellar system and share a common motion. The former is corrected using Eq. (13) of van Leeuwen (2009) and the latter by subtracting the mean proper motion of the system from observed proper motions of each member. Appendix A describes these calculations.

Distance measurements from Gaia DR2 suggest a significant proportion (561/966) of Class III objects are not true members, which is considerably higher than the 29% identified from photometric analysis (R14). While this may reflect the number of true members of the cluster it may also be a symptom of a number of observational biases (e.g. imprecise/too strict distance criteria, unresolved binaries etc.). In addition approximately one-third of Class III objects in our sample lack parallax measurements from which to make a distance-dependant membership determination. As such, it is important to ascertain the effect of a significantly reduced Class III sample size on our results; for example whether the spatial trends found in the Sect. 4.1 and 4.2 hold for these objects identified by the kinematic data. To check, we re-ran our spatial analyses (as outlined in Sect. 4) excluding all Class III sources that did not meet our above discussed DR2 distance criteria and find our conclusions on the spatial behaviour of YSOs in NGC 2264 are unaffected by the exclusion of sources that fail the distance criteria.

Therefore as only 79.0% of Class III, 74.9% of Class II, and 1.8% of Class 0/I sources have reliable DR2 data we used the full sample with Monte Carlo sampling described in Sect. 3.2 for our spatial analyses. For our kinematic analysis we only used sources which met our DR2 distance criteria (Table 1) to ensure the proper motion patterns we observe for Class II and III objects are an accurate reflection of typical member motions.

3. Analysis method

3.1. INDICATE

We analyse the spatial distributions of YSOs using INDICATE (Buckner et al. 2019). A statistical clustering too, INDICATE is used to study the intensity, correlation, and spatial distribution of point processes in 2+D discrete astronomical datasets. It is a local statistic which quantifies the degree of association of each point in a dataset through a comparison to an evenly spaced control field. Advantageously, INDICATE does not make assumptions about (or require a priori knowledge of) the shape of the distribution, nor the presence of any substructure. Extensive statistical testing has shown it to be robust against outliers and edge effects, and independent of the size and number density of a distribution (Buckner et al. 2019).

When applied to a dataset of size n, INDICATE derives an index for each data point j, defined as

$$\begin{array}{c}\hfill I_{5,j}=\frac{N_{\overline{r}}}{5}\end{array}$$(1)

where $N_{\overline{r}}$ is the number of nearest neighbours to data point j within a radius of the mean Euclidean distance, $\overline{r}$, of every data point to its 5th nearest neighbour in the control field. The index is a unit-less ratio with a value in the range $0\le I_{5,j}\le \frac{n-1}{5}$ such that the higher the value, the more spatially clustered a data point. For each dataset the index is calibrated so that significant values can be identified. To do this, 100 realisations of a random distribution of the same size n, and in the same parameter space, as the dataset are generated. Then INDICATE is applied to the random samples to identify the mean index values of randomly distributed data points, ${\overline{I_5}}^{\mathrm{random}}$. Point j is considered spatially clustered if it has an index value above a “significance threshold”, I_sig, of three standard deviations, σ, greater than ${\overline{I_5}}^{\mathrm{random}}$ i.e.

$$\begin{array}{c}\hfill I_{5,j}>I_{\mathrm{sig}},\text{where}{I}_{\mathrm{sig}}={\overline{I_5}}^{\mathrm{random}}+3\sigma \end{array}$$(2)

Table 2 lists the significance thresholds for S₁ to S₄. Appendix A of Buckner et al. (2019) gives an in-depth discussion of the behaviour and properties of the index in random distributions.

3.2. Statistical considerations of field star contamination

As discussed in Sect. 2.2, 29% (345) of the 1189 Class III/F sources in our sample are expected to be field stars. To ensure our spatial analysis results are reflective of the behaviour of Class III YSOs we randomly remove 345 sources flagged as “III/F” from the S₁ and S₄ samples prior to analysis. We limited sources that we removed to those which have not also been identified by K14 as YSO (0/I, II/III, Ambiguous). After analysis the sources are replaced, and the process repeated for a total of 100 iterations. Statistics presented in Sect. 4 for S₁ and S₄ are representative of mean values derived over the 100 samples. The significance thresholds given for S₁ and S₄ in Table 2 were determined for sample sizes of 1450 and 844 respectively. The maximum difference of mean index values for each iteration is $\overline{I_5}<0.1$ for both samples. We find that changing the contamination rate to the higher estimate of 34% (Sect. 2.2) has a negligible impact; the trends found in Sect. 4.1 and 4.2 are unchanged and the difference between mean index values for the two rates is $\overline{I_5}<0.1$.

4. Results

4.1. Distribution of the YSO population

We applied INDICATE to the S₁ sample to investigate the clustering behaviour of YSOs in NGC 2264. As expected the majority of clustered stars are located in the southern region within the star formation active NGC 2264-C and D sub-clusters, whereas clustering in the older northern region is primarily found in the vicinity of S Mon. There is a distinct relationship between evolutionary stage and clustering behaviour of the YSOs. The number of Class 0/I objects with an index above the significance threshold in S₁ is 69.4%, in contrast to 27.9% for Class II, 11.5% for Class TD, and 7.7% of Class III members. Furthermore, there is also a relationship between the degree to which YSOs are clustered (number of neighbours in local neighbourhood) and class; spatially clustered YSOs have median I₅ values of 5.2 for Class 0/I, 4.2 for Class II, 3.2 for Class TD, and 2.8 for Class III. This implies that first, the more evolved an object is the less likely it is to be clustered and second, more evolved objects that are clustered are less concentrated and more dispersed than their less evolved counterparts.

To measure whether these trends are real and significant we compare the index values derived for the different classes using two sample Kolmogorov–Smirnov Tests (2sKSTs) with a strict significance boundary of p <  0.01. The null hypothesis of this test is that differences in the comparative clustering behaviours of two classes are not significant, so their index values have similar empirical cumulative distribution functions (ECDFs). The similarity is quantified by the 2sKST statistic, D, as the distance between two ECDFs (the smaller the statistic, the more similar the distributions). Figure 3 shows the ECDFs of Class 0/I, II, TD, and III objects in S₁ to be dissimilar, and this is confirmed by the 2sKSTs (p ≪ 0.01). We therefore reject the null hypothesis: Class 0/I, II, TD, and III objects have distinct clustering behaviours and our finding that clustering behaviour is a true function of evolutionary stage is both real and significant.

Fig. 3. Empirical cumulative distribution function (ECDF) of index values, I5, calculated for sample S1. The dashed black line denotes the significance threshold of the sample (Table 2). The intercept between the significance threshold and ECDFs is equal to 1−F, where F is the fraction of sources YSOs with an index value greater than this (I5 >  Isig) for each class. As can be seen, the ECDFs of each class are distinct which indicates the differences in their clustering behaviours are significant.

Our assertion is further strengthened by the ECDFs of Class TD and III objects, which are distinct but closely resemble each other (D_TD,III = 0.1). As Class TD objects represent an intermediate evolutionary stage from Class II to III it is reasonable to expect these objects to demonstrate the most similar clustering traits to the Class III objects (the next evolutionary stage).

4.2. Spatial behaviour within classes

We now apply INDICATE to the S₂, S₃, and S₄ samples to evaluate the tendency for objects of the same class to cluster together. Table 3 summarises the statistics of the index values derived for each sample. There is a distinct trend between class and proportion of objects with an index above the significance threshold: 84.7% (Class 0/I), 35.2% (Class II), 2.8% (Class III). In addition, Class 0/I objects are also found to be typically more tightly clustered together than Class II and III objects by a factor of 4.1 and 6.6 respectively. The maximum number of nearest neighbours of the same class decreases with increasing evolutionary stage from 61 (Class 0/I) to just 20 (Class III). This implies that first, the less evolved an object is the more likely it is to be clustered with objects of the same class and second, less evolved objects that are spatially clustered are typically more tightly concentrated together and less dispersed than their more evolved counterparts.

An exception to this trend is in the vicinity of the northern O-type binary, S Mon. In this region we find Class III objects in the local neighbourhood of S Mon to be significantly more self-clustered than the wider NGC 2264 region and exhibit spatial behaviour patterns comparable to that of Class II objects. Sources have higher values than typical within a radius of 0.1° of S Mon, but in particular within 0.05° which contains the most spatially clustered Class III objects in the whole sample. Here 29.1% of objects have an index above the significance threshold, with the median index value of those objects being $\stackrel{\sim }{I_5}=2.0$, which is comparable with Class II’s across the NGC 2264 region (Table 3).

4.3. Kinematic behaviour within classes

We examine the magnitudes of proper motion for our distance selected sample (Sect. 2.4), finding it to have median value of 1.131 mas yr⁻¹ with 1.009 mas yr⁻¹ and 1.192 mas yr⁻¹ for Class II and III sources respectively.

The kinematic distribution of Class II and III sources shown in Fig. 4 are consistent with our findings in Sect. 4.1. Motions of Class III sources are dispersed and randomised with no clear expansion or contraction patterns in the southern region. In the northern region they appear to have a collective outward motion, and there is a grouping in the local neighbourhood of S Mon. While most Class II sources have an outward motion in the northern region, the position and kinematic behaviour of the majority of Class II sources in the NGC 2264-C/D region is consistent with the properties of the “J”,“K”, and “M” sub-clusters identified by Kuhn et al. (2014) using finite mixture models and kinematically characterised by Kuhn et al. (2019) with DR2 data (see Fig. 14 and Table 4 therein).

Fig. 4. Distribution of Gaia proper motions for (top panel) Class II and (bottom panel) Class III objects of our sample overlaid on the Herschel 70 μm map of the region to clearly show the locations of S Mon and sub-clusters C, D for the reader’s reference (outliers not shown). The black solid and dot-dash ellipses around S Mon represent radii of 0.05° and 0.1° respectively.

5. Discussion

We summarise the results of our analysis as follows. There is a difference in spatial behaviour as a function of class in NGC 2264. The youngest, most deeply embedded Class 0/I sources are typically found in strong concentrations with both the general population – and other Class 0/I – sources. While the more evolved Class II and TD sources are also found in such concentrations, the intensity of the concentrations and fraction of the population found in these significantly decreases with increasing evolutionary stage. The trend extends to the Class III sources for which the vast majority are randomly distributed and only a few are found in relatively loose concentrations with the general populous and/or sources of a similar class. This is consistent with previous studies of the region which identified, through qualitative analysis, that Class II objects as being more widely distributed than Class I objects (Sung et al. 2009; Teixeira et al. 2012).

The spatial patterns we find are echoed in the kinematic behaviour of Class II and III sources, which differ considerably. Within the star formation active NGC 2264-C/D regions, Class III sources have predominantly random trajectories and no clear groupings. Objects at the edge of the cluster typically have a larger proper motion than their more central counterparts, which is expected from virial balance as they see a larger enclosed mass. In contrast, Class II sources in this region do not have fully randomised trajectories and demonstrate kinematic behaviour consistent with the known substructure in Kuhn et al. (2019). Although both samples are expected to contain some unresolved binaries, the disparity in their kinematic behaviours suggests that the observed spatial behaviour is age-driven dynamical evolution rather than primordial, this agrees with the work of Venuti et al. (2018) who determined Class III objects to be older than Class II objects and to may have undergone post-birth migration.

With age-driven dynamical evolution sources in the northern region should be significantly less clustered than the southern region, because star formation began there (Sung et al. 2009; Sung & Bessell 2010; Venuti et al. 2017, 2018; González & Alfaro 2017), that is sources have had more time to disperse. Indeed, sources in the north are significantly less clustered than those of the south; 6.9% and 29% have an index above the significance threshold in the north and south respectively. Moreover, there is a correlation between the tightness of clusterings and region, with spatailly clustered sources having median I₅ values of 2.6 (north) and 4.6 (south). The outward motion observed in the kinematics for Class II objects in the north suggests a population at a more advanced stage of dispersal than their southern counterparts.

Interestingly, in the vicinity of the northern O-type binary, S Mon, Class III objects display atypical spatial behaviour. Both Sung et al. (2009) and Venuti et al. (2017) reported a lack of objects with discs within 0.1° of the massive star due to disc disruption caused by stellar feedback. While we found Class III objects in the local neighbourhood of S Mon to be significantly more self-clustered within 0.1° of S Mon than the wider NGC 2264 region, the disparity is more prominent within 0.05°. Here, Class III’s exhibit spatial behaviour patterns comparable to that of Class II sources suggesting disc ablation is causing these objects to appear as more evolved sources.

6. Conclusions

We have characterised the dynamic and spatial distributions of YSOs in the young NGC 2264 cluster. This was achieved through analysis of pre-existing membership catalogues with the new local indicator tool INDICATE and kinematic data from the second instalment of the Gaia catalogue.

In agreement with previous studies, we found the spatial behaviour of objects with and without discs to be distinct, indicating that star formation has been occurring sequentially over a prolonged period. The tool INDICATE has allowed us, for the first time, to quantitatively

1. establish spatial criteria for a source to be considered truly “clustered” or “dispersed” in the region;

2. establish that the proportion of clustered sources decreases with increasing evolutionary stage (Class 0/I, II, TD, III);

3. measure the tightness of these clusterings and establish that this decreases with increasing evolutionary stage;

4. establish that the older northern region has a smaller proportion of clustered sources than the younger southern star formation active region;

5. measure the tightness of these clusterings across the two regions and establish that they are tighter in the south than the north of the cluster;

6. establish that Class IIIs within the local neighbourhood of S Mon exhibit spatial clustering behaviours typical of Class II objects in NGC 2264.

Combining our spatial analysis with kinematic data from Gaia DR2 we derive strong evidence that NGC 2264 is dynamically evolving with stars forming in a centralised, tightly clustered environment, in which they remain for their earliest stage of development before forming part of the dispersed population of NGC 2264. The effect of stellar feedback from S Mon on neighbouring stars is significant, causing these objects to appear as more evolved sources through disc ablation within a radius of 0.1° and particularly within 0.05°.

Thanks to the second data release of Gaia an unprecedented volume of high-precision dynamical data became available for a large number of young clusters. With additional releases planned over the next few years our understanding of star formation and the nature of structures/patterns in these regions is set to increase profoundly. An important consideration going forward therefore is how best to extract, analyse, and interpret these data to produce reliable, robust, and consistent results. In particular, it is important that the community gives careful consideration to terms relating to spatial distribution patterns of sources in these regions, such as “clustered” and “dispersed”, especially in the context of identifying comparative differences. Until now such terms have been frequently used in literature as qualitative descriptors, but when applied subjectively to interpret dynamical behaviours they are at best vague, and at worst could lead to over-interpretation of the data. Building a true picture of star formation history in clusters will therefore require dynamical analyses to be validated by rigorous spatial analysis, in which such terms are clearly, consistently, and quantitatively defined. In this work, we demonstrated with NGC 2264 that the local indicator code INDICATE which quantifies the intensity, correlation, and distribution of stars, can perform this analysis. When combined with Gaia DR2 data can be used to robustly analyse the star formation histories of young clusters.

Acknowledgments

The Star Form Mapper project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 687528.

References

Appendix A: Proper motion perspective corrections

We correct for the perspective contraction of NGC 2264 (caused by radial motions of members) for each source i using Eq. (13) of van Leeuwen (2009) as follows:

$$\begin{array}{c}\hfill {\mu }_{\alpha \ast ,i}^{\mathrm{cor}}=\mathrm{\Delta }{\alpha }_i({\mu }_{{\delta }_0}sin{\delta }_0-\frac{V_{\mathrm{r}0}{\omega }_0}{\kappa }cos{\delta }_0)\end{array}$$(A.1)

$$\begin{array}{c}\hfill {\mu }_{\delta ,i}^{\mathrm{cor}}=-\mathrm{\Delta }{\alpha }_i{\mu }_{{\alpha }_0^{\ast }}sin{\delta }_0-\frac{\mathrm{\Delta }{\delta }_i{V}_{\mathrm{r}0}{\omega }_0}{\kappa }\end{array}$$(A.2)

where ${\mu }_{\alpha *,i}^{\text{cor}}$ and ${\mu }_{\delta ,i}^{\text{cor}}$ are the corrected components of proper motion in right ascension and declination respectively. Table A.1 lists the values of the central coordinates of the cluster (α₀, δ₀), distance unit conversion factor (κ), mean proper motion (${\mu }_{{\alpha }_0^{\ast }},{\mu }_{{\delta }_0}$), parallax (ω₀) and radial velocity (V_r0) used.

For each source we subtract the perspective correction and mean proper motion of the sample to gain the corrected internal proper motion:

$$\begin{array}{c}\hfill {\mu }_{\alpha \ast ,i}^{\mathrm{final}}={\mu }_{\alpha \ast ,i}^{\mathrm{DR}2}-{\mu }_{\alpha \ast ,i}^{\mathrm{cor}}-{\mu }_{{\alpha }_0^{\ast }}\end{array}$$(A.3)

$$\begin{array}{c}\hfill {\mu }_{\delta ,i}^{\mathrm{final}}={\mu }_{\delta ,i}^{\mathrm{DR}2}-{\mu }_{\delta ,i}^{\mathrm{cor}}-{\mu }_{{\delta }_0}\end{array}$$(A.4)

where (${\mu }_{\alpha *,i}^{\text{final}}$, ${\mu }_{\delta ,i}^{\text{final}}$) are the corrected, and (${\mu }_{\alpha *,i}^{\text{DR2}},{\mu }_{\delta ,i}^{\text{DR2}}$) the Gaia DR2, components of proper motion for source i in right ascension and declination respectively.

Appendix B: Inconsistent catalogue classifications

All Tables

All Figures

	Fig. 1. Distribution of our sample overlaid on the Herschel 70 μm map of the region to clearly show the locations of S Mon and sub-clusters C, D. Colours/shapes: Class 0/I (red triangles), Class II (green plus signs), Class TD (purple stars), Class III/F (dark blue dots), Class II/III (grey crosses), Class Ambiguous (pink squares), and non-YSO (light blue pentagons).
In the text

	Fig. 2. Left panel: sample overlaid on a Herschel 250 μm image of NGC 2264. White lines and crosses denote the borders and sources of the sample respectively. Right panel: zoomed-in plot of Fig. 1 in the NGC 2264-C and -D sub-cluster region, overlaid on a corresponding Spitzer 24 μm image. Colours and symbols as defined in Fig. 1.
In the text

Fig. 3. Empirical cumulative distribution function (ECDF) of index values, I5, calculated for sample S1. The dashed black line denotes the significance threshold of the sample (Table 2). The intercept between the significance threshold and ECDFs is equal to 1−F, where F is the fraction of sources YSOs with an index value greater than this (I5 >  Isig) for each class. As can be seen, the ECDFs of each class are distinct which indicates the differences in their clustering behaviours are significant.

In the text

	Fig. 4. Distribution of Gaia proper motions for (top panel) Class II and (bottom panel) Class III objects of our sample overlaid on the Herschel 70 μm map of the region to clearly show the locations of S Mon and sub-clusters C, D for the reader’s reference (outliers not shown). The black solid and dot-dash ellipses around S Mon represent radii of 0.05° and 0.1° respectively.
In the text