# Magnetic Activity-Rotation-Age-Mass Relations in Late Pre-main Sequence
Stars

Konstantin V. Getman Department of Astronomy & Astrophysics
Pennsylvania State University
525 Davey Laboratory
University Park, PA 16802, USA  Eric D. Feigelson Department of Astronomy &
Astrophysics
Pennsylvania State University
525 Davey Laboratory
University Park, PA 16802, USA  Gordon P. Garmire Huntingdon Institute for
X-ray Astronomy
LLC, 10677 Franks Road
Huntingdon, PA 16652, USA

(Accepted for publication in ApJ, May 15, 2023)

###### Abstract

We study the four-dimensional relationships between magnetic activity,
rotation, mass and age for solar-type stars in the age range $5-25$ Myr. This
is the late-pre-main sequence (l-PMS) evolutionary phase when rapid changes in
star’s interior may lead to the changes in magnetic dynamo mechanisms. We
carefully derive rotational periods and spot sizes for 471 members of several
l-PMS open clusters using photometric light curves from the Zwicky Transient
Facility. Magnetic activity was measured in our previous Chandra-based study,
and additional rotational data were obtained from other work. Several results
emerge. Mass-dependent evolution of rotation through the l-PMS phase agrees
with astrophysical models of stellar angular momentum changes, although the
data point to a subpopulation of stars with slower initial rotations than
commonly assumed. There is a hint of the onset of unsaturated tachoclinal
dependency of X-ray activity on rotation, as reported by Argiroffi et al.
(2016), but this result is not confidently confirmed. Both X-ray luminosity
and star spot area decrease approximately as $t^{-1}$ for solar mass stars
suggesting that spot magnetic fields are roughly constant and l-PMS stars
follow the universal solar-scaling law between the X-ray luminosity and
surface magnetic flux. Assuming convective dynamos are dominant, theoretical
magnetic fluxes fail to reveal the universal law for l-PMS stars that enter
late Henyey tracks. Altogether we emerge with a few lines of evidence
suggesting that the transition from the turbulent to solar-type dynamo occurs
at the later stages of l-PMS evolution as stars approach the Zero-Age Main
Sequence.

Pre-main sequence stars (1290) — Stellar rotation (1629) — Single x-ray stars
(1461) — X-ray stars (1823) —- Protoplanetary disks (1300) — Stellar magnetic
fields (1610)

††facilities: CXO, Gaia, ZTF††software: R (R Core Team, 2020)

##  1 Introduction

For main sequence solar-type stars, the relationships between stellar magnetic
activity, rotation and age has been studied for decades (Skumanich, 1972) and
their origins are well-understood (Schrijver & Zwaan, 2000). Stars have a
complicated rotational history during their early stages of evolution
(Bouvier, 2007) but, once settled on the zero-age main sequence (ZAMS), they
undergo a slow, sometimes punctuated spin-down due to the torque of a
magnetized wind. This leads to a slow weakening of their stellar activity:
X-ray flares, ultraviolet emission, chromospheric spectroscopic indicators,
and other manifestations of magnetic fields generated in a $\alpha\Omega$
tachoclinal dynamo. An example of this behavior is the relation between X-ray
emission and Rossby number showing saturated and unsaturated dynamo regimes
for $>40$ Myr old stars in main sequence open clusters (Wright et al., 2011a).

Skumanich-type relationships are not present during the pre-main sequence
(PMS) phase. No relationship between X-ray emission and rotation is seen
during the ‘early-PMS’ (henceforth e-PMS) phase when the stars are fully
convective descending the Hayashi tracks and rotation is regulated by coupling
to the accreting disk (Alexander & Preibisch, 2012; Henderson & Stassun, 2012;
Rebull et al., 2018). Instead, X-ray activity is linked to bulk properties
such as stellar mass, surface area, and volume (Preibisch et al., 2005; Getman
et al., 2022). The intense X-ray luminosities, frequent superflares, large
star spots, and strong surface fields point to a highly efficient distributed
$\alpha^{2}$ magnetic dynamo expected in fully convective stars (Chabrier &
Küker, 2006; Dobler et al., 2006; Käpylä et al., 2008; Christensen et al.,
2009; Browning, 2008; Brown et al., 2011; Yadav et al., 2015). In the solar-
type $\alpha\Omega$ dynamo, a relatively weak magnetic field is generated in
the thin tachocline at the interface between radiative and convective zones.
Higher rotation rates lead to higher surface magnetic fields. Calculations of
$\alpha^{2}$ dynamos in fully convective stars indicate strong fields are
efficiently generated throughout the stellar interior reaching equipartition
strengths. The energy of convective flows rather than rotation sets the
strength of surface magnetic field.

The transition between the e-PMS and ZAMS phases probably occurs during the
‘late-PMS’ (l-PMS) phase that we define to start when the star lies at the
bottom of the Hayashi track on the Hertsprung-Russell diagram and to end when
it begins hydrogen burning on the ZAMS111 Here we assume ’classical’ PMS
evolutionary tracks as initially calculated by Iben (1965) and more recently
by many other researchers (Squicciarini & Bonavita, 2022). There are
substantial differences among the tracks and the evolution may further differ
when the effects of accretion, internal magnetic fields, and surface star
spots are considered (Baraffe & Chabrier, 2010; Feiden & Chaboyer, 2013;
Somers & Pinsonneault, 2015; Kunitomo et al., 2017). . Little is known
observationally about l-PMS surface magnetic activity or changes in internal
magnetic dynamo processes. For a solar mass star, the surface temperature
becomes hotter while the luminosity increases and radius slowly decreases so
that it moves quickly across the Hertzprung-Russell diagram. Internally, the
radiative zone grows, the convective zone size and convective turnover
timescale falls, and brief nuclear burning of 12C may occur (Iben, 1965;
Landin et al., 2022). $ROSAT$ X-ray data show that stars in the IC 2391, IC
2602 and NGC 2547 open clusters exhibit the same X-ray-rotation relation as
seen on the ZAMS (Patten & Simon, 1996; Stauffer et al., 1997; Irwin et al.,
2008b; Wright et al., 2011a). This suggests a transition in dynamo processes
is largely complete at ages $\sim 40-50$ Myr.

Observational evidence for the transition from the e-PMS convective dynamo to
the main sequence tachoclinal dynamo is sparse. The most tantalizing study is
based on the $Chandra$ observatory detection of $\sim 400$ members of the h
Per cluster with age $\sim 13$ Myr (Argiroffi et al., 2016, henceforth A16).
While most stars lie in the (super)saturated regimes of the X-ray-Rossby
number diagram similar to e-PMS stars, several slowly rotating stars with
periods $P\gtrsim 3$ days exhibit X-ray emission declining with Rossby number
suggesting that a tachoclinal dynamo has emerged. This result is best seen in
their $1.2$ M⊙ $<M<1.4$ M⊙ mass stratum.

In the present study and its earlier companion (Getman et al., 2022,
henceforth Paper I), we seek to improve the observational evidence of l-PMS
changes in interior structure and dynamos on activity-rotation-age relations.
Paper I reports X-ray observations of 10 open clusters in the age range $5-25$
Myr at distances 0.3 to 2.5 kpc, including h Per (= NGC 869) studied by A16.
With $Gaia$ observations to assist with selecting cluster members and
$Chandra$ observations to measure magnetic activity in the X-ray band, Paper I
identifies $6,003$ stars with masses above $0.4-1.0$ M⊙ completeness limits.
Most of these stars lie in the l-PMS evolutionary phase. For comparison with
fully convective e-PMS stars, these samples are combined with 40,041 stars
with ages $<$5 Myr from previous $Chandra$ surveys (Feigelson et al., 2013;
Getman et al., 2019).

Paper I found that X-ray luminosities are high and constant during the first
few Myr after which they decay. The activity of higher mass stars that evolve
more rapidly through the l-PMS phase falls more rapidly than in lower mass
stars. This decay in X-ray activity is interpreted as decreasing efficiency of
the $\alpha^{2}$ dynamo as radiative cores grow. But the emergence of an
$\alpha\Omega$ dynamo was not clearly seen.

Several observational challenges must be tackled for a successful
understanding of magnetic activity and dynamos in the l-PMS phase:

  1. 1.

The use of $L_{X}/L_{bol}$ as a measure of X-ray activity that is effective
for revealing dynamo activity for ZAMS and older stars cannot be effectiely
applied to l-PMS stars: the changes in X-ray emission are often overwhelmed by
the radical changes in bolometric luminosity and stellar radius as stars move
across the HR diagram (Paper I). We use $L_{X}$ alone to quantify surface
magnetic activity in the l-PMS era.

  2. 2.

Changes in magnetic activity during the l-PMS should be studied with a range
of masses as higher mass stars evolve from the e-PMS to ZAMS regimes much more
rapidly than lower mass stars. A range of cluster ages must be examined to
capture the l-PMS phase in stars of different masses. While Skumanich-type
relations are three dimensional on the main sequence, they involve four
dimensions during PMS evolution: activity, rotation, age and mass.

  3. 3.

The X-ray brightnesses of an ensemble of stars in a single age and mass
stratum exhibit considerable scatter due, at least in part, to frequent
powerful superflares during the PMS phase (Getman & Feigelson, 2021). For
stars not caught during flares, it is important to incorporate X-ray
nondetections to give unbiased evaluation of magnetic activity using
‘survival’ statistical methods such as the Kaplan-Meier estimator (Feigelson &
Nelson, 1985).

  4. 4.

Unlike e-PMS stars that are concentrated in dense clusters associated with
their natal molecular environment, l-PMS stars are often dispersed. The
closest sample at distance $\sim 120$ pc, the Sco-Cen Association with
outliers like the TW Hya Association and $\beta$ Pic Moving Group, is spread
over large regions of the celestial sphere. For this reason, the l-PMS
clusters in Paper I were selected to be compact. Most of these clusters are at
distances around $1-2.5$ kpc. The greater distance, however, reduces the
fraction of X-ray detections with realistic $Chandra$ exposure times.

  5. 5.

The measurement of stellar rotation from photometric light curves can be
difficult, particularly if a single telescope is used so the star is
observable only for a few hours each day. Alias structure in the periodogram
can lead to ambiguities in the inferred rotation period (VanderPlas, 2018;
Rebull et al., 2022). Careful interpretation of periodograms is necessary.

We tackle these challenges in Paper I and the current study. Items 2$-$5
require that the magnetic activity of large stellar samples with a range of
masses and ages be examined as in Paper I. An additional $\sim 2000$ cluster
members are found in the outskirts of these clusters with reduced Chandra
sensitivity (see Table 4 in Paper I). The X-ray nondetections are reported
based on cluster membership obtained from $Gaia$ satellite astrometry and
photometry. Most of these stars are in the l-PMS phase and, when combined with
$\sim 40,000$ stars in the e-PMS phase from earlier studies (Feigelson et al.,
2013; Getman et al., 2017), large portions of the activity-age-mass space is
covered.

The fifth challenge involving periodogram analysis based on photometric light
curves of the Zwicky Transient Facility (ZTF) is treated in §3 below. We
develop conservative and reliable extraction of rotational periodicities from
Lomb-Scargle periodograms. This results in rotational periods for 566 $[5-25]$
Myr stars with measured X-ray emission. Although it is an increase over the
A16 study of h Per alone, it represents only a small fraction of the full
multi-cluster sample of Paper I.

Our analysis starts with a description of the stellar samples and estimates of
individual stellar masses and X-ray luminosities and other stellar properties
(§2 and Appendix). Section §3 describes methodology for obtaining reliable
rotational periods from ZTF lightcurves. The resulting rotation-mass-age
relations are presented in §4.1 and mass-stratified X-ray activity-rotation
diagrams are presented in §4.2. Section 5 gives some findings relating to the
sizes of star spots, a property indirectly emerging from our search for
rotational periods. For our samples of e-PMS and l-PMS stars from Paper I and
the current study, §6 presents the connection between the inferred X-ray-
rotation quantities and the data- and theory-based stellar surface magnetic
flux estimates. The results are summarized in §7.

##  2 Star Samples and Properties

###  2.1 Late-PMS Stellar Samples With X-ray Measurements

We investigate activity-rotation-age-mass relations for the Gaia-Chandra
stellar members of the six $7-25$ Myr old clusters studied in Paper I that are
accessible to the Palomar 60 inch telescope used for the ZTF survey. These
include the 7 Myr old NGC 1502 (256 members from Table 4 in Paper I), 7 Myr
old NGC 2169 (104 members), 13 Myr old NGC 869 (1370 members) and NGC 884
(1294 members) commonly known as h and $\chi$ Persei, 22 Myr old NGC 1960 (406
members), and 25 Myr old NGC 2232 (117 members)222Cluster ages are obtained in
Paper I by fitting Gaia EDR3 color-magnitude diagrams with isochrones from
PARSEC 1.2S evolutionary models (Bressan et al., 2012). Ages for larger fields
of view are given in §2.2 here.. Out of 3547 members of these clusters, 3410
have ZTF light curves with at least 100 photometric observations that we
examine in this study.

