A Zwicky Transient Facility Look at Optical Variability of Young Stellar Objects in the North America and Pelican Nebulae Complex

We utilize the normalized peak-to-peak variability metric

$$\begin{eqnarray}&&\nu =\displaystyle \frac{{({m}_{i}-{\sigma }_{i})}_{\max }-{({m}_{i}+{\sigma }_{i})}_{\min }}{{({m}_{i}-{\sigma }_{i})}_{\max }+{({m}_{i}+{\sigma }_{i})}_{\min }},\end{eqnarray} \tag{ 1 }$$

featured in Sokolovsky et al. (2017) to quantify variability amplitude for objects in our YSO sample. In this formula, m_i represents a magnitude measurement and σ_i represents the corresponding measurement error, with the minimum and maximum determined from the full light curve. We eliminated possible outliers with the application of a 5σ clip to all light curves, before calculating ν; this is necessary for the metric to be considered a sensitive variability indicator. For each r-band light curve in our sample that meets the measures described in Section 3, we calculate ν. We note that ν correlates well with standard deviation, with values that are about a factor of 6 smaller.

Our final sample for variability analysis includes objects with ν greater than the 15th percentile of ν as a function of magnitude. Although not a rigorously justified cutoff, this criterion ensures that we are studying the fractionally larger amplitude variables at each brightness level. Figure 3 depicts the selection of sources based on ν, and Table 1 shows the effect of the ν cut on our final sample size.

Zoom In Zoom Out Reset image size

We use the Astropy implementation of the Lomb–Scargle periodogram (VanderPlas et al. 2012; VanderPlas & Ivezić 2015) to search for periodic signals in the ∼825 day ZTF data stream. Possible periods between 0.5 and 250 days are considered.

We compute the periodogram for every object in the variable sample, flagging periods between 0.5–0.51, 0.98–1.02, 1.96–2.04, and 26–30 days for further analysis to avoid the most common aliases associated with the solar and sidereal days, as well as the lunar cycle (following Rodriguez et al. 2017; Ansdell et al. 2018). In particular, the window function shows that we have the potential for extremely strong signals at frequencies corresponding to 1.000 and 0.997 day periods. To account for additional potentially aliased periods, we flag periods with half or double multiples that fail at least one of the aforementioned alias checks. We then compute the 90%, 95%, and 99% confidence false alarm probability levels, and consider the period corresponding to the highest power peak greater than the 99% confidence false alarm probability (FAP; < 0.01) to be significant (following Messina et al. 2010). While FAP values may not be as reliable as formal hypothesis tests, they are sufficient for our purpose here, which is to identify potential variability timescales for visual examination and further assessment.

We repeat the above procedure with g-band data for objects lacking a significant period in r. If a significant period is found in g, we return this period for the source.

For each object, the full light curve was phased at the four most significant peaks and visually examined to assess the dispersion across phase. We also marked periods corresponding to beats of the ∼1 day sampling interval to assess the expected Fourier patterns. The periodograms are quite clean in terms of signal-to-noise. In the vast majority of cases, the highest periodogram peak was retained as correct, with recognizable harmonic aliasing and beats having lower power. In some cases, a longer period with a slightly less prominent peak was retained as the most likely true period based on better explanation of the beat patterns, and improved phase dispersion minimization. We call attention to two particularly interesting and nonintuitive examples in Appendix B.

Objects with significant Lomb–Scargle periodogram peaks in either r or g that passed our vetting procedures are reported in Table 2. We have labeled the relevant column as containing both timescales and true periods. In cases where the value of Q (quasiperiodicity metric; see next section) is high, the periodogram peaks are more appropriately thought of as timescales, whereas when Q is low, the light curve is truly periodic in the sense of strictly repeating, with low residual dispersion in the phased light curve. In this latter case, there is also a "P" in the column containing the primary variability classification. A few sources have multiple significant real (nonalias or nonbeat) peaks in the periodogram, and are given "MP" classifications.

