# SPYGLASS. V. Spatially and Temporally Structured Star-Forming Environments
in the Cepheus-Hercules Complex

Ronan Kerr Department of Astronomy, University of Texas at Austin
2515 Speedway, Stop C1400
Austin, Texas, USA 78712-1205
Adam L. Kraus Department of Astronomy, University of Texas at Austin
2515 Speedway, Stop C1400
Austin, Texas, USA 78712-1205
Daniel Krolikowski Steward Observatory, University of Arizona
933 N. Cherry Avenue
Tucson, Arizona, USA, 85721
Luke G. Bouma 51 Pegasi b Fellow  Department of Astronomy, California
Institute of Technology
1216 E. California Boulevard, Pasadena, CA 91125, USA
Juan P. Farias Department of Astronomy, University of Texas at Austin
2515 Speedway, Stop C1400
Austin, Texas, USA 78712-1205

###### Abstract

Young associations record complete histories of star forming events through
their demographics and dynamics, and Gaia has greatly expanded our knowledge
of these associations. We present the first structural and dynamical overview
study of the Cep-Her Complex, which has recently emerged as the largest
stellar population within 500 pc that still lacks substantial follow-up. We
reveal that Cep-Her is not a singular association, but rather an amalgam of
four, consisting of the older ($\tau>100$ Myr) open cluster Roslund 6, in
addition to three dynamically coherent and highly substructured young
associations that we focus on: Orpheus (25-40 Myr), Cinyras (28-43 Myr), and
Cupavo (54-80 Myr). With $9552\pm 960$ stars in Orpheus, $3872\pm 455$ stars
in Cinyras, and $8794\pm 1827$ stars in Cupavo, all three are among the
largest young associations within 500 pc, rivalling major associations like
Sco-Cen. Our novel view of the ages and dynamics of these associations reveals
evidence for sequential star formation in Cinyras, in addition to a multi-
origin pattern of stellar dispersal in Orpheus that may hint to the presence
of multiple clouds at formation. Dynamical simulations indicate that while
some embedded open clusters and central components of these associations are
converging, they likely lack the mass necessary to capture one another and
undergo hierarchical cluster assembly. Our results provide our first view of
the complex star-forming environments that gave rise to the associations of
Cep-Her, which will serve as important laboratories for future studies of star
and planet formation.

Stellar associations (1582); Stellar ages (1581); Star formation(1569) ; Pre-
main sequence stars(1290)

††facilities: Gaia, McDonald Observatory: Tull spectrograph/Harlan J. Smith
Telescope, W. M. Keck Observatory: HIRES Spectrograph/Keck I
Telescope††software: astropy (Astropy Collaboration et al., 2013, 2018), numpy
(Harris et al., 2020), saphires (Tofflemire, 2019), HDBSCAN (McInnes et al.,
2017) AMUSE (Portegies Zwart et al., 2009; Pelupessy et al., 2013; Portegies
Zwart et al., 2013; Portegies Zwart & McMillan, 2018), galpy (Bovy, 2015),
matplotlib (Hunter, 2007)

##  1 Introduction

After formation, most nearby young stars can be found in largely unbound
structures known as associations (de Zeeuw et al., 1999; Lada & Lada, 2003).
These associations leave a record of the star formation event that produced
them, making them essential for revealing the evolution of star formation
processes, especially on longer timescales (Krumholz et al., 2019; Krause et
al., 2020). Within 50 Myr of formation, many associations are still evolving
through a combination of gradual expansion in sparse regions and hierarchical
assembly and relaxation in dense regions (e.g., Lada et al., 1984; Pecaut &
Mamajek, 2016; Vázquez-Semadeni et al., 2017). Studies of young associations
across a range of ages therefore trace the evolving morphologies and dynamical
processes that lead to longer-lived structures such as open clusters
(Dalessandro et al., 2021).

The study of young associations has progressed substantially in the last few
years, accelerated by the monumental contributions of the Gaia spacecraft,
which is providing precise stellar positions and motions for 2 billion nearby
stars in the Milky Way (Gaia Collaboration et al., 2022). Many studies to date
have taken advantage of Gaia’s exquisite photometric and astrometric data,
discovering stellar populations over a wide range of ages and scales (e.g.,
Kounkel & Covey, 2019; Sim et al., 2019; Prisinzano et al., 2022; Kerr et al.,
2021; Hunt & Reffert, 2023). The surveys conducted through the SPYGLASS (Stars
with Photometrically Young Gaia Luminosities Around the Solar System) program
have discovered or more deeply characterized dozens of young ($\tau\lesssim
50$ Myr) associations through Gaia, detecting new associations across a range
of scales and morphologies (Kerr et al., 2021, 2023). Kerr et al. (2021)
(hereafter SPYGLASS-I) reported 27 top-level associations within 333 pc,
nearly half of which remain largely unstudied to date, and Kerr et al. (2023)
(hereafter SPYGLASS-IV) expanded the survey to 1 kpc, revealing 116 young
associations, 30 of which have no direct literature equivalent. The
populations identified by these surveys vary from tenuous, low-density
associations that barely emerge from the sky background, to large and massive
populations with varied densities within, revealing an intricate network of
interconnected star formation.

Among the major emerging stellar complexes identified in SPYGLASS-IV is the
Cep-Her complex, which connects the SPYGLASS-I populations of Cerberus, Lyra,
and Cepheus-Cygnus with more distant open clusters like $\delta$ Lyr and RSG 5
(Stephenson, 1959; Röser et al., 2016). The presence of a larger population in
the area was first noted by Zari et al. (2018), who observed overdensities of
young stars in the direction of Cygnus and Lyra. Kounkel & Covey (2019)
grouped these populations into the broadly-defined Stephenson 1 ($\delta$ Lyr
Cluster) “string”, which contains the populations identified by Zari et al.
(2018), in addition to Cepheus Far North (e.g., Kerr et al., 2022b). SPYGLASS-
IV then provided the first contiguous definition of Cep-Her’s extent in space-
velocity coordinates, along with a comprehensive list of candidate members
present in Gaia DR3. Consisting of 1164 population-defining photometrically
young stars in addition to over $1.5\times 10^{5}$ candidate members spanning
more than 300 pc, SPYGLASS-IV revealed that Cep-Her may be one of the largest
young associations within 500 pc (de Zeeuw et al., 1999; Pecaut & Mamajek,
2016; Cantat-Gaudin et al., 2019a). However, unlike other large associations
such as Sco-Cen and Vela, Cep-Her has never been the target of a dedicated
survey, so its population, structure, dynamics, and star formation history are
almost entirely unknown. Furthermore, while we know that the stars in Cep-Her
are broadly young and co-moving in transverse velocity, we do not know whether
they remain co-moving after the radial velocity component is considered, or
whether they have ages consistent with a single, continuous star-forming
event. Adding these new dimensions of analysis may therefore reveal entirely
new star-forming environments that contribute to the Cep-Her Complex.

The masses and velocities of dynamically coherent sub-populations in Cep-Her
are critical for understanding the history and evolution of the association,
which will inform the future fate of other nearby complexes like Sco-Cen,
Vela, and Orion. The mass and subgroup sizes inform the virial state, which
sets whether stars in a population can diverge freely, or are gravitationally
bound within a dense open cluster (Kuhn et al., 2019). This distinction is
critical during dynamical traceback analyses, as the trajectories of stars in
bound clusters often change dramatically due to two-body encounters and the
overall gravitational field of the cluster. Unbound populations are then well-
suited for traceback, as their motions are largely unchanged since formation.
As has been shown in previous Gaia-based studies (e.g., Miret-Roig et al.,
2020, 2022; Galli et al., 2023) and SPYGLASS publications, the limited
dynamical evolution in young associations after gas dispersal allows for their
point of closest configuration at formation to be established, providing a
straightforward method for computing dynamical ages (Kerr et al., 2022a, b).
The balance between internal motions and gravitational binding can also inform
the future evolution on the scale of an entire association, informing whether
association subgroups will coalesce and hierarchically assemble into a single
clustered population over time, or whether they will drift apart. The
understudied nature of Cep-Her also opens the possibility of dynamically and
demographically distinct structures that have not emerged in the 5D space-
transverse velocity phase space that has been used in past studies.
Identifying any such substructure will better focus future studies by
identifying manageable sets of closely related stars which can be studied
together. By characterizing Cep-Her substructure, we can also take advantage
of Cep-Her’s substantial overlap with the Kepler field by providing ages and
formation environments for numerous young planet hosts, which are critical for
studies of planet formation and evolution (Bouma et al., 2022).

In this publication, we complete the first ever comprehensive demographic and
structural survey of the Cep-Her Complex. We progressively divide the
association using velocities and age to reveal its constituent associations
and subgroups, and revise the association’s stellar membership to estimate the
total mass and population size of Cep-Her and its substructures. In Section 2,
we introduce the data used to complete this survey, which we use to assess
memberships in Section 3. We then analyze Cep-Her’s substructure in Section 4,
before analyzing the demographics of those substructures and the population
more broadly in Section 5. We then discuss the association’s dynamics in
Section 6 and address the implications of Cep-Her’s structure and dynamics in
Section 7, before concluding in Section 8.

##  2 Data

###  2.1 Sample Selection

Our stellar census of Cep-Her is based on the sample that we defined in
SPYGLASS-IV. This publication identified over 4$\times$105 photometrically
young stars in Gaia DR3 using a Bayesian identification framework, which were
clustered into 116 young populations within 1 kpc using the HDBSCAN clustering
algorithm. The Cep-Her Complex (SCYA-96) was one of the largest populations
identified in that work, with its base sample of photometrically young
founding members consisting of 1164 stars. SPYGLASS-IV used that population of
young stars to identify additional stars with positions and velocities
consistent with membership, but uncertain photometric youth assessments. Young
stars in this broader category are generally earlier type stars which have
already merged onto the main sequence at the age of the Cep-Her, or have
fallen close enough to the pre-main sequence to be photometrically
indistinguishable from the field binary sequence. The fraction of the pre-main
sequence that is indistinguishable from main sequence binaries increases as
stars age, resulting in increasing rates of inconclusive photometric
membership until the pre-main sequence and field binary sequence merge
completely for $\tau>50-60$ Myr. To provide a non-photometric method of
assessing stellar youth, SPYGLASS-IV includes an assessment of membership
probability, $P_{\rm spatial}$ (Section 5; Kerr et al., 2023), which compares
the relative abundances of likely young and likely old stars in space-velocity
coordinates. Stars with $P_{\rm spatial}<0.05$ are excluded from the extended
sample of co-moving candidate members in SPYGLASS-IV, and the complete catalog
of potential Cep-Her members contains 152936 stars. The distribution of Cep-
Her founding members and candidates on-sky is shown in Figure 1.

Figure 1: Distribution of Cep-Her photometrically young founding members (dark
crosses) and candidate members with ambiguous photometry that have with
$P_{\rm spatial}>0.2$ (light dots), shown in RA/Dec sky coordinates. We also
label the locations of a few notable open clusters with arrows, and also label
the locations of the SPYGLASS-I groups of Cep-Cyg, Lyra, and Cerberus.

As Cep-Her is a relatively old population (29 Myr on average based on PARSEC
isochrones; Kerr et al., 2023), the SPYGLASS framework is expected to detect
young members at a rate of 10-20%. This suggests that the total population of
the Cep-Her Complex is in the range of 5000 to 10000 members (Kerr et al.,
2021), although Cep-Her is distant enough that not all photometrically young
stars will have high-quality Gaia entries, suggesting the total population may
be substantially higher. Nonetheless, with more than $1.5\times 10^{5}$ stars
in the extended population of candidates, we expect substantial field
contamination rates, likely caused by a combination of its crowded environment
near the galactic plane and spatial and dynamical dispersal of its members
into the field since formation. This contamination is reflected in the
photometric distribution of these candidate members, where the pre-main
sequence does not cleanly separate from a dominant field main sequence. We
therefore restrict the sample to only include stars with $P_{\rm spatial}>0.2$
in our main analysis. This restricts the sample to 34886 candidates, reducing
the population of field stars to allow the young association members to be
separated. This cut certainly removes genuine members, however this is
necessary to ensure manageable contamination rates at the outer edges of our
sample. We discuss the population of genuine members that may have been
removed by this cut in Section 5.8.

###  2.2 Supplemental Data

High-resolution optical spectra provide radial velocities in addition to
lithium 6708 Å (Li) and H$\alpha$ equivalent widths (EWs), which are essential
for membership verification and traceback in addition to dynamical and lithium
depletion age measurements. The Gaia sample provides RVs for approximately 18%
of the 34886 sources in our initial Cep-Her sample. However, Gaia does not
provide measurements of spectroscopic youth indicators, and Gaia RVs have
shown inconsistencies at the 1 km s-1 level among young stars that prevent
accurate dynamical age measurements (e.g., Kerr et al., 2022b). Cep-Her has
only recently been discovered in its entirety (Zari et al., 2018; Kerr et al.,
2023), making its past spectroscopic coverage less complete compared to other
nearby young associations. However, there are still large surveys that have
covered Cep-Her candidate members, as well as some targeted studies that have
focussed on well-known clusters within the larger sample like $\delta$ Lyrae.
The coverage of Cep-Her is nonetheless far from complete, especially in lower-
density subgroups, necessitating that these literature measurements be
supplemented with new observations. We therefore use a combination of
literature spectroscopic measurements and new observations to assemble a
robust spectroscopic sample for Cep-Her.

####  2.2.1 Literature RV and Spectroscopic Line Measurements

To assemble a sample of literature spectroscopic measurements, we compiled a
list of existing sources that provide RVs, Li EWs, or H$\alpha$ EWs. The
search included large surveys like LAMOST, APOGEE, and Gaia (e.g., Abdurro’uf
et al., 2022; Zhang et al., 2021; Gaia Collaboration et al., 2022), as well as
smaller surveys that appear in Vizier and Simbad searches of Cep-Her and other
young associations in the SPYGLASS-IV sample. We also collected additional
sources from Luhman (2022) and Žerjal et al. (2023), who provide their own
compilations of recent spectroscopic studies of young stars. We additionally
consult the RV source compilation in Anderson & Francis (2012). We include
sources with an RV quality grade better than a “C”111A “C” quality source may
have uncertainties that do not account for larger systematic errors, which
avoids sources with substantial uncorrected systematic uncertainties (Anderson
& Francis, 2012). The compilation from Gontcharov (2006) is taken to broadly
represent observations which precede it, with the exception of the northern
hemisphere survey provided by Famaey et al. (2005), which was listed as an
“A”222An “A” grade indicates reliable uncertainties quality source separate
from Gontcharov (2006) in Anderson & Francis (2012). Studies that compile
measurements from multiple sources are otherwise excluded, so our citations
mark the publication where the measurement first appeared.

For targets with multiple observations, precedence was generally given to the
source with the lower uncertainty. However, for sources with $\sigma_{RV}<1$
km s-1, we de-prioritized sources with known systematic uncertainties at the 1
km s-1 level so that high-quality sources are not overwritten in favor of
observations with substantial systematic uncertainties (Anderson & Francis,
2012; Zhang et al., 2021). The sources we deprioritized in this manner were
Gaia, RAVE DR6 (Kunder et al., 2017), and all LAMOST catalogs (Wang et al.,
2022) with the exception of the self-consistent MRS RV measurements described
in Zhang et al. (2021). A complete list of adopted measurements and their
sources is provided in Table 1. Most of the 6518 RVs we adopt are from Gaia
(5456), with the next largest RV sources being LAMOST (506 across LRS and MRS
spectra), APOGEE (148), and (Gontcharov, 2006) (54). We also record flags for
spectroscopic binaries where available, and include those flags in Table 1.