The Chandra X-ray and Gaia-EDR3 data giving cluster memberships are discussed
in Paper I. Inferred cluster properties, including distance from the Sun, age,
visual extinction, numbers of known Gaia-Chandra stellar members, and mass
completeness limits for the Gaia-Chandra stellar samples are listed in Paper
I’s Table 1. Field star contamination depends on cluster distance and field of
view (Paper I, Figure Set 4): clusters with small fields and distances
$\lesssim$1 kpc have less contamination (NGC 1502, NGC 2169, NGC 1960) while
clusters with distances $\simeq$2 kpc suffer higher contamination (NGC 869,
NGC 884, NGC 3293, NGC 4755, NGC 3766). Nearby clusters with larger fields of
view have high contamination as well (IC 2395 and NGC 2232). Gaia color-
magnitude diagrams are shown in their Figure Set 4. The X-ray luminosity
($L_{X}$) measurements and upper limits, as well as stellar properties such as
Gaia-based effective temperature, mass, and bolometric luminosity, are
presented in their Table 4.

Table 1 here gives stellar properties and X-ray luminosities for a subset of
471 stars from Paper I for which we obtain reliable rotational periods from
ZTF light curves (§3). We then add 95 stellar members of the 5 Myr old NGC
2362 open cluster, the oldest region in our earlier MYStIX (Massive Young
Star-Forming Complex Study in Infrared and X-Ray) survey (Feigelson et al.,
2013), for which stellar rotational periods are obtained by Irwin et al.
(2008a). Most X-ray measurements for NGC 2362 members are listed in Table 6 of
Paper I; additional information appears in Appendix A below.

The last column of Table 1 presents estimated convective turnover timescales
obtained by matching stellar ages and masses to the mass-age-turnover
timescale dataset from the stellar evolutionary model of Amard et al. (2019).
In cases of unstable mass regimes when the modeled $\tau_{conv}$ abruptly
drops from high values of $200-300$ days to small values $<10$ days within a
short mass span, $\tau_{conv}$ estimates are omitted.

A straightforward quantitative assessment of the fraction of unresolved
binaries by Gaia among the 566 X-ray stars listed in Table 1 can be conducted
by examining the quantity of Gaia-DR3 re-normalized unit weight error (RUWE).
Out of the 566 stars, only 16 ($<3$%) have a RUWE that exceeds 1.4, which is
considered an indicator of poor Gaia astrometric solutions.

Table 1: $5-25$ Myr PMS Stars With Rotation And X-ray Properties

Src | Cluster | R.A. | Decl. | $\log(T_{eff})$ | $M$ | $\log(L_{bol})$ | $P_{rot}$ | Ampl | $\log(L_{X})$ | $\log(L_{X,up})$ | $\tau_{conv}$
---|---|---|---|---|---|---|---|---|---|---|---
|  | deg | deg | K | $M_{\odot}$ | $L_{\odot}$ | day | mag | erg s-1 | erg s-1 | day
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9) | (10) | (11) | (12)
1 | NGC1502 | 61.971574 | 62.278508 | 3.64 | 0.92 | -0.30 | 11.55 | 0.050 |  $\cdots$ | 28.96 | 268.3
2 | NGC1502 | 61.910148 | 62.277657 | 3.62 | 0.76 | -0.49 | 0.91 | 0.051 |  $\cdots$ | 28.96 | 363.1
3 | NGC1502 | 61.926728 | 62.368586 | 3.55 | 0.56 | -0.76 | 0.30 | 0.131 |  $\cdots$ | 28.97 |  $\cdots$
4 | NGC1502 | 61.920049 | 62.261276 | 4.15 | 3.55 | 2.22 | 3.31 | 0.060 |  $\cdots$ | 28.97 |  $\cdots$
5 | NGC1502 | 61.823610 | 62.298647 | 3.62 | 0.77 | -0.49 | 0.79 | 0.045 |  $\cdots$ | 28.98 | 357.9
6 | NGC1502 | 62.118786 | 62.365436 | 3.62 | 0.78 | -0.48 | 3.70 | 0.061 |  $\cdots$ | 28.99 | 351.9
7 | NGC1502 | 61.862238 | 62.242426 | 3.57 | 0.70 | -0.62 | 5.74 | 0.163 |  $\cdots$ | 29.17 |  $\cdots$
8 | NGC1502 | 62.137218 | 62.374642 | 3.65 | 1.03 | -0.24 | 3.48 | 0.046 |  $\cdots$ | 29.03 | 238.0
9 | NGC1502 | 62.237633 | 62.260949 | 3.63 | 0.80 | -0.40 | 11.00 | 0.056 |  $\cdots$ | 29.38 | 334.1
10 | NGC1502 | 62.228920 | 62.387899 | 3.57 | 0.69 | -0.62 | 10.48 | 0.140 |  $\cdots$ | 29.41 |  $\cdots$

Note. — This table is available in its entirety in the machine-readable form
in the online journal. Out of the 566 table records, the first 471 indicate
stars with ZTF-based rotation periods, which are members of NGC 1502, NGC
2169, NGC 869, NGC 884, NGC 1960, and NGC 2232 (Getman et al., 2022). The
remaining 95 records present stellar members of the 5 Myr old NGC 2362
cluster, with rotation periods taken from Irwin et al. (2008a). Column
descriptions:

Column 1: Sequential source number from 1 to 566. Column 2: Cluster name.
Columns 3-4: Gaia DR3 right ascension and declination (in decimal degrees) for
epoch J2000.0. Columns 5-7: Stellar effective temperature, mass, and
bolometric luminosity derived from the Gaia color-magnitude diagrams. For the
stars $1-471$, these quantities are from Table 4 of Getman et al. (2022). For
the NGC 2362 members, these quantities are re-derived in the current paper.
Column 8: Rotation period. For the stars $1-471$, these are derived here from
ZTF lightcurves; values for NGC 2362 are taken from Irwin et al. (2008a).
Column 9: Peak-to-peak photometric amplitude of the spot rotational modulation
from folded lightcurves. Columns 10-11: $Chandra$ X-ray luminosities or upper-
limits (in the $0.5-8$ keV band) from Tables 4 and 6 of Getman et al. (2022).
For faint X-ray stellar members of NGC 2362, additional Chandra X-ray
luminosities are derived in Appendix A. Column 12: Convective turnover
timescale integrated over the convective envelope based on the stellar mass,
age and rotation using models of Amard et al. (2019).

###  2.2 Additional Late-PMS Stars Without X-ray Measurements

For the purpose of comparison with modern rotation evolution models (§4.1) and
analyses of stellar spot areas (§5), we extend our X-ray-rotation data sample
(§2.1, Table 1) by adding other stars in the $5-25$ Myr age range (mostly
without $Chandra$ data) with previously reported rotation periods. These
additional stars include members of the NGC 2362 and NGC 869 (including many
members outside $Chandra$ fields), Upper Sco (US), Upper Centaurus-Lupus (UCL)
and Lower Centaurus-Crux (LCC) subregions of the Sco-Cen Association, and the
Pleiades. Their membership and rotation period information are taken from
Irwin et al. (2008a), Moraux et al. (2013), Rebull et al. (2018), Rebull et
al. (2022), and Rebull et al. (2016). For US, UCL, and LCC, we only consider
the ‘gold’ members with a single rotation period from Rebull et al. (2018,
2022).

In order to place all these stars on the same age-mass scale as the stars from
our X-ray-rotation sample, we apply PARSEC 1.2S evolutionary models (Bressan
et al., 2012; Chen et al., 2014) to Gaia-DR3 data (Gaia Collaboration et al.,
2022) using the same procedures as with ZTF clusters (Paper I) and NGC 2362
MYStIX stars (Appendix A). Specifically, we match these stars to the Gaia-DR3
catalog using a radius of $1\arcsec$, estimate median cluster distances from
the Sun using only stars with $RUWE<1.4$ and $\sigma_{\omega}<0.1$ mag,
truncate the samples to include only stars with the Gaia color uncertainty of
$\sigma_{BP}<0.1$ mag, create Gaia-DR3 color-magnitude diagrams, fit the
color-magnitude distributions with the reddened PARSEC 1.2S models to derive
cluster extinction/age and stellar masses/temperatures/luminosities.

The Gaia-inferred cluster distance, extinction, and age are as follows: NGC
2362 ($d=1301$ pc, $A_{V}=0.3$ mag, $t=5$ Myr as in Appendix A), US ($d=142$
pc, $A_{V}=0.7$ mag, $t=9$ Myr), NGC 869 ($d=2451$ pc, $A_{V}=1.7$ mag, $t=13$
Myr as in Paper I), UCL ($d=142$ pc, $A_{V}=0.3$ mag, $t=15$ Myr), LCC
($d=114$ pc, $A_{V}=0.1$ mag, $t=15$ Myr), and the Pleiades ($d=$136 pc,
$A_{V}=0.2$ mag, $t=100$ Myr). Appendix B tabulates these additional stars and
shows Gaia color-magnitude diagrams.

###  2.3 X-ray Luminosities

For many thousands of stellar members of very young (generally $t<5$ Myr) star
forming regions, such as the MYStIX regions mentioned in §2.1, Table 6 of
Paper I lists their Chandra-based intrinsic X-ray luminosities in the $0.5-8$
keV band. These quantities were derived using the non-parametric method XPHOT
(Getman et al., 2010) based on a concept similar to the use of color–magnitude
diagrams in optical and infrared astronomy. The XPHOT method requires the
presence of both soft ($<2$ keV) and hard ($>2$ keV) X-ray photons.

But for the six older and less absorbed clusters with ZTF stellar rotation
periods derived here, hard X-ray photons are often lacking. Paper I estimated
stellar intrinsic X-ray luminosities by stacking all individual stellar
Chandra spectra (typically of low S/N) into one high-S/N spectrum, obtaining
intrinsic X-ray luminosity ($L_{X}$) for the merged spectrum. This merged
$L_{X}$ is divided by the sum of apparent X-ray fluxes of the involved stars
to estimate a conversion factor from intrinsic to apparent fluxes which is
then applied to estimate intrinsic X-ray luminosities of individual X-ray
detected cluster members. In addition, for the Gaia cluster members with no
X-ray detections, Paper I estimates upper limits on $L_{X}$ based on spatial
positions in the Chandra ACIS-I detector. These X-ray luminosities and upper
limits for the members of the six ZTF clusters are listed in Table 1.
Treatment of NGC 2362 members is further discussed in Appendix A.

With the inclusion of the Gaia cluster members with X-ray non-detections, the
full Chandra-Gaia stellar samples of the six ZTF clusters listed in Table 4 of
Getman et al. (2022) are complete down to masses of $0.4-0.7$ $M_{\odot}$ for
NGC 1502, NGC 2169, NGC 1960, and NGC 2232, and down to $1$ $M_{\odot}$ for
the more distant NGC 869 and NGC 884 clusters.

The X-ray-rotation stellar samples given in our Table 1 can not be considered
complete for any mass stratum due to the absence of reliable ZTF rotational
periods for most $Chandra-Gaia$ stars. We thus stress that any X-ray-rotation
correlations or slopes observed/measured in the following sections may be
biased with respect to any intrinsic relations.

##  3 Photometric Rotational Periods from ZTF Light Curves

The Zwicky Transient Facility (ZTF) is an automated wide-field survey of $\sim
30,000$ square degrees of the northern sky obtained with the 48-inch Oschin
Schmidt telescope at Palomar Observatory. The ZTF survey is described by Masci
et al. (2019) and Bellm et al. (2019) with access provided by the NASA/IPAC
Infrared Science Archive. Our analysis is based on Data Release 10 (DR10).

After matching the 3547 Chandra-Gaia l-PMS cluster members discussed in §2.1
to the ZTF point source catalog using a matching radius of 2″, we select 3410
ZTF light curves with at least 100 photometric observations in the $r$ or $g$
bands, after flagging for bad photometry. Unreliable data points are further
removed with two rounds of ‘sigma clipping’ using a $4\times IQR$
(InterQuartile Range) criterion. For some analysis steps, the irregular
cadence is converted to a regular cadence with 1.0 day bin width. Note that we
did not attempt to use the more sophisticated custom ZTF data reduction method
by Lu et al. (2022), which involves considering archival ZTF light curves for
multiple stars located near the target star to remove the systematics that are
common to all these stars.

Before proceeding to results, we discuss the difficulties of measuring
rotation periods from irregularly spaced light curves obtained from a single
telescope site (§3.1). We then present our methodology in some detail, as it
differs from other treatments in the field (§3.2).

###  3.1 Obtaining Reliable Photometric Rotation Periods in PMS Stars

Fully convective e-PMS stars often exhibit $\sim 100$ mmag rotational
modulation of large cool star spots, allowing photometric determination of
their rotational periods (e.g. Bouvier et al., 1993; Herbst et al., 2000;
Affer et al., 2013). The amplitude of periodic variability can range over a
factor of 10. The effect declines during the l-PMS phase; in Upper Sco, many
stars have spot amplitudes around $30-80$ mmag (Rebull et al., 2020). Reported
rotation periods mostly fall into the range $\sim 1-10$ day with a few
extending to $15-20$ days during both e-PMS and l-PMS phases. Photometric
measurements of rotation throughout the early phases of stellar evolution,
together with observable bulk stellar properties (mass, radius, age), give
insights into surface magnetic activity and internal dynamo processes.

However, photometric rotation periods of PMS stars are not easily measured.
First, the datasets have significant deficiencies. NASA’s 4-year Kepler
observations did not include the Galactic Plane. The K2 mission covered parts
of the Galactic Plane with $\sim 3$ month light curves but without associated
sensitive X-ray measurements. TESS Full-Frame Image light curves are subject
to crowding and nebulosity effects and have not yet been analyzed for PMS
stars. Periods will be limited to $\lesssim 10$ days due for single-sector
observations. Ground-based surveys such as HATSouth and ZTF have hundreds of
observations spanning several years, but their photometric accuracy is subject
to atmospheric and instrumental limitations.