The distribution of periods for the sources identified as strictly periodic is shown in the lower panel of Figure 4. There are 46 such stars, constituting ∼15% of the variable sample with reported timescales or periods in Table 2. Their period distribution has a bimodal peak, as well as a skew toward longer periods. Both features are consistent with the populations of fast and slow rotators discussed in previous young star literature, e.g., Gallet & Bouvier (2013).

Zoom In Zoom Out Reset image size

Figure 4 also shows, in the upper panel, the full range of variability timescales in the sample. This distribution exhibits a strong peak on the same few week timescales covered by the purely rotationally variable stars in the lower panel. The peak is significantly higher though, indicating a large quasiperiodic and aperiodic population with additional physical processes that are occurring on similar timescales. The timescale distribution is relatively flat out to longer times. The longer variability patterns, from weeks to months, can be seen in the light curves themselves. We emphasize that, even though these timescales were detected and measured using Lomb–Scargle periodogram techniques, they really are approximate timescales for photometric variability, and not true periods. Particularly in cases of high Q values (see below), these timescales are less trustworthy. Other methods for quantifying timescales in aperiodic light curves are discussed in Findeisen et al. (2015).

One particular source with a credible long period or timescale is V1490 Cyg, with a timescale of 31 days. In this case, the phased light curve shows substantially variable depth dips that are more consistent with material orbiting in a disk than with, e.g., a stellar rotation period.

We also performed a Lomb–Scargle period search on the g − r color curves, with colors created from the individual g and r time series as described below in Section 5.2. Few objects were found with significant periodicity in their colors. Among them, FHK 216 and FHK 442 are examples with periodic color variations. While their color curves do not phase particularly well, the color periods are the same as the single-band periods, reported in Table 2. FHK 216 has a color–magnitude slope similar to the reddening vector while FHK 442 is slightly shallower.

Our light-curve classification scheme follows the methodology of using statistics that quantify quasiperiodicity and flux asymmetry, first developed by Cody et al. (2014) and further refined by Cody & Hillenbrand (2018). Quasiperiodicity, Q, is defined as

$$\begin{eqnarray}&&Q=\displaystyle \frac{{\sigma }_{\mathrm{resid}}^{2}-{\sigma }_{\mathrm{phot}}^{2}}{{\sigma }_{m}^{2}-{\sigma }_{\mathrm{phot}}^{2}},\end{eqnarray} \tag{ 2 }$$

where σ_phot is the measurement uncertainty, ${\sigma }_{m}^{2}$ is the variance of the original light curve, and ${\sigma }_{\mathrm{resid}}^{2}$ is the variance of the residual light curve after the smoothed dominant periodic signal has been subtracted. We take σ_phot as the mean photometric error for all observations in an object's light curve, multiplied by a factor of 1.25 to account for an initial compression of Q values between 0 and ∼0.6 that we noticed early in our investigation. We attribute this compression to two effects: (1) the likely underestimate of reported photometric errors, and (2) the cadence dependence of the Q metric, which we discuss in Appendix C. For ${\sigma }_{\mathrm{resid}}^{2}$, we calculate the residual light curve by adopting a modified version of the periodicity routine described in Section 4.2. The light curve is first phased to the best period. Then, following Bredall et al. (2020), we mitigate for edge effects by concatenating three cycles of the period and convolving with a boxcar window of width 25% of the period, effectively smoothing the light curve by a factor of 4. Finally, we subtract the middle cycle of this model light curve from the folded light curve and compute σ_resid as the standard deviation of this residual curve. The resulting Q values are expected to fall between 0 (strong periodicity with low residual noise) and 1 (completely stochastic behavior).

Flux asymmetry, M, is defined as

$$\begin{eqnarray}&&M=\displaystyle \frac{\langle {m}_{10 \% }\rangle -{m}_{\mathrm{med}}}{{\sigma }_{m}}\end{eqnarray} \tag{ 3 }$$

where 〈m_10%〉 is the mean of all magnitude measurements in the top and bottom deciles of the light curve, m_med is the median magnitude measurement, and σ_m is the standard deviation of the light curve (Cody & Hillenbrand 2018). Objects that have M values <0 are predominately brightening, objects with M values >0 are predominantly dimming, and objects with M values near 0 have symmetric light curves. In practice, we follow previous literature and assign M < −0.25 as bursters, M > 0.25 as dippers, and the intermediate M values as relatively symmetric light curves.

