Towards an understanding of YSO variability: A multi-wavelength analysis of bursting, dipping, and symmetrically varying light curves of disc-bearing YSOs

Using simultaneous optical and infrared light curves of disc-bearing young stars in NGC 2264, we perform the first multi-wavelength structure function study of YSOs. We find that dippers have larger variability amplitudes than bursters and symmetric variables at all timescales longer than a few hours. By analysing optical-infrared colour time-series, we also find that the variability in the bursters is systematically less chromatic at all timescales than the other variability types. We propose a model of YSO variability in which symmetric, bursting, and dipping behaviour is observed in systems viewed at low, intermediate, and high inclinations, respectively. We argue that the relatively short thermal timescale for the disc can explain the fact that the infrared light curves for bursters are more symmetric than their optical counterparts.


INTRODUCTION
The photometric variability exhibited by young stellar objects (YSOs) is well established as a key identifying feature (e.g., Joy 1945). For YSOs with circumstellar discs, this variability can arise from processes occurring around the stellar photosphere, or from within the disc. Such varied phenomena can give rise to a variety of light curve morphologies. For example, for some disc bearing YSOs, line-of-sight extinction of the central object by the circumstellar material may manifest itself as short duration dips in the flux. As well as these so-called 'dippers', there are objects whose light curves contain several abrupt increases in flux as a result of an enhanced accretion rate ('bursters'). Sections 3 and 7.1 describes these classes in more detail.
By taking advantage of the unparalleled precision and cadence offered by space-based observatories, recent photometric surveys have been able to examine the light curves of YSOs to an unprecedented level of detail (e.g., Morales-Calderón et al. 2011;Cody et al. 2014;Cody & Hillenbrand 2018). Until recently, however, much of the photometric analysis of YSOs has relied on period searches, e.g. via a periodogram (Lomb 1976;Scargle 1982). Whilst this has been fruitful, only around 30-70 per cent of monitored YSOs are observed to display periodic or quasi-periodic variations (e.g., Herbst et al. 2002;Littlefair et al. 2010;Venuti et al. 2017). Relying too heavily on periodicity searches therefore ignores between one third and two thirds of the disc bearing YSOs. Recent studies (Sergison et al. 2020; ★ E-mail: bsl204@exeter.ac.uk Venuti et al. 2021) tackle this limitation by employing a structure function analysis of YSOs in Cep OB3b and the Lagoon Nebula, respectively. Their analyses have demonstrated that the structure function is a powerful tool to analyse the behaviour of both periodic and aperiodic time series. The difference between this approach and traditional variability analyses can be summarised by noting that the structure function identifies a timescale rather than a coherent period.
In this paper, we use the CoRoT and Spitzer light curves provided as part of the coordinated synoptic investigation of the NGC 2264 star forming region (Cody et al. 2014, see Section 2). We use the morphological classifications provided by Cody et al. (2014). In particular, we focus on objects with CoRoT light curves which demonstrate bursts, dips, and symmetric behaviour (Section 3). For each light curve, we produce a normalised magnitude distribution (Section 4) which shows whether an object is more likely to be found in a brighter or fainter state. Following on from this, we take advantage of the simultaneous optical and infrared light curves to produce single-wavelength (Section 5) and colour (Section 6) structure functions. It is worth noting that the normalised magnitude distributions and structure functions offer complementary analyses. As such, we provide discussion of the analyses from Sections 4, 5, and 6, as well as an overview of the current intepretation of YSO variability, in Section 7 We find results that both are consistent with the current picture of variability of disc-bearing YSOs, and motivate a model of variability with a strong dependence on viewing inclination. Finally, we conclude in Section 8.

DATA AND REGION
NGC 2264 is a young galactic open cluster located in the Mon OB1 association. Located at an approximate distance of ≈ 700 − 800 pc (Park et al. 2000), its rich selection of pre-main-sequence (PMS) objects and relative lack of foreground extinction have made it an attractive target for variability studies (e.g., Alencar et al. 2010;Rebull et al. 2014). Its estimated age of ≈ 3 Myr (Dahm 2008) places NGC 2264 between the sample of ≈ 0.7 Myr-old PMS objects in the Lagoon Nebula (Prisinzano et al. 2019) and the ≈ 6 Myr-old Cep OB3b association (Bell et al. 2013) previously studied with structure functions.
The data used in this paper are taken from the Coordinated Synoptic Investigation of NGC 2264 (CSI 2264) undertaken by Cody et al. (2014). As part of that survey, the authors assign each object a unique identifier ('Mon-XXX') which we will use in place of observatory specific identifiers. A detailed description of data collection and reduction are provided by Cody et al. (2014), which we here summarise.

CoRoT
The optical data used in our analysis were obtained during the fifth short (SRa05) run of the Convection, Rotation, and planetary Transits (CoRoT, Baglin et al. 2006) mission 1 . SRa05 ran from 2011 December 01 to 2012 January 09 and monitored 4235 stars. For this paper, we use the light curves supplied by the CSI 2264 team 2 . These differ from the light curves provided directly from CoRoT in the following ways.
• The WHITEFLUX photometry has been converted to R-band magnitudes by calibrating the logarithmic mean flux to R-band photometry in the literature (Rebull et al. 2002;Lamm et al. 2004;Sung et al. 2008).
• The light curves have been corrected for known systematic effects, most notably the 'hot pixels' which appear in the light curves as sudden discontinuities.
• Some light curves with significant systematics have been entirely omitted from the sample.
The exposure time for most targets was 512 seconds, with a subset of stars with evidence of eclipses or transits being observed with a 32 second cadence mode. This high cadence mode was mainly used for stars with an R magnitude brighter than 14 th . For consistency, the light curves for all of the stars are binned to a 512 second cadence by Cody et al. (2014).
Section 3.2 of Cody et al. (2014) details the light curve production from the basic calibrated data images from the Spitzer Science 1 There are also some optical light curves from observation run SRa01. We choose to omit these data since they lack the simultaneous infrared light curves from which the SRa05 data benefit. 2 https://irsa.ipac.caltech.edu/data/SPITZER/CSI2264/ Centre. The extracted light curves are also made available by the CSI 2264 team 2 . The resulting light curves have ≈ 12 second exposures, separated by ≈ 101 minutes due to the relatively small field of view of the IRAC instrument.