Table 1: Cep-Her data containing both literature values and new observations. The source for each value is included. High-quality literature values supersede our observations in cases of poor RV results or low-signal line measurements. The full table, which contains 6518 entries, is available in the online-only version of this paper.  Gaia DR3 ID | RA | Dec | RV (km s-1) | EWLi (Å) | EWHα (Å)
---|---|---|---|---|---
| (deg) | (deg) | val | err | srcaaThe source for the measurement, either 2.7 Harlan J. Smith Telescope at the McDonald Observatory (HJST), Keck HIRES (HIRES), or an external source. The abbreviations for external sources are: Gaia DR3 (DR3), APOGEE DR17 (Abdurro’uf et al., 2022) (AP17), GALAH DR3 (Buder et al., 2021) (GAL3), LAMOST DR9 MRS spectra (Zhang et al., 2021) (LAM9M), LAMOST DR8 LRS parameters of A, F, G, and K stars (LAM8L), LAMOST DR8 Catalog of M Stars (LAM8LMCAT), LAMOST DR8 Value Added Catalog (Wang et al., 2022) (LAM8LV), Famaey et al. (2005) (F5), Gontcharov (2006) (G6), Grieves et al. (2018) (G18), Klutsch et al. (2020) (K20), Tautvaišienė et al. (2022) (T22)  | val | err | srcaaThe source for the measurement, either 2.7 Harlan J. Smith Telescope at the McDonald Observatory (HJST), Keck HIRES (HIRES), or an external source. The abbreviations for external sources are: Gaia DR3 (DR3), APOGEE DR17 (Abdurro’uf et al., 2022) (AP17), GALAH DR3 (Buder et al., 2021) (GAL3), LAMOST DR9 MRS spectra (Zhang et al., 2021) (LAM9M), LAMOST DR8 LRS parameters of A, F, G, and K stars (LAM8L), LAMOST DR8 Catalog of M Stars (LAM8LMCAT), LAMOST DR8 Value Added Catalog (Wang et al., 2022) (LAM8LV), Famaey et al. (2005) (F5), Gontcharov (2006) (G6), Grieves et al. (2018) (G18), Klutsch et al. (2020) (K20), Tautvaišienė et al. (2022) (T22)  | val | err | srcaaThe source for the measurement, either 2.7 Harlan J. Smith Telescope at the McDonald Observatory (HJST), Keck HIRES (HIRES), or an external source. The abbreviations for external sources are: Gaia DR3 (DR3), APOGEE DR17 (Abdurro’uf et al., 2022) (AP17), GALAH DR3 (Buder et al., 2021) (GAL3), LAMOST DR9 MRS spectra (Zhang et al., 2021) (LAM9M), LAMOST DR8 LRS parameters of A, F, G, and K stars (LAM8L), LAMOST DR8 Catalog of M Stars (LAM8LMCAT), LAMOST DR8 Value Added Catalog (Wang et al., 2022) (LAM8LV), Famaey et al. (2005) (F5), Gontcharov (2006) (G6), Grieves et al. (2018) (G18), Klutsch et al. (2020) (K20), Tautvaišienė et al. (2022) (T22)
2069631150945233280 | 309.7628 | 43.4871 | -19.92 | 0.04 | HJST | 0.138 | 0.011 | HJST | 0.00 | 0.00 | HJST
2069673344703582592 | 309.6291 | 43.8329 | -35.08 | 3.07 | GDR3 |  |  |  |  |  |
2069674512934689536 | 309.6933 | 43.8082 | -14.54 | 4.94 | GDR3 |  |  |  |  |  |
2069755361399374080 | 308.7965 | 43.7338 | -22.32 | 3.02 | GDR3 |  |  |  |  |  |
2069764535449283712 | 308.1051 | 43.3637 | -21.91 | 1.19 | GDR3 |  |  |  |  |  |
2069765287066380416 | 308.3421 | 43.2745 | -24.54 | 4.06 | LAM8LV |  |  |  |  |  |
2069767902703713536 | 308.4111 | 43.4374 | -20.25 | 0.03 | HJST | 0.187 | 0.011 | HJST | 0.00 | 0.00 | HJST
2069836209864147840 | 309.2886 | 43.6723 | -22.28 | 1.88 | LAM8LV |  |  |  | 2.42 | 0.02 | LAM8LMCAT
2069869848049047552 | 309.5119 | 44.1504 | -21.37 | 3.41 | GDR3 |  |  |  |  |  |
2069875757923114240 | 309.5231 | 44.1710 | -35.06 | 6.83 | GDR3 |  |  |  |  |  |
2069876170239978112 | 309.5263 | 44.2266 | -20.05 | 0.04 | HJST | 0.209 | 0.011 | HJST | 0.00 | 0.02 | HJST
2069890154653506432 | 309.5657 | 44.3327 | -18.70 | 3.00 | G6 |  |  |  |  |  |

####  2.2.2 New Observations

To deepen and broaden our spectroscopic coverage of Cep-Her, we performed new
spectroscopic observations with two different telescopes and instruments: the
Tull Coudé Echelle Spectrograph at the McDonald Observatory’s 2.7m Harlan J.
Smith Telescope (HJST; Tull et al., 1995), and the High Resolution Echelle
Spectrometer on the Keck I 10m telescope (HIRES; Vogt et al., 1994). Both
instruments provide high-resolution optical spectra, with the Tull
spectrograph providing non-contiguous coverage from 3400 and 10900 Å at
$R\approx 60{,}000$, and HIRES coverage spanning 3600 to 7900 Å at $R\approx
30{,}000$. The observational configurations were selected to cover the Li 6708
Å and H$\alpha$ lines wherever possible, with resolution sufficient for
resolving the rotational broadening of these lines for each star in our
sample. The high resolution of these instruments combines with their stable
wavelength solutions to enable RV measurements with sub-km s-1 precision.

At the HJST, target selection was optimised to enable dynamical traceback and
provide lithium depletion sequences across a diverse set of locations in Cep-
Her. We first performed preliminary clustering, similar to what is described
in Section 4, identifying stars in regions where subgroup age measurements are
likely possible. We then observed a sample of young stellar populations spread
evenly throughout the preliminary subgroups. Stars in dense regions such as
the $\delta$ Lyrae cluster were de-prioritized due to the higher likelihood of
them being gravitationally bound, which prevents the free expansion of members
after formation that simplifies dynamical analysis. Our limiting magnitude at
the HJST was $m_{G}<14$, and we focused our observations on later-type stars
that have fewer past observations, more informative lithium abundances, and
less rotational broadening, which makes precise radial velocities more
accessible. Our programs observed a total of 351 stars in Cep-Her.

The HIRES data were acquired as part of a targeted survey focused on two of
the Cep-Her subgroups described by Bouma et al. (2022) that contain
statistically-validated transiting planets: RSG-5 and $\delta$ Lyr. The main
purpose of the program was to acquire spectra for stars spanning spectral
types of M1V to G2V in order to improve prospects for lithium-dating the
association. We acquired spectra for 19 stars in $\delta$ Lyr, and 26 in
RSG-5. We selected these stars by beginning with the candidates reported in
Table 2 of Bouma et al. (2022). We selected stars with “weights” exceeding
0.1, and removed any objects with visual, astrometric, spectroscopic, or
photometric indications of binarity based on Gaia DR3. We then required the
stars to have detected rotation periods with either TESS or ZTF, and required
them to be in the range of 12.0$<$V$<$16.2, where the faint limit was set by
time and S/N constraints. The final selected objects were a random draw from
the set of objects that met these constraints. The observations themselves
were executed between April–July of 2023, without iodine, and otherwise using
the standard setup and reduction techniques of the California Planet Survey
(Howard et al., 2010).

The HJST observations were processed and reduced following typical echelle
spectroscopy procedures using a publicly available pipeline tailored to the
Tull spectrograph333https://github.com/dkrolikowski/tull_coude_reduction. We
computed RVs for the HJST targets using spectral line broadening functions
from the saphires package (Tofflemire, 2019). This produced sub-km s-1 radial
velocity measurements for 271 stars in our sample. For HIRES, we measured
radial velocities using SpecMatch-Synth (Petigura et al., 2017), a synthetic
template-matching code built on an interpolation of the model spectra from
Coelho et al. (2005). We also extracted Li and H$\alpha$ equivalent width
measurements from both the HIRES and the HJST measurements by fitting those
lines with Gaussian profiles, following the method described in SPYGLASS-II.
We combine the final results for RV, Li EW, and H$\alpha$ EW from both HIRES
and the HJST with the literature results described in Section 2.2.1, making
final RV and EW selections according to the rules described there.

##  3 Membership

By combining the Cep-Her catalog provided in SPYGLASS-IV with new
spectroscopic measurements, we can assess the membership probability of
candidate members using four different membership markers: the local fractions
of young and old stars, photometry, RVs, and Lithium. Together, these markers
can be used to produce aggregate membership probabilities that take into
account all available data, allowing us to estimate key demographic
properties. Furthermore, we can use membership probabilities to remove
probable non-members, producing a sufficiently pure sample to analyse the
structure of the association.

###  3.1 Spatial-Photometric Membership Probability

SPYGLASS-IV has already produced measurements for probability of youth,
$P_{\rm Age<50Myr}$, which is available for all stars that pass the
photometric and astrometric quality flags presented in that paper. SPYGLASS-IV
also provides a probability of membership calculation, $P_{\rm spatial}$,
which is based on the relative frequencies of high-confidence young and old
objects among adjacent stars in space-velocity coordinates. Since $P_{\rm
Age<50Myr}$ is based on a constant star formation rate over the last 11.2 Gyr
(Binney et al., 2000), it implicitly assumes that $\simeq 0.4$% of the sample
has ages $\lesssim 50$ Myr. However, the stars identified as members of Cep-
Her in SPYGLASS-IV exist at points in parameter space where young stars are
much more dominant, as indicated by their values of $P_{\rm spatial}$. We can
therefore use the ratio between $P_{\rm spatial}$ and the field abundance of
young stars to correct the relative abundances of young and old stars, re-
calibrating $P_{\rm Age<50Myr}$ for the environment in which each star is
found.

We start by defining a corrective factor $C_{y}$ by which the abundance of
young stars in the original model must be inflated to recreate the
environments in Cep-Her where far more young members are present:

$C_{y}=\frac{P_{\rm spatial}}{P_{i}}$ (1)

where $P_{i}$ is the fraction of young stars in the original SPYGLASS model,
or $P_{i}=50/11200\simeq 0.004$. This increase to the fraction of young stars
in our model in turn reduces the fraction $C_{o}$ of older model stars by the
following factor:

$C_{o}=\frac{1-P_{\rm spatial}}{1-P_{i}}$ (2)

which is equal to the local spatial probability that the star is a non-member
divided by the initial model probability that the star is an old non-member.
Using these corrections to the fractions of old and young stars in the model,
we can compute a renormalized spatial-photometric membership probability
$P_{\rm sp}$ which accounts for both the photometric membership probability
and local abundance of young members:

$P_{\rm sp}=\frac{(C_{y}\times P_{\rm Age<50Myr})}{(C_{y}\times P_{\rm
Age<50Myr})+(C_{o}\times(1-P_{\rm Age<50Myr}))}$ (3)

Before calculating $P_{\rm sp}$, we set the maximum value of $P_{\rm spatial}$
to 0.95 to ensure that stars with values of $P_{\rm spatial}$ that are close
to 1 can still be plausibly photometrically removed from the sample. This cap
also reflects results from when values of $P_{\rm spatial}$ were initially
computed. In SPYGLASS-IV, $P_{\rm spatial}$ was generated by computing the
relative abundances of young and old stars across bins defined in clustering
proximity, $D$, a proxy for local density. A Gaussian CDF was fit to the curve
of $D$ vs $P_{\rm spatial}$ computed for that range of $D$, effectively
smoothing the scatter in those bins. However, even for very high values of
$D$, some old contaminants were found to still be present, albeit rare.
Setting the maximum $P_{\rm spatial}$ to 0.95 approximately captures the
maximum scatter in the computed values of $P_{\rm spatial}$ after the fit
results for $P_{\rm spatial}$ exceed 0.99. The resulting values of $P_{\rm
sp}$ provide the probability of youth conditioned on the local abundance of
young members. We provide these results for all stars that pass the
astrometric and photometric cuts that enable the calculation of $P_{\rm
Age<50Myr}$ in SPYGLASS-IV.

###  3.2 Spectroscopic Membership

The new and existing spectroscopic measurements we have collected provide two
additional ways of assessing membership: lithium depletion and radial
velocity. Conclusions of membership or non-membership from velocity or lithium
can therefore dramatically raise or lower the probability of membership for a
spectroscopically observed star relative to its initial spatial-photometric
membership probability. In this section we assess membership using these
spectroscopic markers, taking into account the influence of field
contamination and false negatives that can affect the results from these
methods. We then compute final membership probabilities, which take into
account these spectroscopic membership indicators.

####  3.2.1 RV Selection

Stars in our sample were identified as candidate members using the ubiquitous
on-sky proper motion measurements from Gaia. However, RVs are much less
widespread, limiting their use in dynamical studies, especially in the almost
unstudied Cep-Her Complex. These RVs complete the 3D stellar velocity vector,
establishing whether a star is fully co-moving with the parent association and
greatly strengthening membership assessment.

In Cep-Her, the addition of this new dimension in velocity space reveals a
distinct bimodality in Cep-Her’s RV distribution. In Figure 2, we show the
distribution of RVs in Cep-Her as a function of galactic longitude ($l$). This
figure reveals two distinct radial velocity sequences separated by more than
10 km/s on average. This RV split is present even among populations that
overlap substantially in both space and transverse velocity coordinates. In
limited cases, we find that the space-transverse velocity overlap can be so
great that HDBSCAN clustering in 5D phase-space, as described in Section 4,
can preserve this bimodality within individual subgroups. Therefore, not only
is it necessary to perform the RV cuts on the two samples separately, but the
two populations must be separated entirely to provide internally consistent
clustering results in five dimensions.

We separate the two sequences with a manual cut consisting of two line
segments, as shown in Figure 2. We refer to the larger population which
resides below this line as Cep-Her Main, and stars above this line as Cep-Her
North. The overdensity associated with Cep-Her North occupies a narrower range
in $l$, so we also exclude stars with $l<70$ from that clump. We do not apply
this analysis to stars with poor RV measurements (defined as $\sigma_{RV}>5$
km s-1), so that stars with unreliable radial velocity measurements are not
grouped into one of those sequences according to those velocities. The stars
in these two clumps define regions in 5D space-transverse velocity coordinates
that can be used to group stars without reliable RVs into a parent subregion.

Figure 2: Cep-Her candidate members in galactic $l$ sky coordinates, plotted
against RVs from new observations and the literature. The result shows two
distinct radial velocity sequences. The black lines represent the manual cuts
used to separate Cep-Her Main (red) from Cep-Her North (blue) and the field
(lighter shades of red and blue). Probable field stars are assigned to one of
the two Cep-Her sub-regions according to their position in 5D space-transverse
velocity coordinates for reassessment in the final RV selections (see Fig 3)
and for use in assessments of contamination 3.2.2. The histogram at the right
shows the distribution of RVs for stars with $l>70$, showing the bimodal
distribution of RVs over this range in $l$ where both RV sequences are
present.