We calculate Q and M based on the r-band light-curve data. Figure 5 provides illustrative light curves along the Q and M sequences.

Zoom In Zoom Out Reset image size

We then distinguish YSO light curves as populating nine categories, following Cody et al. (2014). We present these classification results in Table 2. In general, the objects are classified based on their Q and M values and then validated by visual examination, but there are cases where we adjust the classification to defer to the human eye for definitive categorization. This is especially prevalent in cases near the morphological boundaries in the Q − M plane.

Periodic variables (P) are defined as objects with significant periods from Section 4.2 that have −0.25 < M < 0.25 and Q < 0.45. We discuss our rationale for using a higher Q boundary for periodic objects than adopted in previous publications (e.g., Cody & Hillenbrand 2018; Bredall et al. 2020), in Section 7.2. The variability of the strictly periodic objects is most commonly interpreted as driven by rotational modulation due to the presence of star spots, and the periods derived are thus considered measures of stellar rotation (e.g., Rebull et al. 2008). Notable in this regard is that the sources with Q < 0.45 populate only the lowest range of the normalized variability amplitude, with ν < 0.02, while sources with higher Q span the full range of ν. Multiperiodic (MP) objects exhibit more than one distinct period, that is not a beat or alias. One common explanation for these is starspot modulation in binary pairs (e.g., Maxted & Hutcheon 2018; Stauffer et al. 2018).

Objects with Q > 0.45 are designated as quasiperiodic or stochastic. For the quasiperiodics, the variability origin is perhaps related to stellar rotation or the star-disk interaction region, but not dominated by the starspot modulation signal, which can become partially obfuscated in cases of a strong inner disk. There are two flavors of quasiperiodic variables, depending on M.

Quasiperiodic symmetric (QPS) variables are defined as having −0.25 < M < 0.25 and 0.45 < Q < 0.87. The variability of these objects is believed to have two possible sources. The first is a combination of purely periodic spot behavior with longer timescale aperiodic changes (e.g., variable accretion), and the second is a single variability process that is unstable from cycle to cycle (Cody et al. 2014).

Quasiperiodic dipper (QPD) variables are categorized based on M > 0.25 and 0.45 < Q < 0.87. The variability of these objects is believed to stem from variable extinction caused by time-variable occultation by circumstellar material (e.g., Alencar et al. 2010; Morales-Calderón et al. 2011; Ansdell et al. 2016). We do not define a quasiperiodic burster category, as distinct from the stochastic bursters, following previous literature in the field.

Burster (B) variables are defined as having M < −0.25, and correspond to objects with erratic but discrete accretion bursts that are generally short-lived (Cody et al. 2017). These bursts are distinct from larger amplitude and longer timescale outbursts from EX Lup and FU Ori type sources, which are more eruptive (Hartmann & Kenyon 1996; Herbig 2008).

Aperiodic dipper variables (APD) are defined as having M > 0.25 and Q > 0.87. These objects experience variable extinction, and it has been proposed by Turner et al. (2014) that their variability stems from inner disk scale-height changes induced by magnetic turbulence.

Stochastic (S) variables are defined as having −0.25 < M < 0.25 and Q > 0.87. They feature nonrepeating variability patterns with relatively symmetric excursions around a median brightness.

Long-timescale (L) objects exhibit variability on timescales ≳ 100 day timescales. They generally have Q > 0.87. Finally, we assign two objects to the unclassifiable (U) category because we were unable to calculate their Q due to ${\sigma }_{\mathrm{phot}}^{2}$ being larger than ${\sigma }_{\mathrm{resid}}^{2}$ and/or ${\sigma }_{{\rm{m}}}^{2}$.

The results of this classification scheme are reported in Table 2 for individual sources. Additionally, all stars in the variable sample are plotted in the Q − M plane in Figure 6, where they are differentiated in color by the assigned (primary) morphological class. The distribution of primary and any secondary variability classes is summarized in Table 3.