Second, star spot instability can confuse periodic rotational effects for
multi-year light curves. e-PMS star spots may remain large and stationary on
the photosphere for years but may also suddenly change (Vrba et al., 1988;
Grankin & Artemenko, 2009). Spot durations have been estimated for a few l-PMS
stars with durations in the range $\sim 100-400$ days (Maldonado et al.,
2022). But in the Sun, spots often fade in weeks and are replaced by spots at
different longitudes. The resulting changes in phase, even with a constant
period, broadens and weakens the peak in a periodogram.

Third, many stellar light curves from sensitive space-based photometry reveal
autocorrelated stellar variations from magnetic activity, stellar pulsations
or super-granulation motions (Caceres et al., 2019). Autocorrelated variations
from atmospheric conditions must be removed from ground-based photometry using
algorithms like SysRem and the Trend Filtering Algorithm (Tamuz et al., 2005;
Kovács et al., 2005) prior to a search for stellar rotational periodicity.

Fourth, the statistical methodology for obtaining rotational periods from
irregular light curves has considerable limitations. The most common tool is
the Lomb-Scargle periodogram (LSP), the Fourier power spectrum modified for
irregularly spaced cadences (Scargle, 1982). Even in the best of
circumstances, the evaluation of significance of LSP peaks $-$ commonly called
False Alarm Probabilities $-$ is difficult to evaluate (e.g. Koen, 1990;
Baluev, 2008; Süveges, 2014; Delisle et al., 2020).

But the situation is considerably worse for ground-based surveys from a single
observatory where the star is visible for only several hours each day and
several months each year due to solar motion. When the light curves are folded
with periods around 1.0 days or its aliases, sinusoidal fits have unstable
amplitudes leading to false and fluctuating power in the periodogram (Dawson &
Fabrycky, 2010; VanderPlas, 2018). Figure 1 shows an example of the Lomb-
Scargle periodogram for a spotted l-PMS star with rotation period of 11.5 days
observed with ZTF. But strong peaks are also seen at 1.00, 0.50 and 0.25 days
from the diurnal solar cycle. Peaks may also appear around 28 days from the
lunar cycle. These structures in the periodogram can be present even if no
true stellar periodicity is present. Most of these peaks are spurious and
attempts to apply straightforward False Alarm Probabilities to these irregular
cadences will often fail.

As an alternative to the LSP, the nonparametric autocorrelation function (ACF)
is sometimes used for detecting rotational periods. McQuillan et al. (2014)
argue that the first peak of the ACF gives a reliable estimate of the
rotational period. But there is no statistical guidance on how to construct
the ACF: different first peaks will appear with different bin widths or
smoothers used to display the ACF. First peaks may also arise from noise or
aperiodic autoregressive variations rather than meaningful rotational signals.

In the present study, we use both the LSP and ACF methods. We test possible
periodicities with simulated light curves reproducing the peculiar
observational cadence of each light curve. The existence of a possible
rotational signal is then evaluated by conservative subjective judgment:

  1. 1.

An LSP peak stronger than spurious solar and lunar peaks must be present.

  2. 2.

A strong periodic signal must be present in the ACF if the period is longer
than $\sim 3$ days.

  3. 3.

The LSP and ACF periods must agree, understanding that the LSP period is more
accurate than the ACF period.

Figure 1: Results from the rotational period search for 3410 ZTF light curves.
This is the first page from a figure set for 471 stars with reliable ZTF
rotation periods available in the electronic version of the paper. Panel (a)
ZTF lightcurve. The legend gives the number of ZTF observation nights and
magnitude range of the lightcurve. Panel (b) Autocorrelation function after
collecting the light curve into 1 day bins. The log of the Ljung-Box
probability is shown; cases with $\log(P_{LB})<-2$ require autocorrelation
while cases with $\log(P_{LB})>-2$ are consistent with white noise. Panel (c)
Lomb-Scargle periodogram for the unbinned ZTF data. The red plus marks the
chosen best period, and dashed lines show expected diurnal aliases
(VanderPlas, 2018). Panel (d) Folded lightcurve for the best LSP period. The
green curve gives a robust local polynomial fit to the data using R’s loess
algorithm (see text for details). The legend gives the period and amplitude of
a smooth local regression fit. Panel (e) Lomb-Scargle periodogram for a
simulation of the best rotation period and amplitude, plus Gaussian noise,
using the original ZTF observation times. Figure 2: Lomb-Scargle periodograms
for source #260 in NGC 884. Panels a and c show the same LSP for the observed
ZTF lightcurve but assuming two different “true” periods, 1.60 and 0.61 days
(marked by the red pluses). Corresponding expected failure modes are marked by
the blue dashed lines. Panels b and d show the results of simulations for the
trials of the two assumed periods. The simulation trial with the 1.60 day
period more closely resembles the real periodogram. Figure 3: Comparison
between our ZTF periods and previously published rotation periods by Moraux et
al. (2013) for 84 stars in NGC 869 (= h Per). The grey lines mark the unity,
1/2$\times$ and 2$\times$ unity lines.

###  3.2 Periodogram Analysis

The Lomb-Scargle periodogram is calculated using CRAN package $lomb$ in the R
public domain statistical software environment (R Core Team, 2021; Ruf, 1999)
based on the algorithm of VanderPlas & Ivezić (2015). Our LSP examines periods
ranging from 0.2 to 100 days. No attempt is made to estimate statistical
significance of peaks given the difficulties outlined in §3.1 and illustrated
in Figures 1 and 2.

The autocorrelation function is calculated using R’s $acf$ function (Venables
& Ripley, 2002) applied to the 1-day binned light curve. ACFs constructed with
smaller bins were unsuccessful due to the proliferation of missing pairs of
measurements at short lags. As with the Discrete Correlation Function of
Edelson & Krolik (1988), missing contributions to the correlation function due
to cadence gaps are ignored. Confidence intervals corresponding to 95%
probabilities assuming uncorrelated Gaussian noise and no missing values are
shown as dashed lines in Figure 1b. While they underestimate the true
confidence intervals given the erratic number of missing values in each lag,
they are useful to give a rough indication of ACF features with low
statistical significance. More usefully, the probability that no
autocorrelation is present is calculated with the Ljung-Box omnibus hypothesis
test for autocorrelation at any lag time (Ljung & Box, 1978). The Ljung-Box
test generalizes the classical Durbin-Watson test for serial (lag=1 time unit)
autocorrelation. We calculate $P_{LB}$ using R’s $Box.test$ function with the
Ljung-Box test option.

Based on visual inspection of resulting LSPs and ACFs, we select an initial
sub-sample of 657 light curves from the original sample of 3410 ZTF stars with
very low values of $\log(P_{LB})$ and high signal-to-noise LSP periodograms.
These light curves show one or more periodic signals. We find that the ACFs
are insensitive to short rotation periods but are often consistent with one of
the dominant LSP peaks in cases of long rotation periods.

Due to the interaction between observation’s window function and spectral
power, the resulting LSPs exhibit spurious harmonics and aliases that
researchers may mistakenly interpret as true rotation periods (Dawson &
Fabrycky, 2010; VanderPlas, 2018). For the vast majority of our stars, the
ZTF-based LSPs show the presence of strong peaks within both low-period
($P\lesssim 1.5$ day) and high-period ($P\gtrsim 1.5$ day) ranges. Based on
simulations, we often find that one of the peaks, often but not always with
the highest LSP power, is the true period while the others are spurious.

In our simulations, using R’s $sample$ and $sin$ functions, we first scramble
the ZTF photometric measurements to create randomized datasets with erased
periodicity or other autocorrelated structure. Artificial sinusoidal signals
are then injected with LSP trial periods with amplitudes obtained from the
loess fits to the folded ligthcurves333The well-established ‘loess’ (LOcally
Estimated Scatterplot Smoothing) regression model, also known as the Savitsky-
Golay filter, fits quadratic functions in local windows that pass sequentially
with increasing rotational phase (Cleveland & Devlin, 1988). The model
includes robust downweighting of outliers using M estimation with Tukey’s
biweight function. The method is implemented in R function loess.. For each
choice, we mark expected harmonics for $m=1$ and $n=\pm(1,2,3,4)$ according to
equations (45)-(47) from VanderPlas (2018). For all 657 stars with high S/N
periodograms, we ran two trials of simulations. In the first (second) trial,
we assume that the strongest spectral peak at short periods $P<1.5$ day (long
periods $P>1.5$ day) is the true period. The simulation trial where the alias
structure closely matches that of the observed ZTF LSP is considered to be the
true period.

Figure 1 shows an example with a clear long period seen in both the ACF and
LSP. It is a sample page of the electronic Figure Set that presents the ZTF
lightcurve, ACF, LSP, folded lightcurve and simulation trial for the best
period solution (marked with red cross). Figure 2 shows a more challenging
case. All trial periods are too short to be detected with the ACF, and two
trial periods give similar predicted aliases. However, the alias strengths of
the simulated LSPs for $P=1.60$ day more closely match the observed alias
structure than the trial at $P=0.61$ day.

These simulations allowed us to select clear best periods for 471 out of 657
ZTF stars with strong LSP peaks. Only these 471 ZTF stars are tabulated and
used in our X-ray activity-rotation analyses below; true periodic behaviors
are probably present in other stars but we cannot identify the individual
periods in a reliable fashion. The final 471 ZTF-based periods, and associated
amplitudes of the variation in smoothed folded lightcurves, are listed in
Table 1.

The last 95 lines of the table list the NGC 2362 stars in common between the
Chandra X-ray (Broos et al., 2013; Getman et al., 2022) and CTIO
Blanco/Mosaic-II based rotation period (Irwin et al., 2008a) samples. Here the
rotation periods and amplitudes of folded light curves are from Irwin et al.
(2008a); stellar masses, effective temperatures, bolometric and X-ray
luminosities are from §A.

Among the 471 ZTF stars with reliable rotational periods, 154 are members of
the central part of the NGC 869 (= h Per) stellar cluster covered in the
Chandra survey (Getman et al., 2022). Previously, based on their multi-
telescope optical monitoring of a larger area of NGC 869, Moraux et al. (2013)
published rotation periods for 586 candidate cluster members. There are 84
stars that are common between our current ZTF analysis and the Moraux et al.
study. Figure 3 shows that the ZTF periods are generally consistent with the
periods (or half-periods) of Moraux et al. (2013), with the exception of seven
stars, for which our simulations prefer shorter-period solutions. Similar
discrepancies between stellar rotational periods derived from different
datasets analyzed with different methods are shown in Briegal et al. (2022).

Finally we note that some periodic photometric signals that we interpret to be
rotationally modulated spotted single stars may actually be eclipsing binary
systems with a F- or G-type primary and an M-type secondary (Triaud et al.,
2013). We suspect this source of contamination is not common, as our target
stars are restricted to X-ray and kinematic cluster members (Paper I) rather
than drawn from a large population of field stars.

##  4 Skumanich-Like Relations for late-PMS stars

###  4.1  Rotation Dependence on Mass and Age

For the X-ray-rotation sample in Table 1, Figure 4 displays the distribution
of rotation periods as a function of stellar mass in three mass strata $-$
$(0.6-0.9)$, $(0.9-1.2)$, and $(1.2-1.5)$ M⊙ $-$ and four cluster age strata:
5 Myr, 7 Myr, 13 Myr and 22-25 Myr. A handful of stars lie beyond $M<0.6$ M⊙
and $M>1.5$ M⊙. The l-PMS cluster distributions are compared to the well-
studied Pleiades ZAMS-MS stars. At a given cluster age, the mass strata
distinguish between stars of different evolutionary phases, such as Hayashi
versus Henyey tracks at 7 Myr, 13 Myr, and $>20$ Myr. These evolutionary
phases, called $ph_{HRD}$ here, are taken from the PARSEC 1.2S mass tracks by
Chen et al. (2014). Appendix §C provides further details and tabulation of
such phases for all l-PMS stars from Paper I.

Figure 4: Rotation-mass diagrams for seven l-PMS clusters (colored points)
compared to members of the Pleiades cluster (grey points, Rebull et al., 2016,
see Appendix B). The X-ray-rotation clusters are arranged into 4 age strata,
and individual stars are color-coded by their evolutionary phase on the
Hertzprung-Russell diagram: Hayashi track (green; $ph_{HRD}<2$), early Henyey
track (red; $2\leq ph_{HRD}<3$), late Henyey track (orange; $3\leq
ph_{HRD}<4$), and ZAMS-MS (blue; $4\leq ph_{HRD}<6$). The three mass strata
discussed in the text are marked by the grey lines. Figure legends present
numbers of plotted X-ray-rotation stars. Figure 5: Observed and calculated
evolution of stellar rotation. Black (Table 1) and magenta (previously
published) dots show rotational periods with box-and-whisker diagrams for
l-PMS clusters. Colored curves give the rotation evolutionary tracks from
Amard et al. (2019) and Gossage et al. (2021) (see text for details). Figure
legends indicate the cluster names and number of plotted stars.