We assign stars with poor RVs or no RVs to whichever population is closest in
5D space-transverse velocity coordinates. We used the clustering proximity
$D_{N}$, which is defined as the distance to the Nth nearest member, as the
metric to assess closeness. This HDBSCAN-inspired metric has been used in past
SPYGLASS publications, where $D_{N}$ serves as a proxy for the inverse cube-
root of the density surrounding a given star. We defined this distance in 5D
$(X,Y,Z,c*\Delta$v${}_{T,l},c*\Delta$v${}_{T,b})$ space-transverse velocity
anomaly coordinates, where the transverse velocity anomaly $\Delta$vT is
defined as the transverse velocity for the star, minus the expected transverse
velocity at its location given a median UVW 3D velocity vector. The constant
$c$ matches the scales of the space and velocity dimensions, and its value is
chosen to produce similar structure sizes in space and velocity coordinates.
Here we select $c=12$ pc km-1 s, following SPYGLASS-II and -III. The use
transverse velocity anomaly instead of transverse velocity removes some
projection effects, making populations with different velocities and widely-
distributed on-sky locations easier to separate. The value of $N$ in $D_{N}$
serves as a smoothing factor that affects sensitivity to small-scale
structures, and we find that a value of $N=20$ works well to separate these
large-scale clumps. This value of $N$ provides more smoothing compared to the
usage in SPYGLASS-IV, where $N=10$ was employed to assign outlying members to
a range of associations that include smaller-scale features.

We assign each star without a reliable radial velocity (no RV or
$\sigma_{RV}>5$ km s-1) to the group with which it has the lowest $D_{20}$,
measuring $\Delta$vT relative to the median UVW velocity across all of Cep-
Her. We apply the same approach to stars with RVs that do not fall in either
clump to both reassess their membership with more robust group extents, and to
assess field contamination in Section 3.2.2. RV-separated groups that are
otherwise contiguous in 5D space-velocity coordinates will experience some
cross-contamination through this approach, and we discuss potential impacts on
population demographics in Section 5.5.

Figure 3: Results of RV selection across a cross-section of space and velocity
axes in Cep-Her Main (top) and Cep-Her North (bottom). Dark blue diamonds are
accepted as kinematic members, and light blue dots are not. For the plots of
velocity against galactic $l$ coordinates, we use a red curve to show the
median values of U, V, and W as a function of $l$, which were used to define
the velocity loci for the two subtly discontiguous sub-components of Cep-Her.
The shaded red region encloses the maximum and minimum velocities which may be
identified as members in that axis, as defined by the velocity axis, plus or
minus the maximum velocity dispersion.

Once all stars are assigned to a parent population, we make a final 3D
kinematic membership assessment for stars with radial velocities. Both RV-
defined sub-populations in Cep-Her contain two visibly discontiguous
subregions each that emerge in $l$-transverse velocity space (Fig. 3), so we
draw a manual cut between them and consider their velocity distributions
independently. We split each of the four components across Cep-Her Main and
Cep-Her North into between 10 and 25 bins each in $l$, and in each bin we
compute median values of transverse velocity and 3D UVW velocity. We also
compute the maximum variance $\sigma_{v_{max}}$ in transverse velocity
relative to the median, which informs the velocity spread in the plane of the
sky. In regions of parameter space where Cep-Her Main and Cep-Her North
overlap, we use the 95th percentile in transverse velocity rather than the
maximum variance to reduce cross-contamination between the two populations.
The number of bins is chosen by hand to both capture apparent structure and
minimize scatter in the velocity results caused by small sample sizes. We
smooth the results with a Savitsky-Golay filter in $l$ to further reduce the
variations from bin to bin. Taking the transverse velocity extent to be a 2D
projection of a spherical 3D distribution in UVW space, we assume that the
maximum variance in 2D space is also the maximum variance in 3D space. We
therefore identify stars within $\sigma_{v_{max}}$ of the UVW median as
members with a velocity membership flag $V=1$, and identify stars outside this
cut as velocity non-members ($V=-1$). We use an inconclusive membership flag
$V=0$ for stars with $\sigma_{RV}>\sigma_{v_{max}}/2$, as their large
uncertainties provide no strong membership constraints. We show the results of
this selection in Figure 3.

####  3.2.2 RV False Positives from Field Contamination

Much of Cep-Her and its adjacent associations lie close to the galactic
midplane, making them subject to substantial field contamination. As a result,
many field stars have velocities matching those of Cep-Her, resulting in their
incorrect classification as RV members. Pre-main sequence stars and O and B
stars can be reliably identified as members based on photometry alone, however
stars between the pre-main sequence turn-on (PMSTO) and the tip of the main
sequence have no photometric means of membership assessment. Furthermore, many
of the stars in the same range lack conclusive measurements of the Lithium
line, which is a common alternative means of assessing membership. RV
measurements are necessary to assess the membership of these stars, and it is
therefore important to quantify the probability of field contamination. We
account for this possibility by estimating a false positive rate for our RV
membership assessment, and factoring that into our membership probabilities.

To assess this contamination, we must compute a false positive probability
$P_{\rm fp}$, which we define as the conditional probability of RV membership,
given that the object is actually a non-member field star ($P_{\rm
fp}=P(V=1|\textrm{nm})$, where “nm” limits the calculation to non-members
only). In Section 3.2.5, we compute final probabilities of membership
conditioned on the $V$ membership flags, $P_{\rm fin}(\textrm{mem}|V=x)$, so
$P_{\rm fp}$ becomes an important input into Bayes Theorem. To produce these
probabilities, we separately fit the distributions of the field star RVs in
the samples for Cep-Her Main and Cep-Her North with a gaussian, including only
stars with $\sigma_{RV}<5$ km s-1 to remove poor RV measurements, and
excluding stars with $-30<v_{R}<0$ km s-1 to limit the fit to the RV range
where the field dominates. These fits are shown in Figure 4. We can then
compute the false positive rate for a field star as the area under the fit
across the range of RVs where stars are marked as RV members, divided by the
total area under the fit.

The range of RVs where stars are marked as members is proportional to the
radius $\sigma_{v_{max}}$ used to flag velocity membership for stars at a
given position in $l$. However, since we identify stars as RV members using a
spherical cut in velocity space, stars with outlying transverse velocities
have stricter cuts on RV than stars closer to the transverse velocity average.
As a result, not all stars within the radius $\sigma_{v_{max}}$ of the mean RV
are marked as members. We must therefore compute an effective range in RV for
computing $P_{\rm fp}$. The field velocity distribution should be nearly
isotropic over the range where stars are identified as velocity members, so
any shape in 3D velocity space with the same volume should contain
approximately the same number of stars. The initial velocity selection of
potential members in SPYGLASS-IV was performed in transverse velocity space
with an effective 2D search radius of $\sigma_{v_{max}}$. Restricting that
population in the RV dimension turns this 2D bounded circular area into a 3D
cylindrical bounded volume, so we can define the effective RV range
$\Delta_{RV}$ as the length of a cylinder of radius $\sigma_{v_{max}}$ and
length $\Delta_{RV}$ that has the same volume as a sphere of radius
$\sigma_{v_{max}}$. This corresponds to
$\Delta_{RV}=\frac{4}{3}\sigma_{v_{max}}$, and we use this value as the width
of the window where field stars would be identified as members, centered on
the mean RV. In ranges of galactic longitude $l$ where two velocity sequences
are present (see Fig. 3), the velocity range is taken to contain the maximum
extent across both sequences.

Figure 4: Fits to the field contamination in Cep-Her Main and Cep-Her North.
The purple histogram shows the distribution of RVs, while the red curve fits
the sections of the distribution with ($v_{R}<-30$ and $v_{R}>0$) where the
field dominates. We apply dark shading to the Cep-Her-dominated range of RVs
that we do not include in these fits.

The false positive rate for a star with a given value of $l$,
$\sigma_{v_{max}}$, and mean RV is then the integral over the fit to the field
RVs across a window of width $\Delta_{RV}$ centered on the mean RV, divided by
the integral across the entire fit to the field RVs. The resulting values for
$P_{\rm fp}$ range from 12% to 22% in Cep-Her Main and 17% to 21% in Cep-Her
North. Stars on the main sequence, where we largely rely on velocities for
membership, can have values of $P_{\rm sp}\lesssim 0.2$, consistent with a
field population that is dominant over Cep-Her members at a rate of 4:1 or
more. With up to 21% of field stars being assigned as velocity members in
certain locations, this implies that even among stars with $V=1$, outlying
sections of Cep-Her may still host members and field stars at similar rates.
It is therefore important that we do not treat velocity membership as
conclusive, and consider the field contamination when assessing membership.

####  3.2.3 RV False Negatives

Due to a combination of poor RV measurements and apparent motions from the
influence of a binary companion, our RV membership cut may occasionally mark
real members as non-members. We account for this possibility by estimating a
false negative rate $P_{\rm fn}$ for our RV membership assessment, which we
define as the rate at which true members are marked as velocity non-members
($P_{\rm fn}=P(V=-1|\textrm{mem})$). Like for our definition of the false
positive rate in Section 3.2.2, this form is chosen so that it can be directly
inputted into our final membership probability calculations in Section 3.2.5,
$P_{\rm fin}(\textrm{mem}|V=x)$, which are conditioned on the $V$ membership
flags.

Since this probability is conditional on membership, we must estimate it over
an essentially pure sample of members. We first select a sample of stars with
$P_{\rm spatial}=0.95$, which is the maximum value we allow $P_{\rm spatial}$
to assume. In this set, 6.2% of stars with conclusive velocity membership
$V\neq 0$ are classified as velocity non-members. Some of these are likely to
be genuine non-members, as a very weak field main sequence is visible even
among this set. However, of the 5 stars in the section of the CMD where the
pre-main sequence and field main sequence separate, 4 are closer to the pre-
main sequence. This suggests that most of the stars in this set are true false
negatives. We also investigate a set of stars with a spatial-photometric
$P_{\rm sp}>0.99$, in which 18% of stars with $V\neq 0$ are classified as
velocity non-members. This is a much higher value, likely caused by the fact
that it is based almost entirely on a set of photometrically young O and B
stars, which have the highest binarity rates and rapid rotation (Sullivan &
Kraus, 2021), both of which significantly weaken the accuracy of our velocity
membership assessment. Due to the biased selection in the sample chosen by
$P_{\rm sp}$, we take the sample with $P_{\rm spatial}=0.95$ to be a closer
representation of the overall RV false membership statistics, and adopt a
false negative rate of $P_{\rm fn}=0.05$. This roughly reflects the 6.2% value
computed from that sample, with a 1% correction to account for potential field
contamination in that otherwise pure sample. This result is also broadly
consistent with the false negative rates reported in previous studies of young
associations, such as those provided for $\beta$ Pic in Shkolnik et al.
(2017).

####  3.2.4 Lithium Selection

Figure 5: Li selection in both Cep-Her sub-components, showing stars that are
identified as Li members as green stars, Li non-members as red diamonds, and
stars with inconclusive Li membership as black circles. We also include 20
Myr, 30 Myr, 50 Myr, 100 Myr, and 5 Gyr Lithium sequences for reference.

Lithium abundances are a strong indicator of stellar youth, as Lithium
depletes rapidly after star formation for stars less massive than the Sun
(Bodenheimer, 1965). We therefore use the EAGLES package (Jeffries et al.,
2023; Jeffries, 2023) to assess membership for stars with a Lithium equivalent
width (EW) measurement. For a given Li EW and $T_{eff}$, EAGLES computes age
posteriors, which can be converted into a probability of youth for each star.
Jeffries et al. (2023) also includes a formula for converting from
$G_{BP}-G_{RP}$ to $T_{eff}$, making results attainable directly from Gaia
photometry. As SPYGLASS-IV uses the PARSEC 80 Myr isochrone as the maximum age
for inclusion in the SPYGLASS catalog of young associations (Chen et al.,
2015), we identify non-members by requiring that $P_{age<80Myr}<0.005$
according to the EAGLES age posteriors. This strict cut for excluding non-
members reduces the probability of false negatives due to to photometric
contamination from binary companions, which can reduce Li EWs. To further
limit false negatives, we also require that stars that receive a negative
membership assessment also have a value of Gaia’s Renormalized Unit Weight
Error RUWE$>1.2$. This is a common condition used to exclude unresolved
binaries (Bryson et al., 2020), which we discuss further in Section 5.3.1. We
then classify all stars with a lithium age $\tau_{mean}<80$ Myr as Li-verified
members. Unlike for identifying non-members, we do not apply an RUWE
restriction to our identification of members, as any binary contaminant to a
Li-rich star would almost certainly lower the measured Li EW, not increase it.
All other stars receive an ambiguous Lithium membership assessment.

####  3.2.5 Aggregate Membership Probabilities

To produce a final membership probability $P_{\rm fin}$, we must combine our
initial spatial-photometric membership probability estimates $P_{\rm sp}$ with
our spectroscopic membership flags from lithium and velocity. For stars that
lack conclusive spectroscopic membership flags from velocity or Lithium, we
set $P_{\rm fin}$ equal to $P_{\rm sp}$, which takes into account all of our
non-spectroscopic means of membership assessment. Due to the well-established
reliability of Lithium as an indicator of youth (e.g., Binks & Jeffries,
2014), we take Li membership to be conclusive, and set the $P_{\rm fin}$ for
Li members to 1, and $P_{\rm sp}$ for Li non-members to 0. However, as
discussed in Sections 3.2.2 and 3.2.3, velocity false negatives and false
positives are both possible, and it is therefore not appropriate to simply set
the probability to 1 or 0 based on a velocity membership flag.

Figure 6: Color magnitude diagram of Cep-Her Main (left), and Cep-Her North
(right), colored by membership results. Orange diamonds indicate objects with
$P_{\rm fin}>0.5$ that are included in our probable member sample, while blue
squares indicate objects that are not in that sample, but pass the photometric
and astrometric quality cuts. Black triangles fail at least one of the
astrometric and photometric quality cuts.

We can compute the membership probability $P_{\rm fin}$ for a star with a good
RV ($V\neq 0$) as the posterior probability $P(\textrm{mem}|V=x)$, or the
probability of membership given a certain velocity membership flag. We adopt
the spatial-photometric membership probability $P_{\rm sp}$ as the prior
membership probability. Then, for $V=1$, the likelihood is provided by
$P(V=1|\textrm{mem})=1-P_{\rm fn}$, which follows from that fact that $P_{\rm
fn}=P(V=-1|\textrm{mem})$, as described in Section 3.2.2. Bayes theorem for
$P(\textrm{mem}|V=1)$ can then be written as:

$P_{\rm fin}(\textrm{mem}|V=1)=\frac{(1-P_{\rm fn})P_{\rm sp}}{P(V=1)}$

and where the marginal probability $P(V=1)$ can be written as:

$P(V=1)=(1-P_{\rm fn})P_{\rm sp}+P_{\rm fp}(1-P_{\rm sp})$

In the second term, we use the fact that $P_{\rm fp}=P(V=1|\textrm{nm})$, as
defined in Section 3.2.2.

For stars with a $V=-1$, we can apply an analogous calculation for assessing
membership probability given velocity non-membership, which produces the
following:

$P_{\rm fin}(\textrm{mem}|V=-1)=\frac{P_{\rm fn}P_{\rm sp}}{P_{\rm fn}P_{\rm
sp}+(1-P_{\rm fp})(1-P_{\rm sp})}$

This formula also follows Bayes theorem, this time using $P_{\rm
fn}=P(V=-1|\textrm{mem})$ and $1-P_{\rm fp}=P(V=-1|\textrm{nm})$.

For each star with a conclusive velocity membership flag ($V=1$ or $V=-1$), we
set the aggregate probability $P_{\rm fin}$ equal to the output of the
corresponding $P_{\rm fin}$ formula above for its value of $V$. For $V=0$, we
adopt the original value of $P_{\rm sp}$. In cases with a definitive Lithium
flag, we prioritize the lithium result, which returns a $P_{\rm fin}$ value of
either 1 or 0. These membership adjustments for the velocity and lithium flags
complete our results for $P_{\rm fin}$.