Zoom In Zoom Out Reset image size

Objects have secondary classifications in Table 2 when more than one light-curve type is evident in the data. Secondary classifications are typically long-timescale (L), representing objects that exhibit a strong primary behavior such as short timescale periodicity, quasiperiodic behavior, or bursting, along with a pronounced long-term trend such as brightening or dimming, or lasting changes in variability amplitude, etc. We call attention to two objects, V1716 Cyg (FHK 437) and FHK 496 that are designated as periodic in their secondary classification. Both have robust periods that phase very well, but their Q values are high, and the M metric reveals that the full amplitude of the light curve is caused by bursty behavior. As these cases are each somewhat interesting, we detail them here.

For V1716 Cyg, the behavior over most of the light curve is plainly periodic, but for a little over one year of our time series there is superposed continuous bursting at the 1 mag level. Findeisen et al. (2013) also observed bursts in their earlier light curve of this source based on PTF data. The source further distinguishes itself in having—aside from the bursts—a phased light-curve morphology that appears (see Appendix A figure set) similar to a W UMa profile. ⁵ with V-like troughs, and broad inverted-U-like peaks.

In FHK 496, bursting behavior is ubiquitous throughout the light curve, and causes scatter at the 0.4 mag level above the clear periodic signal, which has lower amplitude.

As mentioned above in Section 5.1, there are many different physical mechanisms that can cause photometric variability in young stars. Color time series data can be a helpful tool in distinguishing them. Among others, Carpenter et al. (2001), Günther et al. (2014) and Wolk et al. (2015) showed how near-infrared and mid-infrared multicolor photometry can distinguish extinction-related and accretion-related variability mechanisms. Optical photometry is even more sensitive to these behaviors, due to the steep rise in extinction from the infrared to the optical, and because of the hot nature of the accretion process.

The ZTF g-band and r-band observations are often taken close to one other, within a few hours up to a few days, but they are not simultaneous. Color–magnitude diagrams were thus constructed as follows. First, we partitioned the g-band light curves into segments in which the maximum separation between any two subsequent g observations is less than three days. We adopted this maximum observation separation constraint because linearly interpolating between g observations taken further apart created spurious values in the color–magnitude plane. In practice, within the defined <3.0 day intervals, approximately 20.5% of the g-band observations occur within 0.25 days of the prior g-band observation, 87.0% within 1.25 days, 93.4% within 2.0 days, and 97.5% occur within 2.25 days. For each g-band value that was paired with an r-band observation, we then linearly interpolated between the date bounds of the defined g-band intervals, in order to estimate g − r colors. To account for the propagation of errors throughout this process, we employed the Python package uncertainties (Lebigot 2010), which provided the properly calculated errors on the interpolated g-band magnitude as well as g − r color values.

The time series color–magnitude diagrams exhibit a variety of morphologies, with some sources showing large color and magnitude excursions while others have little color variation, and are essentially constant within the errors. Figures 7 and 8 show examples from among the dipper (M > 0.25) and burster (M < −0.25) categories, which have distinct slope angles. A selection of sources representing different Q values is presented in Appendix A, where a range of CMD slope angles can also be seen. Light curves and CMDs for the full sample are available in an associated online figure set.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For each member of the variable star set, we performed a linear fit to its color–magnitude diagram (CMD). To do so, we applied least-squares linear orthogonal fits using the python package scipy.odr (Boggs & Rogers 1990; Virtanen et al. 2020). This method was chosen following Poppenhaeger et al. (2015), because it can account for the significant and partially correlated errors in both axes. The regression assumes that points lie along a line in the g versus g − r plane, with Gaussian errors perpendicular to the line. While the condition is not strictly true, the assumption is reasonable given that there is no clear independent and dependent variable in this situation. To alleviate the effect of outliers, we performed the fits on the middle 95% of the CMD spans in g and g − r.