CSI 2264 sample selection
The final sample of Cody et al. (2014) contains 162 disc-bearing YSOs, of which 150 have simultaneous light curves from both CoRoT and at least one IRAC band. 140 of these sources have simultaneous light curves from both infrared bands available. Since we wish to investigate the simultaneous multi-wavelength behaviour shown by the objects in our sample, we only retain the 150 objects with at least one infrared light curve. For details on membership criteria, the reader is directed to the appendix of Cody et al. (2014). In addition to the light curves detailed in Sections 2.1 and 2.2, we also have spectral type information available for many of the objects (Walker 1956;Makidon et al. 2004;Dahm & Simon 2005). Cody et al. (2014) built upon long-standing classification schemes (e.g., Parenago 1954;Herbst et al. 1994) by classifying each light curve into one of the following variability types: aperiodic dippers, aperiodic bursters, stochastic, quasi-periodic dippers, quasi-periodic bursters, quasi-periodic symmetric, long timescale variables, eclipsing binaries, periodic variables, multi-periodic variables, or non-variable stars. Objects were initially assigned a class based on a visual inspection of the light curves, and then used to derive a statistical classification by calculating asymmetry and quasi-periodicity statistics to summarise each light curve (Section 3.2). These classifications were done separately for the optical and infrared data. It is worth noting here that in this paper we will not employ the infrared morphology classifications and instead divide the YSOs by their optical behaviour. Optical morphological classifications are based on the SRa05 data (see Section 2.1). This classification scheme was expanded with a study of the Upper Scorpius and Ophiuchus (Cody & Hillenbrand 2018), as part of the second campaign of the K2 (Howell et al. 2014) mission with the re-purposed Kepler spacecraft (Borucki et al. 2010).

VARIABILITY CLASSIFICATION
This study will focus on the YSOs in NGC 2264 identified by Cody et al. (2014) as having bursting, dipping, or symmetrically varying optical light curves. We briefly summarise these variability classes here. Fig. 1 shows example light curves from each class as well as the corresponding structure function for each waveband (see Section 5 for more details).
From the 150 objects with simultaneous CoRoT and Spitzer photometry (see Section 2), we compile a sample of 94 YSOs with optical light curves displaying bursts, dips, or symmetric variations. The morphological classes are described below. The numbers of bursters, dippers, and symmetric variables used in this paper are 21, 35, and 38, respectively. The objects not used in this paper were identified by Cody et al. (2014) as having either purely periodic (or multi-periodic) variations, non-variable light curves, or light curves which are unable to be classified. We note that, unlike Cody et al. (2014), we group the aperiodic and quasi-periodic variables together into the classes outlined in this section. Since we are interested in the simultaneous optical and infrared variability, we do not include objects for which Cody et al. (2014) do not provide Spitzer and CoRoT photometry. The choice to limit our investigation to the morphology classes outlined in this section means that we retain 63 per cent of the entire sample, and 85 per cent of the objects not identified by Cody et al. (2014) as non-variable, or unclassifiable from the optical light curves.

Light curve classes
In this section, we provide an observational overview of the morphological classes of light curves we study in this work. In Section 7.1, we expand on this and summarise the favoured model for the physical processes behind the variability. Cody et al. (2014) define bursters as YSOs whose light curves display rapid and abrupt increases in flux. An important feature that differentiates this class from other brightening events (e.g., stellar flares, FUOr/EXOr outbursts) is that the flux returns to the quiescent value on a similar timescale to the rise (i.e., the bursts are symmetric). It is also noticeable that bursters typically display many bursts over the course of the ∼ 40 day observations of Cody et al. (2014). This contrasts dramatically to the inferred recurrence timescale for large-amplitude, long-lasting bursts in Class II YSOs of 112 kyr (Contreras Peña et al. 2019). We also note that the bursts observed in this sample are generally low amplitude, this is in sharp contrast to the sample of Cody et al. (2017) which can show variability of up to several times the median flux. Using the Stetson index, Cody et al. (2014) identify stars with well correlated optical and infrared light curves and find that the majority of these are identified as optical bursters or stochastic (Section 3.1.3) stars.

Bursters
Bursts are observed to occur aperiodically (with no well defined period), or quasi-periodically (occurring regularly, but with changes to the light curve between bursts). Cody et al. (2014) identified 21 optically bursting light curves, and 8 objects which display bursting in their infrared light curves. Panels a and b in Fig. 1 show examples of the burster class. The sample of bursters used in this paper consist of the 21 optically bursting YSOs identified by Cody et al. (2014).

Dippers
Dipping behaviour is a well documented phenomenon in the optical light curves of YSOs, a notable example being the prototype AA Tau (Bouvier et al. 1999). These light curves are characterised by significant fading events lasting up to a few days. In objects which have simultaneous optical and infrared dips, the dips are generally observed to be more extreme in the optical (e.g., Cody et al. 2014).
In total, Cody et al. (2014) found 35 optical dippers and 7 infrared dippers from the disc-bearing sample. These are a mix of short duration periodic dips (Stauffer et al. 2015, panel c in Fig. 1), AA Tau analogues (long duration periodic dips which have significant structure at their minima, McGinnis et al. 2015, panel d), or aperiodic dippers (panel e in Fig. 1). Even though several of these objects have a detectable period, they are classified as quasi-periodic since the dip profile will change between dips. Our set of dippers consists of the objects identified as optical dippers by Cody et al. (2014).

Symmetric light curves
Objects which display no clear preference for brightening or fading are labelled by Cody et al. (2014) as 'stochastic' for the aperiodic objects and 'quasi-periodic symmetric' for those which display quasi-periodic variations. Those authors suggested that the aperiodic stochastic variables may exhibit complex behaviour which is a combination of bursting and dipping discussed above. The group identified from their optical light curves as quasi-periodic symmetric variables in Cody et al. (2014) encompasses a range of behaviours. A number of the YSOs display almost sinusoidal variations, possibly attributable to features on the stellar surface rotating in and out of view. There are also light curves which display more complex behaviours which may arise from star-disc interactions. Since we are predominantly interested in the complex phenomena arising from the disc and its interaction with the central object, we exclude objects whose light curves are dominated by sinusoidal variations. 3 Owing to the small number of such objects, this does not appreciably change any of our results or conclusions. Panels f and g in Fig. 1 show examples of the quasi-periodic symmetric class, with panel h demonstrating the aperiodic stochastic class.
In total there are 20 aperiodic stochastic variables, and 18 quasi-periodic symmetric variables selected from the optical light curves 4 . In this paper, we refer to all of these objects as symmetric variables. Cody et al. (2014) initially classified objects by visually examining the light curves and assigning them to one of the classes mentioned above. Additionally, they calculate two statistics to quantitatively analyse the light curves. The first is a measure of flux asymmetry, defined as

Statistical classification
where 10 is the mean of the upper and lower 5% of the data, med is the median of all the data, and is the RMS. Since the data supplied are in magnitudes, a large positive corresponds to dipping behaviour, and a large negative value to bursting. The second statistic is a measure of quasi-periodicity, defined as where RMS raw and RMS resid are the RMS of the raw light curve and a phase subtracted light curve, respectively, and is the uncertainty. The phase subtracted light curve is calculated by folding the raw light curve on the dominant period identified by the auto-correlation function and subtracting this phase trend from the raw light curve.