In both Cep-Her North and Cep-Her Main, the sum over $P_{\rm fin}$ is within
5% of the sum over $P_{\rm sp}$, demonstrating that these results do not
deviate significantly from our probability prior of $P_{\rm sp}$ at a
population level. The main impact of including false positive and false
negative rates is to refine membership in cases where contamination and false
negatives are especially common. F and G stars rely almost entirely on
velocities to assess membership, so stars in this regime that have a low
$P_{\rm sp}$ prior may still have relatively ambiguous values of $P_{\rm fin}$
even after being assigned $V=1$, as described in Section 3.2.2. On the O and B
sequence, where $P_{\rm sp}>0.9$ is common, the non-zero value of $P_{\rm fn}$
can become significant, as $P(mem|V=-1)=0.54$ for $P_{\rm sp}=0.95$, $P_{\rm
fp}=0.2$ and $P_{\rm fn}=0.05$. These corrections therefore condition our
membership probabilities for velocity members on the relative dominance of
members and non-members in its section of parameter space, which is a function
of the spatial-photometric membership probability, velocity false negative
rate, and false positive rate from field contamination.

We use the final values of $P_{\rm fin}$ to select a purer probable member
sample that can be used for identifying substructure. We define this sample by
removing all stars that fail the astrometric and photometric quality cuts (and
therefore lack $P_{\rm fin}$), and require $P_{\rm fin}>0.5$. This $P_{\rm
fin}$ cut balances completeness and purity to produce robust substructure-
defining samples for both Cep-Her Main and Cep-Her North. In Cep-Her Main,
7228 stars pass these conditions and the sum of $P_{\rm fin}$ across those
stars is 6268, so we expect that 87% of objects in that sample are members and
13% are field stars. The sum over all $P_{\rm fin}$ in Cep-Her Main is 7351,
so these restrictions accept a 15% loss rate of members. The same calculations
as above predict similar contamination and loss rates in Cep-Her North, where
the probable member sample is expected to contain 89% members and 11% field
stars with an 11% loss rate. Our probable member samples therefore contain the
vast majority of association members with good photometry and astrometry,
while having very limited field contamination.

We summarize the membership selection in Figure 6, showing the distribution of
probable members, probable non-members, and stars that fail our astrometric
and photometric quality cuts and therefore lack a result for $P_{\rm fin}$. We
also compile all membership flags and probabilities for all candidate Cep-Her
members in Table 2, alongside some basic properties for each star.

Table 2: List of Cep-Her members, including basic IDs and properties, substructure assignments, membership probabilities and other flags. The complete version of this table is available in the online-only version, which contains 34886 members Gaia DR3 ID | ASaaAS, SG, AA, and AG indicate membership in parent structures: AS for association-level assignment, SG for the subgroup, AG for alternate subgroup assignment in the case of membership in the opposite top-level group, and AA for the association that contains the assigned alternate subgroup. -1 indicates that it does not have an assigned subgroup. | SGaaAS, SG, AA, and AG indicate membership in parent structures: AS for association-level assignment, SG for the subgroup, AG for alternate subgroup assignment in the case of membership in the opposite top-level group, and AA for the association that contains the assigned alternate subgroup. -1 indicates that it does not have an assigned subgroup. | AAaaAS, SG, AA, and AG indicate membership in parent structures: AS for association-level assignment, SG for the subgroup, AG for alternate subgroup assignment in the case of membership in the opposite top-level group, and AA for the association that contains the assigned alternate subgroup. -1 indicates that it does not have an assigned subgroup. | AGaaAS, SG, AA, and AG indicate membership in parent structures: AS for association-level assignment, SG for the subgroup, AG for alternate subgroup assignment in the case of membership in the opposite top-level group, and AA for the association that contains the assigned alternate subgroup. -1 indicates that it does not have an assigned subgroup. | RA | Dec | d | g | M |  $V$bbVelocity membership assessment. 1 is a kinematic member, 0 is inconclusive, and -1 is a non-member | LiccLithium membership assessment. 1 is a kinematic member, 0 is inconclusive, and -1 is a non-member | FddGeneral flag: 1 indicates that the star has a resolved companion within 10,000 au in the plane of the sky, 2 indicates a bad broadening function solution, 4 indicates a bimodal line profile likely indicative of spectroscopic binarity, 8 indicates an RUWE$>$1.2, indicating likely unresolved binarity, and 16 indicates that the RV recorded was ambiguously attributed to two components of a binary pair. The flags are added in cases where multiple are true; for example, flag 6 indicates both flags 2 and 4. | $P_{\rm Age<50Myr}$ | $P_{\rm spatial}$ | $P_{\rm sp}$ | $P_{\rm fin}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
|  |  |  |  | (deg) | (deg) | (pc) |  | M⊙ |  |  |  |  |  |  |
2125758539488457344 | CUPV | -1 | CINR | 1 | 292.3272 | 42.6479 | 314 | 15.91 | 0.59 | 0 | 0 | 0 | 0.0095 | 0.35 | 0.54 | 0.54
2125766751463294976 | CUPV | -1 |  |  | 292.2994 | 42.8453 | 352 | 20.26 | 0.18 | 0 | 0 | 0 |  | 0.23 |  |
2125769294084167936 | CUPV | -1 |  |  | 292.3366 | 42.9828 | 334 | 19.82 | 0.18 | 0 | 0 | 0 | 0.0000 | 0.24 | 0.00 | 0.00
2125781831091630208 | CUPV | 7 | CINR | 1 | 291.9657 | 42.6111 | 286 | 19.42 | 0.25 | 0 | 0 | 0 | 0.0000 | 0.23 | 0.00 | 0.00
2125788428159994368 | ORPH | -1 |  |  | 291.5863 | 42.5099 | 309 | 11.76 | 1.09 | 1 | 0 | 8 | 0.0033 | 0.27 | 0.21 | 0.61
2125805779829542784 | ORPH | -1 |  |  | 292.0603 | 42.6919 | 297 | 17.83 | 0.33 | 0 | 0 | 0 | 0.0014 | 0.22 | 0.08 | 0.08
2125812381199146112 | ORPH | -1 |  |  | 291.9605 | 42.9027 | 370 | 18.98 | 0.09 | 0 | 0 | 1 |  | 0.24 |  |
2125842033652058368 | ORPH | -1 |  |  | 292.8799 | 42.7379 | 315 | 18.12 | 0.30 | 0 | 0 | 0 | 0.0024 | 0.45 | 0.30 | 0.30
2125856361663343232 | CUPV | 9 | CINR | 1 | 292.7672 | 43.0986 | 345 | 15.76 | 0.60 | 0 | 0 | 0 | 0.0012 | 0.22 | 0.07 | 0.07
2125861717484155008 | ORPH | -1 |  |  | 292.5439 | 43.0932 | 301 | 20.52 | 0.16 | 0 | 0 | 8 |  | 0.26 |  |
2125864161317485824 | ORPH | -1 |  |  | 292.4818 | 43.0602 | 326 | 20.66 | 0.16 | 0 | 0 | 0 |  | 0.27 |  |
2125887663380213504 | ORPH | -1 |  |  | 293.1205 | 43.4306 | 279 | 18.94 | 0.25 | 0 | 0 | 8 |  | 0.21 |  |

##  4 Structure

###  4.1 Clustering

Cep-Her is a large and complicated association that contains diverse
substructure, ranging from dense star clusters to diffuse associations that
blend with both the field and adjacent populations. These populations provide
basic units of co-natal star formation within the broader population that can
trace the structure of the parent cloud, while providing suitable populations
for measuring accurate local ages.

We use our probable member sample described in Section 3.2.5 to identify
subgroups in 5-dimensional $(X,Y,Z,c*\Delta$v${}_{T,l},c*\Delta$v${}_{T,b})$
space-transverse velocity anomaly coordinates, with the constant $c$ set to
$c=6$ pc km-1 s, following the choices for initial group identification in
SPYGLASS-II. For both Cep-Her Main and Cep-Her North, we cluster this 5D
parameter space using the HDBSCAN clustering algorithm with the clustering
parameters min$\\_$samples and min$\\_$cluster$\\_$size set to 8. This
provides clustering that is slightly less sensitive compared to the selection
in SPYGLASS-II, which sets these parameters to 7. This choice reflects the
greater risk of false groups introduced by our less sensitive astrometry at
the distance of Cep-Her, which is approximately twice as distant as Cepheus
Far North (CFN; Kerr et al., 2022b). We perform this clustering in leaf mode,
which identifies the smallest scale populations in HDBSCAN’s clustering tree.
EOM clustering, which identifies the populations that persist over the largest
range of clustering scales, was used in SPYGLASS-I to hint at possible
connections between subgroups, although SPYGLASS-II revealed that true
primordial star-forming sites do not necessarily reflect EOM groups. We
therefore consider connections between leaf subgroups manually.

HDBSCAN is optimised for identifying subtle overdensities within a complex
background, and therefore most objects in a sample are assigned to the
background rather than any group. However, since we are clustering on samples
of candidate members which exclude most contaminants, HDBSCAN’s conservative
group assignments mark many outlying subgroup members as part of the
background, and stars with $P_{\rm fin}<0.5$ are excluded entirely. We can
therefore produce more complete samples for each subgroup by assigning
outlying members and stars that are not in our probable member sample to the
association that agrees most closely in space-velocity coordinates. To assign
these outlying members to a parent subgroup, we compute the clustering
proximity $D_{N}$ with $N=8$ for each star in the outlying sample relative to
the space-velocity distributions of each HDBSCAN cluster in 5D
$(X,Y,Z,c*\Delta$v${}_{T,l},c*\Delta$v${}_{T,b})$ space. Here we use $c=12$ pc
km-1 s, following SPYGLASS-II’s selection for assigning outlying members to a
parent subgroup. Stars are assigned to a given group if the clustering
proximity of that group is 15% lower than that of any other group. This
assigns most stars to a subgroup, while leaving stars that are not closely
connected to any one group as unclustered Cep-Her members. The 15% clustering
proximity cut allows most stars to be assigned to a subgroup, while preventing
cases where substantial portions of Cep-Her’s extended halo are assigned to
outlying subgroups, which can give those groups unphysical shapes. The final
membership assignments are listed in Table 2.

We compile information on the groups identified by HDBSCAN in Table 3. We
identify 50 subgroups in total, including 38 in Cep-Her Main and 12 in Cep-Her
North. Half of the 34886 Cep-Her Candidate members are assigned to a subgroup,
and stars in a subgroup account for 66% of the sum over $P_{\rm fin}$,
indicating that probable members are much more likely to be assigned a
subgroup compared to probable field stars. Before reintroducing outlying stars
and those outside of the probable membership sample ($P_{\rm fin}<0.5$), only
2135 stars were included in these subgroups, demonstrating the importance of
broadening these groups’ membership beyond the initial HDBSCAN selection.