After computing the best-fit slopes in g versus g − r, we define slope angles as the inverse tangent, with the angles in degrees increasing clockwise from 0° (corresponding to color variability with no associated g-band variability) to 90° (colorless g-band variability). To compute errors on the best-fit angles, we adopt the pYSOVAR (Günther et al. 2014; Poppenhaeger et al. 2015) function fit_twocolor_odr and consider slope angles with errors less than 10° to be significant. These values are reported in Table 2. We note that a few sources even exceed 90° slopes: FHK 101 at 93°, 2MASS J21004676+4255265. at 95°, and FHK 473 at 973. Each of these sources is a large-amplitude variable, with a Pythagorean vector length in the range of 2–3, and high Q as well as M < 0.

In Figure 9, all significant CMD slope angles are plotted against their CMD vector lengths, following Poppenhaeger et al. (2015), and differentiated by their flux asymmetry (M) classes. There is a sizable span in the CMD vector lengths, with the longest spans limited to the largest slope angles (indicating the largest amplitude and color changes). A wide range of CMD slope angles is also represented, spanning over 45°. The outlier at small slope angle is FHK 280, which varies over ∼0.2 mag in g − r color with very little variation in g brightness; in the r band it is a low-level burster but we caution that the results could be spurious as the source is near the ZTF bright limit. The outlier with vector length ∼3 is FHK 332, which varies over 2 mag in g − r color and 3 mag in g brightness; it is classified as quasiperiodic symmetric. Finally, the extreme vector length ∼4 source is FHK 238, an aperiodic dipper. We also note FHK 267, the out-of-family burster with vector length ∼3.5; this source was also featured in Figure 5.

Zoom In Zoom Out Reset image size

Small slopes, <35°, are essentially unpopulated, implying a lack of colorless variability. While likely an astrophysical reality, this could also indicate a bias in our slope determination methods. For example, stars with little to no statistical correlation between variability in magnitude, and variability in color, generally have uncertain slope angles. Such sources are not included in the plot. However, the full set of slope and vector length determinations includes only 4 (1.2%) more sources than the set with <10° error that we are plotting.

The entirety of the populated parameter space in Figure 9 is represented by photometric variables that are relatively symmetric in flux (blue). The dipper-type variables (green) tend to dominate the population at the longest vector lengths and highest slope angles; their median slope angle is 738. The burster-type variables (red), conversely, are concentrated toward the shorter vector lengths and populate mainly the lower slope angles, almost exclusively <70°, with median slope angle 477. Beyond just median values, the differences between the morphological dipper and burster classes, and dipper and symmetric classes are very statistically significant based on pairwise K-S tests between their slope angle distributions, with p ≪ 10⁻¹¹. The burster and symmetric difference is also statistically significant, though only with p = 0.013.

We can assess the CMD slope angles in the context of standard reddening due to dust. We adopt the extinction curve of Fitzpatrick (1999), as translated by Schlafly & Finkbeiner (2011) into the PanSTARRS photometric system (to which the ZTF data are calibrated). The expected slopes in g versus g − r CMDs for different reddening laws are: 3.52 for R_V = 3.1, 4.39 for R_V = 4.1, and 5.24 for R_V = 5.1, with steeper slopes corresponding to larger dust grain sizes, as might be expected in molecular clouds and circumstellar disks. Considering a range of possibilities for the extinction law, we expect objects with time-variable CMD behavior dominated by extinction effects to have slope angles between 741 and 792 in the g versus g − r diagram.

Examination of Figure 9 shows, however, that the majority of slopes are much shallower than the expectations from purely reddening variations. We interpret the distribution of slope angles as indicating that much of the observed variability is at least partially related to accretion effects, with extinction effects (likely occurring outside of the accretion zone) perhaps playing some role. We acknowledge that other morphological light-curve types can have slope angles in the CMD consistent with reddening, that is, not all sources that look like they are experiencing variable extinction also exhibit identifiable dipping behavior.

In order to test the quality of our derived periods, we compared them to the periods published in Bhardwaj et al. (2019) and Froebrich et al. (2021). Both papers considered objects in the Pelican Nebula region (see Figure 1).

