NGTS clusters survey $-$ V: Rotation in the Orion Star-forming Complex

We present a study of rotation across 30 square degrees of the Orion Star-forming Complex, following a $\sim$200 d photometric monitoring campaign by the Next Generation Transit Survey (NGTS). From 5749 light curves of Orion members, we report periodic signatures for 2268 objects and analyse rotation period distributions as a function of colour for 1789 stars with spectral types F0$-$M5. We select candidate members of Orion using $\textit{Gaia}$ data and assign our targets to kinematic sub-groups. We correct for interstellar extinction on a star-by-star basis and determine stellar and cluster ages using magnetic and non-magnetic stellar evolutionary models. Rotation periods generally lie in the range 1$-$10 d, with only 1.5 per cent of classical T Tauri stars or Class I/II young stellar objects rotating with periods shorter than 1.8 d, compared with 14 per cent of weak-line T Tauri stars or Class III objects. In period$-$colour space, the rotation period distribution moves towards shorter periods among low-mass (>M2) stars of age 3$-$6 Myr, compared with those at 1$-$3 Myr, with no periods longer than 10 d for stars later than M3.5. This could reflect a mass-dependence for the dispersal of circumstellar discs. Finally, we suggest that the turnover (from increasing to decreasing periods) in the period$-$colour distributions may occur at lower mass for the older-aged population: $\sim$K5 spectral type at 1$-$3 Myr shifting to $\sim$M1 at 3$-$6 Myr.