Table 3: Demographics and mean properties of the subgroups in Cep-Her, identified by their parent association and ID. ASSOC | ID | N | M |  $\sigma_{dem}$aa$\sigma_{dem}$ is the fractional uncertainty in both the number of members $N$ and the mass $M$. | RA | Dec | l | b | d | $\overline{v_{T,l}}$ | $\overline{v_{T,b}}$ | $\overline{v_{r}}$ | $R_{hm}$ | $\sigma_{1D}$ | $\sigma_{vir}$ | Vir. Rat | Age
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
|  |  | (M⊙) |  | (deg) | (deg) | (pc) | (km s-1) | (pc) | (km s-1) |  | (Myr)
CUPV | 1 | 129 | 50.4 | 0.379 | 292.6 | 31.0 | 64.7 | 6.1 | 347.0 | -6.8 | -5.4 | -18.2 | 14.1 | 0.96 | 0.06 | 12.2 $\pm$ 6.1 | 68.2 $\pm$ 9.2
CUPV | 2 | 187 | 68.6 | 0.364 | 291.4 | 30.5 | 63.8 | 6.8 | 309.0 | -5.6 | -2.5 | -19.4 | 17.5 | 0.52 | 0.06 | 6.3 $\pm$ 5.2 | 62.6 $\pm$ 11.8
ORPH | 1 | 185 | 76.3 | 0.100 | 286.5 | 19.0 | 51.5 | 5.5 | 359.0 | -8.9 | -5.2 | -14.3 | 11.8 | 0.94 | 0.07 | 8.9 $\pm$ 2.6 | 27.9 $\pm$ 2.4
ORPH | 2 | 88 | 35.2 | 0.100 | 284.8 | 18.7 | 50.5 | 6.8 | 361.0 | -10.3 | -8.9 | -14.2 | 11.9 | 0.59 | 0.05 | 8.3 $\pm$ 5.4 | 26.5 $\pm$ 2.8
ORPH | 3 | 272 | 102.9 | 0.140 | 295.5 | 48.8 | 81.8 | 12.6 | 294.0 | 3.5 | -4.7 | -21.9 | 22.4 | 1.74 | 0.07 | 18.9 $\pm$ 5.3 | 39.0 $\pm$ 1.8
CUPV | 3 | 1435 | 561.3 | 0.224 | 320.4 | 51.2 | 93.1 | 0.9 | 362.0 | 12.1 | -5.6 | -19.3 | 3.0 | 0.31 | 0.42 | 0.5 $\pm$ 0.2 | 55.3 $\pm$ 2.3
CUPV | 4 | 401 | 137.0 | 0.219 | 287.0 | 24.9 | 57.1 | 7.7 | 332.0 | -5.8 | -3.5 | -14.6 | 15.2 | 1.23 | 0.09 | 9.9 $\pm$ 3.0 | 68.4 $\pm$ 8.8
CUPV | 5 | 140 | 51.3 | 0.257 | 315.9 | 46.4 | 87.6 | -0.2 | 295.0 | 11.6 | -6.3 | -18.4 | 12.8 | 0.48 | 0.06 | 5.8 $\pm$ 2.2 | 71.0 $\pm$ 6.5
CUPV | 6 | 371 | 158.1 | 0.291 | 315.2 | 54.1 | 93.2 | 5.2 | 307.0 | 13.1 | -4.0 | -19.5 | 17.3 | 0.75 | 0.09 | 5.8 $\pm$ 2.1 | 68.1 $\pm$ 6.7
CUPV | 7 | 391 | 143.7 | 0.325 | 302.9 | 44.3 | 80.5 | 5.7 | 295.0 | 6.1 | -2.6 | -20.4 | 16.5 | 0.55 | 0.09 | 4.6 $\pm$ 2.2 |  $>$ 80.0
ORPH | 4 | 180 | 73.7 | 0.100 | 276.9 | 15.4 | 44.2 | 12.1 | 363.0 | -13.2 | -2.8 | -14.8 | 13.0 | 0.32 | 0.07 | 3.3 $\pm$ 1.5 | 34.5 $\pm$ 2.5
ORPH | 5 | 155 | 64.7 | 0.100 | 281.1 | 22.6 | 52.6 | 11.6 | 375.0 | -10.5 | -2.8 | -15.3 | 11.8 | 0.52 | 0.07 | 5.3 $\pm$ 1.9 | 30.0 $\pm$ 1.9
CUPV | 8 | 887 | 397.2 | 0.258 | 309.3 | 43.3 | 82.3 | 1.4 | 341.0 | 8.1 | -6.6 | -20.1 | 17.1 | 0.69 | 0.15 | 3.4 $\pm$ 1.0 | 55.5 $\pm$ 4.1
CUPV | 9 | 892 | 360.9 | 0.243 | 300.9 | 39.3 | 75.5 | 4.3 | 354.0 | 3.3 | -4.1 | -20.7 | 16.5 | 0.78 | 0.14 | 3.9 $\pm$ 1.2 | 53.8 $\pm$ 4.7
CUPV | 10 | 559 | 250.1 | 0.242 | 303.4 | 34.4 | 72.4 | 0.0 | 325.0 | 4.3 | -6.7 | -18.5 | 7.1 | 0.21 | 0.18 | 0.8 $\pm$ 0.3 | 58.4 $\pm$ 3.0
ORPH | 6 | 109 | 65.2 | 0.100 | 272.1 | 18.0 | 44.7 | 17.4 | 344.0 | -12.1 | -0.7 | -15.8 | 12.1 | 0.53 | 0.07 | 5.5 $\pm$ 2.9 | 33.7 $\pm$ 3.8
ORPH | 7 | 164 | 65.5 | 0.100 | 276.5 | 30.1 | 58.1 | 18.3 | 307.0 | -10.2 | -1.7 | -20.2 | 14.2 | 0.49 | 0.06 | 5.5 $\pm$ 2.8 | 39.8 $\pm$ 3.3
ORPH | 8 | 197 | 67.2 | 0.100 | 288.6 | 32.4 | 64.5 | 9.8 | 324.0 | -3.5 | -6.4 | -16.8 | 17.1 | 0.54 | 0.06 | 6.5 $\pm$ 2.8 | 27.1 $\pm$ 1.3
ORPH | 9 | 229 | 107.0 | 0.100 | 275.9 | 27.2 | 55.0 | 17.7 | 364.0 | -11.2 | -0.8 | -18.1 | 10.2 | 0.80 | 0.09 | 5.9 $\pm$ 1.6 | 32.8 $\pm$ 2.8
ORPH | 10 | 199 | 70.7 | 0.100 | 274.8 | 23.6 | 51.1 | 17.3 | 386.0 | -9.6 | -4.3 | -18.4 | 15.4 | 0.57 | 0.06 | 6.4 $\pm$ 2.1 | 24.9 $\pm$ 1.9
ORPH | 11 | 268 | 110.6 | 0.100 | 290.8 | 35.5 | 68.1 | 9.5 | 372.0 | -2.5 | -4.7 | -19.0 | 13.7 | 0.39 | 0.08 | 3.3 $\pm$ 1.5 | 32.4 $\pm$ 2.0
ORPH | 12 | 151 | 73.9 | 0.100 | 276.4 | 20.9 | 49.1 | 14.9 | 304.0 | -10.5 | -2.9 | -17.3 | 8.7 | 0.22 | 0.09 | 1.8 $\pm$ 1.0 | 27.7 $\pm$ 1.6
ORPH | 13 | 143 | 55.4 | 0.100 | 275.4 | 20.0 | 47.9 | 15.4 | 345.0 | -11.7 | -5.4 | -16.8 | 9.8 | 0.38 | 0.07 | 3.8 $\pm$ 1.3 | 26.5 $\pm$ 3.8
ORPH | 14 | 108 | 34.5 | 0.100 | 275.7 | 26.6 | 54.3 | 17.7 | 323.0 | -7.6 | -3.4 | -18.6 | 5.0 | 0.14 | 0.08 | 1.3 $\pm$ 1.3 | 26.0 $\pm$ 2.3
ORPH | 15 | 172 | 76.4 | 0.100 | 277.1 | 26.2 | 54.5 | 16.3 | 347.0 | -7.7 | -4.7 | -18.5 | 2.2 | 0.19 | 0.17 | 0.8 $\pm$ 0.4 | 26.5 $\pm$ 2.1
ORPH | 16 | 461 | 185.5 | 0.100 | 282.7 | 33.0 | 62.9 | 14.5 | 282.0 | -3.6 | -4.4 | -17.8 | 7.1 | 0.18 | 0.15 | 0.8 $\pm$ 0.2 | 34.5 $\pm$ 1.6
ORPH | 17 | 162 | 58.2 | 0.100 | 273.0 | 26.7 | 53.5 | 19.9 | 282.0 | -7.2 | -0.6 | -18.2 | 10.6 | 0.58 | 0.07 | 6.0 $\pm$ 2.2 | 27.9 $\pm$ 2.8
ORPH | 18 | 119 | 40.8 | 0.100 | 271.8 | 19.5 | 46.0 | 18.2 | 329.0 | -11.4 | -3.6 | -16.1 | 9.1 | 0.34 | 0.06 | 3.9 $\pm$ 1.4 | 30.0 $\pm$ 2.5
ORPH | 19 | 202 | 68.9 | 0.100 | 266.9 | 16.0 | 40.7 | 21.1 | 307.0 | -13.5 | -3.0 | -16.0 | 7.6 | 0.25 | 0.09 | 2.0 $\pm$ 0.8 | 27.3 $\pm$ 1.3
ORPH | 20 | 79 | 24.5 | 0.100 | 271.0 | 18.9 | 45.1 | 18.7 | 307.0 | -11.4 | -3.0 | -16.3 | 4.4 | 0.33 | 0.07 | 3.4 $\pm$ 1.2 | 32.6 $\pm$ 2.0
ORPH | 21 | 141 | 49.4 | 0.100 | 269.4 | 20.9 | 46.4 | 20.8 | 286.0 | -10.5 | -2.3 | -17.5 | 8.1 | 0.29 | 0.07 | 2.8 $\pm$ 1.3 | 31.4 $\pm$ 1.9
ORPH | 22 | 189 | 71.3 | 0.100 | 279.9 | 31.1 | 60.2 | 16.0 | 361.0 | -6.1 | -2.8 | -16.1 | 11.0 | 0.26 | 0.07 | 2.5 $\pm$ 1.0 | 31.7 $\pm$ 2.7
ORPH | 23 | 342 | 134.3 | 0.100 | 285.0 | 31.6 | 62.4 | 12.1 | 363.0 | -4.0 | -6.1 | -15.9 | 10.0 | 0.21 | 0.11 | 1.4 $\pm$ 0.4 | 30.1 $\pm$ 1.6
ORPH | 24 | 277 | 116.8 | 0.100 | 281.3 | 34.2 | 63.6 | 16.0 | 391.0 | -5.0 | -3.0 | -17.2 | 10.6 | 0.29 | 0.10 | 2.1 $\pm$ 0.7 | 27.4 $\pm$ 1.9
ORPH | 25 | 300 | 115.7 | 0.100 | 286.3 | 41.7 | 72.4 | 15.3 | 328.0 | -0.7 | -3.1 | -19.4 | 9.4 | 0.33 | 0.10 | 2.3 $\pm$ 0.5 | 34.9 $\pm$ 2.6
ORPH | 26 | 237 | 101.5 | 0.100 | 285.2 | 39.7 | 70.1 | 15.2 | 353.0 | -1.7 | -3.3 | -17.9 | 9.3 | 0.47 | 0.10 | 3.4 $\pm$ 0.9 | 28.0 $\pm$ 2.7
ORPH | 27 | 117 | 45.3 | 0.100 | 282.8 | 34.6 | 64.5 | 15.0 | 361.0 | -4.2 | -3.4 | -17.8 | 6.3 | 0.28 | 0.08 | 2.5 $\pm$ 1.1 | 28.7 $\pm$ 1.8
ORPH | 28 | 615 | 247.4 | 0.100 | 283.6 | 36.9 | 66.9 | 15.3 | 352.0 | -4.0 | -3.8 | -18.4 | 1.3 | 0.27 | 0.41 | 0.5 $\pm$ 0.1 | 32.1 $\pm$ 1.2
ROS6 | 1 | 1030 | 467.8 | 0.230 | 306.9 | 40.0 | 78.5 | 0.9 | 346.0 | 8.8 | -5.8 | -7.9 | 6.6 | 0.27 | 0.25 | 0.8 $\pm$ 0.2 |  $>$ 80.0
CINR | 1 | 139 | 63.4 | 0.346 | 295.7 | 51.4 | 84.2 | 13.6 | 276.0 | 5.0 | -2.6 | -10.1 | 16.3 | 0.41 | 0.05 | 5.6 $\pm$ 3.6 | 27.9 $\pm$ 1.9
CINR | 2 | 151 | 79.0 | 0.212 | 322.2 | 50.2 | 93.3 | -0.6 | 355.0 | 8.6 | -6.4 | -6.9 | 17.1 | 0.98 | 0.06 | 11.5 $\pm$ 3.5 | 37.1 $\pm$ 4.5
ROS6 | 2 | 246 | 100.4 | 0.385 | 318.1 | 49.5 | 90.9 | 0.8 | 301.0 | 7.9 | -6.3 | -7.6 | 12.3 | 0.63 | 0.08 | 5.6 $\pm$ 3.4 | 71.2 $\pm$ 6.4
CINR | 3 | 96 | 42.3 | 0.281 | 314.1 | 47.8 | 87.9 | 1.6 | 373.0 | 6.2 | -5.5 | -6.1 | 8.1 | 0.35 | 0.07 | 3.6 $\pm$ 1.9 | 29.3 $\pm$ 2.6
CINR | 4 | 135 | 55.6 | 0.311 | 310.3 | 43.3 | 82.7 | 0.8 | 354.0 | 4.2 | -5.5 | -6.7 | 12.3 | 0.31 | 0.07 | 3.2 $\pm$ 1.8 | 38.5 $\pm$ 4.8
CINR | 5 | 475 | 177.9 | 0.109 | 312.7 | 55.5 | 93.3 | 7.2 | 330.0 | 8.5 | -2.4 | -7.1 | 9.0 | 0.27 | 0.13 | 1.5 $\pm$ 0.5 | 31.4 $\pm$ 1.4
CINR | 6 | 375 | 142.0 | 0.111 | 318.7 | 61.0 | 99.5 | 8.5 | 314.0 | 10.6 | -2.7 | -6.1 | 11.3 | 0.29 | 0.10 | 2.0 $\pm$ 0.4 | 29.9 $\pm$ 1.6
CINR | 7 | 606 | 225.6 | 0.126 | 305.8 | 49.3 | 85.8 | 6.9 | 354.0 | 5.6 | -2.0 | -7.2 | 8.1 | 0.26 | 0.16 | 1.2 $\pm$ 0.3 | 32.9 $\pm$ 1.2
CINR | 8 | 149 | 52.7 | 0.268 | 306.0 | 49.3 | 85.8 | 6.7 | 314.0 | 6.4 | -3.1 | -8.5 | 9.0 | 0.25 | 0.07 | 2.4 $\pm$ 1.1 | 39.3 $\pm$ 4.5
CINR | 9 | 93 | 35.3 | 0.537 | 302.7 | 44.3 | 80.4 | 5.9 | 312.0 | 4.3 | -3.3 | -7.0 | 10.8 | 0.28 | 0.06 | 3.2 $\pm$ 2.0 | 43.1 $\pm$ 3.5
CINR | 10 | 486 | 213.4 | 0.112 | 303.1 | 45.0 | 81.1 | 6.0 | 339.0 | 5.0 | -3.6 | -7.0 | 3.6 | 0.19 | 0.23 | 0.6 $\pm$ 0.2 | 35.7 $\pm$ 1.5

###  4.2 Ages

Figure 7: Age solution for Orpheus subgroup 22, best fit indicated by the
thick red curve. We mark stars used in the fit as black diamonds, and include
isochrones of 10 Myr, 20 Myr, 40 Myr, 80 Myr, and 1 Gyr, from top to bottom,
for reference. The age fits for all other Cep-Her subgroups are provided in
the online-only version of this figure.

Broad age categories are a useful diagnostic for verifying common origins.
While complete dynamical and lithium depletion boundary ages are beyond the
scope of this paper and will be published in an upcoming publication (Kerr, in
prep), isochronal ages are readily available using Gaia photometry and
astrometry. These ages provide a reliable relative age scale, which supports
the identification of broader formation trends. These ages can therefore
identify stars that formed in the same generation, helping us to identify co-
natal structures within the larger Cep-Her Complex.

Figure 8: Distribution of ages in Cep-Her Main by subgroup, revealing a
bimodality in age that we use to split Cep-Her Main into the older Cupavo
association and younger Orpheus association. Stars associated with Orpheus are
marked with diamonds, and stars associated with Cupavo are marked with
squares. Figure 9: Distribution of ages in Cep-Her North by subgroup. Like in
Cep-Her Main, Cep-Her North also shows an age bimodality between the older
Roslund 6 and younger Cinyras association. Stars associated with Roslund 6 are
marked with dots, and stars associated with Cinyras are marked with triangles.

We compute age by fitting the absolute magnitudes of stars in each subgroup
against PARSEC v1.2S isochrone grids (Chen et al., 2015). We use Gaia
photometry in $M_{G}$ /$G_{BP}-G_{RP}$ space for these fits, focusing on the
pre-main sequence with $1.8<G_{BP}-G_{RP}<4$, and correcting for distance and
reddening using Gaia distances (Bailer-Jones et al., 2018) and Lallement et
al. (2019) reddening maps. We limit the stellar selection to stars that
satisfy the astrometric and photometric quality cuts in SPYGLASS-IV. We also
exclude stars with $RUWE>1.2$ to removing likely unresolved binaries, as the
additional flux from the companion makes these objects appear brighter and
therefore photometrically younger than single stars. We further restrict the
sample to RUWE $<1.1$ in subgroups where at least 15 stars would remain after
applying the cut. This is an even more restrictive binarity cut that removes
nearly all stars with astrometrically detectable influence from binaries
(Bryson et al., 2020), provided that the population has enough stars to
support a strong fit after the imposition of this harsh restriction. Finally,
we require that stars included in the fit lie above the 80 Myr PARSEC
isochrone, which is well below the sequence for Cep-Her groups in most cases.
We then employ a bootstrapping routine to both limit the impact of outlying
field and binary stars and to improve our handling of the systematic
uncertainties in our age solutions. For each population and model grid, we
select half of the stars in the sample at random with a minimum selection size
of 8, with selection probability weighted by their $P_{\rm spatial}$. We do
this 10000 times for each population, defining the 2-$\sigma$ clipped mean as
the age solution and the corresponding standard deviation as the uncertainty.

In most cases, the solutions fit the densest part of the pre-main sequence,
with little impact from apparent spectroscopic binaries that lie above the
pre-main sequence, or probable field contaminants below. However, there is a
limited set of populations with age solutions exceeding 40 Myr where the cut
on the 80 Myr isochrone may remove a substantial number of photometrically
older good members, artificially lowering the age solution. For groups with
age $\tau>40$ Myr, we therefore remove this cut and re-calculate ages using a
different set of restrictions designed to avoid introducing age biases. We
reduce contamination by requiring that all stars used in the fit lie above a 1
Gyr PARSEC isochrone (which corresponds to a typical field sequence), and also
require that $P_{\rm spatial}>0.5$, provided that the result has at least 20
stars. This latter restriction is meant to counteract the increased risk of
field contamination with the looser isochrone cut. We find this choice reduces
contamination well where possible, while still producing healthy sample sizes
for bootstrapping. Ages $\tau>80$ Myr lie in a range where the pre-main
sequence becomes an increasingly poor age indicator, so we cap solutions at
this value and record them as having ages older than 80 Myr. We present the
age solutions for all Cep-Her subgroups in Figure 7. The results closely
follow the visually identifiable young sequences. While these ages are
effective for comparisons between Cep-Her subgroups, we caution against the
use of these solutions in absolute terms, as isochronal ages are highly model-
sensitive, especially in this age range.