Of the 40 periodic variables identified in Froebrich et al. (2021), 32 are within 2'' of a source in our sample, of which 20 we classify as periodic in Q − M space. All 20 have periods within 5% of our derived periods. Additionally, among the 56 objects found to have strong periodic variations in Bhardwaj et al. (2019), 27 variables match within 2'' of a source in our sample and we classify four as periodic in Q − M space. In this subsample, we recover three periods within 5% of their Bhardwaj et al. (2019) values. We compared the phase-folded light curve of the remaining object with the discrepant period and found that forcing it to the period in Bhardwaj et al. (2019) introduces pronounced additional scatter. Among the additional 23 stars in common, 18 have periods in agreement within 5%, though we have labeled these periodogram peaks as timescales in Table 2 rather than periods, due to poor phasing and/or high Q values.

Overall, between these two samples our period recovery rate for strictly periodic objects is 95.8%. These period comparisons are presented in Figure 12.

Zoom In Zoom Out Reset image size

The quasiperiodicity metric Q has proven to be an extremely powerful tool for classifying the variable behavior of YSOs (Cody et al. 2014; Cody & Hillenbrand 2018). In this work, we have shown that the metric is sufficiently applicable to ZTF data, affirming its utility for ground-based as well as space-based photometric data sets. However, recent studies have revealed issues in translating Q results between various sources of data (Cody & Hillenbrand 2018; Bredall et al. 2020).

When comparing objects studied with different space-based telescopes, systematic effects present in the precision photometry from different observatories, or those caused by slight wavelength dependencies of the astrophysical phenomena, could produce systematically different Q values for sources with similar intrinsic behavior (Cody et al. 2014; Cody & Hillenbrand 2018). Notably, the morphological boundaries in Q were adjusted going from analysis of Convection, Rotation, and Transit (CoRoT) satellite data in Cody et al. (2014) to K2 data in Cody & Hillenbrand (2018). For the CoRoT data, Cody et al. (2014) set Q boundaries at 0.11 and 0.61, with M boundaries at ±0.25. For K2 data, in Cody & Hillenbrand (2018) the adopted Q boundaries were 0.15 and 0.85, with M remaining at ±0.25.

In contrast to this adjusting of the Q metric, Bredall et al. (2020) did not adjust the morphological boundaries of Q in their ground-based study based on ASAS-SN data, instead opting to simply apply the morphological boundaries defined for K2 data in Cody & Hillenbrand (2018). These authors did note the discrepancies that appeared in their classification scheme, such as objects having been classified in previous variability studies as QPDs, assessed in their numerical analysis as strictly periodic objects.

Indeed, prior to tailoring our Q routine to the ZTF data set, we noticed a similar compression and shift of our objects toward the periodic end (low Q) of the quasiperiodicity scale, with an absence of objects having Q > 0.65. This would be a very different distribution from the results of Cody et al. (2014) and Cody & Hillenbrand (2018), as in both of these studies, the range from Q = 0 to Q = 1 is densely populated (in fact, exceeding these bounds in Cody et al. 2014).

Given the substantially larger errors and the lower observational cadence associated with ground-based data, significant discrepancies in Q values might be expected when comparing objects studied with space-based versus ground-based telescopes. In Appendix C we investigate in detail the effects that cadence and photometric uncertainty have on Q.

For our analysis, we opted to follow the methodology of Cody & Hillenbrand (2018), tailoring both our morphological boundaries as well as our procedure for calculating Q to the nuances of our particular ZTF data set. As stated in Section 5.1, our adopted Q boundaries are 0.45 and 0.87. We selected these particular values after visually inspecting all variable object light curves, both unfolded and phase-folded, and comparing them to light curves from both Cody et al. (2014) and Cody & Hillenbrand (2018). We assigned the Q boundaries to occur at the most obvious behavioral transitions we noticed in the light curves when ordered by Q, and to minimize the number of edge-case objects that we end up reclassifying by eye. After visual inspection, approximately 11% of objects were reclassified.

We believe that such a flexible treatment of the morphological boundaries indicated by Q is the best way to ensure that variables classified in the Q − M parameter space are behaviorally representative of the original classes as defined by Cody et al. (2014). In other words, there is an unquantified error inherent in the calculation of the Q and M metrics, and it is more useful for furthering our understanding of the origins of young star variability, and for comparing sources observed with different instruments, if the quantitative procedures and morphological boundaries are tailored to the instruments used in each analysis. While objects of the same intrinsic morphological class could appear in different Q ranges in different studies, they should be properly grouped with their light-curve peers.