INTRODUCTION
The evolution of a star depends primarily upon its initial mass, metallicity and angular momentum. These properties help to determine the nuclear reaction rates and processes, energy transport modes, pressure support mechanisms and mass loss rates which lead to distinct evolutionary pathways. While high-mass stars without convection zones do not easily spin down, and so spend their lives as rapid rotators, low-mass stars ( 1.3 M ) evolve appreciably, affected along the way by both internal and external factors. During formation, stellar phases in the evolution of angular momentum in low-mass stars: the first few million years of the PMS phase consists of a largely constant surface rotation rate, followed by an abrupt increase towards the zeroage main sequence (ZAMS), and a steady decline on the MS (Gallet & Bouvier 2015). This picture relies on observations of coeval populations of stars in open clusters and interpreting them in a coherent theoretical framework, such that we can attempt to understand the sequence of evolutionary steps and the relative ages at which they take place. Well-populated open clusters are ideal sources, because they contain stars which span a large range of masses, but which are essentially of the same age and composition. For field stars whose properties are otherwise stable (rendering dating by other means such as via isochrones difficult), the empirical MS spin-down can provide a valuable age predictor. It is the basis of gyrochronology (Barnes 2003), potentially opening the door to stellar ages for large datasets (e.g. McQuillan et al. 2013McQuillan et al. , 2014Angus et al. 2015;Davenport 2017;Davenport & Covey 2018;Lu et al. 2021), and a better characterisation of exoplanets and their host stars (e.g. Gallet & Delorme 2019;Zhou et al. 2021). Models (e.g. Gallet & Bouvier 2013Amard et al. 2016Amard et al. , 2019 have demonstrated that the main PMS and MS evolutionary trends, as seen in observations of rotation in lowmass stars, can be reproduced, although the picture is not complete (e.g. Roquette et al. 2021;Godoy-Rivera et al. 2021).
The census of observed rotation periods for low-mass stars is quite well-populated at many stages of evolution, from star-forming regions to mature clusters, e.g. at 1 Myr in the Orion Nebula Cluster (ONC) (Herbst et al. 2002;Rodríguez-Ledesma et al. 2009) and Oph (Rebull et al. 2018); 3 Myr in Taurus (Rebull et al. 2020) and NGC 2264 (Lamm et al. 2004;Venuti et al. 2017); 8 Myr in Upper Sco (Rebull et al. 2018); 13 Myr in h Per (Moraux et al. 2013); 35 Myr in NGC 2547 (Irwin et al. 2008); 110-120 Myr in the Pleiades (Rebull et al. 2016) and Blanco 1 (Gillen et al. 2020); 150 Myr in NGC 2516 (Bouma et al. 2021); 700-800 Myr in Praesepe and the Hyades (Brandt & Huang 2015a,b); 1 Gyr in NGC 6811 (Meibom et al. 2011); 3 Gyr in Ruprecht 147 (Gruner & Barnes 2020); and 4 Gyr in M67 (Esselstein et al. 2018). However, due to variable mass ranges, time coverage, data quality and field contamination between surveys, new contributions and updates remain valuable. This is particularly so in the era of Gaia (Gaia Collaboration et al. 2016Collaboration et al. , 2021, with its unmatched astrometry and uniform broadband photometry across the sky. Rotation periods for low-mass stars in the youngest-observed clusters (1-3 Myr) have often revealed a broad distribution across the observed mass range, with typical periods of 1-10 d (e.g. Bouvier et al. 2014;Rebull et al. 2018Rebull et al. , 2020. The cause of the initial dispersion of rotation rates at ages younger than 1 Myr is not well understood, but it may be linked to star-disc interactions in the embedded protostellar phase, with more massive discs being more efficient at preventing protostars from spinning up (Gallet & Bouvier 2013).
For young stars, there is significant observational evidence suggesting that, at a given age, stars with a circumstellar disc are, on average, slower rotators than those without a disc (e.g. Herbst et al. 2002;Rodríguez-Ledesma et al. 2009;Affer et al. 2013;Serna et al. 2021). Angular momentum evolution models commonly invoke a disc-locking mechanism, caused by magnetic interactions between the young star and its accretion disc, which maintains a constant stellar angular velocity for the lifetime of the inner disc. There are, however, a number of studies which find no significant differences in rotation properties between systems with and without an accretion disc (e.g. Stassun et al. 1999;Nguyen et al. 2009;Le Blanc et al. 2011;Karim et al. 2016). During the early PMS (< 10 Myr), the lowest-mass stars appear to preferentially spin up, leaving a dearth of slow rotators (Bouvier et al. 2014;Roquette et al. 2021). Such observations are often correlated with the presence or absence of excess near-mid infrared emission, indicative of an accretion disc (Herbst et al. 2002;Rodríguez-Ledesma et al. 2009;Affer et al. 2013;Venuti et al. 2017;Rebull et al. 2018). Hence, they can be explained by different disc-locking timescales, on the understanding that slow rotators are prevented from spinning up due to ongoing star-disc interactions, whilst the discs of fast rotators have been dissipated, allowing the young stars to spin up as they contract towards the ZAMS, which they reach at ages of approximately 22, 33, 66, and 100 Myr for masses of 1.2, 1.0, 0.7, and 0.5 M , respectively (Moraux et al. 2013). The mass-dependent disc-locking timescales may be influenced by the local environment, i.e. via external photoevaporation driven by the far-ultraviolet emission of massive stars, which disperses discs around very low-mass stars more quickly than those around higher mass stars (Roquette et al. 2021). At higher masses ( 0.3 M ), bimodal distributions, which may also be a consequence of variations in disc longevity, have been reported, e.g. at ∼1 Myr in the ONC (Herbst et al. 2002;Rodríguez-Ledesma et al. 2009) and ∼3 Myr in NGC 2264 (Venuti et al. 2017), as well as at post-accretion ages (> 10 Myr), e.g. at ∼13 Myr in h Per (Moraux et al. 2013)there interpreted as evidence for core-envelope decoupling -and at ∼120 Myr in the Pleiades (Rebull et al. 2016).
The Orion Star-forming Complex is one of the largest and most active regions of nearby star formation, comprising numerous wellstudied clusters with ages up to ∼10 Myr. It is located at an average distance of ∼400 pc, toward the Galactic anticentre (Kounkel et al. 2018, hereafter K18). On the sky, it spans approximately 75 to 90 • in right ascension and −10 to 13 • in declination. The current epoch of star formation is confined to the Orion A and B molecular clouds, home to familiar clusters such as the Orion Nebula Cluster (ONC), NGC 2024, and NGC 2068, where typical ages are 1-3 Myr (K18). Gaia astrometry has enabled more extensive and accurate membership lists to be compiled, coupled with a wealth of other photometric and spectroscopic data that have made study of the Orion Complex more tractable.
Godoy-Rivera et al. (2021) conducted a systematic revision of open cluster sequences based on Gaia DR2 and noted that rotation sequences measured from the ground can be as informative for stellar rotation studies as those from space. Furthermore, while observing from space brings many benefits, cutting-edge ground-based facilities can comfortably measure rotation periods from time-series photometry, and can often do so for longer and in more crowded environments than some space-based instruments. This paper presents a study of rotation in Orion using ∼200 d of ground-based data from the Next Generation Transit Survey (NGTS; Chazelas et al. 2012;Wheatley et al. 2018). The observations were taken as part of the NGTS Clusters Survey (Gillen et al. 2020;Jackman et al. 2020;Smith et al. 2021;Moulton et al. 2023).

Selection of candidate cluster members
K18 performed a kinematic analysis of the Orion Complex using spectroscopic and astrometric data from APOGEE-2 and Gaia DR2. They applied a hierarchical clustering algorithm in five (sometimes six) dimensions to identify distinct groups of young stellar objects (YSOs). Here, we take astrometric and photometric data from the Gaia EDR3 1 data release and create a new candidate membership list, using the K18 members 2 to set bounds on the astrometric parameters of potential members in each field. We performed an EDR3-DR2 crossmatch on the K18 members using the gaiaedr3.dr2_neighbourhood query tool 1 When this work was done, the full DR3 release was not yet available, so it was the EDR3 release from which the Gaia measurements were extracted. We note, however, that the astrometry and photometry used are identical between DR3 and EDR3. 2 K18 members being those stars assigned to a named group in their paper. (https://gaia.aip.de/metadata/gaiaedr3/dr2_neighbourhood/), taking the smallest angular separation match if multiple EDR3 sources matched a DR2 ID. K18 objects were discarded if parallax ( ), proper motion ( , ), or photometric data was absent, or if the astrometric precision was such that / > 0.1, > 0.2 mas yr −1 , or > 0.2 mas yr −1 was satisfied. After clipping single outliers in three of the four fields, bounds for the new candidates were set to the K18 members' minimum and maximum values of , and for each field, with the exception of field NG0535, which contained a large number of outliers in parallax; here, the bounds were set to the mean ±4 standard deviations. All EDR3 sources were then extracted, subject to the K18-based bounds on parallax and proper motion, and / < 0.1, < 0.2 mas yr −1 , and < 0.2 mas yr −1 as requirements on precision. The EDR3 parallaxes were corrected for the zero-point bias using the expression given in Lindegren et al. (2021), and the -band magnitudes for sources with six-parameter astrometric solutions were corrected in accordance with Riello et al. (2021). The distributions of the resulting EDR3 candidates and the K18 members are shown in Figure 2, and colour-magnitude diagrams (CMDs) are shown in Figure 3.
The approach of using astrometric cuts based upon the K18 members' properties leads to naturally similar, but not identical, distributions. If the goal had been replication of the K18 distributions, their clustering algorithm would need to be applied, along with radial velocity data from APOGEE. A significant difference also lies in the use of Gaia EDR3 vs DR2, where the former is an expanded catalogue, with a level of precision that brings more objects within the boundaries set; this is particularly evident when comparing the -magnitude histograms of Figure 2. The priority in this work was, rather, to try and capture as many members as possible, maximising the yield of YSO rotation periods. The elevated false positive fraction ought to be mitigated somewhat in the periodic sample, following inspection of light curves and the identification of those objects with rotational modulation patterns characteristic of young stars.  . CMDs of the candidate Orion members for the four NGTS fields studied. Gaia magnitudes have been converted to absolute magnitudes via their parallaxes. No extinction correction has been applied at this stage.

PERIOD DETECTION PIPELINE
In this section, we explain the processing of light curves, identify spurious signals, distinguish rotation signatures from other variability, and make a comparison to literature measurements.

Light curve pre-processing
Data points flagged by the NGTS pipeline (e.g. from pixel saturation, blooming spikes, cosmic rays, laser crossing events) and 7-sigma outliers were masked (light curves >80 per cent masked were removed). The light curves were binned in time to 20 min, and those with a single gap greater than half the baseline of the observations, or, with three or more gaps greater than 30 d, were dropped.

Periodic signals
Period measurements were made by calculating the Lomb-Scargle periodogram (Lomb 1976;Scargle 1982) for each candidate member using the (Astropy Collaboration et al. 2013) package, over a search grid covering frequencies 0.001 to 24 d −1 . Periods corresponding to the highest peaks in the periodograms were taken as provisional rotation periods. 3 While statistical uncertainties on Lomb-Scargle periodogram measurements do not capture the real uncertainties inherent to the technique, e.g. inaccuracies associated with false peaks and aliases, the effects of long-term trends and spot evolution, we measure the half 3 False-alarm probabilities associated with the highest peaks were all effectively zero, with an extreme outlier maximum value of 10 −6 and a median of 10 −190 .  = 1789; see section 5.1). Right: percentage relative error as a function of period. The grey dashed line is a precision boundary used during filtering (section 3.2). We note that the slightly separated population below the main grouping consists of targets observed in field NG0523.
width at half maximum (HWHM) of the periodogram peaks in frequency space to estimate the precision of our measurements (on the assumption that the correct peak has been selected). Figure 4 summarises these uncertainties 4 by displaying the relative error as a histogram and as a function of period.
NGTS light curves have been found to sometimes retain the imprint of flux from the moon for fainter targets. A typical moon-affected light curve exhibits periodic dips in flux in phase with the lunar cycle, due to an over-correction of the sky background. In order to detect potential rotation signals present in these light curves, a simple trend removal step was incorporated for objects with >14 and an initial period between 27 and 30 d. A Savitzky-Golay (SG) filter was applied to the light curves (phase-folded on the detected moon period), followed by a convolution, with the target light curve being detrended by the result of these two steps (SG filter + convolution). An equivalent detrending was applied to light curves with initial periods greater than half the baseline of the observations.
The NGTS pipeline includes a calculation of the dilution affecting each target, with stars within 7 pixels of the target and brighter than magnitude 16 in the TESS band contributing. In this work, targets were dropped from the period analysis if the summed flux of the contaminating stars exceeded the target flux. Additionally, objects separated by less than 20 arcsec (a distance at which the flux contribution of an average source falls close to zero) were dropped if the percentage difference in their periods was below the precision boundary line plotted in the right-hand panel of Figure 4 (at the relevant period), and the measured amplitudes of the signal could not identify the source 5 . Objects with otherwise suspect periods were provisionally rejected, subject to inspection. The periods in question were (1) those likely to be caused by the diurnal pattern of observations, i.e. the one-day signal and its aliases, (2) those ( > 14) which could be an alias of the lunar period (or half the lunar period), and (3) those ( > 14) of ∼half the lunar period (13.5-15 d).
A detected period was classified as an alias if it fell within calculated boundaries of the expected alias periods described by which comfortably enclosed the corresponding peaks in histograms of the detected periods. The remaining periods were provisionally accepted, subject to inspection.

Injection-recovery tests
In an attempt to evaluate the ability of the period detection pipeline to recover real rotation signals, injection-recovery tests were conducted on each field. For a given stellar magnitude and period, the results take the form of a distribution of threshold amplitudes, above which signal recovery was successful. A percentile score was assigned to each detection in the main sample, based on its amplitude among the test distributions (at the corresponding magnitude and period). We refer the reader to appendix A and Figure A1 for the full details.

Rotation periods
Pipeline output and light curves (in time and in phase) were inspected, filtered and labelled based on their likelihood of representing stellar rotation. Objects not in one of the three spurious categories previously described were provisionally accepted if either their percentile score from the injection-recovery tests was above 80, or the detection had been labelled as 'clean' -an attempt to identify objects where the detected signal is unique and unambiguous. Following Xiao et al. (2012) and Covey et al. (2016), the designation is given to objects whose periodogram contains no secondary peaks exceeding 60 per cent of the height of the primary peak, aside from beat periods between the primary peak and the window function (the one-day sampling period). If, upon inspection of the data, a provisionallyaccepted signal appeared suspect, that object was then rejected. Objects initially classified as either the one-day signal, an alias of the moon, or half the lunar period, were accepted following inspection in 0.4, 3 and 12 per cent of cases, respectively. This initial stage of inspection left a sample of detections believed to be of astrophysical -but not necessarily rotational -origin. In order to identify the signals which most-likely represent stellar rotation periods, each remaining detection was given a period quality label of 1, 2, 3 or 4. '1' indicates a signal believed to be a clear stellar rotation period (although aliases cannot be ruled out); '2' indicates a signal which is also believed to be stellar rotation, but where the detected signal is relatively weak; '3' indicates a signal which could possibly be stellar rotation, but which could easily be attributed to other forms of variability; and, '4' indicates a signal which, whilst likely to be real, is almost definitely not rotational. Objects in category 3 or 4 tend to have light curves without the typical smooth, starspotinduced modulation patterns for which stellar rotation is a good explanation. They are more stochastic, sometimes displaying signs of accretion bursts or the presence of additional material in the system -57 have literature designations of Type I/II YSO or classical T Tauri star (CTTS), compared with 6 of Type III YSO or weak-line T Tauri star (WTTS). They may also display variations on multiple timescales, which makes the identification of a rotation signal hard to pinpoint (even if present), particularly for a rigid, single-component model like that used in Lomb-Scargle. The forthcoming analysis of rotation period distributions is restricted to category 1 and 2 objects.
Label 1 was assigned to all objects whose percentile score exceeded 95 and which held the 'clean' designation. A second round of inspection identified class 3 and 4 objects, with the remainder being assigned to class 2. Out of 5749 stars with NGTS light curves and 4964 stars for which period measurements were attempted, 2268 periods were retained, with 2179 of those being assigned to class 1 or 2. All period measurements are available in Table 2, along with supplementary data on the targets. Figure 5 shows some example light curves, periodograms and phase folds.

NG0523/NG0533 duplicates
There is a region of overlap between fields NG0523 and NG0533 (see Figure 1), resulting in 165 objects with two light curves, 84 of which had period detections in at least one field surviving the abovedescribed filtering. Duplicate objects were removed as follows, based on period detections where available. The most convincing detection was manually selected in seven cases where period estimates disagreed (generally, the differences are attributable to beat periods or harmonics). If an object had a valid period measurement from just one field (10 objects in NG0523 and 11 objects in NG0533), the corresponding data was retained. The data from field NG0533 was preferred in all other cases, due to the large gap in observations for field NG0523. Excluding the seven objects where different Lomb-Scargle peaks were preferred between fields, approximately 100 (90) per cent of objects with detections in both fields agree to within 3(1) per cent. Table 3 gives an indication of the completeness of the NGTS sample and the periodic sample as a fraction of the full Orion candidate members list, described in section 2.1. Approximately 80 per cent of the candidate members have NGTS light curves, and we obtain period measurements for ∼32 per cent (the accepted periods described above). In the most crowded region -around the centre of the Trapezium cluster (RA = 83.82, Dec = −5.39) -dilution restricts the number of successfully retrieved periods; although NGTS has light curves for 87 per cent of candidate members in this inner region of the ONC, we retrieve periods for only 15 per cent. As a comparison, observations for the classic Herbst et al. (2002) study of rotation in this same region, were made using 0.24-arcsec pixels, compared with the 5-arcsec pixels of NGTS, which results in a small overlap.

Completeness
From the candidate members with NGTS light curves, we recover periods for ∼40 per cent, whilst from the K18 members with light curves (which constitute approximately half of the candidate members with light curves) we recover ∼56 per cent 6 . This difference is partly attributable to the fainter stars incorporated in this work, but may also reflect a smaller false positive fraction of members in the K18 sample.

Comparison with literature rotation periods
We compare our rotation periods to literature values from Stassun et al. (1999) A selection of NGTS light curves and output plots from the periodic detection pipeline. Each object appears in a row, with its light curve (binned to 20 min), followed by a Lomb-Scargle periodogram and light curve phase-folded on the selected period. In each periodogram, a red vertical line locates the periodogram peak of the adopted period, while blue dashed lines locate some of the beat periods resulting from the 1-d sampling. The phase-fold plots cycle through a colourmap with observation time: beginning (blue) to end (yellow). The period and period quality designation (on a scale of 1-4; section 3.4) are shown above. The second and third objects from the bottom are examples where the primary periodogram peak was not selected as the most-likely period. The third-from-bottom star exhibits significant structure in its phase-folded light curve, most likely due to dust enshrouding the system (Stauffer et al. 2017;Zhan et al. 2019;Günther et al. 2022), while the system below is an eclipsing binary, where the most-likely period was selected based on the out-of-eclipse variability and differing eclipse depths. The final example is an object where the variability was thought less likely to reflect rotation, with a quality 3 designation given.  of those objects with periods differing by more than 5 per cent are explainable as being either beat periods related to the 1-d sampling of the observations, or as 2:1 or 1:2 harmonics of the periods identified in this work.

STELLAR AND CLUSTER PARAMETERS
In what follows, we begin by explaining our procedure for estimating interstellar extinction and for obtaining effective temperatures, before assigning stars to kinematic groups and deriving individual and cluster ages.

Extinction from broadband photometry
We estimated extinction on a star-by-star basis by comparing the observed BP − RP , − RP and − colours affected by reddening, with a table of standard colours (SC table hereafter), e.g.
where BP and RP are the extinctions in the Gaia BP and RP photometric bands in this case. The standard colours came from Luhman (2022). We did not use colours involving WISE or -band photometry, so as to mitigate the worst effects of infrared excess from circumstellar discs. Additionally, we dropped the − colour from the fit on occasions when the 2MASS source was matched to multiple Gaia objects (a consequence of Gaia's higher angular resolution), as determined by the 'number_of_mates' parameter in the gaiadr3.tmass_psc_xsc_best_neighbour table from the Gaia documentation.
The SC table incorporates the colours and spectral types from Table 4 of Luhman (2022), and the eff values corresponding to spectral types F0-M4 from   eff values for spectral types later than M4 were taken from Table  5 of Herczeg & Hillenbrand (2014), following Fang et al. (2017) and Fang et al. (2021). Linear interpolation in eff -colour space was applied to obtain a particular intrinsic colour prediction for a given effective temperature. The extinction values were obtained via the reddening law from Fitzpatrick et al. (2019) and synthetically reddened PHOENIX spectra (Husser et al. 2013 spectroscopically-derived TIC-8 temperatures, a correction only applied to field NG0535, which was the only field showing significant discrepancies. Our adopted uncertainties generally increase with stellar mass and (in the case of spectroscopically-derived temperatures) are typically ∼250 K for M spectral types, increasing to ∼500 K for spectral types K and G, before increasing steeply for spectral types earlier than mid-F. We refer the reader to appendix B for details concerning the sourcing of effective temperatures, the corrections, and the derivation of uncertainties.

MCMC
With a small number of exceptions, to be explained in section 4.1.3, the MCMC runs were initialised with the derived eff values and Gaussian priors described in appendix B, and a uniform prior on the extinction parameter, , in the range 0-15 (initial values from the 3D dust maps of Green et al. (2019)). 100 'walkers' explored the posterior parameter space for 5000 steps. The first 3000 steps were discarded as 'burn-in', and the values corresponding to the 50th, 84th−50th and 50th−16th percentiles from the marginalized distributions over the remaining 2000 steps constitute the final adopted values and 1-sigma errors for eff and . For each step in the MCMC, the current value of was used to redden a PHOENIX spectrum best-matched with the current value of eff and a fixed value of log 8 . Reddening of the PHOENIX spectra was applied using the P Dust Extinction package. The filter response functions were obtained for each photometric band from the Filter Profile Service, and the effective wavelength for each filter (λ) was calculated using the filter transmission ( λ ) and the stellar flux ( λ ) in the respective bandpass. The reddening law (Fitzpatrick et al. 2019) was finally interpolated to the effective wavelengths for each bandpass to give values for λ / . This process was pre-computed for all bandpasses, for all available PHOENIX spectra, for a range of in increments of 0.2. with linear interpolation of all parameters to produce a finer grid.

Extinction constraints
For the objects without any temperature constraint (2 per cent of the sample with rotation periods, and only 0.3 per cent when considering the best-quality sample used for much of the forthcoming analysis), a Gaussian prior was placed on . For objects in fields NG0523, NG0533 and NG0531, we used the 16th, 50th and 84th percentiles of the reddening predictions from the 3D dust maps of Green et al. (2019), converting to assuming = 3.1 and using the coefficient from Table 6 of Schlafly & Finkbeiner (2011). For stars in the ONC-centred field, NG0535, dust map predictions appear to substantially over-predict the extinction, when compared with spectroscopic estimates from K18 9 . In these cases, the approach adopted was to take the values computed from the targets with eff constraints 8 Log , not being well constrained by broadband photometry, was fixed.
It was taken from APOGEE Net where available (41 per cent of periodic sample), but in other cases the MIST v1.2 (Dotter 2016;Choi et al. 2016) stellar models were used to predict its value by interpolation to the star's absolute -band magnitude, derived from 2MASS photometry and Gaia parallax at the median age of the stars in the field (based on K18 HR ages). 9 We note that the use of similar dust maps (Green et al. 2018) in the dereddening procedure applied in the derivation of photometry-based effective temperatures in the TIC-8 catalogue, could explain the systematically high values when compared with spectroscopically-derived temperatures for field NG0535 (see Figure B2). from literature spectral types, APOGEE Net, or the spectroscopic TIC-8 sample, and to interpolate in RA-Dec space to the position of each star requiring a constraint on , i.e. to targets with no available temperature. In addition, due to very large eff uncertainties for the objects with non-spectroscopic TIC-8 temperatures (see Figure  B3), a constraint on was added to these stars too, in an attempt to break the eff -degeneracy. The Gaussian prior was centred on the interpolated value, with a width equal to the standard deviation of the values of the 10 nearest stars to the target. The interpolation was implemented using the S G routine; Figure 7 illustrates the process. Linear interpolation was used, except for the seven stars outside of the interpolation limits; in those cases, the value of the nearest neighbour was adopted. To supplement the stars available for interpolation, we included stars not in the membership list, but which still lie within the parallax bounds and meet the requirements on parallax precision. In order to filter out stars with poor quality photometry from the interpolation, we excluded objects by way of a cut on the corrected Gaia BP and RP flux excess factor, the cut being made at the 5 level relative to the Stetson and Ivezic standards sample (see equation 18 and section 9.4 of Riello et al. 2021).

Binary identification
We identify binary and higher-order systems by two methods. Firstly, using the Renormalised Unit Weight Error (RUWE) goodness-of-fit statistic reported in Gaia EDR3, which is expected to be around 1.0 for sources where the single-star model provides a good fit to the astrometric observations. We consider objects with a RUWE > 1.4 to be likely binary or higher-order multiple star systems (see e.g. Stassun & Torres 2021). Secondly, we draw on the spectroscopic analysis of Kounkel & Covey (2019) -a study of high-resolution APOGEE spectra of nearby star-forming regions, searching for double-lined spectroscopic binaries (SB2s). From their catalogue, we take binary candidates to be objects where multiple components were identified in the cross-correlation functions (CCFs), and also objects labelled as 'inconclusive SB2/Spotted star pair'. The latter group contains stars where the authors are unsure whether the structure in the CCFs is attributable to multiple stellar components or the impact of star spots, i.e. spots can affect the shape of spectral lines, and hence the CCF profile, as the flux deficit they impart moves across the stellar disc, sometimes resembling a spectroscopically unresolved SB2.

Star clusters in Orion
The Orion Complex is of significant volume and is home to a large number of stellar associations reflecting its star-formation history. These associations, or clusters, represent groups of stars, presumably of very similar age. Hence, we attempted to identify our target stars with their parent cluster. We cross-matched the kinematic groups from K18 with Gaia EDR3, cutting those outside of the parallax bounds previously described, and those with RUWE > 1.4. Target stars were then assigned to the best-matching group, i.e. the group for which was minimised, where , and represent the target value, target error and cluster mean for parameter ∈ {RA, Dec, parallax, pmRA, pmDec}. A further stipulation was that each of the target's astrometric parameters was within the minimum and maximum bounds of the group 10 . Targets in common with the K18 objects were automatically given the K18 designation. The K18 groups are specified as sub-groups of parent clusters, e.g. 'onc-1', ' Ori-3' etc. Collecting these sub-groups into their parent clusters and plotting in 2D and 3D space yields Figure 8. The ONC has been split into inner and outer regions, with the inner region being 2000 × 2000 arcsec, centred on the Trapezium cluster 11 .

Individual stellar ages
We derive model-dependent ages (and masses) by linearly interpolating the MIST v1.2 (Dotter 2016;Choi et al. 2016) stellar evolution models in the HRD (log L vs log eff ) and CMD (M G vs M BP − M RP ). For the HRD, we used the eff posterior distributions from the MCMC output, and calculated the stellar luminosities from the absolute extinction-corrected -band magnitudes, using the Gaia parallax and the bolometric corrections for PMS stars given by Pecaut & Mamajek (2013). For stars of spectral type earlier than F, which are absent from the Pecaut & Mamajek (2013) PMS table, the luminosity was calculated based on the Gaia -band absolute magnitudes, using the DR3 bolometric correction tool (Creevey et al. 2022). Distributions of luminosities and magnitudes were calculated using the posterior distributions, enabling age and mass distributions to be calculated for each target in both the HRD and CMD. We tabulate 50th, 84th−50th and 50th−16th percentiles of the age and mass distributions, as well as the 1.4826 × MAD error estimate 12 . In addition, we calculate HRD and CMD quantities based upon the individual median eff and values from their respective posterior distributions, interpolate to the corresponding single points in the HRD and 10 The minimum and maximum values were replaced by the group mean ± two standard deviations in cases where the latter constituted wider bounds. 11 This matches the location studied by Herbst et al. (2002) in their work on stellar rotation in the ONC. 12 i.e. the median absolute deviation scaled to estimate the standard deviation, assuming normally distributed data.

Figure 8.
Left: RA-Dec distribution of the candidate Orion member stars with NGTS light curves. Each star is coloured according to its assigned parent cluster. Right: 3D distribution formed by the inclusion of distance. We take distance to be the reciprocal of the Gaia parallax, which is reasonable given the median and maximum relative parallax errors of 2 and 9 per cent, respectively. The colour coding is the same as in the left-hand panel.
CMD, and quote these values as our best estimates of age and mass. In the vast majority of cases the values are very similar to taking the median of the age and mass distributions, but when a significant fraction of the points in the HRD or CMD fall outside of the model bounds, the distributions are shifted and are no longer centred on the best estimates of the individual parameters. For that reason, and for reasons of consistency, in what follows we use the age and mass values that are derived from the quoted values of eff and . The MIST models were sampled using isochrones between 0.1 and 100 Myr. All post-MS data and parameter space belonging to stars with eff > 13 000 K ( 3.5 M ; spectral type B7.5) was removed. Stars situated in the region corresponding to a younger age than the minimum model age were assigned the minimum age of 0.1 Myr.

Cluster ages
Cluster ages were taken to be the median age of the corresponding member stars, with upper (84th−50th) and lower (50th−16th) percentile uncertainties. The member stars were first filtered based on: (i) the corrected Gaia BP and RP flux excess factor being below the 5 level (as previously described); (ii) the fraction of each distribution in the HRD and CMD lying within the model bounds being greater than 0.5 13 (which was true and ∼1.0 in 91 per cent of cases); (iii) the target not being an identified or candidate binary (see section 4.2); (iv) the photometry being free from blending issues 14 ; (v) eff < 7280 K. Extreme (7 ) outliers were also removed (unless the number of cluster members making it through the above cuts was less than five). This process was applied to the parent clusters and to the subclusters, for both HRD-and CMD-derived ages. The stars used in the derivation of cluster ages are shown in the HRD and CMD in Figure  9. HRD-based age estimates are generally found to be younger than those derived from the CMD. In this work we find that the CMD cluster ages are, on average, a factor of 1.2 older than those from the HRD. Individual HRD stellar ages (for the same stars), as derived from the MIST models, appear as histograms at the top of Figure 10, and as a function of colour in the left panel of Figure 11. From Figure 11, it is apparent that the older stars are, in general, bluer. It is possible that drawbacks to do with the de-reddening method employed, or the more-rapid evolution of higher-mass stars in the HRD, which could 13 With the exception of stars whose points fell below the minimum model age. 14 A blended source was taken to be a star with a stellar companion within 1.75 arcsec (limit from Riello et al. (2021)), where the companion was less than three magnitudes fainter in any Gaia bandpass. In practice, this cut removed very few additional objects on top of the flux excess filter. make interpolation to model grids more sensitive to any inaccuracies in either observation or theory, could be a factor 15 . However, we note that trends of increasing stellar age with increasing stellar mass in model predictions have been seen many times before, e.g. Hillenbrand (1997); Hillenbrand et al. (2008); Herczeg & Hillenbrand (2015); Feiden (2016). It should also be noted that, as a consequence of the location of the NGTS fields, we sample only a small fraction of the members of some clusters. Table 2 includes our derived stellar and cluster ages. Feiden (2016) find that their evolutionary models, which incorporate magnetic inhibition of convection, are able to ameliorate the age discrepancy between high-and low-mass stars in the HR diagram for Upper Scorpius. Inhibition of convection produces lower effective temperatures, slowing the contraction rate of young stars. Therefore, stars have a larger radius and a higher luminosity at a given age. The effect is more dramatic for cool, low-mass stars, having relatively little influence on high-mass stars. Hence, a 10 Myr isochrone from the magnetic Feiden models, looks like a 5 Myr isochrone from nonmagnetic models for stars with eff < 5000 K. To see if the magnetic Feiden models fix the age discrepancy which we observe with our 15 We tested restricting the stars used in deriving cluster ages to those with log eff < 3.75, but the differences were minimal. data and the MIST models, we calculate HRD ages with said models, and plot the results in the bottom half of Figure 10 and in the right-hand panel of Figure 11 16 . It is clear, from Figure 11, that the trend of increasing age with mass remains. It is also evident that, as expected, the ages derived from the magnetic models are older than those from the non-magnetic MIST models. For our purposes, the absolute ages are less important, as we wish to test, primarily, whether there is any noticeable evolution in the period-colour relation. Hence, the sequence of ages is what matters. We see that three pairs of neighbouring clusters (in Figure 10) switch places, but that the youngest five clusters are identical for both the MIST-and Feidenderived ages. This means that a division at age ≤3 Myr, based upon the MIST models (as is adopted in the subsequent analysis) would be equivalent to a division at age ≤6 Myr using the Feiden models.
We have adopted HRD (rather than CMD) ages in Figures 10-15 and in the forthcoming analysis for two main reasons: (1) for a clearer comparison between magnetic and non-magnetic evolutionary models (the magnetic models not being available in the colour-magnitude plane); and (2) the age at which to divide samples between old and young being more obvious using the HRD ages. We note, however, that the general age order of clusters is preserved with HRD or CMD ages, bar a few clusters shifting by one or two places, and so the main trends we will highlight in period-colour space (more faster rotators at the blue and red ends and slowest rotators shifting to lower mass for the older populations) are present even if the exact make-up of the young-and old-aged samples changes slightly when adopting CMD instead of HRD ages. Both sets of ages can be found in Table 2. With stellar properties, cluster membership and cluster ages determined, we can now investigate the rotation period distribution in Orion.

Period-colour distribution
The rotation period distribution of each cluster as a function of ( BP − RP ) 0 colour is shown in Figure 12, ordered by age, as derived from the HRD (see section 4.5). The filled circles represent stars which met all of the criteria stated in section 4.5 and have period quality designation 1 or 2. The open circles also have period quality designation 1 or 2, but are objects which did not meet all of the criteria.
The top row of Figure 13 shows the rotation period distribution as a function of ( BP − RP ) 0 colour for each of two age groupings: 1-3 Myr on the left and 3-6 Myr on the right. Overlaying the periodcolour plots are lines representing the 10th, 50th and 90th percentiles of the rotation period distributions. The black circular markers again represent stars meeting all of the criteria from section 4.5, while the red triangles locate stars which were identified as candidate binaries (see section 4.2), but which otherwise meet the criteria. The binary candidates were incorporated into the percentile calculations, which are shown for both ages together in the bottom-left of the figure. 16 We found that Feiden isochrones for ages younger than 1 Myr cross over those of older ages in a way that makes interpolation problematic. Hence, the minimum-age isochrone used was 1 Myr. Any objects in a region of parameter space corresponding to younger ages were assigned an age of 1 Myr. In addition, there are regions of parameter space (high mass, young age) to which the Feiden models do not extend, but the MIST models do. We find ∼40 objects in this category, but do not expect this to affect any of our conclusions. Additionally, we show the percentiles with all binary candidates excluded in the neighbouring panel.
In order to assess uncertainties on the percentiles, we generate percentile distributions based on the extinction samples from the MCMC data. We take all samples from the ±1 region of the corresponding distribution of ( BP − RP ) 0 for each star, shuffling the order, giving × samples. We then calculate the 10th, 50th and 90th percentiles of the rotation period distributions for the 1-3 Myr and 3-6 Myr populations times, i.e. each sample is used once. We incorporate period uncertainties by taking each period to be a random draw from the Gaussian distribution derived by fitting the relevant periodogram peak in frequency space 17 . We repeat the process for five separate shuffles of the ( BP − RP ) 0 samples and then plot the full extent of the resulting percentile distributions (bottom-right panel of Figure 13). 18 One of the most prominent differences between the rotation distributions at 1-3 and 3-6 Myr is at the red end, where ( BP − RP ) 0 2.25 (M2 spectral type). Here, we see a distribution shifted towards shorter rotation periods in the older-age population, with no periods longer than 10 d for ( BP − RP ) 0 > 2.65 (M3.5). The median rotation period at 1-3 Myr decreases from 4 d to 2 d in the range 2.25 < ( BP − RP ) 0 < 3, but the equivalent decrease at 3-6 Myr is from 4 d to 0.9 d, with additional faster rotators at redder colours. One possible explanation is that the circumstellar discs of very low-mass stars could be more readily dispersed, facilitating an earlier spin-up as they continue their contraction towards the main sequence (e.g. Roquette et al. 2021). We also see a population of high-mass, fast rotators in the older-aged group, not present in the younger ensemble.
In order to test the shift towards shorter rotation periods for the older stars, we ran permutation tests comparing the 10th, 50th and 90th percentiles of the young and old populations for ( BP − RP ) 0 > 2.25. The tests were run times: once for each of the ( BP − RP ) 0 sample sets described above. 99 per cent of the resulting p values lay below 0.001 and 0.004 for the 10th and 50th percentile tests, respectively, with ∼85 per cent of the p values from the 90th percentile tests below 0.05. Hence, these results favour the alternative hypothesis that the rotation periods are longer for the younger population.
Another feature of Figure 13 is that the turnover from increasing to decreasing periods is located at lower mass for the older-aged ensemble: ( BP − RP ) 0 ≈ 2 (M1 spectral type) at 3-6 Myr, compared with ( BP − RP ) 0 ≈ 1.5 (K5) at 1-3 Myr. From the percentile distributions described above, the turnover (as assessed by the 50th percentile) is found at lower mass for the older-aged population 60 per cent of the time, at higher mass 3 per cent of the time, and at approximately the same mass for the remainder. We note that the percentile bin size is similar to the average uncertainty on ( BP − RP ) 0 , hence there is an issue of resolution. It would be interesting for future rotation studies of young clusters to further investigate this feature. Figure 13 (top row) also depicts a decreasing fraction of binary candidates moving towards redder colours, which marries with previous findings that bluer, more massive stars are more likely to have companions than redder, less massive ones (Raghavan et al. 2010;Duchêne & Kraus 2013;Lee et al. 2020;Belokurov et al. 2020). However, we also expect binaries to be more difficult to detect in observations of fainter targets, where the signal-to-noise ratio is less favourable. The most prominent difference between the candidate binary fractions of the two age groups is the spike at ( BP − RP ) 0 1 in the 3-6 Myr sample. However, the impact on the percentiles is small, noticeably affecting only the 10th percentile of the 3-6 Myr group for 1 ( BP − RP ) 0 1.5. Figure 14 shows the same period-colour distributions (minus the binary candidates), but this time overlaying a density map to more clearly highlight the relative concentration of points across periodcolour space. From the density distributions, the steeper slope and later turn-over at the red end (in the older population) are emphasised.
17 Uncertainty on colour due to the uncertainty on extinction is, however, the dominant uncertainty. Typical uncertainties on ( BP − RP ) 0 are 0.2-0.4. 18 All of the percentiles in Figure 13 were calculated with a rolling window of width 0.5 in ( BP − RP ) 0 colour. Centre-of-bin plotting was applied, except at the extremes of the distributions (shown with dotted lines), where there is left-side-of-bin plotting at the blue end and right-side-of-bin plotting at the red end, with bin widths of 0.25.
We also see the density distribution for the older-aged population extend to shorter periods at the blue and red end. Figure 15 displays the rotation periods as a function of cluster age 19 . Red, green and blue coloured circles mark the 10th, 50th and 90th percentiles of the period distribution for each cluster. The lower envelope of the distributions, as described by the 10th percentiles, transitions to slightly shorter periods after 3 Myr, in line with the observations in period-colour space:

Period-age distribution
<3 Myr median ≈ 1 d and >3 Myr median ≈ 0.7 d. We note that in plotting all cluster members at a single age, the uncertainty and spread in ages is not represented. The figure is nonetheless informative, so long as the sequence of cluster ages is accurate.

Amplitude
In Figure 16, we plot amplitude (90th−10th percentiles of the flux, converted to a percentage) as a function of colour for the same selection as Figure 14. We see the smallest amplitudes appearing among the bluest stars, which may reflect the smaller spot coverage expected to be present, but we see some high-amplitudes present as well. These high-amplitude signals do, however, correspond with large values of − 2 colour, indicative of a circumstellar disc, where the large flux variations likely originate from accretion bursters or dippers. We note that the increasing − 2 colour trend from left to right is attributable to differing stellar photosphere shapes for different stellar masses, rather than being due to extinction by additional material in the system. Conversely, changes in − 2 colour in the vertical direction, at a particular ( BP − RP ) 0 colour, are indeed likely to be caused by material external to the photosphere.

Disc-rotation relation
Excess emission at infrared wavelengths -thought to originate in the warm dust heated by irradiation from the central star -is often used as an indicator for the presence of a circumstellar disc. In light of this, Figure 17 shows the rotation periods as a function of − 2 colour for stars with identifications found in the literature indicative of the presence or absence of a circumstellar disc. The blue markers locate objects classified as either Class III YSOs (blue circles) or WTTS (blue crosses), and the orange markers locate objects classified as either Class I or II YSOs (orange circles) or CTTS (orange crosses). Additionally, stars found to belong to the category of variable stars known as 'dippers' -objects displaying transient, aperiodic or quasiperiodic dimming events, possibly caused by a warped or clumpy inner-disc as seen from a nearly edge-on viewpoint (Cody et al. 2014) -are highlighted with open black circles (Moulton et al. 2023).
The YSO classifications are based on photometry and were extracted from Hernández et al. (2007), Megeath et al. (2012) and Marton et al. (2016). The T Tauri classifications on the other hand (sourced from Briceño et al. (2019) and Serna et al. (2021)) are derived spectroscopically, based on the relation between the equivalent widths of the H line and spectral types. CTTS and WTTS labels distinguish stars which show or lack evidence of active accretion, respectively. While some WTTS may retain a passive, non-accreting circumstellar disc, the expectation is that there is a high degree of correlation between accretion and the presence of an inner circumstellar disc, as indicated by the Class I or II YSO designation, e.g. Nguyen et al. (2009) find accretion signatures based on H equivalent widths to be highly correlated with 8 m excess in 63/67 cases in their study of T Tauri stars in the young (∼2 Myr old) Chamaeleon I and Taurus-Auriga star-forming regions.
In Figure 17, we observe that the subset of stars displaced to the right, i.e. the population with significant infrared excess, is made up  per cent of WTTS or Class III YSOs. The paucity of short-period rotators with significant infrared excess is consistent with the idea that disc braking plays an important role in the evolution of angular momentum in YSOs.

CONCLUSIONS
We conducted a ∼200-d monitoring campaign across 30 square degrees of the Orion Star-forming Complex. We determined probable members using astrometry from Gaia and corrected for extinction on a star-by-star basis. We reported periodicity 2268 out of 5749 stars and analysed rotation period distributions for 1789 stars with spectral types F0-M5. We assigned stars to clusters within Orion and determined their ages using MIST v.1.2 and Feiden magnetic evolutionary models. The vast majority of rotation periods lie in the range 1-10 d. We observe some evolution in period-colour space between younger and older populations. For older (3-6 Myr) clusters, we notice a shift towards shorter rotation periods for low-mass (>M2) stars, with no periods longer than 10 d among stars later than M3.5. This could indicate a mass-dependence in the dispersal of circumstellar discs. The turnover of the period-colour distribution also occurs at lower mass for the older-aged ensemble, e.g. we see the slow (90th percentile) rotators ( rot ≈ 10 d) shift from ∼K5 (1-3 Myr) to ∼M1 spectral type (3-6 Myr). The fraction of binary candidates decreases towards redder colours in both young and old populations.
Finally, we find that only 4 per cent (1.5 per cent) of CTTS and Class I/II YSOs rotate with periods shorter than 2 d (1.8 d), compared with 17 per cent (14 per cent) for WTTS and Class III YSOs.
ing this for a range of injected periods, would produce (recovered) amplitude distributions as a function of magnitude and period. The tests were conducted as follows. 35 periods were selected spanning 0.05 d to half the baseline of the observations, with a small random jitter added to each. 12 evenly-spaced samples were taken from phase space, again with random jitter. Then, for each star, for each period, for each phase, sinusoidal signals of increasing amplitude were injected, until the injected period was recovered in the Lomb-Scargle detection pipeline. The criterion for detection was that the recovered period fell within bounds based on the injected period, the baseline of observations and the sampling of the Lomb-Scargle algorithm: where is the Lomb-Scargle sampling in period space around the injected period. Figure A1 (top row) displays two views on the injection-recovery results for field NG0535 at a coarse level. The left-hand plot shows the cumulative distribution function for the recovered signal amplitudes, grouped into bins of size two stellar magnitudes, while the righthand plot shows detected amplitude as a function of injected period. The most important variable is stellar magnitude. The Lomb-Scargle periodogram can be thought of in terms of least-squares fits around a constant reference model and a periodic model at each frequency, with best-fit sums of residuals 2 ref and 2 , i.e.
Hence, the periodogram peak height relative to the background noise depends primarily on the signal-to-noise of the data (i.e. the stellar magnitude) and the number of data points (VanderPlas 2018). The final stage of the process was to compare each periodic detection in the main sample to the injection-recovery results. For each object in the main sample, the detected amplitude of its best-fitting sinusoidal signal was compared with the amplitudes recovered in the injection-recovery tests, for stars of similar magnitude and for signals of similar period. The percentile of the detected amplitude relative to the injection-recovery amplitudes (at the corresponding magnitude and period) was then recorded as a score. The implementation was as follows. For each detection, take the injection-recovery results corresponding to the five nearest periods. Then, using a sliding window across the magnitude range of the injection-recovery sample, record the percentile of the target amplitude among the injection-recovery amplitudes for each step of the sliding window. The recorded percentile values were smoothed using a rolling mean filter, before the final percentile score for the target magnitude was obtained via linear interpolation. Figure A1 (bottom row) shows an example of the results of a sliding window calculation for a particular amplitude and period. converted to eff by linear interpolation using the SC table described in section 4.1. Spectral types were converted to integers for this process, i.e. 0-59 for classes B, A, F, G, K, M and their 10 subclasses. Where more than one spectral type was available for a source, the mean was used, avoiding duplicate values from compilation catalogues. Accompanying uncertainties in the spectral types were taken to be two subclasses for stars earlier than M0 and one subclass for M0 and later, which approximates the reported errors for YSOs in the Young Stellar Object Corral (YSOC; Hillenbrand 2021) (see Figure  8 in Cao et al. 2022).

B2 APOGEE Net temperatures
APOGEE Net is a deep convolutional neural network designed to predict eff , log and Fe/H for stars with APOGEE spectra (Abolfathi et al. 2018). Version 1, described in Olney et al. (2020), built on the data-driven approach of Ting et al. (2019), which was trained on Kurucz atmospheric models, to incorporate training labels for PMS and low-mass MS stars based on empirical photometric relations and theoretical isochrones. It yielded properties for stars with eff < 6700 K in the DR14 APOGEE data release. Sprague et al. (2022) extended APOGEE Net to create a pipeline for estimating the parameters of stars across the full mass range in a self-consistent manner, applying it to DR17.
In their study of λ Orionis, Cao et al. (2022) noticed a trend in temperature in the cross-sample of sources with spectral types from the literature and APOGEE Net stars, which they attributed to the fact that the APOGEE Net PMS temperatures are generated from synthetic stars drawn from PARSEC isochrones. In our cross-sample of sources with both APOGEE Net and spectral type temperatures 20 , which covers a much wider range of temperatures than the Cao et al. (2022) sample, we too find disagreement in eff between the two sources, but by way of two separate trends for low and high mass stars. The top-left plot in Figure B1 displays two linear fits in logarithmic temperature space using orthogonal distance regression, with the division being set at eff,ApNet = 4730 K. The lower-left plot shows the cross-sample after correcting for these trends. So as to bring APOGEE Net temperatures in line with temperatures derived from spectral types, this correction was applied to all APOGEE Net stars used. Where both an APOGEE Net and a spectral type temperature existed, the spectral type temperature was adopted. The reason for the greater divergence from a 1:1 relation in the low-temperature domain is uncertain, but one potential contributing factor is the use, by APOGEE Net, of training labels made from synthetic photometry reliant on PARSEC v1.2S stellar models (Chen et al. 2014). PAR-SEC v1.2S models included a shift in the temperature-Rosseland mean optical depth relation, − , in order to reproduce the observed mass-radius relation for low-mass dwarf stars. However, such a correction may not simultaneously be a good recipe for contracting PMS stars in Orion. The shift was applied from 4730 K, increasing towards lower temperatures, which is our reason for placing the division at eff,ApNet = 4730 K. Also in Figure B1, is an illustration of how the residual scatter in the eff,ApNet − eff,SpT relation was combined with the spectral type errors previously described. The residuals were fit with rolling 16th and 84th percentile filters across log eff space, with the mean of the (absolute) percentile values being added in quadrature with the spectral type errors, and then smoothed, to yield eff as a function of eff . The eff values were then used as constraints on eff in the MCMC (section 4.1.2) and were applied to both APOGEE Net and spectral type temperatures. That is, Gaussian priors were placed on eff , with mean values set to the spectral type or (corrected) APOGEE Net temperatures and standard deviations equal to eff .

B3 TIC-8 temperatures
In order to fit the stars without a sourced spectral type temperature or APOGEE Net temperature, we used the values from the TIC-8 catalogue. Effective temperatures in the TIC-8 catalogue are derived in three different ways for the sources in this work: from external spectroscopic catalogues, from the Cool Dwarf List (a carefully vetted list of stellar parameters for K-and M-dwarf stars with eff < 4000 K), or from photometric colours via empirical relations and a dereddening procedure (see Stassun et al. (2019) for details). Figure B2 shows how these temperatures compare with the available cross-sample of spectral type temperatures sourced from the literature. What is clear is that, whilst an approximate 1:1 relation is apparent in three out of the four fields, there is considerable disagreement for the ONCcentred field, NG0535, which is likely attributable to the high levels of extinction affecting observations of these stars. Because of this, we opted to treat TIC-8 temperatures for all fields except NG0535 in the same way as spectral type temperatures from the literature, with errors calculated as described above and displayed in Figure B1. The same approach was also applied to TIC-8 spectroscopic temperatures for objects in field NG0535, i.e. to objects with spectroscopic temperatures in the TIC-8 catalogue, but where a spectral type temperature 20 To increase the cross-sample size, we use all sources within the cluster parallax bounds and precision previously described, i.e. not all objects in the cross-sample are in the cluster members list. from the literature had not been sourced and no APOGEE Net temperature existed. For the other TIC-8 stars in field NG0535, with eff values from the Cool Dwarf List or the standard TIC-8 photometric relations, we attempted an equivalent procedure to that which was applied to the APOGEE Net temperatures (i.e. a linear correction), the results of which appear in Figure B3. We limited the correction to stars with TIC-8 temperatures below 7280 K (∼F0 spectral type Pecaut & Mamajek (2013)), where the scatter is reduced, and above which rotation by the detection of spot-modulation patterns in light curves is not expected. Two extreme outliers were also removed prior to the fit. Large amounts of scatter remain post correction, which is reflected in the final eff error estimates used as constraints in the MCMC and displayed in the bottom-right plot of the figure. This paper has been typeset from a T E X/L A T E X file prepared by the author.  Centre right: eff residuals as a function of log eff along with rolling 16th and 84th percentiles. We take the maximum of the absolute values of the 16th and 84th percentiles of the plotted residuals, rather than the mean, as the error contribution here. Bottom right: Combining (in quadrature) the scatter in the log eff,TIC − eff,SpT residuals with the spectral type errors to give eff error as a function of log eff , used as a constraint in the MCMC (section 4.1.2).