We show the distribution of ages in Cep-Her Main and Cep-Her North in Figures
8 and 9, respectively. There we show a regional bimodality in ages in both
components of Cep-Her. These regions with common age correspond roughly to the
subtly separated sequences shown in our velocity selection fits (Figure 3),
indicating that these different ages also correspond to physical structures.
We discuss the structures indicated by this division in Section 4.3.

###  4.3 Association-Level Structures

Through the wide gap between Cep-Her Main and Cep-Her North in RV space and
the age and kinematic discrepancies between different regions of those
populations, our analysis suggests the presence of four distinct regions of
star formation within the broader Cep-Her Complex. In Cep-Her North, we see a
clear age bimodality, with the larger, young component having ages between 28
and 43 Myr, and the older component having ages near our 80 Myr limit for
analysis. The larger of the two subgroups in the older component represents
the known cluster Roslund 6 (Roslund, 1960), which is thought to have a much
older age (470 Myr; Kharchenko et al., 2013). Roslund 6 is also the sole
representative of the smaller Cep-Her North velocity sequence that we fit in
Figure 3, so its distinction from the rest of the association appears to be
supported through dynamics as well as age. The other old subgroup in Cep-Her
North does not show the same distinct dynamics, and its wide range of stellar
magnitudes suggests it may contain some stars from a younger component.
However, due to the apparent dominance of the older stars, we group it and
Roslund 6 into a population distinct from the rest of Cep-Her North. The
additional old structure may be related to the tidal tails that have already
been proposed around Roslund 6 (Bhattacharya et al., 2022), although it is far
enough from the Roslund 6 core that it could be just a component of a larger
dissolving association.

At the age of Roslund 6, SPYGLASS-IV was only able to detect populations with
a strong binary sequence. These binary-dominated populations were flagged as
old, and Roslund 6 would likely have been included in this category if
SPYGLASS-IV identified it separately from the rest of Cep-Her. Older ages
reduce the accuracy of both pre-main sequence isochronal ages and traceback,
substantially reducing our ability to reconstruct its star formation history.
We therefore exclude Roslund 6 and its companion from the discussion of
association dynamics in Section 6.

In Cep-Her Main, the sub-populations divide into two broad components. The
younger component is centered on $\delta$ Lyr, and has relatively consistent
ages between 25 and 40 Myr with a median of $\tau=30$ Myr. The older component
has a broader age range from 54 Myr to our 80 Myr mass limit. This division is
most evident in Figure 8, where two regions of differing ages emerge in
spatial coordinates. Most of this age division also reflects the split in the
space-velocity sequences shown in Figure 3. We therefore use this wide
separation in age as the basis to divide Cep-Her Main into a younger group
($\tau<45$ Myr) and an older group ($\tau>45$ Myr). This simple age cut alone
divides the association into two largely contiguous populations in space and
velocity coordinates.

The subgroups of the older component that are furthest to the galactic west
($l\sim 60\arcdeg$, see Fig. 8) are potential space-velocity outliers within
that component, and like the companion group to Roslund 6, these populations
appear to contain some young stars in a region overall dominated by older
objects. Further spectroscopic follow-up may therefore reassign some objects
to the young component. However, for the purpose of this analysis, we group
them into the old component.

The gaps in age and velocity between the four regions identified here suggest
that these populations, not the full Cep-Her complex, are closest to the
standard definition of an association as a largely unbound population of
common origin. We therefore refer to the four populations that emerge from
this analysis as “associations”, while referring to Cep-Her as a “complex” of
intertwined populations. The definition of Cep-Her nonetheless remains useful
due to the difficulty in distinguishing between members of each association
contained within. To reflect their distinct nature, we introduce new
nomenclature to refer to the four association-level populations. The
population connected to Roslund 6 is dominated by that cluster, so we refer to
it using the name of that cluster. The other three associations have extents
unlike any presented in the literature, so we assign new names to each. We
refer to the young component of Cep-Her Main that contains $\delta$ Lyr as the
Orpheus Association444This name refers to Greek hero that lulls the multi-
headed dog Cerberus to sleep with a lyre, referencing the two regions covered
by this association - Lyra, the constellation of the lyre, and the obsolete
constellation Cerberus, which is found in SE Hercules., the old component of
Cep-Her Main as the Cupavo Association, and the remaining component of Cep-Her
North as the Cinyras Association555Cinyras and Cupavo are the children of
Cycnus, the king most often associated with the story of Cygnus in Greek
mythology, and due to the youth of these associations the title of “children
of Cygnus” is fitting. The use the story of King Cycnus also makes reference
to the portion of these associations in Cepheus, the constellation of the
king..

To enable more complete demographic studies for these newly-defined
associations, we assign all unclustered stars to whichever association has the
highest local density as defined using the same clustering proximity-based
metric that was used to assign outlying members in Section 4. We provide
summary plots for each association in Figures 10, 11, 12, and 13. These show
the extents of association subgroups in on-sky space-velocity coordinates,
galactic XY coordinates, and transverse velocity anomaly relative to the mean
UVW vector of each association. Online-only, we also provide an interactive
view of the Cep-Her Complex as part of Figure 10, which allows for the
exploration of the region in sky-plane $l$/$b$ coordinates, 3D galactic XYZ
coordinates, and transverse velocity anomaly, with the ability to
independently show or hide each association. These figures therefore show the
overall structure of these associations, while also allowing the comparison of
trends in velocity with the spatial positions of the association, allowing
dynamical trends to be unveiled, which we explore further in Section 6.

Figure 10: Overview plot for the Orpheus Association, showing the distribution
of groups in $l$/$b$ sky coordinates (top), X/Y galactic coordinates (bottom
left), and $\Delta v_{T,l}$/$\Delta v_{T,b}$ transverse velocity anomaly
(bottom right), measured relative to the mean UVW velocity vector. The
transverse velocity anomaly plot shows average values to improve legibility,
as the scatter of the groups overlap heavily. The mean velocity standard
deviation for a subgroup is marked with the black error bars in the lower
right of that panel. We annotate the subgroup IDs in the top panel. Members of
bound open clusters or groups with a bound core as assessed in Section 5.7 are
marked with open icons. We assign similar colors to qualitatively similar
populations in spatial coordinates and age to aid visual comprehension. The
grey curves are lines of equal distance drawn at 50 pc intervals, which
indicate both the location of the sun and the distance spread of the
association. We also annotate the location of known bound clusters. Online-
only, we include an interactive figure in galactic XYZ cartesian space that
covers Orpheus as well as the other associations covered in Figures 11, 12,
and 13. In that figure, associations can be shown and hidden using the legend
and the display can be rotated or zoomed as desired.  Figure 11: Same as
Figure 10, but for Cinyras. The subgroup IDs are annotated in the XY galactic
coordinates panel, as the groups in Cinyras are better separated in that
parameter space. Figure 12: Same as Figure 10, but for Cupavo. Group IDs are
annotated in the XY galactic coordinates panel.  Figure 13: Overview plot for
Roslund 6. The axes and layout are the same as Figure 10, however since there
are only two subgroups we show the entire velocity distribution.

##  5 Demographics

Our substructure identification and association membership assignment enables
the first estimates of Cep-Her’s demographics for both the constituent
associations and the subgroups within. Of particular interest is the mass,
which helps to inform to virial state of Cep-Her substructures. Masses of
individual stars can be derived from stellar models, however corrections must
be made to account for missing stars, such as unresolved binary companions and
stars beyond the sensitivity limit of Gaia. In this Section, we compute
individual stellar masses, and then correct for factors that influence these
calculations to produce population-wide results for mass and stellar
population size. We then use the masses and velocity dispersions of these
groups to perform virial analyses, assessing whether any of our newly-
identified substructures can be classified as open clusters.

###  5.1 Stellar Masses

We compute stellar masses using solar-metallicity PARSEC v1.2S isochrones
(Chen et al., 2015). Cep-Her’s large overall mass makes the presence of high-
mass stars likely, so we assemble a new isomass grid from the PARSEC v1.2S
isochrones which are spaced every 0.005 M⊙ between 0.09 and 1 M⊙, every 0.01
M⊙ between 1 and 2 M⊙, every 0.02 M⊙ between 2 and 4 M⊙, and every 0.05 M⊙ for
between 4 and 20 M⊙. We assign stars the mass associated with the nearest
isomass track in Gaia de-reddened $M_{G}$ /$G_{BP}-G_{RP}$ color-absolute
magnitude space.

The highest mass we report for any star in Cep-Her is 9.55 $M_{\odot}$ for 102
Herculis in Orpheus, which, based on those same PARSEC models, places an upper
age limit for the youngest stars in the population at $\sim 24$ Myr (Chen et
al., 2015). This matches closely with our ages computed using the pre-main
sequence in Orpheus (Section 4.2). As expected, the older populations of
Cupavo and Roslund 6 both have lower maximum masses at $7.15M_{\odot}$ and
$7.6M_{\odot}$, respectively, while Cinyras, which has a younger average age
similar to Orpheus, has a larger maximum mass of $8.15M_{\odot}$. These
maximum masses correspond to an upper age limit for the youngest generation of
45 Myr in Roslund 6, 51 Myr in Cupavo, and 40 Myr for Cinyras. Of these, all
but the result for Roslund 6 are consistent with their corresponding pre-main
sequence isochronal age solutions. The unexpectedly massive star in Roslund 6
is highly outlying within that subgroup and a plausible member of Cinyras,
which would better explain its mass. The next most massive star in Roslund 6
is only 4.65 M⊙, which would imply a more plausible upper age limit at 140
Myr.

The stellar sample in Cep-Her contains 204 stars that have masses consistent
with O and B stars ($M>2.7$ M⊙; Pecaut & Mamajek, 2013). This is substantially
more than the 150 O and B stars that were used to define Sco-Cen (de Zeeuw et
al., 1999), suggesting that the total population of Cep-Her at least exceeds
the $10^{4}$ members that are estimated to reside that population (Pecaut &
Mamajek, 2016). Of these high-mass stars in Cep-Her, 101 are in Orpheus, 56
are in Cupavo, 35 are in Cinyras, and 12 are in Roslund 6.

###  5.2 Older Subgroups in Cep-Her

Our age calculations in Section 4.2 revealed 12 subgroups across Cupavo and
Roslund 6 with ages exceeding 50 Myr. The youth probability calculations in
SPYGLASS are based on the probability of an age younger than 50 Myr ($P_{\rm
Age<50Myr}$), and as a result, stars older than this limit can fall completely
outside of the range of photometry where SPYGLASS identifies stars as young.
$P_{\rm spatial}$ is also affected by this issue, as that calculation is based
on a ratio of near-certain old stars to near-certain young stars. In these
older populations, members occasionally fall in that near-certain old
category, which lowers our results for $P_{\rm spatial}$. The result is that
very low photometric and spatial-photometric membership probabilities are
possible among genuine members of these older populations. The straightforward
use of $P_{\rm fin}$ for assessing membership probability will underestimate
the mass in these older populations, necessitating an alternate means of
assessing membership. We therefore compute subgroup-wide values of membership
probability $P_{bulk}$ for populations with $\tau>50$ Myr that account for the
minimal separation between the field main sequence and pre-main sequence in
this age regime.

To do this, we select stars in each subgroup with
$0.25M_{\odot}<M<0.5M_{\odot}$, which is roughly the range where the pre-main
sequence of these older populations has some separation from the field
sequence, while avoiding the dimmest stars, which have higher photometric
uncertainties that can result in higher rates of misclassification. We exclude
unresolved binaries from these samples by applying the maximally strict cut
$RUWE<1.1$, which prevents binaries above the main sequence from breaking the
symmetry of the field stars around the 1 Gyr isochrone. We then mark all stars
located below a 1 Gyr PARSEC isochrone in this mass range as field stars.
Taking the 1 Gyr isochrone to be roughly representative of the field main
sequence, we assume that single field stars are distributed symmetrically on
either size of the field sequence. We therefore take the number of field stars
($N_{field}$) in the sample to be double the number of stars found below the 1
Gyr isochrone. The mean membership probability in the population $P_{bulk}$ is
then 1 - ($N_{field}$/$N_{tot}$), where $N_{tot}$ is the total number of
potential association members in the chosen mass range.

This result for $P_{bulk}$ assumes that stars of a given mass are represented
equally in the association and the field, however this is not true for stars
more massive than $\sim 1$M⊙, where field stars may be old enough to evolve
into compact objects and thus leave our sample. We must therefore down-weight
the expected field contribution in this mass range. Following Binney et al.
(2000), we assume a constant star formation in the field over the last 11.2
Gyr. We take a sequence of PARSEC solar-metallicity isochrones from 1 Myr to
11.2 Gyr, and record the maximum stellar mass in each isochrone. This provides
a lifespan for a star of a given mass. If that lifespan is less than 11.2 Gyr,
a star of that mass may have evolved into a compact object. We therefore
define a corrective factor $f_{alive}$ as the lifespan of a star at a given
mass, divided by 11.2 Gyr. This value approximates the fraction of field stars
of a given mass that have yet to become compact objects, assuming a constant
star formation rate over the age of the solar neighborhood (Binney et al.,
2000). The inclusion of this corrective factor changes our bulk membership
probability to $P_{bulk}=1-\frac{f_{alive}N_{field}}{N_{tot}}$. Uncertainties
in these membership probabilities are dominated by the binomial uncertainty in
the number of field stars $N_{field}$, and we propagate these uncertainties
through to our mass and population results.

We calculate $P_{bulk}$ for all subgroups in the older associations of Cupavo
and Roslund 6. We also calculate $P_{bulk}$ for the outlying unclustered
regions separately in each of these associations. These calculations are used
to produce approximate overall membership rates in these older populations,
which we use to assess association mass and population size in Section 5.6.
These values therefore replace $P_{\rm fin}$ when calculating cluster
demographics, however we do not change the values of $P_{\rm fin}$ in Table 2.

Populations older than 50 Myr are only marginally detectable through SPYGLASS,
and are typically identified through their binary sequence alone. These
populations therefore tend to be quite dominant in their local parameter space
if detectable through SPYGLASS. As such, the bulk membership probabilities in
most old populations are high, with stars in all subgroups being assigned
$P_{bulk}>0.7$. Outlying components of Cupavo and Roslund 6 receive
$P_{bulk}=0.72$ and $P_{bulk}=0.78$ for low-mass stars, respectively, both of
which are quite high for an unclustered halo. This may suggest that there are
substantial lower-density populations in Cupavo beyond the SPYGLASS
sensitivity limit for groups in this age range.

###  5.3 Binaries

Binaries are a major source of uncertainty in the mass estimates and overall
stellar demographics of a young association. Unresolved companions require
dedicated spectroscopic measurements for their full characterization,
resulting in their mass being difficult to ascertain and often overlooked. We
therefore identify probable binaries in Cep-Her, and use these to correct our
demographic calculations.

####  5.3.1 Identifying Binaries

To identify resolved binaries, we search an on-sky radius of $10^{4}$ AU
around each candidate member in Gaia, marking all stars with $\Delta v_{T}<3$
km s-1 and $\frac{\Delta\pi}{\pi}<0.2$ as a probable binary companion. This is
similar to the selection in SPYGLASS-II, only changing the velocity cut to
adjust for the distance of Cep-Her. 3 km s-1 corresponds roughly to the
circular velocity of a companion at 1 arcsecond separation from a G or F star
at the near edge of Cep-Her, with the addition of a 1 km s-1 buffer to account
for the typical Gaia measurement uncertainties for dim objects at the distance
of Cep-Her. The 1 arcsecond value corresponds to the smallest separations
resolvable through Gaia that have the highest orbital speeds, so this 3 km s-1
velocity limit represents the maximum velocity difference expected for a
resolved stellar companion (Rizzuto et al., 2018). Changing the selection
limit to 2 km s-1 or 4 km s-1 only changes the number of binaries by
$\sim$20%, verifying that most of the binaries are found within 2 km s-1, and
that 3 km s-1 provides a reasonable buffer for high-uncertainty objects
without introducing substantial contamination.