An initial result from Figure 4 is that while some are more rapidly rotating
than Pleiades stars, others have quickly spun down to the $P_{rot}-M$ locus
seen during ZAMS stage. This is seen on both the Hayashi and Henyey tracks. In
these slower rotating stars, the changes in internal structure during the
l-PMS phase do not cause strong spin-up at the surface. This result was found
by Moraux et al. (2013) who state: “at least a fraction of slow rotators are
prevented from spinning up” during the l-PMS phase. We caution that our
methodology restricts period determination from ZTF lightcurves to only
reliable cases (§3), it is quite likely that some fast rotators are missing
from our sample. For example, rapid l-PMS rotators may produce multi-modal
surface field morphologies with smaller star spots rather than one large
hemispheric star spot needed for reliable photometric detection.

We note that the slower rotators ($2<P_{rot}<10$ days) in our l-PMS sample
show greater scatter than the narrow locus of slow rotators in the Pleiades in
the $P_{rot}-M$ diagrams. This is not due to observational uncertainties, as
$P_{rot}$ values are measured with high accuracy from our periodograms. This
suggests that the angular momenta of an ensemble of l-PMS stars have not yet
settled down to the well-defined rotational patterns seen on the ZAMS.

The $P_{rot}$ distribution can be compared to modern stellar evolutionary
models with self-consistent treatment of rotation internal angular momentum
transport and rotational mixing processes (Amard et al., 2019; Gossage et al.,
2021). Gossage et al. present models with two magnetic effects: traditional
magnetic braking of Matt et al. (2015) and consideration of the morphology of
stellar surface magnetic fields by Garraffo et al. (2018).

Figure 5 compares the temporal evolution of stellar rotation predicted by
these models with our rotational data. Individual stars from our X-ray-
rotation ZTF sample (Table 1) and their associated box-and-whisker diagrams
are shown in black while stars with previously published rotation periods
(Table 3) are shown in magenta. The colored curves give the rotation
evolutionary tracks. Models from Amard et al. (2019) are in red for three
different initial rotation rates (slow, median, fast) and three different
masses per figure panel ($0.6,0.8,0.9$ M⊙, $0.9,1,1.2$ M⊙, and $1.2,1.3,1.5$
M⊙). Models from Gossage et al. (2021) are shown for magnetic braking
prescriptions of Matt et al. (2015, in blue) and Garraffo et al. (2018, in
green) for a single mass per figure panel ($0.8,1,1.3$ M⊙) and for two initial
rotation rates (fast and slow). During the e-PMS phase at ages $\lesssim 4$
Myr, the models assume stellar rotation rates are constant due to disk
coupling. Figure legends indicate numbers of plotted stars.

Figure 5 shows broad agreement between observations and theory. The median and
quartile rotation rates increases following disk decoupling during $7-25$ Myr
are due to contraction along the Hayashi and Henyey tracks. The spin-up is
particularly notable in our NGC 1960 and NGC 2232 stars at ages $22-25$ Myr.
Contraction slows during $30-100$ Myr and rotation slows on the main sequence
due to centrifugal momentum loss by a magnetic wind. The predicted U-shaped
pattern of rotation periods is seen in all mass strata with the fastest
rotations seen in NGC 1960 and NGC 2232.

Some discrepancy, however, can be noted. Particularly for the lower $0.6-0.9$
M⊙ mass stratum, a quarter or more of the observed rotations are slower than
the rotation model tracks: the models are missing a very slow rotator
population. This could be be attributed to disk locking timescales
substantially longer than $\sim 4$ Myr in many lower-mass stars (Ribas et al.,
2015; Richert et al., 2018). It thus does not represent a fundamental
disagreement between data and models.

The bifurcation in $P_{rot}$ distribution in the NGC 869 cluster reported by
Moraux et al. (2013) is confirmed using our independent rotation data. It is
also present in the rich NGC 884 cluster; it appears as a deficit of stars
with $P_{rot}\simeq 1$ day (Figures 4 and 5).

For the $0.9-1.2$ M⊙ mass stratum we compare the distributions of rotation
periods for the rich samples of similar aged UCL/LCC and NGC 869 clusters
taken from space-based data of Rebull et al. (2022) and ground-based data of
Moraux et al. (2013), respectively. Using the Anderson-Darling 2-sample test,
the period distribution of Moraux et al. (2013) is significantly different
from UCL ($p<0.001$) and is probably different from LCC ($p=0.02$). UCL and
LCC have indistinguishable period distributions. The NGC 869 and UCL/LCC
differences are mainly due to the difference in spread rather than in median.
The nonparametric Ansari-Bradley test, where the alternative hypothesis
assumes same medians but different spreads, give $p<0.001$ that they have the
same spread. Compared to the space-based UCL/LCC data, the NGC 869 data
exhibit an excess of rapidly rotating stars.

Moraux et al. suggest that the excess may be due to a population of tidally
locked close binary systems identified by their photometric locus in a color-
magnitude diagram. The bifurcation in rotational periods thus may not reflect
a stage in the rotational evolution of single stars. Alternatively, the
existence of very rapid rotators in l-PMS clusters with $P_{rot}$ as short as
8 hours may be related to the truncation of protoplanetary disks in wide
binary stars (Meibom et al., 2007).

However, while Rebull et al. (2022) derive space-based rotation periods for
the vast majority of the UCL/LCC cluster members, the NGC 869 data of Moraux
et al. (2013) and our own NGC 869+884 ZTF data are highly incomplete in the
determination of the rotation periods. We thus can not evaluate the true
fraction of rapid rotators in these clusters as well as their nature.

###  4.2 X-ray Dependence on Mass, Rotation and Age

We use the stellar sample from Table 1 to study relationships between X-ray
luminosity and stellar rotation, particularly to examine whether
supersaturated, saturated and unsaturated regimes are present in $5-25$ Myr
l-PMS stars as they are in ZAMS stars.

Figure 6 shows the l-PMS mass-stratified $L_{X}-P_{rot}$ distributions. As the
X-ray data consist of a mixture of detections and upper limits, we use methods
of survival analysis to evaluate the presence of (anti)correlation (Feigelson
& Babu, 2012). We calculate Helsel’s generalized Kendall’s tau $p$-value for
null hypothesis of no correlation, and the Akritas-Theil-Sen linear regression
slope (Helsel, 2012; Akritas et al., 1995). These are implemented with
functions cenken and akr in CRAN package NADA within the R statistical
software environment (Lee, 2020). We looked for correlations in mass-
stratified $L_{X}-P_{rot}$ diagrams for individual clusters as well as cluster
ensembles, for all rotation periods as well rotations constrained to rapid
($P_{rot}\lesssim 1$ day) and slow ($P_{rot}\gtrsim 1$ day) ranges with
different boundaries.

In the $0.6-0.9$ M⊙ stratum (dominated by Hayashi-track stars), stellar
samples for individual clusters do not show any $L_{X}-P_{rot}$ correlations,
but the cluster ensemble sample exhibits a weak positive correlation with
slope $L_{X}\propto P_{rot}^{0.2}$ (black line in Figure 6a). The generalized
Kendall’s $\tau$ p-value is $p\simeq 0.0001$. This positive correlation in
l-PMS stars is opposite to the well-known anticorrelation between $L_{X}$ and
$P_{rot}$ for unsaturated ZAMS stars. But the l-PMS effect may be
coincidental. $L_{X}$ and $P_{rot}$ decay with age for different reasons:
X-ray luminosity falls as the distributed dynamo weakens (Getman et al., 2022)
while rotation period declines as contracting stars spin up following disk
decoupling (Figure 5).

In the $0.9-1.2$ and $1.2-1.5$ M⊙ mass strata, our X-ray-rotation samples are
dominated by the Henyey-track stellar members of $13$ Myr old, rich NGC $869$
and NGC $884$ clusters (Figure 6b and c). Other individual clusters are
represented by sparse samples with no $L_{X}-P_{rot}$ correlations. For the
entire rotation range, using either the entire sample of all open clusters (as
shown in figure panels b and c) or just a sub-sample of the dominant $13$ Myr
old clusters, $L_{X}$ appears independent of $P_{rot}$ (Kendall’s $\tau$
p-value $p>0.03$). But one can construct a sub-sample of rapidly rotating
l-PMS stars where $L_{X}$ increases with $P_{rot}$ (magenta lines, $p\simeq
0.005$). This is consistent with the supersaturation result of A16 for
$1.0-1.4$ M⊙ stars in NGC 869.

In Figure 6 panel d where only X-ray detections are considered (to mimic the
X-ray sensitivity level of the shallower Chandra observation from A16), we
find relations $L_{X}\propto P_{rot}^{0.6}$ for $P_{rot}<1$ day ($p\simeq
0.05$, which is not significant) and $L_{X}\propto P_{rot}^{-0.4}$ ($p\simeq
0.005$, which is significant) for $P_{rot}>1$ day. The latter is similar to
the non-saturation result of A16.

A16 interpret these results as the onset of the supersaturated and unsaturated
tachoclinal dynamo regimes that dominate the ZAMS phase.

Figure 6: X-ray luminosity as a function of rotation period for the l-PMS
stars in Table 1. Panels (a,b,c) show $L_{X}-P_{rot}$ distributions for
different mass strata. X-ray detections (circles) and X-ray non-detections
(triangles) are color-coded according to stellar age: $5$ Myr (pink), $7$ Myr
(red), $13$ Myr (blue), and $>20$ Myr (green). Black and magenta lines are
calculated using the Akritas-Thiel-Sen statistical procedure; see text for
detail. Panel (d) shows $L_{X}-P_{rot}$ distribution for 13 Myr old stars that
were detected in X-rays in the $0.9-1.5$ M⊙ range.

But we should also emphasize strong selection effects in our sample of l-PMS
stars in Figure 6. While the underlying Gaia-$Chandra$ samples of cluster
members are large and reasonably complete for solar-mass stars, only $\simeq
13$% have reliable rotational periods from our ZTF periodogram analysis. Only
stars with large hemisphere-size star spots have reliable rotational periods.
The ZTF-Gaia-Chandra sub-sample lacks weaker X-ray stars, so the distribution
of X-ray luminosities differs considerably from the full Gaia-$Chandra$ sample
for NGC 869$+$884 (Paper I). It is further possible that our rotational
periods for NGC 869 and NGC 884 stars are missing long-period stars that were
found by Moraux et al. (2013) (Figure 5) and is thereby unable to reveal the
unsaturated anticorrelation between $L_{X}$ and $P_{rot}$. We thus conclude
that our findings in Figure 6 are consistent with, but not a convincing
independent demonstration of results reported by A16.

We note that the relation $L_{X}\propto P_{rot}^{-0.4}$ for $P_{rot}>1$ day
shown in Figure 6d is much shallower than the slope of $\simeq-2$ measured for
main sequence stars (Pizzolato et al., 2003; Wright et al., 2011b).
Furthermore, the longest rotation periods of $\sim 10$ days for our l-PMS
stars correspond to Rossby number of $Ro=P/\tau_{conv}\sim 0.1$ (Table 1)
whereas the unsaturated $L_{X}$-rotation regime in older stars starts at
$Ro\gtrsim 0.13$. It may be more realistic to view the $L_{X}\propto
P_{rot}^{-0.4}$ effect found here and by A16 as a mild (factor of 2-3)
decrease of X-ray activity within a saturated dynamo regime.

##  5 Late Pre-Main Sequence Star Spot Sizes

A byproduct of photometric searches for rotational periods is a measurement of
the fractional surface area occupied by large star spots. If the stellar
radius is known from Gaia multiband photometry and parallax measurements, then
the star’s spot area can also be calculated. This quantity is not evaluated
from the periodogram where the LSP peak power is difficult to interpret, but
from the light curve folded modulo the period selected by the process
described in §3. We discuss here star spot sizes as a measure of surface
activity in addition to X-ray luminosity. A forthcoming study will examine the
utility of star spot sizes to estimate ages for PMS and ZAMS stars.

Star spot fractional coverages are obtained from $loess$ local regression fit
to folded ZTF light curves, as shown by the green curve in panel d of Figure
1. The spot coverage fraction is the peak-to-peak amplitude of this fitted
function in magnitudes. The measured amplitudes for the ZTF lightcurves are
listed in column 9 of Table 1. Previously reported amplitudes by Irwin et al.
(2008a) and Moraux et al. (2013) for NGC 2362 and NGC 869 stars with Gaia-
based masses are listed in Table 3.

Spot areas are calculated assuming that the observed periodic photometric is
caused by rotational modulation of a single dominant spot (or a conglomerate
of nearby smaller spots). Spot sizes are estimated following Herbst et al.
(2021). Spot area $A_{spot}$ is proportional to the observed drop in
brightness relative to the average stellar flux ($\Delta F/F_{av}$), projected
stellar hemisphere area ($A_{star}=\pi R_{star}^{2}$), and stellar ($T_{eff}$)
and spot ($T_{spot}$) temperatures:

$A_{spot}=A_{star}\Delta F/F_{av}(1-(T_{spot}/T_{eff})^{4})^{-1}$ (1)

where the stellar and spot temperatures are related as $T_{spot}=-3.58\times
10^{-5}T_{eff}^{2}+1.0188T_{eff}-239.3$. $\Delta F/F_{av}$ is obtained from
the peak-to-peak magnitude amplitude ($Ampl$) as $\Delta
F/F_{av}=10^{0.2Ampl}-10^{-0.2Ampl}\approx Ampl$. See Kounkel et al. (2022)
for a related calculation.