NORMALISED MAGNITUDE DISTRIBUTIONS
As a measure of variability in the light curves, we use half the 16-84 inter-percentile range (see Sergison et al. 2020). For a normally distributed signal, this is equivalent to the RMS, and is an estimate of variability amplitude which is robust against outliers and reduces the dependence on the number of observations. Fig. 2 shows the cumulative distribution function of the variability amplitude for the bursters, dippers, and symmetric variables. The distributions of 68 show that the overall variability level is larger for the dippers than for the bursters, in all wavelengths. Fig. 3 shows that this is not simply due to the objects occupying a different range of magnitudes. From this half 16-84 inter-percentile range, 68 , we calculate the distributions of normalised magnitudes, defined by Sergison et al. (2020) where is an individual magnitude and is the mean magnitude of the light curve. This is constructed in such a way as to allow us  to compare magnitude distributions between objects (see Fig. 4). A positive value of normalised magnitude indicates a star being in a brighter state than its mean level and vice versa. This is related to the asymmetry metric devised by Cody et al. (2014). We find a very strong anti-correlation between and the Fisher skew of for each light curve 5 . This is unsurprising when one considers that the definition of is a modified skew calculation. We choose to investigate the full distribution as it allows for a more general investigation into the shape of the light curve distribution.
To compare the classes identified by Cody et al. (2014), we combine the distributions of for each source, divided by optical morphology 6 . The top left panel of Fig. 4 shows that the bursters and dippers have rather different distributions of . Although both are asymmetric (which is unsurprising when one considers that these are objects chosen for their light curve asymmetry), the bursters have a significant bright and faint tail, whereas the dippers have a very sharp drop-off for a normalised magnitude greater than about two. The symmetric group predictably has a more symmetric distribution than either the bursters or dippers, but the distribution is generally more similar to the distribution for the bursters.
The lower two panels of Fig. 4 show the same distribution for the 3.6μm and 4.5μm data, respectively. We group the light curves by their optical morphology, with the same colours as in the first panel.
In both distributions, the optical burster Mon-804 has been removed since the infrared light curves demonstrate anomalous behaviour which skews the study of the bursters as a whole (see Appendix A). We can see that both the dippers and bursters have normalised magnitude distributions comparable to the symmetric class. To verify that this is not an effect of the different sampling of the light curves, we resampled the R-band light curves to the IRAC cadence. The top right panel of Fig. 4 shows the resulting normalised magnitude distribution. The distributions for the dippers and bursters keep their characteristic shapes after downsampling, demonstrating that the difference between the optical and infrared distributions is a genuine effect of wavelength, not simply a difference in sampling. For clarity, any further references in this work to the CoRoT data are to the full, non-downsampled, data. We will discuss these results along with those from Sections 5 and 6 in Section 7.

STRUCTURE FUNCTIONS
Whilst the normalised magnitude distributions detailed in the previous section demonstrate the range of magnitudes YSOs may explore, it does not tell us anything about the timescale of this variability. For this we wish to analyse the correlation within the light curve. To this end, we follow the prescriptions of Sergison et al. (2020) and Venuti et al. (2021) and use the structure function, since many of the objects in our sample show aperiodic variability. Other studies which have used similar analyses to investigate aperiodic behaviour of YSOs include the 'pooled sigma' (Rigon et al. 2017) and Δm-Δt (Findeisen et al. 2015) studies, however we aim to compare results with Sergison et al. (2020) and Venuti et al. (2021), and so opt for the structure function.
Whilst the structure function is a longstanding analysis tool in the field of extra-galactic astrophysics (e.g., Simonetti et al. 1985;Hughes et al. 1992;de Vries et al. 2003), its use in the study of YSOs is relatively new. The first extensive study of YSOs with structure functions was undertaken by Sergison et al. (2020) who used -band photometry of over 800 YSOs in the Cep OB3b association at timescales ranging from one minute to ten years. The choice of the structure function as an analysis tool was motivated in part by their unevenly sampled observations, and the need to be able to analyse aperiodic signals. Amongst their results, they found that the -band variability at all timescales is systematically lower for later YSO evolutionary class. In addition, the authors found that the observed structure functions for Class IIs could be explained if all Class IIs exhibit a roughly power-law variability spectrum driven by accretion rate changes with some Class IIs exhibiting an additional periodic component arising from the rotation of the star.
Later, Venuti et al. (2021) used structure functions to analyse YSOs in the Lagoon nebula. Unlike Sergison et al. (2020), those authors benefit from evenly sampled data from the ninth campaign of the K2 (Howell et al. 2014) mission. Those authors divide the sample into the variability types introduced in Cody et al. (2014, see Section 3.1). We will compare our findings with those of Venuti et al. (2021) in Section 7.3 The CSI 2264 data set is unique when compared to these previous studies, in having simultaneous time series in optical and infrared wavelengths. This allows us not only to examine the optical or infrared variability spectrum, but also to examine the colour variability (see Section 6). Additionally, like the data set of Venuti et al. (2021), we also have information about the spectral types of the objects, with the sample covering the G to M star range (Walker 1956;Makidon et al. 2004;Dahm & Simon 2005).