Regarding the full range spanned by our derived Q values (about 0.11–0.96), we believe there are reasonable explanations why the full zero-to-one parameter space is not fully covered. On the low side, none of our objects have Q values between 0.0 and 0.1 partially due to our choice to represent the estimated photometric uncertainty with the mean photometric error for each object, and our further inflation of the photometric errors during the Q calculation process. This imposes an artificial minimum to Q. There is also the likely reality that no objects in our sample have negligible phase dispersion (including photometric uncertainty) relative to the amplitude of a hypothetical periodic oscillation. On the high side of Q, the maximum at 0.96 instead of 1.00 may indicate simply that we have not inflated our errors enough to produce a truly stochastic Q metric. On the other hand, there is no rationale for fine tuning so as to have any particular number of objects with the theoretical maximum value of Q = 1.00.

Regarding the range spanned by the M values (about −1.2–1.1) for our aperiodic objects, there is again no expectation that any given variable star sample should span, or exist constrained within, the theoretical range from −1 to +1; although we observe our results seem to appropriately populate the parameter space in almost all cases. We note that, as the procedure for calculating the flux asymmetry M metric is more agnostic to the instrumental effects than the procedure for calculating Q, we were comfortable adopting the same flux asymmetry boundaries used in previous studies to designate dippers, symmetrics, and bursters. We note that the periodic objects at Q < 0.45 occupy a more restricted range in M, closer to the symmetry line of M = 0.

Summarizing, as these metrics are further applied, we recommend the continued utilization of visual inspection. We further suggest that Q values be scaled between approximately [0–1] in future works, in order to standardize the use of this metric across data sets.

In Table 3, we summarized the fraction of objects in the member sample we have assembled for the NAP nebulae region, in the different variability classes. We can compare our relative distribution to those resulting from previous uses of the Q − M parameter space to study young stars.

To do so, we refer to Table 3 in Cody et al. (2022), which includes the Taurus, Oph, and NGC 2264 regions studied with K2 and CoRoT data taken from space-based platforms. We focus our comparisons on the young regions with comparable ages to the NAP. The only other region with such information, Upper Sco, has variability properties that are distributed somewhat differently, perhaps because Upper Sco is somewhat older. We do not consider the results of Bredall et al. (2020) from ground-based data in Lupus, due to the apparent distortion in their Q distribution toward low values; this may be a sign of the issue we faced in this study (Section 5.1) where inflating our photometric errors by 25% recovered more of the expected range in Q for a young star variable sample. There is also the issue of whether ground-based and space-based data sets are comparable. As mentioned earlier, we assess the effects of data precision and data cadence on Q in Appendix C.

Relative to the previously studied star-forming regions, we find a similar fraction of periodic sources, including the periodic (P) and multiperiodic (MP) classes, at around 15% of the sample. This is interesting given that our sample is designed to include all members and not just disk-bearing sources that comprised the samples studied previously. The result can be understood if we accept that our sample is more biased toward including periodic sources, but at the same time we are less capable of detecting them as such, given the lower quality of our data set. Another commonality is that the quasiperiodic symmetric (QPS) category is our largest, and similar to the other star-forming regions at about 28%. Our stochastic fraction at 11% also matches well to the previous results. The quasiperiodic and aperiodic dipper (QPD and APD) groups are also similar, cumulatively around 29% of our sample. There is a range in these fractions among the other star-forming regions over a factor of 2, with our value in the middle. Our burster fraction is 14%, on the low side of the other regions (which range from about 13%–21%). An explanation here is that much of the burster activity identified from K2 is short-lived (e.g., Cody et al. 2017) and stands out relative to what could be detected from the ground.

Summarizing, our ground-based data set has matched the variability properties of space-based studies of comparably aged stars rather well, with any differences plausibly attributed to our lower sensitivity from the ground to narrow bursts and accretion flaring activity.