Figure 7 shows that X-ray luminosity and star spot area are strongly
correlated as $L_{X}\propto A_{spot}^{0.4}$ using the Akritas-Thiel-Sen
procedure. NGC 869 and NGC 884 stars, dominating our X-ray-rotation sample,
are omitted from this plot since their limited $L_{X}$ and $A_{spot}$ dynamic
ranges and low sensitivity towards $<0.9$ M⊙ stars weakens the correlation. If
these clusters are included, the dependence is still statistically significant
but is weaker, $L_{X}\propto A_{spot}^{0.1}$. The power law index ranges from
0.3 to 0.6 when other mass ranges or subsamples are considered.

We believe this is the first demonstration that X-ray emission and spot sizes
are correlated measures of surface magnetic activity in a stellar population.
However, the derived slope should not be interpreted astrophysically because
our X-ray-rotation stellar samples from Table 1 are badly incomplete due to
the absence of reliable ZTF rotational periods for most Chandra-Gaia stars
from Paper I. The spot areas are only measurable in ZTF periodograms when the
star spot area is both large and distributed asymmetrically across the surface
so that one hemisphere is substantially darkened compared to the other
hemisphere.

Figure 7: X-ray luminosity as a function of surface spot area for l-PMS stars.
The sample in this plot is restricted to stars with ages $t=[5-7]$ Myr and
$t>20$ Myr (omitting NGC 869 and NGC 884 stars at 13 Myr stars) and masses
$0.6-1.5$ M⊙. X-ray detections and upper limits are marked by the solid
circles and triangles, respectively. The red line shows the Akritas-Thiel-Sen
slope.

More insight might emerge from independent examination of how $L_{X}$ and
$A_{spot}$ evolve during the l-PMS phase in Figure 8. The $L_{X}-t$ plot in
panel a is taken from Paper I where the decay slope during the l-PMS phase was
derived employing all known Chandra-Gaia stellar members of open clusters that
have mass completeness limits $\lesssim 0.7$ M⊙. The evolution of spot area in
panels b and c are based on all stars with known peak-to-peak amplitudes given
in Tables 1 and 3.

Figure 8: Temporal evolution of X-ray luminosity and spot area in late PMS
stars. (a) Evolution of X-ray luminosity through the e-PMS and l-PMS phases,
adapted from Paper I (Getman et al., 2022, see their Figure 8 for details).
(b) Empirical cumulative distribution function of spot area for four age
strata: 5 Myr (pink), 7 Myr (red), 13 Myr (dark blue), and $22-25$ Myr
(green). The 13 Myr sample from Moraux et al. (2013) is shown for comparison
(light blue). (c) For the distributions in panel (b), median values with
bootstrapped 68% confidence intervals with weighted linear fit. The samples
for panels (b,c) are drawn from Tables 1 and 3 of the current paper. Table 2:
Temporal Evolution of X-ray Luminosity and Spot Area

Mass | $N(L_{X})$ | $b$ |  $p$-val(b) | $N(A_{spot})$ | $a$ |  $p$-val(a) | $b/a$
---|---|---|---|---|---|---|---
$M_{\odot}$ |  |  |  |  |  |  |
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8)
$0.7-0.9$ | 556 | $-0.52\pm 0.22$ | 0.07 | 178 | $-0.70\pm 0.33$ | 0.12 | $0.74\pm 0.47$
$0.9-1.2$ | 845 | $-0.61\pm 0.16$ | 0.01 | 370 | $-0.88\pm 0.10$ | 0.00 | $0.69\pm 0.20$
$0.7-1.2$ | 1769 | $-0.83\pm 0.30$ | 0.03 | 548 | $-0.77\pm 0.17$ | 0.02 | $1.08\pm 0.46$
$0.9-1.5$ | 1335 | $-1.16\pm 0.15$ | 0.00 | 636 | $-0.99\pm 0.17$ | 0.01 | $1.17\pm 0.25$
$0.7-1.5$ | 2259 | $-1.10\pm 0.24$ | 0.00 | 814 | $-0.85\pm 0.16$ | 0.01 | $1.29\pm 0.37$

Note. — Column 1: Mass range. Column 2: Number of PMS stars evolved in the
analysis of the temporal evolution of $L_{X}$ using data and methods from
Paper I. Columns 3-4: Results from the linear regression fits for the relation
$L_{X}\sim t^{b}$. The results include the slope ($b$) with 68% error, and
p-value for the hypothesis of zero slope. Column 5: Number of PMS stars
evolved in the analysis of the temporal evolution of spot area ($A_{spot}$).
Columns 6-7: Results from the linear regression fits for the relation
$A_{spot}\sim t^{a}$. The results include the slope ($a$) with 68% error, and
p-value for the hypothesis of zero slope. Column 8: Inferred slopes and 68%
errors for the relation $L_{X}\sim A_{spot}^{b/a}$.

Figure 8 demonstrates that spot area, like X-ray luminosity, drops several-
fold during the l-PMS phase. The plots show that $A_{spot}\propto t^{a}$ with
$a\simeq-0.9$ and $L_{X}\propto t^{b}$ with $b\simeq-1.1$ for the $0.7-1.5$ M⊙
mass range. Evolutionary trends for other mass ranges are listed in Table 2.
These results are consistent among different rotation datasets: for ages $5-7$
Myr, the spot areas are similar for our ZTF (red point in Figure 8c) and CTIO-
based data/analyses of Irwin et al. (2008a, pink point); and for age $13$ Myr,
our ZTF spot areas (blue point) agree with Moraux et al. (2013, cyan point).

Figure 8c dramatically shows that the decline of the median stellar spot area
with age follows a power law relation.

Figure 9: Same as in Figure 8c but for the amplitude of the photometric
variations.

Some studies of the evolution of star spots (e.g. Rebull et al., 2018) use the
amplitude of the photometric variations rather than the spot area. In our
sample of $5-25$ Myr stars when stellar radii can vary substantially, we find
that the terms $A_{star}$ and $Ampl$ in equation 1 also decay with time, but
slower than $A_{spot}$. For instance, Figure 9 shows that for the $(0.7-1.5)$
M⊙ stellar sample from Figure 8b,c the amplitude decreases with time as
$Ampl=t^{-0.53\pm 0.15}$. The stellar hemisphere area decreases approximately
as $A_{star}\sim t^{-0.3}$, while $A_{spot}\propto t^{-0.9}$ (Table 2). This
temporal evolution of $A_{spot}$ supports the magnetic activity paradigm (see
below). Spot area is thus an effective measure of star spot evolution in this
l-PMS regime.

Hemispheric spot area fraction ($f_{spot}=A_{spot}/A_{star}$) decays
approximately as $f_{spot}\sim t^{-0.5}$. This is consistent with the Kepler
and TESS study of six clusters with ages ranging from $10$ Myr to 4 Gyr
(Morris, 2020), for which the authors report $f_{spot}\sim t^{c}$ with
$c=-0.37\pm 0.16$.

Table 2 shows that both $A_{spot}$ and $L_{X}$ decrease approximately as
$t^{-1.0}$ for solar mass stars with a somewhat shallower drop around
$t^{-0.7}$ for sub-solar mass stars.

We thus find that the time-averaged X-ray luminosity and surface spot area of
l-PMS stars are nearly proportional to each other with $b/a\simeq 1$. A
plausible explanation is that local spot magnetic field strengths are roughly
constant but the areas of spots vary and evolve among l-PMS stars. This is
motivated by Zeeman broadening measurements in solar-type stars from PMS to
main sequence. They show a constant spot magnetic field at $B\sim 3$ kG but
magnetic filling factor ($f$), and thus disk-averaged magnetic field strength
$<B>=Bf$, decreasing with age (Kochukhov et al., 2020; Sokal et al., 2020).
For our l-PMS stars, assuming that the spot magnetic field strength is also
constant, the $L_{X}\propto A_{spot}^{1}$ relationship (Figure 8 and Table 2)
implies that the time-averaged X-ray luminosity would scale linearly with
disk-averaged stellar magnetic flux $L_{X}\propto\Phi\sim B\times
f_{spot}\times R^{2}\sim B\times A_{spot}$. The implied $L_{X}\propto\Phi$
scaling appears to be universal in the Sun and stars (§6). This scaling would
imply that $L_{X}\propto A_{spot}^{1}$ would be seen in l-PMS stars, or any
sample of normal stars.

Figure 10: The Pevtsov relation linking X-ray luminosity to magnetic flux.
Panel (a) is reproduced from Getman et al. (2021). The solid and dashed black
lines are the relations $L_{X}\sim\Phi^{m}$, where $m=1.15$ (Pevtsov et al.,
2003) and $m=1.48$ (Kirichenko & Bogachev, 2017), respectively. Magnetic
structures include the quiet Sun (brown), solar X-ray-bright points (light
green), solar active regions (cyan), solar disk averages (gray), and old G, K,
and M dwarfs (magenta). The solar microflares (golden oval) are from
Kirichenko & Bogachev (2017) and solar M-class flares (red oval) are from Su
et al. (2007). Individual solar-type dwarfs (blue squares) are from Kochukhov
et al. (2020). Individual Orion and Taurus e-PMS stars (dark green squares)
are from Sokal et al. (2020). Panels (b), (c) and (d) expand the high-
luminosity regime with the Pevtsov et al. (2003) and Kirichenko & Bogachev
(2017) regression lines shown. Panel (b) shows e-PMS stellar members of Orion
Nebula, IC 348, and NGC 2264 from Table 6 of Paper I (light green) with
kinetic magnetic flux estimated according to equation 2. All these stars are
on Hayashi tracks ($ph_{HRD}<2$). The purple curves are spline fits to the
25%, 50%, and 75% quartiles of the $L_{X}-\Phi$ distribution (see text for
details). Panels (c) and (d) show l-PMS stars from Paper I with stars on early
Henyey tracks in red ($2.0\leq ph_{HRD}<2.7$) and stars on later Henyey and
ZAMS-MS tracks in orange ($2.7\leq ph_{HRD}<6.0$). The purple curves show
smoothed 50% and 75% Kaplan-Meier quartiles of the $L_{X}-\Phi$ distributions
(see text for details).

##  6 X-ray Emission and Surface Magnetic Flux

Two decades ago, Pevtsov et al. (2003) found a remarkable power law
relationship that links solar and stellar magnetic events; a version is shown
in Figure 10a. Here the characteristic level of X-ray luminosity $L_{X}$
scales with the magnetic flux integrated over the responsible surface region
of area $A$, $\Phi=B\times A$, according to $L_{X}\propto\Phi^{m}$ with
$m=1.0-1.5$ (Pevtsov et al., 2003; Kirichenko & Bogachev, 2017)444 The
individual early PMS (dark green squares) and older solar-type stars (blue
squares) shown in Figure 10a have magnetic fluxes based on Zeeman broadening
measurements. Zeeman Doppler imaging measurements produce a slope $m\gtrsim
1.5$ while Zeeman broadening follow the standard slope (Vidotto et al., 2014;
Reiners et al., 2022). It is unclear whether the slope difference is arises
from methodology or is intrinsic to the stars (Kochukhov et al., 2020;
Kochukhov, 2021). . It points to a universal mechanism of coronal plasma
heating from a supply of magnetic flux to the stellar surface, regardless of
the nature of underlying magnetic dynamos. The relationship is present from
the smallest solar structures to powerful solar and stellar flares.

There are indications that the relationship changes for highly magnetically
active conditions. For instance, the X-ray super-flares that frequently occur
in e-PMS stars may require a steeper slope with $m\simeq 3$ and unusually
strong local surface fields of up to 20 kG (Getman et al., 2021; Zhuleku et
al., 2021).

For rapidly rotating stars with convection-driven dynamos such as e-PMS or
mid-M dwarfs, the available kinetic energy of convective flows may set the
strength of the surface magnetic field (Christensen et al., 2009; Reiners et
al., 2009; Reiners & Christensen, 2010; Reiners et al., 2022). These
researchers provide a simple analytic formula for the average strength of
surface magnetic field, where magnetic pressure ($\propto B^{2}$) depends on
the stellar density $\rho^{1/3}$ (where $\rho\sim M/R^{3}$) and stellar
bolometric flux $q^{2/3}$ (where $q\sim L/R^{2}$). Normalized to the empirical
measurements of magnetic field strength for various stars and planets in units
of $G$, the expression is

$B_{kin}=4800(ML^{2}/R^{7})^{1/6}$ (2)

where $M$, $L$, and $R$ are stellar mass, luminosity, and radius in solar
units.

We applied this formula to fully and partially convective l-PMS Gaia-Chandra
stars listed in Tables 6 and 4 of Paper I. For comparison to e-PMS stars, we
choose a sub-sample of 1038 lightly absorbed ($A_{V}<5$ mag) members of the
nearby IC 348, Orion Nebula, and NGC 2264 star forming regions with masses
$0.2-2$ M⊙. These are the nearest, rich regions for which we have most
sensitive Chandra X-ray data. These e-PMS stars are shown in Figure 10b (light
green) with quantile regression curves (purple)555 These curves are generated
with the Constrained B-Splines method that combines spline and quantile
regression fits. The 50% (median) curve is shown with a solid purple curve and
75% and 25% quartiles are marked by dashed curves. The procedure uses
quadratic splines with knot selection based on the Akaike Information
Criterion. The method is implemented in the R CRAN cobs package (Ng &
Maechler, 2007, 2020). . The fits show that e-PMS Chandra X-ray luminosities
and kinetic magnetic surface fluxes ($B_{kin}\times 4\pi R^{2}$ from equation
2) roughly follow the universal relationship of Pevtsov et al. (2003) and
Kirichenko & Bogachev (2017) though with X-ray luminosities reduced by factors
of $2-3$. A linear slope of $m=1.3$ resulting from an independent linear
regression fit to these e-PMS data (not shown) supports this finding.