Definition
For a light curve with magnitudes ( ) observed at times , we define the structure function as where the average is taken over all pairs of observations separated by a timescale . In practice, this was calculated for logarithmically spaced bins in . ( where the sum is taken over the ( 1 , 2 ) pairs of observations that satisfy Throughout the entirety of this paper, we require each structure function bin to have ( 1 , 2 ) ≥ 10 to mitigate the effect of small number statistics.
Note that we follow the example of Hawkins (2002) and de Vries et al. (2003) by calculating the structure function in magnitudes, rather than flux (e.g., Simonetti et al. 1985;Hughes et al. 1992;Sergison et al. 2020). We do this because we believe that the ratio of fluxes is a better metric of variation than the difference. That is to say that we believe a 90 per cent rise in flux is a less significant change than a 90 per cent drop in flux since the latter is a much larger fractional change. Calculating the structure function in flux would treat these two changes equally. Of course, this may be more naturally achieved by calculating the structure function in log (Flux) or ln (Flux), but we opt to use magnitudes for historical reasons.  The square root of the structure function at a timescale , √︁ ( ) is the RMS of all observations separated in time by . A small value of the structure function indicates that observations separated by the corresponding timescales exhibit little variation; if we know the magnitude at a time , it is easier to predict the magnitude at a time + . Conversely a large value of the structure function implies there is little correlation between points separated by the corresponding timescale.
A typical schematic of a structure function is shown in fig. 7 of Sergison et al. (2020). At short timescales, the structure function is dominated by the measurement uncertainties and point-to-point scatter of the light curve. There is then a region of increasing variability with timescale as we begin to probe the timescales for relevant physics. This tends to follow a roughly power law relation. Finally, when we are probing timescales longer than the longest timescale of the variability, the structure function will no longer be increasing, reaching a plateau. This usually corresponds to a 'knee' in the structure function at the characteristic timescale (see Section 5.3.1).
A general feature of structure functions is that (for most types of variability) ( 1 ) ( 2 ) provided 1 > 2 . A notable exception to this rule is for a purely periodic signal since observations separated by multiples of the period would have zero variability. For example, the structure function corresponding to a signal given by ( ) = sin ( ) is ( ) = 1 − cos ( ) (Findeisen 2015). This gives rise to a characteristic shape in log-log space (see curve 1 in Fig. 5). Fig. 9 in Sergison et al. (2020) demonstrates the effect of finite sampling on the shape of the structure function. Table 6 in Sergison et al. (2020) compares the power law index for some well characterised noise profiles and we expand on this discussion in Appendix C.

Comparison to PeakFind
Cody et al. (2014) employ a conceptually similar analysis to calculate a variability timescale for aperiodic light curves. Their approach is, for a given variability amplitude , to select maxima and minima in the light curve that vary by at least . The corresponding timescale for this amplitude is the median time between the consecutive minima and maxima. This 'PeakFind' algorithm can then be used to produce an amplitude-timescale (or timescale-amplitude) plot, shown in fig. 32 in Cody et al. (2014). This amplitude-timescale plot is then used to derive a timescale (corresponding to the 70% of maximum amplitude).
Since PeakFind only considers the relative locations of minima and maxima, any signals which have coinciding extrema will appear identical to PeakFind, regardless of the behaviour of the light curve between the extrema. This is demonstrated in Fig. 5. The first panel shows three basic 'light curves'. We have included a sinusoid (curve 1), periodic Gaussian dips from a constant continuum (curve 2), and alternating positive and negative -like functions 7 (curve 3). These are constructed to be dramatically different signals whilst having identical minima and maxima. This ensures that PeakFind treats them identically despite their different shapes. Panel 2 demonstrates that the structure function is able to distinguish between these cases. Notably, we highlight that the 'knee' timescales identified in curves 1 and 2 correspond to the 'width' of the sinudoidal variation and Gaussian dips, respectively. By contrast, both structure functions show a dip at the true period of ≈ 7 days. For a more thorough description of the relative benefits of PeakFind analysis, the reader is directed to Findeisen et al. (2015).

Results
We calculated structure functions from the CoRoT and Spitzer light curves for each object in our sample. Fig. 6 shows the structure functions calculated from the CoRoT light curves and Figs. 7 and 8 show the same for the IRAC 3.6 μm and 4.5 μm data, respectively. As in Section 4, we divide both the optical and infrared structure functions by the optical morphology. In each plot, the structure functions are colour coded by spectral type, and shown with solid or dashed lines to indicate whether the object varies aperiodically or 7 Strictly this is a top hat function with a single observations in the brighter or fainter state. (i.e., into aperiodic and quasiperiodic variables) statistics. Note that in both cases, the division is done on the optical light curves. Columns indicate the waveband or colour (see Section 6) the structure functions are calculated from. Objects with fewer than 150 observations in a particular light curve are omitted from the relevant count.
In Fig. 6 the dipping objects clearly exhibit stronger variability than the bursters, especially at longer timescales, in agreement with Fig. 2. This indicates that the dips are generally larger amplitude than the bursts. The systematic difference in variability amplitudes we demonstrate here is not an artefact of the selection performed by Cody et al. (2014) as that selection is made only on the light curve asymmetry and is therefore independent of the overall variability level. We therefore believe this to reflect a genuine difference between the variability classes.

Knees in the structure functions
We have identified structure functions with clear knees in their structure functions (i.e., structure functions which do not show increasing variability at the longest timescales accessible from these data). Table 1 presents the results from this. It it clear that the fraction of structure functions with knees compared to those without varies between the variability types, as well as observational waveband. The majority of IRAC structure functions demonstrate variability that is still growing at longer timescales. It is worth noting that the CoRoT light curves have a longer baseline and higher cadence, so we are naturally sensitive to a larger range of timescales in the CoRoT structure function. Even considering this, it appears that the infrared variability generally has a longer maximum timescale. Also of note is that some stars whose structure functions appear to still be increasing at the longest timescales also show the oscillating behaviour attributed to periodic signals and there are some structure functions which demonstrate a second rise after the plateau. This corresponds to longer term, approximately linear, trends (e.g., Mon-1187) or longer period variations (e.g., Mon-525) being imposed on the continuum flux. This demonstrates that the structure function is sensitive to multiple variability modes.
Following Sergison et al. (2020), we divide the light curves by whether they contain a knee. Those authors find that the structure functions with knees display more variability than the structure functions without knees, despite the selection being independent of overall variability level. We do not find this to be the case for all morphology classes. Specifically, we only find a significant difference for the dippers. This can be seen clearly in Fig. 9 which  shows that the median structure functions with knees and without knees are very similar for the bursters and symmetric variables, especially at longer timescales. The black lines group the structure functions as in Sergison et al. (2020) (i.e., by presence or absence of a knee in the structure function, with no knowledge of light curve morphology). In doing so, we recover a similar result to that obtained by Sergison et al. (2020). It is possible therefore that Sergison et al. (2020) see a difference between the amplitude of the structure functions with and without knees because their sample is dominated by dippers. This dominance of dippers was also suggested by Sergison et al. (2020) to explain the shape of the Class II normalised magnitude distribution. We also considered the possibility that the difference between objects with and without detectable knees were due to the lower-amplitude variables being drowned out by noise, and thus not having sufficient signal-to-noise to probe astrophysical effects. This posed two difficulties however; firstly, it is not clear why such an amplitude dependence should affect the dipper population more than the bursters and symmetric variables. Secondly, if it were intrinsically harder to identify physical timescales for the lower-amplitude variables, we might expect that the aperiodic sample of Cody et al. (2014) to have lower overall variability than the quasi-periodic sample. Comparing the 68 distributions of these two samples, we find no such difference.
It is also interesting that of the 71 YSOs whose optical structure functions have visible knees, 34 are identified as aperiodic by Cody et al. (2014). Similarly, a large fraction of the infrared structure functions with knees are identified as aperiodic in Cody et al. (2014) (44 per cent for IRAC1 and 45 per cent for IRAC2). This illustrates that the the structure function is sensitive to the variability timescale, rather than simply identifying periods. To highlight the ability of the structure function to identify timescales in aperiodic objects, we generate five artificial light curves with Gaussian dips (Fig. 10). The locations and widths of the Gaussian dips are chosen at random 8 , thus ensuring aperiodicity. Despite this, each corresponding structure function clearly has a knee. It is worth highlighting that the location of the knees are not the same in each structure function, despite the average time between dips being the same. This is to due to the light curves having dips of different durations, and thus different timescales. The difference between a structure function timescale and a period can also be seen in Fig. 5, where the structure functions have knees at different timescales despite being calculated from synthetic light curves that were constructed to have identical periods. Hence we highlight that the structure function can identify timescales in signals, even when a period search may fail, but that the timescale identified is not necessarily a period.

Timescales
Following Sergison et al. (2020) and Venuti et al. (2021), we extract the longest timescale of variability by fitting the following power law to all structure functions with knees, where = brk / brk is defined to ensure continuity. Optimal values of brk , brk , and are found by a least squares fitting to the logarithm 8 The depths of the Gaussian dips are fixed to be proportional to the widths. . For YSOs without knees, we simply fit a power law relation to the structure function (i.e., the brk → ∞, brk → ∞ limit of Equation 7). Figs. 11 and 12 show the distributions of brk and derived from the CoRoT and Spitzer structure functions, respectively. Table 2 shows the median and RMS values of for the CoRoT and Spitzer wavelength structure functions. Due to the limited number of observations for some of the IRAC light curves, we only extract the above parameters from IRAC structure functions corresponding to light curves with at least 150 observations. This leaves 69 objects for which we fit Equation 7 to the 3.6μm structure functions, and 76 objects for which we do the same for the 4.5μm structure functions. In all bands, the median value of brk for the dippers is less than that for the bursters and symmetric variables, which are similar to each other. We find a median optical timescale of 2.1, 1.5, and 2.3 days for the bursters, dippers, and symmetric variables, respectively.
In Fig. 13 we show how the optical and infrared timescales are related. We find an overall correlation, with most of the object having a near 1:1 correspondence. Cody et al. (2014) perform a similar analysis and find that when periods are detected in the optical and infrared light curves, that they are similar. It is important to note the difference between our Fig. 13 and fig. 34 in Cody et al. (2014), namely that we only calculate timescales for objects with knees in both the optical and infrared structure functions, whereas Cody et al. (2014) plot timescales for aperiodic objects. The 1:1 correlation we find suggests a common physical process for the optical and infrared behaviour.

MULTI-WAVELENGTH STRUCTURE FUNCTIONS
Using the multi-wavelength data set at our disposal, we were able to produce time series of the CoRoT-Spitzer colours (see Appendix B). From these colour curves, we calculate the structure function in the same way as for the single-band light curves. Whilst Venuti et al. (2021) have simultaneous K2 and VST/OmegaCAM (Kuĳken 2011) observations, the latter is too sparsely sampled (with only 17 observations over the ∼ 70 day campaign) to enable a full multi-wavelength structure function analysis.
The choice to calculate the structure function in magnitude space (see Section 5.1) as opposed to flux space is further justified by the fact that the colour structure function we calculate here is independent of our choice of colour (i.e., choosing to calculate the colour as Spitzer-CoRoT rather than CoRoT-Spitzer would yield identical results). As with the monochromatic structure functions, the value of the colour has no effect, only the relative change (i.e., we cannot tell how red or blue an object is, only how its colour varies at different timescales). Fig. 14 shows the colour time series and colour structure functions for the example YSOs in Fig. 1. For comparison, we also include the respective single-waveband structure functions (right panels of Fig. 1) as faint dashed lines.  Figures 15 and 16 show the structure functions of the colour. It can be seen that the G stars, which show somewhat lower variability than their later-type counterparts in the first panel of Fig. 6, display very similar colour behaviour to the K and M type dippers. In addition, the colour series for the dippers are more variable than those for the other variability types. This implies that the optical and infrared light curves are less correlated for dippers than they are for bursters. This is especially noticeable in the R-[4.5] colour. As in Section 5.3.1, we visually examine the structure functions to identify those objects with visible knees (see Table 1). Similarly to the single wavelength structure functions, we use Equation 7 to extract the best fit parameters which we plot in Figs. 17 & 18. The colour timescales in both cases seem reasonably similar between the morphology classes, though our inferences are limited by the relatively small number of structure functions with knees. We also see that the dipper colour variability is much steeper than that of the bursters or symmetric variables.

Current model
Before discussing the observational findings of this paper, we wish to first summarise the current favoured model for YSO variability in order to better place our findings into context.

Dippers
Optical dipping behaviour is common in YSOs, accounting for around 30 per cent of all disc bearing stars in a given region (e.g., Stauffer et al. 2015;Cody & Hillenbrand 2018;Roggero et al. 2021). It is believed that the cause of this variability is extinction by circumstellar material crossing our line of sight to the central object. There are a number of specific models to explain various light curves, for example a warped inner disc wall, caused by a tilted stellar dipole in the prototype AA Tau (Bouvier et al. 1999), dust trapped in accretion funnels rising above the plane of the disc (Stauffer et al. 2015), or clumps of dust being produced by instabilities in the disc and levitated into our line of sight (Turner et al. 2010). A key observational prediction of these models is that dips should be fainter in MIR wavelengths than in optical wavelengths due to the ∼ μm sized dust particles in the circumstellar material (e.g., Cody et al. 2014).

Bursters
Bursting behaviour is attributed to flux from hot spots on the stellar surface caused by accretion shocks (Herbst et al. 1994). These bursts differentiate themselves from flares caused by stellar activity by their roughly symmetric shape (i.e., the rise and decay timescales of the burst are roughly equal). For a more comprehensive discussion of current explanations for bursting behaviour, the reader is directed to section 7 of Cody et al. (2017). Cody et al. (2014) find that the majority of the objects with strong correlation between the optical and infrared light curves are classified as bursters or symmetric variables. This is an important prediction of the accretion shock model, as both optical and infrared wavebands are in the Rayleigh-Jeans tail of the ∼ 10 6 K accretion material (Calvet & Gullbring 1998). We argue in Section 7.2 that rotational effects also strongly contribute to the bursting phenomena.

Comparison with observation
We present two pieces of evidence which support the above picture. We have shown in Section 6 that the variability in dippers is more chromatic than in the bursters, as expected from the variable extinction model of dipping variability. Furthermore, the shapes of the normalised magnitude distributions presented in the first panel of Fig. 4 show that the dippers have a much better defined maximum flux than the bursters do. We argue that this is also consistent with the current explanations. The reasoning is as follows. If the dippers are caused by occultation of circumstellar material, and the variation in continuum flux is relatively slow, there is a well defined maximum in brightness corresponding to the light from the unobscured star. By contrast, variability caused by accretion rate would presumably have a more significant bright and faint tail given that the accretion rate could increase or decrease.

Light curve morphology as an effect of inclination
It is commonly reasoned that dipping systems must be observed at high inclinations since the central object must be variably occulted by asymmetries in the innermost disc regions or the stellar magnetosphere. In this section, we will explore how a similar inclination argument might explain the difference between the bursters and symmetric variables. In particular, we argue that systems observed at low inclinations are more likely to display symmetric variations, and systems at intermediate inclinations are more likely to be classified as bursters.  Robinson et al. (2021) use a composite model consisting of a one-dimensional hydrodynamic simulation of the accretion column (Robinson et al. 2017), non-local thermal equilibrium accretion shock (Calvet & Gullbring 1998;Robinson & Espaillat 2019), and rotational modulation to produce synthetic optical light curves of YSOs with accretion modulated variability. They assume aligned rotational and magnetic axes, and that the in-falling material is uniformly distributed over a single hot-spot. Robinson et al. (2021) showed that by varying magnetic field properties and viewing inclination of accretion dominated YSOs, the light curves could occupy the symmetric or even dipping regions of the -plot (see Section 3.2). Fig. 9 of Robinson et al. (2021) demonstrates the general trend for the model light curves as the inclination is varied from 90 • to 0 • : light curves vary from periodic to aperiodic and from purely bursting to symmetric or dipping behaviours, as characterised by the asymmetry statistic of Cody et al. (2014)  (see Section 3.2). The dipping behaviour is especially interesting considering the models of Robinson et al. (2021) do not include the effect of occulting disc material and serves to highlight the stochastic nature of accretion.

Accretion variability
The explanation of why changing the viewing inclination can change the morphology of the light curve is as follows. The light curve consists of a contribution from a reasonably constant photosphere, and a bright hot spot whose flux varies with the accretion rate. As the spot rotates into view, the total flux increases (and varies stochastically). So for low inclination systems, where the spot is in view for longer, the light curve is dominated by the stochastic variations driven by the variable accretion rate, with potential dips as the spot rotates out of view and we are left with just the photospheric flux. By contrast, the high inclination systems have a much smaller fraction of the rotational phase at the brighter flux level associated with the spot, and so we would describe these light curves as bursting more. Additionally, there are a number of models which would be classified as dippers despite being very low inclination models, and thus having the hot-spot visible for the entire rotational phase, highlighting the effect of the stochastic driving function. Fig. 19 shows a simplified illustration of how inclination can affect light curve asymmetry. Unlike the more models of Robinson et al. (2021), we do not offer a full treatment of accretion variability. Instead we consider a constant photosphere with a single, infinitesimally small hot spot which shows stochastic variability. We place this hot spot at a latitude of 75 • and calculate the resultant light curve when viewed from a range of inclinations (note we allow > 90 to account for the case where we are viewing the opposite side of the star to the spot). We also add a small level of white noise to simulate observational uncertainties. To illustrate the effect of inclination alone, the accretion and white noise have identical time series for each model. The red points indicate the observations when the spot is in view and the black points indicate we only see the photosphere. Clearly the different viewing angles give rise to a range of light curve behaviours, despite the simulated accretion noise and white noise being identical between each instance. Note that, as in the simulations of Robinson et al. (2021), we do not include the effect of circumstellar discs. We wish to emphasise that we are not intending this to be taken as a physical model of accretion. Instead, Fig. 19 should serve simply as an illustration of how varying inclination can manifest itself as different light curve behaviour.

Effects of the circumstellar disc
As previously stated, neither the schematic diagram of Fig. 19 nor the more sophisticated treatment of accretion variability in Robinson et al. (2021) include the effects of the circumstellar disc. In this section, we briefly discuss how we expect the inclusion of the disc to affect the results of the model.
The first, most obvious, effect occurs when the disc intersects our line of sight to the star. In this case, we would expect the optical light curve to show dips of the form discussed in Section 7.1.1. As we have previously discussed, these dips will be less prominent in the mid-infrared wavelengths of the Spitzer observations. A second, more subtle, effect is that as the disc absorbs light from all directions from the star, it reaches thermal equilibrium and re-emits the light in longer wavelengths, thus acting to provide a measure of the azimuth-integrated flux from the star. Therefore by observing in infrared wavelengths, which more efficiently probe the disc flux, we can get information about the total light emitted by the star, largely unaffected by the rotational phase of the star. In principle, the signal originating in the disc should lag the stellar signal due to a combination of light-travel time and the disc thermal timescale. Radiative transfer simulations (e.g., Harries 2011) and photo-reverberation mapping studies (Meng et al. 2016) of circumstellar discs show that this lag occurs at a shorter timescale than our cadence allows us to sample.
We demonstrate in Section 4 that the infrared normalised magnitudes for all variability types demonstrate similar distributions. This is perhaps unsurprising considering the large difference in wavelength. Nevertheless it seems contrary to the observations of Cody et al. (2014), who find that a large number of bursters and symmetric variables show strong correlation in the optical and infrared light curves. This might imply that the shape of the magnitude distributions should be reasonably independent of wavelength.
We can explain this by noting that in the infrared, we are more sensitive to the light from the disc than in the optical. This implies that in the Spitzer photometry, we are seeing the accretion variability without the rotational modulation highlighted in Fig. 19, similar to the optical light curves of the symmetric stars. Hence the bursters and symmetric variables have similar magnitude distributions in the mid-infrared wavelengths we study. Additionally, we can also see in Fig. 2 that the overall variability of the symmetric variables and the bursters are more consistent in the infrared than in the optical. Finally, the similarity between the classes can also be seen quantitatively by performing two-sample Kolmogorov-Smirnoff tests between each pair of distributions. In doing this, we find that the distributions of power law index between the bursters and symmetric variables are considerably more similar in the infrared wavebands than the optical. It is worth noting that these three pieces of evidence are complementary; the similarity of the magnitude distributions demonstrates that the infrared light curves explore a similar range of magnitudes, the half 16-86 inter-percentile range demonstrates that the overall level of variability is similar, and the power law indices show that the light curves have similar levels of structure.
Following this logic, the infrared structure functions should probe the accretion variability. Taking this to be the case, we find the accretion to follow a power law with = 0.55 ± 0.29 from the [3.6] structure functions, and = 0.70±0.30 from the [4.5] structure functions 9 . Note that the range quoted here is the RMS of the individual power law indices as opposed to an uncertainty for . Although this is a reasonably large spread in , it is important to note that this is the structure function power law index, not the more conventional Fourier spectral index (see table 6 for comparison). The spread may also represent an intrinsic spread in the accretion profiles between different objects. The values of structure function power law indices observed here correspond to Fourier power spectra roughly comparable to a random walk (see Appendix C for a discussion of this comparison). This is also consistent with the 'long' subclass of Sergison et al. (2020), which have values of varying between ≈ 0.3 and ≈ 1.2. This indicates that the power law relation we observe at timescales of roughly hours to weeks, may extend to timescales of years as in Sergison et al. (2020). Whilst we include the values of of the entire sample, rather than just the objects without a knee in the structure function, we argue that this is a valid comparison in two ways. Firstly, there are relatively few objects with knees in the infrared structure functions, and so their effect on the overall results will be minimal. Secondly, the near-infrared data of Sergison et al. (2020) will be more sensitive to the rotation modulation from the central object than the mid-infrared data used here, requiring those authors to apply a selection to identify light curves without strong rotational effects. This is not necessary in our sample due to the longer-wavelength photometry.

Threshold inclinations
Since we suspect that dipping systems are observed edge on, owing to the requirement that the optically thick circumstellar material intersect our line of sight, we expand on the predictions of Robinson et al. (2021) and suggest a simple picture of the variability types, where dipping, bursting, and symmetric behaviour arises from systems observed at high, intermediate, and low inclinations, respectively. To estimate the threshold inclinations for these behaviour types, we assume that the systems in our sample have randomly distributed orientations. Therefore the fraction of systems with inclination less  Table 3. Population-derived threshold inclinations for observing bursting and dipping behaviour in the light curves of disc-bearing YSOs in a number of clusters.
than a particular 0 is given by We can therefore use the fraction of the sample exhibiting dipping, bursting, or symmetric behaviour to calculate the corresponding inclinations. Table 3 shows the threshold inclinations for light curves to demonstrate bursting ( th, B ) and dipping ( th, D ) behaviour in NGC 2264, Oph, and Upper Sco, calculated from the Cody et al. (2014) and Cody & Hillenbrand (2018) variability fractions, respectively. We also include the threshold inclinations derived from the discbearing young stars in the Lagoon Nebula (Venuti et al. 2021). We note that, with the exception of the Lagoon Nebula sample, there is remarkable consistency between the different samples, with threshold inclinations varying by only a few degrees between the different star forming regions. McGinnis et al. (2015) investigate objects in this sample which show behaviour similar to the prototype AA Tau. The inclinations derived from the extinction models in that work are all consistent with the dipping threshold in Table 3.

Caveats
There are a couple of important caveats to this picture and the inclination values quoted in Section 7.2.3. Firstly, in quoting these population averages, we have assumed that all the dippers are observed at the highest inclinations, whereas the simulations of Robinson et al. (2021) demonstrate dipping phenomena arising from low-inclination systems. Since the majority of the simulations in that paper result in either bursting or symmetric behaviour, it is reasonable to assume that most of the observed dippers are caused by line-of-sight extinction.
Secondly, the threshold inclinations separating the morphology types we calculate in Section 7.2.3 are calculated from the entire population of variable Class II YSOs in each region. The simulations of Robinson et al. (2021), whilst predominantly showing the same general trend, have different threshold inclinations separating the bursters and symmetric stars, depending on the parameters of each model. Similarly, the minimum viewing angle to observe dipping will strongly depend on the geometry of the disc since we require at least part of the disc to intersect our line of sight. For example, a star with a misaligned disc (e.g., Davies 2019) would presumably have a different range of viewing angles which would produce dipping behaviour. We note that we have implicitly assumed that all the stars in the sample have the same threshold inclinations, and therefore that the stars have similar magnetic and disc properties. We must therefore emphasise that the values quoted in Table 3 are population-level statistics, rather than definite thresholds for each individual object.
A further complication is highlighted in the light curves of Mon-119 and Mon-577. These stars have symmetric light curves which show narrow dips attributed to dust in accretion columns rising above the plane of the disc. Stauffer et al. (2015Stauffer et al. ( , 2016 estimate that these objects are observed at inclinations 60 • owing to their lack of deep dips. This picture implies that the objects categorised as dippers may also exhibit bursting phenomena, given that they lie at larger inclinations than our supposed threshold for bursting behaviour to manifest itself. As shown in Sections 4 and 5, the dipper light curves are generally more variable than their bursting counterparts. Therefore it is possible that, while both dipping and bursting phenomena are present, the latter is drowned out by the former. Further, there is evidence for variable accretion in many dippers (Bouvier et al. 2003(Bouvier et al. , 2007McGinnis et al. 2015), indicating that even in the picture described in Section 7.1, there is a significant overlap between the bursting and dipping objects.

Summary
We suggest a model for TTS variability dependent on viewing angle to the star. Recent simulations by Robinson et al. (2021) have demonstrated that inclination can have a significant effect on the optical behaviour of the system by introducing a rotational modulation. By additionally considering the effect of a disc, we can also describe the observed infrared behaviour by arguing that the disc reprocesses light from all rotational phases, leading to much more similar behaviours for bursters and dippers. Consequently, we argue that the variability in the infrared light curves of the symmetric variables and bursters directly probes the accretion variability. Following this logic, we find that the majority of burster and symmetric variables display accretion signatures roughly consistent with a random walk process.

Comparison with Venuti et al. (2021)
A seemingly natural comparison to draw is between this work and the work of Venuti et al. (2021), owing to the similarity in the data, and the analysis methods. However, there are significant differences between the ages and spectral types of the stars studied in these two works, and we are limited in the conclusions we can draw from any comparison of them. Nevertheless, we present here a brief comparison of the differences between our findings and those of Venuti et al. (2021).
A key result of Venuti et al. (2021) is the trend for stars of earlier spectral type to show systematically less variability than the later type stars. Additionally, there is a over-representation of non-variable stars in the early type sample, compared with the late type sample. Neither of these effects are observed in the sample studied here. Figs. 6, 7, and 8 do not show the trend of increasing variability with later spectral type found by Venuti et al. (2021). Indeed, perhaps with the exception of the CoRoT light curves of the dippers, the spectral type of the star appears to bear no relation to the structure function. Additionally, the over representation of non-variable stars seen in the early type objects in Venuti et al. (2021) is absent from the CSI 2264 data set, where the non-variable fraction of the G, K, and M stars is 13, 4, and 25 per cent, respectively. We find in section 5.3 that the dippers exhibit more variability than the bursters at all but the shortest timescales. We note that this is contrary to the findings of Venuti et al. (2021) who find that the maximum intrinsic variation for bursters is greater than that for dippers or symmetrically varying stars.
Finally, in section 5.3.2, we find that the bursters and symmetric variables have longer median timescales than the dippers. In contrast, Venuti et al. (2021) find that the dippers have the longest timescale ( brk ∼ 3 d), with the bursters and symmetric stars having brk ∼ 1 − 2 d.
As earlier highlighted, it is perhaps unsurprising that the objects in this sample and those studied by Venuti et al. (2021) show different properties, given the difference in region properties. Nevertheless we wish to draw attention to the differences in this section.

CONCLUSIONS
In this paper we have used the simultaneous optical and infrared photometry of 94 disc-bearing YSOs made available as part of the CSI 2264 project (Cody et al. 2014) to examine the variability properties of bursting, dipping, and symmetrically varying stars. We have analysed YSOs with variability signatures, classified by their optical light curve morphology as in Cody et al. (2014). Our primary conclusions are detailed below.

Observational findings
We find in Section 4 that the magnitude distributions for the optical light curves of dippers and bursters are not simply mirror images of one another but show qualitatively different properties. By contrast, the infrared flux histograms are much more symmetric for each morphology type.
In Section 5, we use structure functions to investigate the variability of the YSO light curves at a range of timescales. We find that the majority of the infrared light curves show increasing variability at the longest timescales accessible by the observations. This contrasts starkly with the optical light curves which generally show a plateau (or 'knee') in their structure function, corresponding to the maximum timescale for variability. We emphasise that this is a timescale, not necessarily a period. This is notable since many objects classified as aperiodic variables by Cody et al. (2014) show knees in their structure functions. We also find that, similar to the results of Section 4, the infrared structure functions for the bursters and symmetric variables show considerable similarity, with similar distributions for the structure function power law index.
Using the simultaneous optical and infrared data at our disposal, we perform, to the best knowledge of the authors, the first structure function analysis of colour time series data (see Section 6). In doing so, we see that the variations in the bursters are significantly less chromatic than the variations in the dippers. We argue that this is consistent with the current picture of dipping and bursting behaviour, being caused by variable extinction and accretion, respectively. We note that the symmetric variables demonstrate colour behaviour which lies broadly in the middle of these two classes.
In Section 4 we demonstrate that the dippers in this sample show a larger variability amplitude than the bursters and symmetric variables. In Section 5, we show this effect to be persistent at all timescales longer than a few hours. Interestingly, when Venuti et al. (2021) perform a similar analysis of disc-bearing YSOs in the Lagoon nebula, they find that the bursters have a larger intrinsic variability than the dippers, and that both are more variable than the symmetric class. We wish to highlight the fact that the significant differences between the regions limits any conclusions we can draw from this comparison.

Interpretation
In this section, we aim to place the findings of Section 8.1 into context by providing the interpretation of the authors. We firstly note that the findings of Section 6, that the dippers have variability which is more chromatic, is consistent with the findings of Cody et al. (2014) the current accepted model of YSO variability (Section 7.1). To this end, we also argue that the difference in optical magnitude distributions in Section 4 is also evidence for this model, as the dippers have a much smaller 'bright tail' than the 'faint tail' of the bursters. This is, we argue, consistent with the idea that the dips are caused by extinction of circumstellar material and that the variation in the bursters are caused in part by a change of accretion rate (Section 7.1).
Additionally, we suggest a model of YSO variability dependent on viewing angle to the star (see Section 7.2), motivated in part by recent simulations by Robinson et al. (2021). We argue that symmetric, bursting, and dipping light curves will generally arise from systems with low, intermediate, and high inclinations, respectively. This model suggests that the primary difference between the bursters and symmetric variables is the presence or absence of rotational modulation of the accretion hot spots.
Finally, we point out the importance of the fact that the disc reprocesses light from all directions of the star. Thus by observing the system in mid-infrared wavelengths which more effectively probe the discs, one is able to negate most of the rotational modulation which separates the symmetric and bursting classes in the optical wavelengths. Invoking this model can explain why, in the mid-infrared observations, the structure functions and light curve histograms of the bursters and symmetric variables behave similarly. It also suggests that the mid-infrared variability is a direct probe of the accretion rate changes, largely free from the effects of rotation. Assuming this to be the case, we find that the power law relation for the accretion variability time series has index = 0.55 ± 0.29 and = 0.70 ± 0.30, derived from the [3.6] and [4.5] data, respectively. This is roughly consistent with the power law indices found by Sergison et al. (2020) in the 'long' variables, which those authors also believe to represent accretion variability. Comparison between the Fourier power spectrum and structure function of correlated noise profiles implies that disc-bearing YSOs have accretion variability which roughly follows a random walk.

APPENDIX A: THE CASE OF MON-804
Fig . A1 shows the normalised magnitude distribution for all the objects in the sample, including Mon-804. There is a notable feature in both IRAC channels at a normalised magnitude of around -1.5 to -2.0. This has the effect of making the distribution appear significantly more skewed than if Mon-804 were excluded (see Fig. 4). We therefore choose to exclude Mon-804 in the main work to ensure the conclusions are more representative. Fig. A2 shows the individual normalised magnitude distributions for each burster. The red distribution corresponds to Mon-804 and highlights how it can affect the results for the entire sample. Fig. A3 shows the optical and infrared light curves of Mon-804.

APPENDIX B: COLOUR TIME SERIES
The R-[3.6] and R-[4.5] time series were produced by linearly interpolating the CoRoT magnitudes to the timestamps of the IRAC observations. The colour time series is produced by taking the difference between the IRAC and interpolated CoRoT measurements.
We followed Cody et al. (2014) and interpolated the denser grid onto the more sparsely measured data since the effect of interpolation will be negligible. In effect, this means we interpolate between the two CoRoT magnitudes either side of Spitzer observation we are considering. We demonstrate in Fig. B1 that the median offset in sampling between the IRAC observations and the corresponding CoRoT observations is around 2 × 10 −3 days.

APPENDIX C: STRUCTURE FUNCTION POWER LAW
For some particular Fourier power spectra, there exist analytic corresponding structure function power laws (see Findeisen 2015). To provide a more granular comparison between Fourier power spectra and structure functions, we simulated noise profiles using the algorithm described in Timmer & Koenig (1995). We then fitted a power law to both the Fourier power spectrum and to the structure function of each noise time series. Fig. C1 shows the empirical relationship between Fourier power index and structure function index. We found a best fit relationship (shown in red in Fig. C1) between the power law indices for the structure function and Fourier power spectrum ( SF and F , respectively) of SF = 0.99 − 1.06 × tanh [0.94 ( F + 1.98)] . (C1) We do however note that there is a spread in the relation, and there is not a one-to-one correspondence between the Fourier power and structure function. We can see this in Fig. 5, where the structure functions for the sinusoid and the Gaussian dips have structure functions with the same power law, despite having very different Fourier spectra. This paper has been typeset from a T E X/L A T E X file prepared by the author.     Figure B1. The distribution of time offset when finding the nearest points between bands. Blue is the offset between each CoRoT data point and its nearest IRAC counterpart, and orange is the offset between each IRAC data point and its nearest CoRoT counterpart. The blue and orange vertical lines mark the median offset for each scenario.  Figure C1. Relationship between the Fourier power spectrum power law index and the structure function power law index for 4000 instances of correlated noise. The best fit relation is shown in red. The solid, dashed, and dot-dashed vertical lines show the value of F for white noise, flicker noise, and a random walk, respectively.