We avoid duplicate systems by removing searches if a brighter Cep-Her member
is found within our search radius, only recording the result from the search
around the brighter object, which we mark as the primary. Due to subtly
different distances to different stars in a system, it is possible for the
projected on-sky distance between two stars to disagree on whether a star is
within the 10000 AU limit depending on which distance is adopted. In these
cases, the distance to the brighter star is adopted, and cases where the
primary is already included in another system are removed. In cases where one
companion is assigned to multiple systems but the primary of one system is not
assigned as a companion of the other, the primaries of both systems are marked
as primaries, and common companions are listed under both possible parent
systems. We record probable resolved binary or multiple systems in Appendix A.

For completeness, we also mark unresolved binaries, which are identifiable
using other available metrics. As mentioned in Section 3.2.4, the Renormalized
Unit Weight Error (RUWE), which is a metric based on the goodness of fit from
Gaia’s photometric solution, is a strong indicator of an unresolved companion.
We therefore mark stars with RUWE$>1.2$ as probable unresolved binaries,
following Bryson et al. (2020). We also mark stars as unresolved binaries if
they have spectroscopic line profiles consistent with a double-lined
spectroscopic binaries (SB2), either through the visual identification of a
bimodal line profile in our observations, or through SB2 flags provided in
literature RV catalogs. Stars with evidence of spectroscopic binarity are
marked with a flag in Table 2.

####  5.3.2 Mass in Binaries

To compute the mass in stellar companions, we follow the methods employed in
SPYGLASS-II, first stripping each probable binary system identified in Section
5.3.1 to only its primary. This ensures that no systems are double-counted in
our statistics, as the population-wide contribution of binaries, resolved and
unresolved, can be inferred using the combination of known companion
demographics and primary mass. The presence of an unresolved companion moves
stars in a direction near-perpendicular to the pre-main sequence in the HR
diagram, which has a relatively limited effect on the photometric mass
estimate of the primary. These mass estimates therefore tend to accurately
reflect the mass of the primary regardless of the presence of a companion,
allowing us to use these masses to set the binary demographics regardless of
unresolved binarity.

For each primary star, we use existing measurements of binarity rate and mean
mass ratio to compute an expected number of companions and total companion
mass, following the methods of SPYGLASS-II. We use the binarity rates as a
function of mass from Sullivan & Kraus (2021), interpolating a binarity rate
for each star off of that curve. We then use the power law indices for mass
ratio to compute an expected mass ratio for a hypothetical companion to each
star in the sample, again interpolating the values of the expected mass ratio
between mass bins (Sullivan & Kraus, 2021). For stars with $M>2$M⊙, we adopt
the mean mass ratio of the highest mass bin from Sullivan & Kraus (2021). We
also include the Duchêne & Kraus (2013) companion fraction for stars with
$1.5M_{\odot}<M<5M_{\odot}$ and adopt that result for a mean mass of 3.25
$M_{\odot}$. We interpolate binarity rate and expected mass ratio for a star
of a given mass using the resulting curves. The expected number of binaries is
then provided by the binarity rate, while the expected mass in binaries is the
binarity rate multiplied by the expected mass ratio and the mass of the
primary. On average, the inclusion of stellar companions increases the number
of stars by a factor of 1.42, and increases the mass by a factor of 1.24.

###  5.4 Correction for Low-Mass Stars

Due to the distance of Cep-Her, the lowest-mass stars lack quality Gaia
entries. Many of these stars have marginal detections, likely forming a
substantial fraction of the widely-scattered Gaia measurements in Figure 6
that fail our astrometric and photometric quality cuts. Other dim stars in
Cep-Her lack $G_{BP}$ and/or $G_{RP}$ photometry entirely, especially M dwarfs
where $G_{BP}$ falls below the Gaia flux limit. In these cases, it is
impossible to assess their membership through Gaia photometry. In Figure 14,
we show a mass histogram for all stars in Cep-Her that pass our astrometric
and photometric quality cuts. This histogram provides the number of stars per
mass bin, weighted by the membership probability $P_{\rm fin}$. We compare
that result to the individual object Initial Mass Function (IMF) from Chabrier
(2005), which we scale to best fit to the mass histogram using least-squares
optimization.

The form of our mass distribution for Cep-Her is broadly similar to the
canonical IMF, however also has a notable deficit of stars at $M\sim
0.5M_{\odot}$ and an overabundance around $M\sim 0.25M_{\odot}$. This pattern
has been seen in other young populations, and is likely caused by slight model
inaccuracies in the evolution of stars in this section of the pre-main
sequence (Kraus et al., 2014). We also see a deficit of stars more massive
than $\sim 2M_{\odot}$, although we caution against the overinterpretation of
this result, as both the saturation of stars with apparent magnitude $m_{G}<6$
in Gaia and the overabundance of binaries in this mass range can contribute to
a lower rate of quality Gaia astrometry (Lindegren et al., 2018). Cep-Her may
therefore benefit from a follow-up survey targeting bright stars that fail the
astrometric and photometric quality cuts to confirm whether this deficit of
massive stars is real. However, the most significant deficit of stars that we
see is for masses $M\lesssim 0.2M_{\odot}$, where the mass histogram drops
well below the canonical IMF, and no stars have $M<0.11M_{\odot}$.

Figure 14: Our correction for stars below Gaia’s sensitivity limit. The
histogram shows the distribution of stars in Cep-Her, which we show alongside
a smoothed version of that histogram and the Chabrier (2005) individual star
IMF. The shaded area between the curves for low-mass objects indicates the
missing mass which we correct for. The deficit of stars at $M\sim 0.5$ and
overabundance at $M\sim 0.25$ is a known issue, likely caused by model
inaccuracies in stellar evolution times that displace stars to lower masses
(Kraus et al., 2014).

This deficit of low-mass stars relates directly to the sensitivity limits of
Gaia, and we must therefore apply a correction to account for the mass and
population of Cep-Her members in this mass range. To do this, we split our
stellar sample into 4 distance bins, and for stars in each of those bins we
create a separate weighted mass histogram. All of these have the same form as
the histogram shown in Figure 14, and we smooth each result to suppress
unphysical variations using a Savitsky-Golay filter with a window size of 11
and a polynomial of order 2. We assign an IMF to each result, which is scaled
to produce the best fit with the histogram for $M>0.25$M⊙ via least-squares
optimization. We then set the lowest-mass crossover point between the scaled
IMF and the smoothed mass histogram as the point at which mass starts being
lost. The region of parameter space where stars are lost to the Gaia
sensitivity limit is then provided by the region between the IMF and the
smoothed mass histogram for stars less massive than this crossover point, as
illustrated in Figure 14. The corrective factor for the number of stars is
then the sum across the entire stellar histogram, plus the sum over the region
where stars are missing (red shaded area in Fig. 14), all divided by sum of
the stellar histogram. We limit this calculation to $M>0.09M_{\odot}$, which
is the mass limit in the PARSEC isochrones. We compute the missing mass in a
similar way, simply multiplying the histograms and corresponding IMFs by the
mass at each bin element and following the same calculation as for the number
of stars. We do not limit this calculation in mass, as mass in objects below
our $M>0.09M_{\odot}$ cut can still contribute to the boundedness of clusters.

Upon computing the missing number and missing mass fractions for each distance
bin, we see a gradual linear trend decreasing completeness and increasing
corrective factors with distance, as expected for a deficit caused by Gaia
sensitivity limits. For each star in our sample, we assign a correction by
interpolating on the number of missing stars and the missing mass between the
distance bin centers. This produces corrective factors that raise the number
of members by a factor of between 1.35 and 1.41, and raise the total mass by a
factor of between 1.12 and 1.13. The sum of spatial-only membership
probabilities $P_{\rm spatial}$ for objects that fail our astrometric and
photometric quality cuts implies that 4725 of them are members, compared to
the 6798 stars added by this correction. This suggests that more than two-
thirds of the stars added through this correction may have Gaia entries, but
lack the quality photometry necessary to assess their membership.

###  5.5 Association-Level Cross-Contamination

A consequence of the substantial overlap between the sky-plane velocities and
positions of Cep-Her Main and Cep-Her North is that the membership in one over
the other cannot always be confidently established without RVs. In Section 3.2
we separated the two populations in RV space, assigning stars without RVs to
the closest population in 5D space-velocity coordinates. While this procedure
assigns all stars to their most plausible group, there remains a risk of
cross-contamination, or the assignment of a star to the wrong parent group due
to the lack of an RV measurement.

To assess the impact of this cross-contamination on our measurements of
population size and mass, we investigate the rate at which we assign stars to
one population versus another with two-dimensional velocities and compare it
to the same rate among RV-confirmed members, producing a corrective factor
that we can apply to our mass and population size estimates. We first pass
each star in both Cep-Her Main and Cep-Her North through the procedure used to
assign outlying population members to a parent group in Section 4, this time
finding an alternate subgroup assignment, or the most likely parent subgroup
in the RV-defined region that each star was not initially assigned to. We
assign stars to an alternate subgroup if they pass the clustering proximity
dominance metric used in Section 4, provided that the $D_{8}$ for that star
within its best-fit subgroup is less than the maximum value of $D_{8}$ for the
assigned subgroup members. For stars that receive RV follow-up, these
alternate group assignments provide a straightforward way to update
membership, as their assigned parent population can be easily updated from
their current assigned group to their alternate group if membership in the
other RV-defined region is confirmed. We therefore record these alternate
group assignments in Table 2.

The combination of the original assigned members and alternate members provide
a complete list of stars that can be plausibly connected to a given subgroup.
The ratio between the number of stars assigned to the subgroup and the
combined number of assigned and alternate stars therefore provides the
fraction of stars with consistent space-velocity positions that are included
in that subgroup’s stellar sample. We can compute the same fraction using only
stars with reliable RVs ($\sigma_{RV}<5$), which provides the fraction of
stars with plausible positions in 5D space-velocity coordinates that remain
consistent with the assigned subgroup when an RV is available. We take this
latter ratio to represent the true membership rate in the subgroup, and
therefore compute the cross-contamination corrective factor as the fraction of
members among stars with good RVs, divided by the fraction of all plausible
members that are assigned to the subgroup rather than its alternate.

We assign the cross-contamination factor of the parent subgroup to each
stellar primary. For stars without a parent subgroup, we re-compute the value
of cross-contamination across the entire association, and assign the
unclustered stars this averaged value. The use of these fractions introduces a
binomial uncertainty into these corrections, which we propagate through to our
mass and population size measurements in Section 5.6.

In this correction, values less than one imply that more stars are assigned to
the subgroup than is merited by the membership rates among stars with good
RVs. Our cross contamination corrective factors are usually at or near 1,
however they can be as low as 0.72 and as high as 1.23. The corrective factor
of 1.23 is assigned to CINR-1, where there are far more RVs consistent with
the population than its relatively small assigned membership would predict.
The corrective factor of 0.72 is assigned to CINR-9, which overlaps with
several RV-dominant populations in central Cupavo. Most of the corrective
factors that are not near 1 are found in Cinyras and Cupavo, as these
associations overlap the most in 2D velocity space.

###  5.6 Expected Masses and Population Sizes

Expected masses and population sizes are computed by combining our stellar
populations and masses with the expected binary contribution, losses due to
stars without adequate Gaia coverage, and corrections for cross contamination.
The expected number of stars present for a stellar single or primary of a
given mass is provided by the following formula:

$N_{x}=P_{\rm fin}*(1+f_{bin})*c_{N,lm}*c_{cc}$ (4)

where $f_{bin}$ is the binarity rate for a star given its mass, $c_{N,lm}$ is
the correction for the number of stars below the mass detection limit, and
$c_{cc}$ is the cross contamination correction for the star’s parent region.
The expected mass is then computed as:

$M_{x}=M_{*}*P_{\rm fin}*(1+(f_{bin}*\overline{r_{M}}))*c_{M,lm}*c_{cc}$ (5)

where $\overline{r_{M}}$ is the expected mass ratio of a star given its mass,
and $c_{M,lm}$ is the mass corrective factor for objects below our detection
limit. The total number of stars $N_{tot}$ and mass $M_{tot}$ in an
association or subregion can then be computed by summing $N_{x}$ and $M_{x}$,
respectively, across all singles or system primaries in the population. We
report the expected masses and population sizes for each subgroup in Table 3,
and use the population masses to assess the virial states of subgroups in
Section 5.7.

We also provide fractional uncertainties in Table 3, which combine the
binomial uncertainties in the cross-contamination correction and membership
probability calculation for older subgroups with systematic uncertainties. The
mass models, IMF correction for low-mass stars, and the binarity correction
all contribute to the systematic uncertainties, and we adopt a value of 10% in
most cases, which combines the approximate maximum discrepancies between the
IMF and our distribution of stellar masses, and typical literature variation
in binary statistics (Duchêne & Kraus, 2013; Sullivan & Kraus, 2021). We raise
the systematic uncertainty to 20% for the old populations where we use
$P_{bulk}$ in place of $P_{\rm fin}$, as the use of a limited stellar mass
range in the $P_{bulk}$ calculation and potential for biases from binary
contaminants introduce additional systematic uncertainties. The values of
fractional uncertainty we report average the upper and lower uncertainty
intervals in cases where the uncertainty is asymmetric. In cases of young
subgroups without plausible cross-contaminants, only the 10% systematic
uncertainties affect the reported values.

The binomial uncertainties in the cross contamination decrease substantially
at the association level, due to the large number of stars on that scale. We
therefore use the binomial uncertainty in the association-wide cross
contamination calculation as the statistical uncertainty in the cross
contamination rate for an association, calculated using the same methods
employed in Section 5.5. We also compute association-level uncertainties in
$P_{bulk}$ for the older associations, computing binomial uncertainties based
on the number of field stars across all subgroups in both Roslund 6 and
Cupavo. These two sources of binomial uncertainties are added in quadrature
with our systematic uncertainties to produce statistics at the association
level.

After all corrections, we find that Orpheus is the largest association in the
Cep-Her complex with $9552\pm 960$ members totalling $3881\pm 390$ M⊙. This
gives it a membership rivalling the $10^{4}$-star population of Sco-Cen
(Pecaut & Mamajek, 2016), which is substantially larger than the $\sim 2000$
M⊙ Orion Nebula Cluster (Hillenbrand, 1997). Cupavo is next largest, with
$8794\pm 1827$ members containing $3620\pm 752$ M⊙, followed by Cinyras with
$3872\pm 455$ members and $1604\pm 188$ M⊙, and finally Roslund 6, with
$1507\pm 347$ members and $700\pm 161$ M⊙. With nearly $10^{4}$ M⊙ and over
$2.3\times 10^{4}$ stars across all four component associations, Cep-Her as a
whole may be larger than the entire Orion Complex, which is thought to contain
less than $2\times 10^{4}$ stars (Bally, 2008).

###  5.7 Boundedness

We assess the boundedness of our subgroups following the methods used in
SPYGLASS-II, with some optimizations to better automate the process and to
minimize the effects of contamination. Bound clusters satisfy
$\sigma_{1D}<\sqrt{2}\sigma_{virial}$, where $\sigma_{1D}$ is the square root
of the average variance in a multi-dimensional velocity-space, and the virial
velocity, $\sigma_{virial}$, is defined as follows (Portegies Zwart et al.,
2010; Kuhn et al., 2019):