We similarly apply equation (2) to the l-PMS members of open clusters in Table
4 of Paper I. Here the Gaia-Chandra sub-samples used are complete down to
$\simeq 0.7$ M⊙, not $\simeq 0.2$ M⊙ as in the e-PMS case. Like their younger
e-PMS siblings, numerous fully-convective l-PMS stars lying on Hayashi tracks
($ph_{HRD}<2.0$) follow the universal law (figure is not shown). So do the
l-PMS stars that lie on the early segments of Henyey tracks ($2.0\leq
ph_{HRD}<2.7$). These are shown as red dots in Figure 10c with spline fits666
The l-PMS samples in Figure 10c and d include a substantial number of X-ray
upper limits (shown as open trianges) as well as detections (shown as filled
circles). To treat the limits in a statistically consistent fashion, we
calculate the nonparametric maximum likelihood Kaplan-Meier estimator (KM;
Kaplan & Meier, 1958) for the $L_{X}$ distributions in sliding bins of $\Phi$.
Local quadratic fits are then calculated to the 50% and 75% quartiles of the
KM estimators to generate smooth curves. Due to high fractions of X-ray
nondetections in many $\Phi$ bins, the unreliable KM 25% quartiles are
omitted. Software calculations are made in R with the survival (Therneau,
2020) and loess (Cleveland et al., 1992) functions. .

However, more massive l-PMS stars on the more evolved Henyey track and on the
ZAMS-MS show a different behavior. Stars on later segments of the Henyey
tracks have evolutionary phases of $2.7<ph_{HRD}<4.0$ in Appendix C, and ZAMS-
MS stars have $4<ph_{HRD}<6$. These stars, shown as orange dots in the
$L_{X}-\Phi$ diagram of Figure 10d with spline fits, exhibit an $anti-
correlation$ between $L_{X}$ and $\Phi$, much different from the universal
$correlation$ found by Pevtsov et al. (2003); Kirichenko & Bogachev (2017) and
seen in our younger l-PMS and e-PMS stars in Figure 10b and c.

As the universal $L_{X}-\Phi$ relationship is present for both e-PMS and main
sequence stars, we do not believe it is absent for l-PMS stars. Our finding in
§5 that $L_{X}\propto A_{spot}^{1}$ (Table 2) further supports a correlation
between X-ray emission and surface magnetic flux levels. Rather, the anti-
correlation seen in Figure 10d is arguably due to a failure of equation (2) in
estimating magnetic field strengths based on a convective dynamo. This
provides indirect evidence that the dynamos during the later stages of l-PMS
and ZAMS evolution are switching from $\alpha^{2}$ to $\alpha\Omega$
processes.

This change in $L_{X}-\Phi$ behavior occurs approximately two-thirds of the
duration between the base of the Henyey track and the time of maximum
bolometric luminosity on the Henyey track before descent to the ZAMS. Recall
that these results are based on PARSEC 1.2S evolutionary models.

This result agrees with different arguments in Paper I that pointed to the
$\alpha\Omega$ dynamo as dominant in these more massive and older l-PMS stars.
Paper I found that the X-ray luminosity of $1-3.5$ M⊙ stars (most of which lie
on late Henyey or ZAMS tracks) decays as $L_{X}\propto t^{-1.8}$. This decay
rate is consistent with that of main sequence stars, whose magnetic activity
is driven by the tachoclinal $\alpha\Omega$ dynamo (Maggio et al., 1987; Güdel
et al., 1997; Feigelson et al., 2004).

##  7 Conclusions

###  7.1 Observational Accomplishments

From an observational viewpoint, this study together with Paper I has greatly
increased our knowledge of stellar activity and rotation for solar-type stars
in the elusive late-pre-main sequence phase of stellar evolution. We generate
samples of several thousand stars within open clusters with ages $5-25$ Myr.
Combining Gaia and $Chandra$ observations allows cluster membership to be
established with high reliability. Paper I presents a catalog of $\sim 8000$
l-PMS stars ($\sim 6000$ and $\sim 2000$ lying in the cluster centers and
outskirts, respectively), most of which have X-ray detections.

In this study, we examined photometric light curves for 3410 of these l-PMS
using the northern hemisphere multi-epoch survey of the Zwicky Transient
Facility. We establish a rigorous procedure for identifying reliable
photometric periodicities due to rotationally modulated star spots using light
curves obtained at a single telescope site. This analysis gives rotation
periods for 471 Gaia-Chandra-ZTF stars ranging from around 0.5 to 8 days (§3).
Several dozen stellar members of NGC 2362 with Chandra-CTIO measurements were
added to this X-ray rotation sample (Appendix A and Table 1). For some
analyses, we also consider previously published rotation periods for
additional $>3500$ stars, mostly without Chandra data (Appendix B).

The resulting mass-dependent evolution of rotation through the l-PMS phase
agrees well with astrophysical models, although our sample points to a sub-
population of stars with slower initial rotations than commonly assumed
(Figure 5).

###  7.2 Switching from Convective to Tachoclinal Dynamos

The primary scientific question of the study is raised in the Introduction:
When during the complex changes in stellar structure of the l-PMS evolutionary
phase does the magnetic dynamo change from the $\alpha^{2}$ distributed
convective process to the solar-type $\alpha\Omega$ tachoclinal process? The
X-ray/rotation study of the l-PMS cluster h per by Argiroffi et al. (2016)
gave evidence, based on a few slowly rotating stars, that the transition to an
$\alpha\Omega$ dynamo was beginning for more massive stars at age $\simeq 13$
Myr. Our $L_{X}-P_{rot}$ diagram for stars in a wider range of cluster ages
from $\simeq 5-25$ Myr (Figure 6) is consistent with this result but does not
provide convincing independent support. The distributions are also consistent
with the hypothesis that X-ray luminosities have no relationship to rotational
period in l-PMS stars.

Our ZTF-based periods are limited to periods shorter than $\sim 10$ days where
the stars are still in a magnetically saturated regime with Rossby number
being below $\sim 0.1$. Definitive results may require measuring rotation
rates spectroscopically for mass-stratified samples of l-PMS covering a range
of X-ray luminosities.

However, this study together with Getman et al. (2022) provide two other lines
of evidence for dynamo changes. First, Paper I showed strong declines in X-ray
luminosity in l-PMS stars. The decline is dramatic in higher-mass $1-3.5$ M⊙
stars and slower in lower-mass $0.75-1$ M⊙ stars. The $L_{X}\propto t^{-1.8}$
decline in $1-3.5$ M⊙ l-PMS stars is consistent with that of young main
sequence stars where surface activity is driven by the $\alpha\Omega$ dynamo.

Second, a subtle effect is seen at the later stages of l-PMS evolution as
stars approach the ZAMS. There is clear evidence that the universal relation
between $L_{X}$ and stellar magnetic flux $\Phi$ applies at all evolutionary
phases from e-PMS through l-PMS, ZAMS and older main sequence stars. Our
empirical results support this with the correlation $L_{X}\propto
A_{spot}^{1}$ (§5, Figure 8, Table 2). Further assuming convective dynamos are
dominant, the $L_{X}\propto\Phi^{1}$ relation, where $\Phi$ is a theory based
magnetic flux, holds for stars in the e-PMS and beginning of l-PMS stages (§6,
Figure 10b,c). But the convective dynamo flux generation equation (2) fails to
reveal the expected $L_{X}-\Phi$ correlation at later l-PMS stages (Figure
10d). This suggests that a tachoclinal $\alpha\Omega$ processes emerge as the
dominant dynamo mechanism as stars move along the late part of Henyey track in
the Hertzsprung-Russell diagram and approach the ZAMS.

##  8 Acknowledgments

We offer this study in memory of our late friend and colleague Leisa K.
Townsley (Penn State). We are grateful to the anonymous referee for thoughtful
and helpful comments that improved the manuscript. This project is supported
by the Chandra grant GO9-20011X (K. Getman, Principal Investigator) and the
Chandra ACIS Team contract SV4-74018 (G. Garmire & E. Feigelson, Principal
Investigators), issued by the Chandra X-ray Center, which is operated by the
Smithsonian Astrophysical Observatory for and on behalf of NASA under contract
NAS8-03060. The Chandra Guaranteed Time Observations (GTO) data used here were
selected by the ACIS Instrument Principal Investigator, Gordon P. Garmire, of
the Huntingdon Institute for X-ray Astronomy, LLC, which is under contract to
the Smithsonian Astrophysical Observatory; contract SV2-82024.

This research has also made use of the NASA/IPAC Infrared Science Archive,
which is funded by the National Aeronautics and Space Administration and
operated by the California Institute of Technology. Based on observations
obtained with the Samuel Oschin Telescope 48-inch and the 60-inch Telescope at
the Palomar Observatory as part of the Zwicky Transient Facility project. ZTF
is supported by the National Science Foundation under Grants No. AST-1440341
and AST-2034437 and a collaboration including current partners Caltech, IPAC,
the Weizmann Institute for Science, the Oskar Klein Center at Stockholm
University, the University of Maryland, Deutsches Elektronen-Synchrotron and
Humboldt University, the TANGO Consortium of Taiwan, the University of
Wisconsin at Milwaukee, Trinity College Dublin, Lawrence Livermore National
Laboratories, IN2P3, University of Warwick, Ruhr University Bochum,
Northwestern University and former partners the University of Washington, Los
Alamos National Laboratories, and Lawrence Berkeley National Laboratories.
Operations are conducted by COO, IPAC, and UW.

This work has made use of data from the European Space Agency (ESA) mission
Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing
and Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement.

##  Appendix A NGC 2362 Properties

$Chandra$-ACIS-I 100 ks exposure of this cluster was analysed by our team as
part of the 20-region Massive Young Star-forming complex study in Infrared and
X-ray (MYStIX) survey (Feigelson et al., 2013). The X-ray study emerged with
491 X-ray selected cluster members tabulated by Broos et al. (2013). Stellar
X-ray luminosities and ages based on the $Age_{JX}$ chronometer (Getman et
al., 2014) are listed in Table 6 of Paper I.

Here we re-examine the earlier distance and age estimates for the cluster
using Gaia DR3 data (Gaia Collaboration et al., 2016, 2022) and methods
discussed in Paper I. Four hundred forty six Chandra cluster members from
Broos et al. (2013) have Gaia DR3 counterparts within $1\arcsec$ matching
radius. The median distance from the Sun (with 68% bootstrap error) based on
186 stars with Gaia re-normalized unit weight errors (RUWE) values $<1.4$ and
parallax accuracy $\sigma_{\omega}<0.1$ mas, are $1301\pm 9$ pc. This is
consistent with, and more accurate than, the previous Gaia DR2-based distance
of $1332\pm 70$ pc (Kuhn et al., 2019).

Fitting redden theoretical PARSEC 1.2S isochrones (Bressan et al., 2012; Chen
et al., 2014) to the Gaia DR3 color-magnitude distribution of 351 Chandra
members with the error on the Gaia $G_{BP}-G_{RP}$ color of $<0.1$ mag, gives
a cluster extinction and age of $A_{V}=0.3$ mag and $t=5$ Myr, respectively.
The extinction coefficients from Luhman & Esplin (2020) were used to redden
the theoretical PARSEC isochrones. This isochronal fit is shown in the Gaia
color-magnitude diagram in Appendix B. Following Paper I, effective
temperatures, bolometric luminosities, and masses for individual cluster
members are obtained through a comparison of the observed $G$-band absolute
magnitudes with theoretical predictions for the best-fit ($5$ Myr) reddened
PARSEC 1.2S isochrone.

Ninety-five X-ray cluster members are located in the central part of the
cluster covered by the $\it Chandra$ observation; have Gaia DR3-based masses;
and have rotation periods measured by Irwin et al. (2008a) as part of the
Monitor project. The cluster was observed with a 6 minute cadence in the $i$
band for several hours on 8 nearly-contiguous nights with the CTIO Blanco 4-m
telescope. Typical photometric precision is $2-3$ mmag at $i=15$ mag and $5-6$
mmag at $i=18$ mag. Rotational periods are obtained from least-squares fits to
sinusoidal functions, a threshold on the reduced chi-squared value, and visual
vetting. Reliability is established with simulations. The last 95 records in
Table 1 here present these stars and rotational periods.

The X-ray luminosities for NGC 2362 members are calculated using the XPHOT
procedure (Getman et al., 2010). This requires the presence of both soft ($<2$
keV) and hard ($>2$ keV) X-ray photons, but many NGC 2362 members lack hard
X-ray photons. Table 6 of Getman et al. (2022) provides XPHOT-based intrinsic
X-ray luminosities ($L_{X}$) for only 207 (out of 491) brighter X-ray members
of NGC 2362. Using the 207 bright X-ray stars, we calculate a single apparent-
to-intrinsic X-ray flux conversion factor by averaging the ratios of $L_{X}$s
to apparent X-ray photometric fluxes $Photon\\_Flux\\_t$, given in Table 2 of
Broos et al. (2013). This factor is applied to the apparent Chandra fluxes of
the remaining weaker stars to determine their intrinsic X-ray luminosities.
For all 95 NGC 2362 X-ray members with available rotation periods from Irwin
et al. (2008a), their inferred intrinsic X-ray luminosities are listed in
Table 1.

No additional non-X-ray Gaia members were added to this NGC 2362 stellar
sample, because the X-ray sample by itself is already complete down to $M=0.6$
$M_{\odot}$, a low mass limit of interest in this X-ray-rotation study, thanks
to the relatively long Chandra observation exposure time for this unobscured
cluster. Thus all X-ray upper limit records in Column 11 of Table 1 for this
cluster are marked as “…” (no data).

##  Appendix B Gaia-DR3 Masses For Previously Published Stars With Rotation
Periods

This section presents Table 3, listing Gaia-based masses for stars most of
which lack Chandra-X-ray data but have previously published rotation periods,
as well as Figure 11 showing Gaia-DR3 color-magnitude diagrams for such stars.

##  Appendix C Updated HRD Phases for l-PMS Stars from Paper I

Table 4 in Paper I lists the properties of 7924 Chandra-Gaia stellar members
of 10 open clusters. The Column 11 of that table gives an integer flag
indicating the evolutionary status of the star on the Hertzsprung-Russell
diagram based on the predictions of the PARSEC 1.2S evolutionary models.
Original HRD statuses in the PARSEC 1.2S mass-track files of Chen et al.
(2014) are given as floating-point real numbers, with the integer parts
indicating critical evolutionary points and the fractional parts indicating
fractional time durations between the critical points. The critical points for
young stars include: $=1$ \- the onset of Hayashi track; $=2$ \- onset of
early Henyey track; $=3$ \- onset of late Henyey track (when $L_{bol}$ reaches
maximum on Henyey track); $=4$ \- onset of ZAMS; $=5$ \- onset of phase when
hydrogen burning is fully active. According to Chen et al. (2014), fractional
parts are “proportional to the fractional time duration between the current
point and the beginning of that phase”, as
$frac=(t-t_{beg,i})/(t_{end,i}-t_{beg,i})$ and $i=\\{1,2,3,4,5\\}$ for young
stars. Table 4 provides more accurate (including the fractional parts) HRD
statuses for the 7924 l-PMS stars of Paper I. These new stellar HRD statuses
are further employed in Figures 4 and 10.

Figure 11: Gaia-DR3 color-magnitude diagrams for stars with previously
published rotation periods. For NGC 2362, the shown sample is the MYStIX X-ray
sample discussed in §A. For the remaining clusters, the shown stars are from
Table 3. The red curves represent the best-fit reddened PARSEC 1.2S
evolutionary models. Figure titles list the inferred best fit cluster
extinction and age. Corresponding stellar masses are indicated by orange
markers. Figure legends provide numbers of plotted stars. Table 3: Gaia-DR3
Masses For Stars With Previously Published Rotation Periods

Cluster | R.A. | Decl. | $G$ | $\sigma_{G}$ | $B-R$ | $\sigma_{BR}$ | $\log(T_{eff})$ | $M$ | $\log(L_{bol})$ | $P_{rot}$ | $Ampl$
---|---|---|---|---|---|---|---|---|---|---|---
| deg | deg | mag | mag | mag | mag | K | $M_{\odot}$ | $L_{\odot}$ | day | mag
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9) | (10) | (11) | (12)
NGC2362 | 109.395103 | -25.168670 | 16.816 | 0.003 | 1.379 | 0.007 | 3.62 | 0.79 | -0.38 | 7.34 | 0.010
NGC2362 | 109.412882 | -25.192423 | 16.510 | 0.003 | 1.396 | 0.007 | 3.63 | 0.85 | -0.29 | 7.77 | 0.006
NGC2362 | 109.500341 | -25.137874 | 19.014 | 0.003 | 3.090 | 0.065 | 3.51 | 0.39 | -0.93 | 1.14 | 0.030
NGC2362 | 109.518947 | -25.240841 | 18.059 | 0.003 | 2.118 | 0.020 | 3.55 | 0.60 | -0.69 | 9.34 | 0.018
NGC2362 | 109.554830 | -25.236658 | 20.009 | 0.004 | 2.270 | 0.093 | 3.47 | 0.23 | -1.23 | 10.80 | 0.050
NGC2362 | 109.564310 | -25.247717 | 18.288 | 0.003 | 2.496 | 0.032 | 3.54 | 0.49 | -0.76 | 5.46 | 0.040
NGC2362 | 109.596027 | -25.194042 | 17.907 | 0.003 | 1.887 | 0.014 | 3.56 | 0.63 | -0.64 | 18.17 | 0.010
NGC2362 | 109.592932 | -25.261247 | 17.210 | 0.006 | 1.706 | 0.040 | 3.59 | 0.70 | -0.48 | 2.80 | 0.194
NGC2362 | 109.601028 | -25.141686 | 16.454 | 0.007 | 1.672 | 0.030 | 3.64 | 0.94 | -0.26 | 12.02 | 0.052
NGC2362 | 109.638649 | -25.191095 | 18.971 | 0.003 | 2.805 | 0.061 | 3.51 | 0.43 | -0.93 | 3.61 | 0.022

Note. — This table is available in its entirety (3879 stars with Gaia
properties and rotation periods) in the machine-readable form in the online
journal. Gaia-based stellar properties (effective temperatures, masses, and
bolometric luminosities) are estimated in the current paper. Rotation periods
are from Irwin et al. (2008a, NGC 2362), Rebull et al. (2018, Upper Sco),
Moraux et al. (2013, NGC 869), Rebull et al. (2022, UCL & LCC), and Rebull et
al. (2016, Pleiades). For NGC 2362 and NGC 869, there is overlap between the
current sample and the X-ray-rotation sample given in Table 1. Column 1:
Cluster name. Columns 2-3: Gaia DR3 right ascension and declination (in
decimal degrees) for epoch J2000.0. Columns 4-7: Gaia DR3 magnitude and color
and their uncertainties. Only stars with $\sigma_{BR}<0.1$ mag are included in
this table and corresponding Figure 11. Columns 8-10: Stellar effective
temperature, mass, and bolometric luminosity derived from the Gaia color-
magnitude diagrams. Columns 11-12: Rotation periods and peak-to-peak
amplitudes from Irwin et al. (2008a); Moraux et al. (2013); Rebull et al.
(2016, 2018, 2022). The amplitudes are only available for NGC 2362 and NGC
869.

Table 4: New HRD Statuses For l-PMS Stars From Paper I

Cluster | R.A. | Decl. | $ph_{HRD}$
---|---|---|---
| deg | deg |
(1) | (2) | (3) | (4)
NGC1502 | 62.257371 | 62.385897 | 1.11
NGC1502 | 61.944486 | 62.153143 | 5.00
NGC1502 | 61.577576 | 62.266908 | 1.25
NGC1502 | 61.904523 | 62.509605 | 1.20
NGC1502 | 62.374648 | 62.351225 | 1.11
NGC1502 | 61.665573 | 62.324773 | 2.05
NGC1502 | 61.679861 | 62.340361 | 2.23
NGC1502 | 61.758327 | 62.315094 | 1.71
NGC1502 | 61.764918 | 62.335727 | 1.11
NGC1502 | 61.769816 | 62.227697 | 2.44

Note. — This table is available in its entirety (7924 stars) in the machine-
readable form in the online journal. Additional stellar properties for these
l-PMS stars including the X-ray and bolometric luminosities, stellar masses,
and effective temperatures can be found in Table 4 of Paper I. Column 1:
Cluster name. Columns 2-3: Gaia DR3 right ascension and declination (in
decimal degrees) for epoch J2000.0. Columns 4: Updated status of the stellar
evolutionary phase. See text for details.