$\sigma_{virial}=\left(\frac{GM}{r_{hm}\eta}\right)^{-\frac{1}{2}}$

where $M$ is the subgroup mass, $r_{hm}$ is the half-mass radius, and $\eta$
is a value that represents the cluster’s mass profile. For a Plummer profile,
$\eta=10$, however young clusters often have lower values of $\eta$ due to
having broader density profiles (Kuhn et al., 2019). We take $\eta=5$ for our
calculations, typical of these broader profiles. We define $r_{hm}$ in 2D
$l$-$b$ coordinates to remove the substantial uncertainties in the distance
dimension, and use two different methods to compute $r_{hm}$. First, we
identify $r_{hm}$ as the radius out from the subgroup center that contains
half of the mass, and use the full mass of the subgroup as $M$. Second, we
compute the density profile of the subgroup using a KDE, fit a bivariate
gaussian to the resulting profile, and compute $r_{hm}$ using the mass
expected by the fit, rather than the measured mass. This secondary method
accounts for dense populations with wide halos, where a core is bound but the
outer reaches are not. At the age of Cep-Her, clusters are already dissolving,
so the presence of a plausibly bound region indicates that the subgroup as a
whole was likely fully bound at formation, with all of the mass associated
with a cluster’s widely-scattered halo in the central core. Our quoted
$\sigma_{virial}$ uses the smaller of the two $r_{hm}$ values. The uncertainty
in $\sigma_{virial}$ then combines the fractional uncertainty in mass with a
10% systematic uncertainty for $r_{hm}$.

For stars within the selected value of $r_{hm}$, we compute $\sigma_{1D}$ in
$\Delta$vT coordinates, excluding probable binaries, and list the results in
Table 3. The restriction to single stars within $r_{hm}$ ensures a maximally
pure sample of association members, while minimizing any internal sources of
velocity uncertainty. In cases where the binarity restrictions would leave
fewer than 20 stars, we remove the binarity cuts, which prevents unphysical
results for velocity dispersion that can emerge through small number
statistics. For the chosen subset of stars on which to compute $\sigma_{1D}$,
we then bootstrap, selecting half of the stars at random to a minimum of 10
and a maximum of 60. For each bootstrapped sample, we compute a standard
deviation, and subtract the median velocity uncertainty from the result in
quadrature to remove the contribution from uncertainty to the velocity
dispersion. To further limit the effect of outliers, we sigma-clip each result
at $2\sigma$. The value for $\sigma_{1D}$ is then the median of the sigma-
clipped standard deviation across 1000 bootstrapped samples, and its
uncertainty is the corresponding standard deviation.

Figure 15: The virial velocities and 1D velocity dispersions for Cep-Her
subgroups. The curve shows $\sigma_{1D}=\sqrt{2}\sigma_{virial}$, which is the
boundary between bound and unbound populations. Subgroups above this curve are
open clusters, marked with orange hexagons, and stars below the line are
unbound substructures, marked with purple triangles. Marker sizes are scaled
with the subgroup mass. The bound substructure found in ORPH-23 is marked with
a red X.

From the framework in Kuhn et al. (2019), populations with
$\sigma_{1D}<\sqrt{2}\sigma_{virial}$ are bound, so we can define a virial
ratio as $\sigma_{1D}/\sqrt{2}\sigma_{virial}$, with bound populations having
a virial ratio less than 1. For groups with an initial virial ratio less than
2, we inspect the population by hand to check whether a small high-density
core is present that may be smaller than the kernel size for the fit to
$r_{hm}$. In seven cases where the kernel appears to be too large for the core
that we find, we reduce the kernel width and re-calculate. This is mainly
necessary in cases like CUPV-3, which has a dense core that is only about 3 pc
across (0.5∘ on-sky), however some populations with slightly looser cores such
as ROS6-1 are also included. In cases with a low virial ratio but a kernel
width that appears suitable for the cluster, we experiment further with
changes to the kernel size, and find that our choice of kernel does not
substantially change our virial results.

The final virial ratios are presented in Table 3. We present the distribution
of $\sigma_{virial}$ vs $\sigma_{1D}$ in Figure 15, coloring each population
according to their most probable virial state. There we demonstrate the
presence of two distinct sequences - one containing virialized clusters
roughly consistent with with a linear trend elevated above the boundedness
limit, and a sequence of unbound subgroups, which have similar values of
$\sigma_{virial}$ but varied values of $\sigma_{1D}$. The virial ratio of 1
that we include as a black curve in Figure 15 roughly separates the two
sequences, enforcing this limit as a transition between two different
dynamical states.

Three populations, CUPV-3, ORPH-28, and CINR-10, are high-confidence open
clusters, meaning that they have a low virial ratio that is not within
uncertainties of 1. Four more, ROS6-1, CUPV-10, ORPH-15, and ORPH-16, are
probable clusters, as they have virial ratios less than 1, but have
uncertainties high enough that they are consistent with an unbound state.
Finally, four populations, CINR-7, ORPH-12, ORPH-14, and ORPH-23 appear
unbound in our analysis, however their uncertainties do not rule out the
possibility of them being bound. ORPH-23 is a unique case, as it appears to
contain an embedded high-density substructure. Since this overdensity is near
the edge of OPRH-23, it does not appear to represent a central core, and we
therefore do not take the virial state of this region to represent the entire
subgroup. However, its size and shape is similar to a series of emerging low-
mass and extremely compact clusters that have been discovered recently, such
as EE Dra in CFN and $\eta$ Cha in Sco-Cen (Mamajek et al., 1999; Kerr et al.,
2022b). By limiting our $\sigma_{virial}$ and $\sigma_{1D}$ calculations to
only this core, we found a virial ratio of $0.54\pm 0.28$, here, indicating
that this substructure is bound. ORPH-23 therefore contains a embedded open
cluster, which may disrupt stellar motions. We include this bound substructure
separately from ORPH-23 in Figure 15.

Six of the probable bound structures we identify in Cep-Her have already been
recognized in published open cluster catalogs (Cantat-Gaudin & Anders, 2020,
e.g.,). The largest of these is the dense open cluster NGC 7058, which is
found in our catalog as CUPV-3. The other two high-confidence open clusters,
ORPH-28 and CINR-10, are also known structures, corresponding to the $\delta$
Lyrae cluster and RSG-5, respectively (Stephenson, 1959; Röser et al., 2016).
Of the probable bound structures, ROS6-1 represents the core of the open
cluster Roslund 6 (Roslund, 1960), and ORPH-15 corresponds to UBC 9 (Castro-
Ginard et al., 2018). The other two probable bound structures, CUPV-10 and
ORPH-16, have not yet been recognized as open clusters. Finally, the dense
core we identify in ORPH-23 appears to correspond to the small open cluster
ASCC 100 (Piskunov et al., 2006).

###  5.8 Completeness

While our expected mass and population estimates account for most sources of
missing mass within the boundaries of Cep-Her, there are certain to be
additional members located beyond the limits to $P_{\rm spatial}$ which we
have imposed. Stars in these outer reaches are likely to have been ejected
soon after or even during formation, limiting their effect on the long-term
virial state of Cep-Her substructures. However, the number of members in these
outer reaches may constitute a substantial fraction of Cep-Her’s total
population (e.g., Fujii & Portegies Zwart, 2016).

The sum of $P_{\rm spatial}$ across all stars beyond our $P_{\rm spatial}>0.2$
limit is over $1.5\times 10^{4}$, suggesting that as much as half of the mass
of Cep-Her may be outside of our $P_{\rm spatial}$ limit. However, we also
notice that 63% of the stars in these outer reaches have poor photometric
quality or high RUWEs, compared to just 45% in the core. This suggests that
some of these stars may be poor-quality Cep-Her members with sufficiently high
uncertainties that they are not readily associated with Cep-Her’s core
populations. These are already partially accounted for by our correction for
stars without quality Gaia entries, although only $\sim 2000$ stars predicted
by this correction are unaccounted for based on $P_{\rm spatial}$ for stars in
our core Cep-Her sample (see Sec. 5.4). This would leave approximately
$1.3\times 10^{4}$ stars in this outer halo as plausible additional members of
Cep-Her, which would raise the total number of stars in the complex by over
50%.

The outlying stars likely disproportionately originate in Cupavo and Roslund
6, as their older ages give stars more time to reach these widely-distributed
locations. Stars in these older populations are also more likely to be missed
in the initial clustering in SPYGLASS-IV, so the stars in this outlying region
could form entire coherent sub-populations that SPYGLASS is unable to detect.
Better constraints on the sizes of these outlying populations will require
traceback to reveal the origins of these outlying objects, as well as the
modelling of any hidden coherent associations found outside the extent of Cep-
Her we explore here. Updates to the Gaia astrometric solution are also likely
to help, as they will improve the overabundant low-quality Gaia solutions that
dominate this extended sample. However, true in-depth studies of these
populations are likely to require entirely new methods that are more suitable
for isolating these older populations. One possible approach may be to select
stars based on the Gaia excess variability metric, which has recently been
shown to predict age within 10-20% for stars younger than 2.5 Gyr (Barber &
Mann, 2023).

##  6 Population Dynamics

###  6.1 Expansion Velocities

Young associations gradually expand after their gas cloud disperses, as they
typically have densities far too low to overcome the initial internal velocity
dispersions of stars, except in relatively localized overdensities like
embedded open clusters (e.g., de Zeeuw et al., 1999). Dynamically coherent
young associations are therefore expected to have stars with velocity vectors
pointed radially outwards from the association center, with higher velocities
for more outlying members. Deviations from a divergent trend can indicate the
presence of expanding substructures within the broader population, or
contaminant populations with velocities that do not place them close to the
rest of the association at formation.

Figure 16: Expansion trends in Orpheus, Cinyras, and Cupavo, shown as distance
along an axis versus velocity along that axis. The top row shows distance vs
RV, and the bottom and middle rows show the corrected transverse velocity in
$l$ and $b$, respectively versus the linear distance in the $l$ and $b$
direction relative to the mean at the distance of each star. We include
uncertainties for the distance axes, and error bars for the $l$ and $b$ axes
are typically smaller than the data points.

By Wright & Mamajek (2018), we can identify expansion through a positive trend
in velocity along a given axis. Along the distance/RV axis, closer members of
an expanding group have an RV component pointing away from the cluster center
and towards the observer, while stars on the far edge have a velocity
component away from the observer (Wright & Mamajek, 2018). Along this axis and
any other linear axis in 3D space, a line can be fit to the distribution of
velocities along that axis, producing a rate of expansion, $\kappa$. In on-sky
galactic $l$ and $b$ coordinates, this calculation is more complicated, as the
sky-plane component of stellar motion is subject to a perspective effect in
which stars moving towards the observer show virtual expansion and vice versa
(Kuhn et al., 2019). Since we have well-constrained RVs for all three
associations in our sample, we can correct for this effect using the following
formulae, following Brown et al. (1997):

$\displaystyle v^{\prime}_{T,l}$ $\displaystyle=v_{T,l}-\overline{v_{R}}*\cos
b_{c}\sin(l_{c}-l)$ (6) $\displaystyle v^{\prime}_{T,b}$
$\displaystyle=v_{T,b}-\overline{v_{R}}(\sin b_{c}\cos b-\sin b\cos
b_{c}\cos(l_{c}-l))$ (7)

Where $\overline{v_{R}}$ is the mean radial velocity, and $l_{c}$ and $b_{c}$
are the mean position of the association in galactic sky coordinates. This
allows us to calculate the on-sky expansion rate in the $l$ and $b$
directions.

Figure 16 shows the expansion trends along the distance axis and the on-sky
galactic $l$/$b$ axes. For the distance axis, we plot distance against RV,
while for the $l$ and $b$ axes we plot corrected transverse velocity against
$\Delta d_{l}$ and $\Delta d_{b}$, the on-sky distance in pc between the mean
value of $l$ and $b$ for the cluster and a given star at its distance. To each
trend, we fit a line, which provides the expansion rate along that axis. For
our fits to all axes, we remove binaries, as well as stars with $P_{\rm
fin}<0.8$ for Cinyras and Orpheus and $P_{\rm fin}<0.5$ for Cupavo. These
restrictions limit the presence of interlopers and apparent motions from
stellar companions, with the looser cut in Cupavo justified by the
underestimation of $P_{\rm fin}$ in this older population. For the distance/RV
axes, we additionally require that $\sigma_{RV}<1$ km s-1 to ensure that the
stars we include have adequate RV data. All fits assume equal uncertainties,
as most of the scatter seen is likely to be physical, rather than due to
measurement uncertainties.

The results show that 8 of the 9 subplots (3 axes across 3 associations) show
statistically significant expansion at more than a 2-$\sigma$ level. The rest
show a general trend of expansion, although often with very different rates
along the three different axes. These varied expansion rates by axis are a
common feature of studies like this, and may reflect spatially varied stellar
feedback after the first stellar generation or other forms of substructure
(e.g., Wright & Mamajek, 2018; Cantat-Gaudin et al., 2019b; Pang et al.,
2021). Galactic orbital effects are another potential cause, especially in the
$b$ axis. As stars orbit the galaxy, the rise and fall through the disk. This
results in stars reaching zero vertical velocity at the peaks of their
epicycles before falling back towards the galactic plane. The result is that
for an association with a wide range in galactic $Z$ and consistent epicyclic
amplitudes for its members, stars will converge in $Z$ and $b$ space before
diverging as the last stars begin heading back to the galactic midplane. This
is the likely cause for the relatively weak divergence trend in $b$ that we
display in Orpheus, which, with its high position in galactic $b$ and low
$v_{T,b}$, is near the peak of its epicycle. A similar effect is visible in
Cinyras, where we see an overall positive divergence trend in $b$ that turns
negative for stars with the highest values of $b$.

Figure 17: Vector plots for the corrected transverse velocities of stars in
Orpheus, Cinyras, and Cupavo, spatially averaged. A 5 km s-1 velocity vector
is shown in the bottom-right of each figure for reference. All three broadly
show expansion from a common origin along at least one axis, with some
evidence for further dynamical substructure. The arrow widths are proportional
to the square root of the number of stars included in the velocity averages
for each grid cell to indicate where substantial stellar overdensities lie.
The dark dots and arrows in the Orpheus panel show the locations and velocity
vectors of the three fully-bound subgroups (ORPH-15, 16, and 28).

###  6.2 Dynamical Coherence

We have confirmed that all three associations in Cep-Her are expanding,
however deviations within these overall velocity trends may indicate the
presence of additional distinct structures at formation. We map these
expansion trends in Figure 17, which shows the spatial distribution of
corrected on-sky transverse velocity vectors in Orpheus, Cinyras, and Cupavo,
averaged by bin in $l$ and $b$. The conditions for inclusion in these average
velocities are the same as those we use for our fits to $\Delta d$ and
$v^{\prime}_{T}$ (Figure 16). Velocities are shown relative to the group
average. The result provides a visual representation of expansion that reveals
regional patterns that are averaged over during bulk measurements of
expansion. All three populations show visual indications of expansion from a
common origin in at least some form, however they also show some anisotropies
that reflect the anomalies visible in Figure 16.

Cinyras shows especially consistent expansion in all directions on-sky, with
most vectors pointing radial outward from an apparent origin. Stars with high
$b$ are an exception to this trend, however, as the $+b$ component of the
corrected transverse velocity decreases as it approaches $b=10$ before turning
over entirely. This is the same pattern shown in Figure 16, and is likely
attributable to motion near the peak of an epicycle.

Cupavo also has relatively consistent expansion across most of its extent,
however Figure 17 shows substantially higher rates of divergence in the
galactic west ($-l$) of the population near Orpheus. These high velocity