## References

  * Affer et al. (2013) Affer, L., Micela, G., Favata, F., Flaccomio, E., & Bouvier, J. 2013, MNRAS, 430, 1433, doi: 10.1093/mnras/stt003
  * Akritas et al. (1995) Akritas, M. G., Murphy, S. A., & Lavalley, M. P. 1995, Journal of the American Statistical Association, 90, 170, doi: 10.1080/01621459.1995.10476499
  * Alexander & Preibisch (2012) Alexander, F., & Preibisch, T. 2012, A&A, 539, A64, doi: 10.1051/0004-6361/201118100
  * Amard et al. (2019) Amard, L., Palacios, A., Charbonnel, C., et al. 2019, A&A, 631, A77, doi: 10.1051/0004-6361/201935160
  * Argiroffi et al. (2016) Argiroffi, C., Caramazza, M., Micela, G., et al. 2016, A&A, 589, A113, doi: 10.1051/0004-6361/201526539
  * Baluev (2008) Baluev, R. V. 2008, MNRAS, 385, 1279, doi: 10.1111/j.1365-2966.2008.12689.x
  * Baraffe & Chabrier (2010) Baraffe, I., & Chabrier, G. 2010, A&A, 521, A44, doi: 10.1051/0004-6361/201014979
  * Bellm et al. (2019) Bellm, E. C., Kulkarni, S. R., Graham, M. J., et al. 2019, PASP, 131, 018002, doi: 10.1088/1538-3873/aaecbe
  * Bouvier (2007) Bouvier, J. 2007, in Star-Disk Interaction in Young Stars, ed. J. Bouvier & I. Appenzeller, Vol. 243, 231–240, doi: 10.1017/S1743921307009593
  * Bouvier et al. (1993) Bouvier, J., Cabrit, S., Fernandez, M., Martin, E. L., & Matthews, J. M. 1993, A&A, 272, 176
  * Bressan et al. (2012) Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127, doi: 10.1111/j.1365-2966.2012.21948.x
  * Briegal et al. (2022) Briegal, J. T., Gillen, E., Queloz, D., et al. 2022, MNRAS, 513, 420, doi: 10.1093/mnras/stac898
  * Broos et al. (2013) Broos, P. S., Getman, K. V., Povich, M. S., et al. 2013, ApJS, 209, 32, doi: 10.1088/0067-0049/209/2/32
  * Brown et al. (2011) Brown, B. P., Miesch, M. S., Browning, M. K., Brun, A. S., & Toomre, J. 2011, ApJ, 731, 69, doi: 10.1088/0004-637X/731/1/69
  * Browning (2008) Browning, M. K. 2008, ApJ, 676, 1262, doi: 10.1086/527432
  * Caceres et al. (2019) Caceres, G. A., Feigelson, E. D., Jogesh Babu, G., et al. 2019, AJ, 158, 58, doi: 10.3847/1538-3881/ab26ba
  * Chabrier & Küker (2006) Chabrier, G., & Küker, M. 2006, A&A, 446, 1027, doi: 10.1051/0004-6361:20042475
  * Chen et al. (2014) Chen, Y., Girardi, L., Bressan, A., et al. 2014, MNRAS, 444, 2525, doi: 10.1093/mnras/stu1605
  * Christensen et al. (2009) Christensen, U. R., Holzwarth, V., & Reiners, A. 2009, Nature, 457, 167, doi: 10.1038/nature07626
  * Cleveland & Devlin (1988) Cleveland, W. S., & Devlin, S. J. 1988, Journal of the American Statistical Association, 83, 596–610, doi: 10.2307/2289282
  * Cleveland et al. (1992) Cleveland, W. S., Grosse, E., & Shyu, W. M. 1992, in Statistical Models in S, ed. J. Chambers & T. Hastie (Wadsworth & Brooks/Cole).  https://www.bibsonomy.org/bibtex/249dc786266d86d0cc54ed6ccb2fe9894/gron
  * Dawson & Fabrycky (2010) Dawson, R. I., & Fabrycky, D. C. 2010, ApJ, 722, 937, doi: 10.1088/0004-637X/722/1/937
  * Delisle et al. (2020) Delisle, J. B., Hara, N., & Ségransan, D. 2020, A&A, 635, A83, doi: 10.1051/0004-6361/201936905
  * Dobler et al. (2006) Dobler, W., Stix, M., & Brandenburg, A. 2006, ApJ, 638, 336, doi: 10.1086/498634
  * Edelson & Krolik (1988) Edelson, R. A., & Krolik, J. H. 1988, ApJ, 333, 646, doi: 10.1086/166773
  * Feiden & Chaboyer (2013) Feiden, G. A., & Chaboyer, B. 2013, ApJ, 779, 183, doi: 10.1088/0004-637X/779/2/183
  * Feigelson & Babu (2012) Feigelson, E. D., & Babu, G. J. 2012, Modern Statistical Methods for Astronomy with R Applications (Cambridge University Press)
  * Feigelson & Nelson (1985) Feigelson, E. D., & Nelson, P. I. 1985, ApJ, 293, 192, doi: 10.1086/163225
  * Feigelson et al. (2004) Feigelson, E. D., Hornschemeier, A. E., Micela, G., et al. 2004, ApJ, 611, 1107, doi: 10.1086/422248
  * Feigelson et al. (2013) Feigelson, E. D., Townsley, L. K., Broos, P. S., et al. 2013, ApJS, 209, 26, doi: 10.1088/0067-0049/209/2/26
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272
  * Gaia Collaboration et al. (2022) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2022, arXiv e-prints, arXiv:2208.00211.  https://arxiv.org/abs/2208.00211
  * Garraffo et al. (2018) Garraffo, C., Drake, J. J., Dotter, A., et al. 2018, ApJ, 862, 90, doi: 10.3847/1538-4357/aace5d
  * Getman et al. (2017) Getman, K. V., Broos, P. S., Kuhn, M. A., et al. 2017, ApJS, 229, 28, doi: 10.3847/1538-4365/229/2/28
  * Getman & Feigelson (2021) Getman, K. V., & Feigelson, E. D. 2021, ApJ, 916, 32, doi: 10.3847/1538-4357/ac00be
  * Getman et al. (2010) Getman, K. V., Feigelson, E. D., Broos, P. S., Townsley, L. K., & Garmire, G. P. 2010, ApJ, 708, 1760, doi: 10.1088/0004-637X/708/2/1760
  * Getman et al. (2021) Getman, K. V., Feigelson, E. D., & Garmire, G. P. 2021, ApJ, 920, 154, doi: 10.3847/1538-4357/ac1746
  * Getman et al. (2022) Getman, K. V., Feigelson, E. D., Garmire, G. P., et al. 2022, ApJ, 935, 43, doi: 10.3847/1538-4357/ac7c69
  * Getman et al. (2019) Getman, K. V., Feigelson, E. D., Kuhn, M. A., & Garmire, G. P. 2019, MNRAS, 487, 2977, doi: 10.1093/mnras/stz1457
  * Getman et al. (2014) Getman, K. V., Feigelson, E. D., Kuhn, M. A., et al. 2014, ApJ, 787, 108, doi: 10.1088/0004-637X/787/2/108
  * Gossage et al. (2021) Gossage, S., Dotter, A., Garraffo, C., et al. 2021, ApJ, 912, 65, doi: 10.3847/1538-4357/abebdf
  * Grankin & Artemenko (2009) Grankin, K. N., & Artemenko, S. A. 2009, Information Bulletin on Variable Stars, 5907, 1
  * Güdel et al. (1997) Güdel, M., Guinan, E. F., & Skinner, S. L. 1997, ApJ, 483, 947, doi: 10.1086/304264
  * Helsel (2012) Helsel, D. R. 2012, Statistics for censored environmental data using minitab and R, 2nd edn. (Hoboken, NJ: Wiley-Blackwell)
  * Henderson & Stassun (2012) Henderson, C. B., & Stassun, K. G. 2012, ApJ, 747, 51, doi: 10.1088/0004-637X/747/1/51
  * Herbst et al. (2021) Herbst, K., Papaioannou, A., Airapetian, V. S., & Atri, D. 2021, ApJ, 907, 89, doi: 10.3847/1538-4357/abcc04
  * Herbst et al. (2000) Herbst, W., Rhode, K. L., Hillenbrand, L. A., & Curran, G. 2000, AJ, 119, 261, doi: 10.1086/301175
  * Iben (1965) Iben, Icko, J. 1965, ApJ, 141, 993, doi: 10.1086/148193
  * Irwin et al. (2008a) Irwin, J., Hodgkin, S., Aigrain, S., et al. 2008a, MNRAS, 384, 675, doi: 10.1111/j.1365-2966.2007.12725.x
  * Irwin et al. (2008b) —. 2008b, MNRAS, 383, 1588, doi: 10.1111/j.1365-2966.2007.12669.x
  * Kaplan & Meier (1958) Kaplan, E. L., & Meier, P. 1958, Journal of the American Statistical Association, 53, 457, doi: 10.1080/01621459.1958.10501452
  * Käpylä et al. (2008) Käpylä, P. J., Korpi, M. J., & Brandenburg, A. 2008, A&A, 491, 353, doi: 10.1051/0004-6361:200810307
  * Kirichenko & Bogachev (2017) Kirichenko, A. S., & Bogachev, S. A. 2017, Sol. Phys., 292, 120, doi: 10.1007/s11207-017-1146-8
  * Kochukhov (2021) Kochukhov, O. 2021, A&A Rev., 29, 1, doi: 10.1007/s00159-020-00130-3
  * Kochukhov et al. (2020) Kochukhov, O., Hackman, T., Lehtinen, J. J., & Wehrhahn, A. 2020, A&A, 635, A142, doi: 10.1051/0004-6361/201937185
  * Koen (1990) Koen, C. 1990, ApJ, 348, 700, doi: 10.1086/168277
  * Kounkel et al. (2022) Kounkel, M., Stassun, K. G., Bouma, L. G., et al. 2022, AJ, 164, 137, doi: 10.3847/1538-3881/ac866d
  * Kovács et al. (2005) Kovács, G., Bakos, G., & Noyes, R. W. 2005, MNRAS, 356, 557, doi: 10.1111/j.1365-2966.2004.08479.x
  * Kuhn et al. (2019) Kuhn, M. A., Hillenbrand, L. A., Sills, A., Feigelson, E. D., & Getman, K. V. 2019, ApJ, 870, 32, doi: 10.3847/1538-4357/aaef8c
  * Kunitomo et al. (2017) Kunitomo, M., Guillot, T., Takeuchi, T., & Ida, S. 2017, A&A, 599, A49, doi: 10.1051/0004-6361/201628260
  * Landin et al. (2022) Landin, N. R., Mendes, L. T. S., Vaz, L. P. R., & Alencar, S. H. P. 2022, arXiv e-prints, arXiv:2212.13958.  https://arxiv.org/abs/2212.13958
  * Lee (2020) Lee, L. 2020, NADA: Nondetects and Data Analysis for Environmental Data.  https://CRAN.R-project.org/package=NADA
  * Ljung & Box (1978) Ljung, G. M., & Box, G. E. P. 1978, Biometrika, 65, 297, doi: 10.1093/biomet/65.2.297
  * Lu et al. (2022) Lu, Y. L., Curtis, J. L., Angus, R., David, T. J., & Hattori, S. 2022, AJ, 164, 251, doi: 10.3847/1538-3881/ac9bee
  * Luhman & Esplin (2020) Luhman, K. L., & Esplin, T. L. 2020, AJ, 160, 44, doi: 10.3847/1538-3881/ab9599
  * Maggio et al. (1987) Maggio, A., Sciortino, S., Vaiana, G. S., et al. 1987, ApJ, 315, 687, doi: 10.1086/165170
  * Maldonado et al. (2022) Maldonado, J., Colombo, S., Petralia, A., et al. 2022, arXiv e-prints, arXiv:2204.12206.  https://arxiv.org/abs/2204.12206
  * Masci et al. (2019) Masci, F. J., Laher, R. R., Rusholme, B., et al. 2019, PASP, 131, 018003, doi: 10.1088/1538-3873/aae8ac
  * Matt et al. (2015) Matt, S. P., Brun, A. S., Baraffe, I., Bouvier, J., & Chabrier, G. 2015, ApJ, 799, L23, doi: 10.1088/2041-8205/799/2/L23
  * McQuillan et al. (2014) McQuillan, A., Mazeh, T., & Aigrain, S. 2014, ApJS, 211, 24, doi: 10.1088/0067-0049/211/2/24
  * Meibom et al. (2007) Meibom, S., Mathieu, R. D., & Stassun, K. G. 2007, ApJ, 665, L155, doi: 10.1086/521437
  * Moraux et al. (2013) Moraux, E., Artemenko, S., Bouvier, J., et al. 2013, A&A, 560, A13, doi: 10.1051/0004-6361/201321508
  * Morris (2020) Morris, B. M. 2020, ApJ, 893, 67, doi: 10.3847/1538-4357/ab79a0
  * Ng & Maechler (2007) Ng, P., & Maechler, M. 2007, Statistical Modelling, 7, 315.  http://smj.sagepub.com/content/7/4/315.abstract
  * Ng & Maechler (2020) Ng, P. T., & Maechler, M. 2020, COBS – Constrained B-splines (Sparse matrix based).  https://CRAN.R-project.org/package=cobs
  * Patten & Simon (1996) Patten, B. M., & Simon, T. 1996, ApJS, 106, 489, doi: 10.1086/192346
  * Pevtsov et al. (2003) Pevtsov, A. A., Fisher, G. H., Acton, L. W., et al. 2003, ApJ, 598, 1387, doi: 10.1086/378944
  * Pizzolato et al. (2003) Pizzolato, N., Maggio, A., Micela, G., Sciortino, S., & Ventura, P. 2003, A&A, 397, 147, doi: 10.1051/0004-6361:20021560
  * Preibisch et al. (2005) Preibisch, T., Kim, Y.-C., Favata, F., et al. 2005, ApJS, 160, 401, doi: 10.1086/432891
  * R Core Team (2020) R Core Team. 2020, R: A language and environment for statistical computing, R Foundation for Statistical Computing, Vienna, Austria.  https://www.R-project.org
  * R Core Team (2021) —. 2021, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria.  https://www.R-project.org/
  * Rebull et al. (2020) Rebull, L. M., Stauffer, J. R., Cody, A. M., et al. 2020, AJ, 159, 273, doi: 10.3847/1538-3881/ab893c
  * Rebull et al. (2018) —. 2018, AJ, 155, 196, doi: 10.3847/1538-3881/aab605
  * Rebull et al. (2022) Rebull, L. M., Stauffer, J. R., Hillenbrand, L. A., et al. 2022, AJ, 164, 80, doi: 10.3847/1538-3881/ac75f1
  * Rebull et al. (2016) Rebull, L. M., Stauffer, J. R., Bouvier, J., et al. 2016, AJ, 152, 113, doi: 10.3847/0004-6256/152/5/113
  * Reiners et al. (2009) Reiners, A., Basri, G., & Christensen, U. R. 2009, ApJ, 697, 373, doi: 10.1088/0004-637X/697/1/373
  * Reiners & Christensen (2010) Reiners, A., & Christensen, U. R. 2010, A&A, 522, A13, doi: 10.1051/0004-6361/201014251
  * Reiners et al. (2022) Reiners, A., Shulyak, D., Käpylä, P. J., et al. 2022, A&A, 662, A41, doi: 10.1051/0004-6361/202243251
  * Ribas et al. (2015) Ribas, Á., Bouy, H., & Merín, B. 2015, A&A, 576, A52, doi: 10.1051/0004-6361/201424846
  * Richert et al. (2018) Richert, A. J. W., Getman, K. V., Feigelson, E. D., et al. 2018, MNRAS, 477, 5191, doi: 10.1093/mnras/sty949
  * Ruf (1999) Ruf, T. 1999, Biological Rhythm Research, 30, 178
  * Scargle (1982) Scargle, J. D. 1982, ApJ, 263, 835, doi: 10.1086/160554
  * Schrijver & Zwaan (2000) Schrijver, C. J., & Zwaan, C. 2000, Solar and Stellar Magnetic Activity (Cambridge University Press)
  * Skumanich (1972) Skumanich, A. 1972, ApJ, 171, 565, doi: 10.1086/151310
  * Sokal et al. (2020) Sokal, K. R., Johns-Krull, C. M., Mace, G. N., et al. 2020, ApJ, 888, 116, doi: 10.3847/1538-4357/ab59d8
  * Somers & Pinsonneault (2015) Somers, G., & Pinsonneault, M. H. 2015, ApJ, 807, 174, doi: 10.1088/0004-637X/807/2/174
  * Squicciarini & Bonavita (2022) Squicciarini, V., & Bonavita, M. 2022, A&A, 666, A15, doi: 10.1051/0004-6361/202244193
  * Stauffer et al. (1997) Stauffer, J. R., Hartmann, L. W., Prosser, C. F., et al. 1997, ApJ, 479, 776, doi: 10.1086/303930
  * Su et al. (2007) Su, Y., Van Ballegooijen, A., McCaughey, J., et al. 2007, ApJ, 665, 1448, doi: 10.1086/519679
  * Süveges (2014) Süveges, M. 2014, MNRAS, 440, 2099, doi: 10.1093/mnras/stu372
  * Tamuz et al. (2005) Tamuz, O., Mazeh, T., & Zucker, S. 2005, MNRAS, 356, 1466, doi: 10.1111/j.1365-2966.2004.08585.x
  * Therneau (2020) Therneau, T. M. 2020, A Package for Survival Analysis in R.  https://CRAN.R-project.org/package=survival
  * Triaud et al. (2013) Triaud, A. H. M. J., Hebb, L., Anderson, D. R., et al. 2013, A&A, 549, A18, doi: 10.1051/0004-6361/201219643
  * VanderPlas (2018) VanderPlas, J. T. 2018, ApJS, 236, 16, doi: 10.3847/1538-4365/aab766
  * VanderPlas & Ivezić (2015) VanderPlas, J. T., & Ivezić, Ž. 2015, ApJ, 812, 18, doi: 10.1088/0004-637X/812/1/18
  * Venables & Ripley (2002) Venables, W. N., & Ripley, B. D. 2002, Modern Applied Statistics with S, 4th edn. (New York: Springer).  http://www.stats.ox.ac.uk/pub/MASS4
  * Vidotto et al. (2014) Vidotto, A. A., Gregory, S. G., Jardine, M., et al. 2014, MNRAS, 441, 2361, doi: 10.1093/mnras/stu728
  * Vrba et al. (1988) Vrba, F. J., Herbst, W., & Booth, J. F. 1988, AJ, 96, 1032, doi: 10.1086/114863
  * Wright et al. (2011a) Wright, N. J., Drake, J. J., Mamajek, E. E., & Henry, G. W. 2011a, ApJ, 743, 48, doi: 10.1088/0004-637X/743/1/48
  * Wright et al. (2011b) —. 2011b, ApJ, 743, 48, doi: 10.1088/0004-637X/743/1/48
  * Yadav et al. (2015) Yadav, R. K., Christensen, U. R., Morin, J., et al. 2015, ApJ, 813, L31, doi: 10.1088/2041-8205/813/2/L31
  * Zhuleku et al. (2021) Zhuleku, J., Warnecke, J., & Peter, H. 2021, arXiv e-prints, arXiv:2102.00982.  https://arxiv.org/abs/2102.00982