# Theia 456: Tidally Shredding an Open Cluster

Kyle R. Tregoning Department of Physics, University of Florida, 2001 Museum
Rd, Gainesville, FL 32611, USA
Jeff J. Andrews Department of Physics, University of Florida, 2001 Museum Rd,
Gainesville, FL 32611, USA
Institute for Fundamental Theory, University of Florida, 2001 Museum Rd,
Gainesville, FL 32611, USA  Marcel A. Agüeros Department of Astronomy,
Columbia University, 550 West 120th Street, New York, NY 10027, USA
Laboratoire d’astrophysiquede Bordeaux, Univ. Bordeaux, CNRS, B18N, Allée
Geoffroy Saint-Hilaire, 33615 Pessac, France  Phillip A. Cargile Center for
Astrophysics — Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138,
USA  Julio Chanamé Instituto de Astrofísica, Pontificia Universidad Católicade
Chile, Av. Vicuña Mackenna 4860, 782-0436 Macul, Santiago, Chile  Jason L.
Curtis Department of Astronomy, Columbia University, 550 West 120th Street,
New York, NY 10027, USA  Simon C. Schuler University of Tampa, Department of
Physics and Astronomy, Tampa, FL 33606, USA

###### Abstract

The application of clustering algorithms to the Gaia astrometric catalog has
revolutionized our census of stellar populations in the Milky Way, including
the discovery of many new, dispersed structures. We focus on one such
structure, Theia 456 (COIN-Gaia-13), a loosely bound collection of $\sim 320$
stars spanning $\sim 120$ pc that has previously been shown to exhibit
kinematic, chemical, and gyrochronal coherency, indicating a common origin. We
obtain follow-up radial velocities and supplement these with Gaia astrometry
to perform an in-depth dynamical analysis of Theia 456. By integrating stellar
orbits through a Milky Way potential, we find the currently dispersed
structure coalesced into a small cluster in the past. Via Bayesian modeling,
we derive a kinematic age of $245\pm 3$ Myr (statistical), a half-mass radius
of $9\pm 2$ pc, and an initial velocity dispersion of $0.14\pm 0.02$ km s-1.
Our results are entirely independent of model isochrones, details of stellar
evolution, and internal cluster dynamics, and the statistical precision in our
age derivation rivals that of the most precise age-dating techniques known
today, though our imperfect knowledge of the Milky Way potential and simple
spherical model for Theia 456 at birth add additional uncertainties. Using
posterior predictive checking, we confirm these results are robust under
reasonable variations to the Milky Way potential. Such low density structures
that are disrupted by the Galactic tides before virializing may be ubiquitous,
signifying that Theia 456 is a valuable benchmark for studying the dynamical
history of stellar populations in the Milky Way.

††software: AstroPy (Astropy Collaboration et al., 2022), gala (Price-Whelan,
2017), MINESweeper (Cargile et al., 2020), NumPy (Harris et al., 2020),
Matplotlib (Hunter, 2007), SciPy (Virtanen et al., 2020), Seaborn (Waskom,
2021), Pandas (pandas development team, 2020), Emcee (Foreman-Mackey et al.,
2013), Corner (Foreman-Mackey, 2016), IMF (https://github.com/keflavich/imf),
pytreegrav (https://github.com/mikegrudic/pytreegrav)

##  1 Introduction

The Gaia catalog, with over 1.8 billion sources in its third data release, has
transformed the study of stars and stellar populations within our own Milky
Way and its surrounding satellite galaxies (Gaia Collaboration et al., 2016,
2023; Babusiaux et al., 2023). Of these 1.8 billion sources, roughly 80% have
full five-dimensional (5D) astrometric data with unparalleled precision. In
particular, this wealth of high quality data has aided in uncovering numerous
coeval populations, such as open clusters (Castro-Ginard et al., 2022; Hao et
al., 2022), moving groups and young associations (Faherty et al., 2018; Ujjwal
et al., 2020), extended coronae and tails connected to known clusters
(Meingast et al., 2021), as well as linking newfound structures to well-
studied open clusters and moving groups (Gagné et al., 2021). Along with
revealing new structures, Gaia data has been instrumental in understanding
star cluster kinematics (Kuhn et al., 2019; Vasiliev & Baumgardt, 2021), Milky
Way substructure (Castro-Ginard et al., 2021) and star formation history (Kerr
et al., 2021; Briceño-Morales & Chanamé, 2023), and refining censuses for open
clusters and associations (Cantat-Gaudin et al., 2018; Luhman & Esplin, 2020).

Star clusters are extensively used as benchmarks to study stellar evolution
(e.g., Kalirai & Richer, 2010), reliant on the self-consistent age and
metallicity of their members. The need for accurate and precise age dating
methods has been discussed extensively in the literature, each with its own
benefits and drawbacks (for a review, see Soderblom, 2010). Model dependent
methods, like isochrone fitting, can reach uncertainties as large as $50\%$,
and rely heavily on accurate de-reddening corrections. Astroseismology, which
utilizes models of stellar pulsations, can provide ages accurate to $\sim
10\%$, but relies heavily on high-resolution photometry and spectroscopy (see
Li et al., 2022, for an example). Skumanich (1972) uncovered an empirical
relation between rotation period, age, and stellar activity by comparing
chromospheric Ca II H & K emission for stars in the Pleiades, Ursa Major,
Hyades, and the Sun. This discovery laid the foundation for gyrochronology
(Barnes, 2003), which relies on measuring the rotation period of a star,
$P_{\rm rot}$, which is now relatively easy to observe and can age-date single
stars. However, a poor understanding of the details of spin-down (e.g., the
temporary phase of stalling; Curtis et al., 2020) makes it especially
difficult to age-date lower mass K- and M-type stars. Dynamical age-dating,
utilizing the basic assumption that stars were likely born when they were most
compact in position space, is typically only applied to relatively young open
clusters (e.g., Jilinski et al., 2008). For these reasons, nearby common
origin structures, where all of the techniques detailed above and more can be
applicable, are the predominant sources for age-dated stars in our Milky Way.

It is widely believed that the bulk majority of star formation occurs in these
clustered environments through the gravitational collapse of cold molecular
gas (Lada & Lada, 2003). However, only a small fraction of these systems
survive as bound star clusters (e.g., Bressert et al., 2010; Kruijssen, 2012).
The majority of these stars are destined to disperse into the Milky Way field
(Goodwin & Bastian, 2006; Coronado et al., 2020). The evolution of stars from
their clustered birth-places to their eventual position within the Milky Way
field is not yet fully understood, though significant steps have been made
(see Helmi, 2020, for a review). For example, Kamdar et al. (2019) reproduce
the age–velocity relation, surface-density profile, and metallicity
distribution function in the Solar Neighborhood via dynamical modeling of
clustered star formation.

Ever increasing astrometric precision from Gaia has spawned the discovery of
many exciting new stellar structures that may probe this transition from
cluster to field star: stellar streams found within the Milky Way disk. Prior
to the advent of Gaia, only stellar streams within the Milky Way halo were
discussed extensively in the literature. Within the disk, disentangling
kinematically coherent streams from foreground and background stars has been
made possible via the combination of machine learning and clustering
algorithms with highly precise astrometry from Gaia. Kounkel & Covey (2019)
and Kounkel et al. (2020) apply HDBSCAN (Campello et al., 2013; McInnes et
al., 2017) to the Gaia DR2 catalog (Gaia Collaboration et al., 2018) to locate
phase space over-densities in search of disperse stellar structures; their
Theia catalog contains 8292 structures, often extending hundreds of parsecs in
length with hundreds or thousands of stars, all within 3 kpc of the sun. Hunt
& Reffert (2023) similarly apply HDBSCAN to the Gaia DR3 catalog to improve
the open cluster census, finding 739 high-probability open clusters previously
unmentioned in the literature. Numerous similar searches for open clusters in
the solar neighborhood and beyond have been performed recently in the
literature (e.g., He et al., 2022; Chi et al., 2023).

One such object from the Kounkel & Covey (2019) catalog, Theia 456, extends
nearly $120$ pc in length and 20∘ across the sky in the Northern Hemisphere.
This structure, which was originally discovered by Cantat-Gaudin et al. (2019)
and named COIN-Gaia-13, was shown to be a kinematically, chemically, and
gyrochronally coherent, thin-disk stream by Andrews et al. (2022).
Additionally, Bai et al. (2022) use the unsupervised clustering algorithm
pyUPMASK (Pera et al., 2021) to determine astrometric membership of the
extended structure, and find that it is likely in the process of dissolution
via differential Galactic rotation. Despite these follow-up studies, imprecise
radial velocity measurements among the subset of members observed by
additional instruments, $\sim 5\;\rm km~{}s^{-1}$ precision for Gaia (Gaia
Collaboration et al., 2018, 2023) and $\sim 6\;\rm km~{}s^{-1}$ precision for
LAMOST (Cui et al., 2012; Luo et al., 2015), precluded an in-depth dynamical
analysis.

In this work, we re-visit Theia 456, taking a much more detailed look at the
dynamics of the system. In Section 2 we present new radial velocities from
MMT’s Hectochelle, a precise multiobject echelle spectrograph, and discuss the
quality cuts we make to obtain our dynamical sample. We utilize this sample to
provide a qualitative dynamical history of Theia 456 by integrating orbits in
a Milky Way potential in Section 3, and quantitatively estimate Theia 456’s
age, birth scale, initial velocity dispersion, and initial phase space
position via a Bayesian approach in Section 4. In Section 5, we comment on
relevant assumptions made related to the Milky Way potential and Theia 456’s
self-gravity, and propose potential applications of stellar streams as probes
of star forming regions. We provide some concluding thoughts in Section 6.

Figure 1: An example spectrum from MMT-Hectochelle for a member of Theia 456
(Gaia DR3 191357324369704448) in black, and the corresponding MINESweeper fit
to the spectrum in purple. Residuals are shown in the top panel. Figure 2:
Measurement uncertainties for stellar radial velocities determined from
Hectochelle spectra. Top: Radial velocity precision plotted versus observed
Gaia DR3 $G$ magnitude. Bottom: Similarly, radial velocity precision shown
against the median signal-to-noise ratio of the Hectochelle spectrum. Figure
3: Left: Proper motion space of Theia 456 stars from the A22 catalog, color-
coded by right ascension, which is a proxy for position along the stream.
Circled stars (labelled as tail stars) appear systematically offset from the
observed positional gradient in proper motion space. Right: Parallax vs. right
ascension - tail stars have systematically lower parallaxes than other Theia
456 stars at the same right ascension. We cull these tail stars from our
catalog.

##  2 Spectroscopic Data

Data are obtained for members of Theia 456 with the Hectochelle multi-object
spectrograph (Szentgyorgyi et al., 2011) on the MMT. We use the RV31 filter
setup resulting in spectra with a resolution of $R\sim 32,000$ over a single
wavelength range of 5150–5300 $\rm\AA$. This instrument setup and wavelength
range has been shown to provide reliable stellar parameters (Cargile et al.,
2020). The spectra were processed using the standard reduction pipeline for
Hectochelle data (HSRED v2.1111https://bitbucket.org/saotdc/hsred).

We observe six fields along Theia 456, placing Hectochelle fibers on
kinematically selected cluster members. Each field was carefully selected to
simultaneously maximize the number of high-confidence, bright cluster members
(down to $G$ $\sim$ 18 mag), while providing information about stars across
its extent. We derive stellar parameters by modeling the spectra using the
MINESweeper code (Cargile et al., 2020), which has been used extensively to
model similar Hectochelle spectra (e.g., Conroy et al., 2019). Figure 1
presents an example spectrum and its corresponding MINESweeper fit.

###  2.1 Radial Velocities

Radial velocities are determined from spectral fitting for all Hectochelle
data as output from the MINESweeper stellar parameter pipeline. Barycentric
corrections were calculated using Astropy (Astropy Collaboration et al., 2022)
and applied for each observed MMT hectochelle observation. Our measured
velocities have a median precision of $\sim$0.1 km s-1, with precision only
rising above $\sim$1 km s-1 for the faintest, lowest signal-to-noise ratio
spectra (i.e., $G$ $>$ 17, $S/N<$ 10 per pixel; see Figure 2).

###  2.2 Our Sample

In order to obtain a useful dynamical sample, we first refine the catalog of
Theia 456 members reported in Andrews et al. (2022) (hereafter A22) in order
to accurately determine Theia 456’s median astrometric quantities (i.e.
$\mu_{\alpha^{\star}}\equiv\mu_{\alpha}\cos\delta$, $\mu_{\delta}$, $\varpi$)
for comparison with the stars with Hectochelle radial velocities. We omit
members that appear as a tail in proper motion space in the A22 catalog, with
systematically larger values of $\mu_{\alpha^{\star}}$ compared with other
stars within the same region of Theia 456, visible in Figure 3. These stars
also feature systematically lower parallaxes when compared with the rest of
the stream; on average, these tail stars sit $0.3$ mas below the median stream
parallax of $\varpi\approx 1.9\pm 0.1$ mas reported in A22. We cull these
discrepant members from our catalog, reducing its size from 362 to 345 stars.

We obtained a sample of 270 stars with highly precise radial velocities from
MMT’s Hectochelle spectrograph from six fields, shown in the left panel of
Figure 4. These fields contain members of Theia 456, suspected members,
potentially erroneous members, and background stars. We wish to define an
additional catalog of stars containing both Gaia astrometric data and
Hectochelle radial velocities; we undergo the following procedure to do so:

  1. 1.

Cull stars in the MMT-Hectochelle sample with large uncertainties in radial
velocity ($>$ 5 km s-1), $\mathrm{[Fe/H]}$ ($>$ 0.05 dex), or low signal-to-
noise ratios ($<$ 5). This reduces our sample from 270 to 92 stars, visible in
the right panel of Figure 4.

  2. 2.

Separate stars in our catalog into two separate lobes in R.A. and decl. with a
disperse bridge joining them. We designate stars with R.A.
$>\mathrm{5^{h}20^{m}}$ as one lobe (L1), and stars with R.A.
$<\mathrm{4^{h}52^{m}}$ as another (L2). Stars in between are labelled as
bridge stars. This qualitatively matches the observed structure of Theia 456.
The left panel of Figure 4 shows that four MMT pointings fit firmly into L1,
with the other two pointings in L2.

  3. 3.

From our catalog of 345 stars with Gaia data, we compute the median and
dispersion of $\mu_{\alpha^{\star}}$, $\mu_{\delta}$, and $\varpi$. We do this
by lobe in proper motion space, and across the whole structure in parallax,
because clear positional gradients exist in proper motion space, but not in
parallax:

     * •

$\mu_{\alpha^{\star}}^{\rm L1}$: $-3.8\pm 0.2$ mas yr-1

     * •

$\mu_{\alpha^{\star}}^{\rm L2}$: $-4.9\pm 0.2$ mas yr-1

     * •

$\mu_{\delta}^{\rm L1}$: $-1.7\pm 0.2$ mas yr-1

     * •

$\mu_{\delta}^{\rm L2}$: $-1.0\pm 0.2$ mas yr-1

     * •

$\varpi$: $2.0\pm 0.1$ mas

We restrict our sample of stars with Hectochelle radial velocities to stars
with proper motions no further than $0.6$ mas yr-1 ($\approx$ 1.5 km s-1 at
the distance of Theia 456) from the median quantity of its respective lobe,
and within $0.3$ mas ($\approx$ 90 pc at the distance of Theia 456) of the
median parallax. We show our astrometric selection criteria as boxes in Figure
5; proper motion in the leftmost panel and parallax in the middle panel.

  4. 4.

Restrict our sample to stars with metallicities no further than 0.36 dex away
from the median ([Fe/H]: $-0.07\pm 0.12$ dex, reported in A22 via analysis of
LAMOST DR5 metallicities). This selection criterion is shown in the rightmost
panel of Figure 5.

  5. 5.

Compute the median Hectochelle-measured radial velocity of the remaining stars
($-13.6\;\rm km~{}s^{-1}$). At this point, many non-members have been filtered
out, making the median radial velocity a reliable estimate of Theia 456’s bulk
radial velocity. We then cut stars further than $6\;\rm km~{}s^{-1}$ from the
median, similar to the radial velocity dispersion found in Tian 1 (Tian,
2020), a ‘stellar snake’ with similar extent across the sky and distance as
Theia 456. The middle panel of Figure 5 shows this selection criterion.

Figure 4: Left: MMT telescope pointings used to obtain Hectochelle radial
velocities, overlaid on our Theia 456 catalog. Right: 92 stars from the
Hectochelle sample with sufficiently precise radial velocities and
metallicities and high signal-to-noise ratios, compared to all Theia 456
members in proper motion space. Circle markers denote stars from our catalog,
while star markers denote stars from the Hectochelle sample with precise
radial velocities. Black star markers come from the MMT pointings in the
southern lobe (L1) of Theia 456, while blue star markers come from the
northern lobe (L2). Figure 5: Plot of our sample and culled stars in various
parameter spaces, with selection criteria shown as rectangles. Bold triangles
show our sample, while crosses represent stars culled via our refinement
process. In proper motion space (the left-most plot), we color-code stars and
selection criteria by lobe; black for L1 and blue for L2, and overplot our
sample on top of the entire Theia 456 catalog. Crosses that appear within the
selection criteria were culled via some criterion not shown on that plot
(i.e., crosses within the rectangle in radial velocity vs. parallax space were
culled via proper motion or metallicity criteria). The red triangle shows a
single star that passed all selection criteria, but was culled manually due to
an inconsistent parallax and large parallax uncertainty
($\sigma_{\varpi}=0.073$ mas, roughly three times larger than the median
parallax uncertainty in our sample).

Following this procedure, our dynamical sample contains 44 stars with
Hectochelle radial velocities, from 270 stars observed. We cut one star from
L1 manually, which has a particularly large uncertainty in parallax ($\sim$ 3
times larger than the typical parallax uncertainty in our sample), shown as a
red triangle in Figure 5. In total, this procedure provides us a sample
comprised of 43 qualifying stars with Hectochelle radial velocities, 36 in L1
and 7 in L2. We also remove 25 stars from the original A22 catalog that did
not pass our radial velocity or metallicity cuts. Three relatively faint stars
from the A22 catalog with MMT-Hectochelle radial velocities appear fully
consistent with the cluster, but were removed from the dynamical sample due to
poor spectral quality; they are included in our final catalog but their radial
velocities are not used in our dynamical analysis. Our final sample of stars
with precisely measured radial velocities is shown in the top panel of Figure
6, over-plotted on our Theia 456 catalog (referred to as T24 catalog
hereafter). Throughout this work, we will separately use both the sample of 43
members with radial velocities as well as all 321 identified members (one of
the stars in our sample observed by Hectochelle was not originally in the A22
catalog but nevertheless fit all selection criteria, bringing our Theia 456
catalog to 321 members). An abbreviated catalog of members, including stellar
parameters derived via MINESweeper for stars with precise Hectochelle radial
velocities, is presented in the appendix in Table 2.

##  3 The Two Lobes of Theia 456

Having selected our sample with Hectochelle-derived radial velocities, we
first analyze the compactness of Theia 456 at earlier times by using
astrometric positions, velocities, and radial velocities to integrate stars
backwards in time. We utilize gala, a python package dedicated to integrating
galactic orbits (Price-Whelan, 2017). With the default Astropy v4.0 (Astropy
Collaboration et al., 2022) parameters, we use the built-in MilkyWayPotential,
which combines an NFW halo potential (Navarro et al., 1996), a spherical bulge
and nucleus, and a multi-component disk model described in Bovy (2015), for
all integrations. Note that we consider variations to this model in Section
5.1.

To obtain a qualitative understanding of the evolution of the two lobes, we
define a metric $L$ which denotes the 3D spatial distance between the mean
$X$, $Y$, and $Z$ positions of L1 and L2 stars. The top panel of Figure 6
shows a schematic of our metric $L$ between the two lobes, and the bottom
panel shows $L$ as a function of backwards integration time. The two lobes
start roughly 100 pc apart and come closer in space until reaching a local
minima 50 Myr ago. They then expand away from each other, reaching a distance
of 200 pc at 175 Myr, before coming very close together ($\simeq$ 20 pc)
$\simeq$ 250 Myr ago. This latter minima at 250 Myr is consistent with the
isochrone and gyrochronology ages reported in A22.

Figure 6:  Top: Our sample of stars used for dynamical analysis (black
points), over-plotted on the entire Theia 456 catalog. Stars are segmented
into lobes—stars encircled in the southern portion are considered lobe 1 (L1)
stars, stars encircled in the northern portion are considered lobe 2 (L2)
stars. The dashed line shows the 3-dimensional distance $L$ between the two
lobes, taken as the distance between the mean X, Y, and Z Galactocentric
positions of L1 and L2 stars. Bottom: The L1–L2 distance $L$ as a function of
backwards integration time within a Milky Way potential. The two lobes nearly
overlap in space $\simeq$ 250 Myr ago.

Figure 7:  Slices in the $X$–$Y$ ($1$ kpc $\times$ $1$ kpc) and $X$–$Z$ ($1$
kpc $\times$ $1$ kpc) Galactocentric planes of the positions of Theia 456
stars, backwards integrated at various times under the MilkyWayPotential in
gala. Each column is a time-slice of backwards integrated orbits. Arrows in
the $X$–$Y$ plane point toward the Galactic center. Closed black circles show
L1 stars, open blue circles show L2 stars. The two lobes overlap 250 Myr ago.

In Figure 7, we show the positions of L1 and L2 stars in the Galactocentric
plane at various backwards integrated time-slices. Arrows in the $X$–$Y$ plane
point towards the Milky Way center. Here, the overlap between L1 and L2 stars
is visible at $50$ and $250$ Myr, and the expansion and contraction of our
metric $L$ can be seen as a function of time. In the $250$ Myr panel, we note
significant spread along the y-axis. Some fraction of this spread is due to
measurement uncertainty in Gaia parallax, and some is due to the natural
dynamics of Milky Way orbits; we address the shape drawn out by the error
ellipses in parallax and radial velocity in Section 4.2.

Theia 456 is an excellent candidate for the deployment of this method due to
the high quality kinematic data at both ends of the structure. By considering
the bulk motion of both lobes and not focusing on the individual positions of
stars integrated backwards, we are able to reduce the effect of observational
uncertainties from the data. Typical astrometric uncertainties for individual
Theia 456 stars in $\mu_{\alpha^{\star}}$ ($\simeq$ $0.13\;\rm km~{}s^{-1}$),
$\mu_{\delta}$ ($\simeq$ $0.07\;\rm km~{}s^{-1}$), $\varpi$ ($\simeq$ $0.023$
mas), and radial velocity ($\simeq$ $0.1\;\rm km~{}s^{-1}$) translate to a
positional uncertainty of roughly 50 pc when integrated backwards in time by
250 Myr, though parallax uncertainties dominate this spread. Theia 456’s
structure is particularly advantageous; its two distinct lobes, which are
quite far apart today ($\simeq$ 120 pc), provide a great deal of leverage when
compared with more compact open clusters.

##  4 Quantitative Dynamical History

###  4.1 Statistical Model

To perform a more robust constraint on the origin of Theia 456, we formulate a
Bayesian statistical model that defines the distribution of Theia 456 at birth
in an effort to best generate the distribution of 5D Gaia astrometry and
Hectochelle radial velocities seen today. We define our model parameters as:

$\psi_{0}:\\{X,Y,Z,V_{X},V_{Y},V_{Z},\sigma_{v},R,T\\}$

where $X,Y,Z$ and $V_{X},V_{Y},V_{Z}$ refer to the cluster’s mean position and
bulk velocity in a Galactocentric coordinate frame at birth, $\sigma_{v}$ is
an isotropic random velocity component given to each model star, $R$ is the
initial scale parameter, and $T$ is the cluster age. The model is populated as
a 3D Gaussian, such that $R$ refers to the standard deviation in position in
all three dimensions, and $\sigma_{v}$ refers to the standard deviation about
the bulk velocity in all three dimensions.

Figure 8: Posterior distributions from our MCMC model describing Theia 456’s
birth conditions ($\psi_{0}$), derived via maximization of equation 7.

Given measured observables
$x_{f,k}^{\,\prime}=(\alpha,\delta,\varpi,\mu_{\alpha^{\star}},\mu_{\delta},V_{\rm
rad})$ for the $k$-th star in Theia 456 in its final observed phase space
position, we seek to find the following posterior distribution:

$P(\psi_{0}|\\{x_{f}^{\,\prime}\\})\propto
P(\\{x_{f}^{\,\prime}\\}|\psi_{0})P(\psi_{0})$ (1)

where $\\{x_{f}^{\,\prime}\\}$ refers to the set of all $x_{f,k}^{\,\prime}$.
Since stars are observed independently, we can split our conditional
likelihood on the r.h.s.:

$P(\\{x_{f}^{\;\prime}\\}|\psi_{0})=\prod_{k}P(x_{f,k}^{\,\prime}|\psi_{0})$
(2)

Our model ($\psi_{0}$) defines a distribution of initial positions ($x_{i}$),
that are then uniquely mapped to final positions ($x_{f}$) via our chosen
Milky Way potential, such that $x_{f}=f(x_{i},\,T)$. We marginalize over the
distribution of initial and final positions:

$P(x_{f,k}^{\,\prime}|\psi_{0})=\int
P(x_{i},x_{f},x_{f,k}^{\,\prime}|\psi_{0})\,dx_{i}\,dx_{f}\equiv\int
I\,dx_{i}\,dx_{f}$ (3)

where the integrand simplifies as:

$\displaystyle I$ $\displaystyle=$ $\displaystyle
P(x_{f,k}^{\,\prime}|x_{i},x_{f},\psi_{0})P(x_{f}|x_{i},\psi_{0})P(x_{i}|\psi_{0})$
(4) $\displaystyle=$ $\displaystyle
P(x_{f,k}^{\,\prime}|x_{f})P(x_{f}|x_{i},T)P(x_{i}|\psi_{0})$

where we utilize the fact that $x_{f,k}^{\,\prime}$ depends only on $x_{f}$,
and $x_{f}$ depends only on $x_{i}$ and $T$. Since we uniquely map the
distribution of initial positions to final positions, we have
$P(x_{f}|x_{i},T)=\delta(x_{f}-f(x_{i},T))$. This simplifies equation 3 as:

$\begin{array}[]{l}P(x_{f,k}^{\,\prime}|\psi_{0})=\int
P(x_{f,k}^{\,\prime}|x_{f}^{\star})P(x_{i}|\psi_{0})\,dx_{i}\\\
P(x_{f,k}^{\,\prime}|\psi_{0})=\int\mathcal{N}(x_{f}^{\star};x_{f,k}^{\,\prime},\sigma^{\prime})P(x_{i}|\psi_{0})\,dx_{i}\end{array}$
(5)

where $x_{f}^{\star}=f(x_{i},T)$ in ICRS coordinates, and $\sigma^{\prime}$
refers to the observational uncertainties on $x_{f,k}^{\,\prime}$, which we
keep consistent across all $K=321$ stars. We discuss the reasoning for this
later in this section. Note that we include radial velocities for the 43 stars
in our sample outlined in Section 2.2, and only use Gaia astrometry for the
remaining 278 stars.

We approximate the integral in equation 5 as a sum such that:

$P(x_{f,k}^{\,\prime}|\psi_{0})\approx\frac{1}{N}\sum_{n}\mathcal{N}(x_{f,n}^{\star};x_{f,k}^{\,\prime},\sigma^{\prime})$
(6)

where $x_{f,n}^{\star}=f(x_{i,n},T)$ and $x_{i,n}$ is drawn from the
distribution $P({x_{i}|\psi_{0}})$.

Equation 1 simplifies, in $\log$-space to:

$\log
P(\psi_{0}|\\{x_{f}^{\,\prime}\\})\propto\sum_{k}\log\bigg{[}\frac{1}{N}\sum_{n}\mathcal{N}(x_{f,k}^{\,\prime};x_{f,n}^{\star},\sigma^{\prime})\bigg{]}+\log
P(\psi_{0})$ (7)

where we apply flat priors on all model parameters.

Figure 9: Theia 456, backwards integrated by $T=250$ Myr. Black contours show
the size and shape of the error ellipses of one star with typical
uncertainties, with full propagation of 5D astrometric and radial velocity
uncertainties. Green contours show error ellipses for the same star, but only
considering uncertainties in parallax. Red contours only consider
uncertainties in proper motion (both $\mu_{\alpha^{\star}}$ and
$\mu_{\delta}$). Typical radial velocity uncertainties produce negligibly
small error ellipses. Contours are drawn at the $2\sigma$ level containing 95%
of the data.

To find optimum model parameters, we employ $N=10,000$ model stars, and use
emcee (Foreman-Mackey et al., 2013) to maximize our posterior via a Markov
Chain Monte Carlo (MCMC) simulation with 64 walkers for 75,000 steps. To
minimize effects of stochasticity, we take the following steps. First, we
increase the uncertainties in our 5D Gaia astrometric data and Hectochelle
radial velocities, and apply these uncertainties consistently to all $K=321$
stars in our sample. We choose uncertainties of:

$\sigma^{\prime}=(\sigma_{\alpha},\sigma_{\delta},\sigma_{\varpi},\sigma_{\mu_{\alpha^{\star}}},\sigma_{\mu_{\delta}},\sigma_{V_{\rm
rad}})$
$\sigma^{\prime}=(1^{\circ},\,1^{\circ},\,0.1\;\rm{mas},\,0.15\;\rm{mas\;yr^{-1}},\,0.15\;\rm{mas\;yr^{-1}},\,1\;\rm{km~{}s^{-1}})$

so that the observable data points have overlapping error bars, but do not
significantly overshoot the width of the distribution in phase space. We also
ignore covariances in equation 7 because our newly defined uncertainties
dominate. This is a reasonable step to take because we are not interested in
precisely reproducing individual stars’ observables; rather we are interested
in finding model parameters that reproduce the size, shape, and location of
the distribution of stars in Theia 456. Attempting to fit model stars to the
extremely precise Gaia and Hectochelle uncertainties leads to numerical
convergence only with an unreasonably large value of $N$. To maintain detailed
balance, we draw one set of $\\{x_{i,n}\\}$ at the beginning of our
simulation, and shift and re-scale this distribution according to the new
model parameters at each step. This ensures that returning to the same model
parameters at a later step reproduces the same likelihood. To confirm that our
chosen value of $N$ is sufficiently large, we perform convergence tests and
find that our results are robust to different initial random draws with our
choice of $\sigma^{\prime}$ and $N$.

###  4.2 Velocity Dispersion, Birth Size, and Age of Theia 456

Figure 8 shows our resulting model parameters and their covariances. Of
particular importance to us, we find:

$\sigma_{v_{\rm opt}}=0.14\pm 0.02\;\;\rm km~{}s^{-1}$ $R_{\rm opt}=6\pm
1\;\;\rm pc$ $T_{\rm opt}=245\pm 3\;\;\rm Myr$

where all quoted uncertainties are statistical. Since we generate positions by
populating a Gaussian in 3D, the half-mass radius, assuming that mass is
randomly distributed at birth, is $R_{h,\rm opt}\approx 1.5\,R_{\rm
opt}\approx 9$ pc. We note that our model does not take any consideration of
natal gas, and that derived model parameters likely represent Theia 456 in its
post-embedded phase, potentially making the true age of the structure
marginally older (for an analysis on the timescale at which natal gas is
expelled, see Dinnbier & Walch, 2020). Our derived age is consistent with both
the isochrone and rotation period age derived in A22, as well as the ages of
individual stars derived via spectral fitting reported in Table 2.

In Figure 9, we show backwards integrated positions of Theia 456 at $T=250$
Myr, this is essentially a zoom-in of the 250 Myr panels in Figure 7. We note
significant spread ($\sim$ 500 pc) along the y-axis. Contours represent the
size and shapes of 2$\sigma$ error ellipses of a single star with typical
astrometric and radial velocity uncertainties for our sample. The black
contours, which span roughly half of the entire spread along the y-axis, take
into account full error propagation of the 5 astrometric properties and radial
velocity. Green and red contours show error ellipses due only to parallax and
proper motion, respectively. Typical uncertainties in Hectochelle-derived
radial velocity sweep out negligibly small error ellipses. The error ellipses
swept out by parallax uncertainties cover almost the entire length of the full
error ellipse, indicating that these uncertainties are the predominant cause
of the observed spread about the y-axis.

Figure 9 further showcases the difficulties of backwards modeling a cluster of
this age, as the propagation of errors leads to an unavoidable spread in birth
positions. With currently available data, deriving precise constraints about
the birth position and velocity of any individual star is fruitless.
Nevertheless, the close match between the error ellipses and Theia 456’s
elongated distribution, suggests the derived birth positions are an artifact
of the propagation of observational uncertainties under our Milky Way model;
Theia 456’s shape at birth is consistent with that of a sphere as described by
our generative forward model.

###  4.3 Posterior Predictive Checking

Figure 10: Theia 456 observables (black points), compared with forward
integrated model stars generated using the optimal model parameters shown in
Figure 8 (red points). In position and proper motion space, we show all 321
stars in Theia 456, and show only stars with Hectochelle-derived radial
velocities in parallax-radial velocity space.

While our model clearly favors an age of $T_{\rm opt}=245$ Myr, size of
$R_{\rm opt}=6$ pc ($R_{\rm h,opt}=9$ pc), and velocity dispersion of
$\sigma_{v_{\rm opt}}=0.14\;\rm km~{}s^{-1}$, we wish to check if these
parameters provide a reasonable representation of the observables when forward
integrated. To qualitatively assess the validity of these parameters, we
generate 10,000 model stars from the derived birth parameters, and integrate
them forward for a time $T_{\rm opt}$. These forward integrated model stars
are compared with observables in 6D phase space in Figure 10. We see good
agreement in all phases of 6D parameter space, and are able to reproduce the
bulk observables of Theia 456. Particularly, the observed scale and direction
of Theia 456’s morphology in position and proper motion space corresponds with
the predicted dispersion of a loosely bound open cluster via Galactic tides.
However, our model does not perfectly reproduce what we see on the sky today.
Theia 456 appears on the sky as two overdense lobes connected by a diffuse
bridge, whereas our generative model produces an elongated structure with a
monomodal distribution. This indicates that a spherically symmetric ball is
not a perfect physical description of Theia 456 at birth, and that initial
substructure may play a role; it may have been a more elongated structure, or
a more dumbbell-like structure to produce the disperse bridge seen. We leave a
more dedicated study, with a focus on initial substructure, to reproduce these
features for future work. We also note that we are able to accurately
reproduce the spread in all phase space dimensions but parallax, which may
indicate that uncertainties on Gaia parallaxes are somewhat underestimated for
our sample.

To check the validity of $\sigma_{v_{\rm opt}}$, $R_{\rm opt}$, and $T_{\rm
opt}$ individually we perform the following procedure. First, we integrate our
sample backwards by a time $T$ (which may be different than $T_{\rm opt}$)
which forms the initial position for our trial solution. Next, we generate a
cluster of 10,000 model stars with $\sigma_{\rm pos}=R$ at the mean birth
position of our backwards integrated sample and initialize their velocities
with the cluster bulk velocity and isotropic random component $\sigma_{v}$.
Finally, we integrate the stars drawn from this trial solution forward in time
for three different velocity dispersions, ages, and birth sizes.

Figure 11: Comparison of Theia 456’s phase space at the present epoch with
that of a sphere forward-integrated under the MilkyWayPotential, centered at
the mean position given by backwards integrating dynamical members to time
$T_{\rm opt}$, with a radius of $R_{\rm opt}$. Each contour contains 95% of
model stars forward integrated, drawn from a Maxwell–Boltzmann distribution at
different values of $\sigma_{v}$. Black and grey points denote our sample of
stars with Hectochelle radial velocities and the entire Theia 456 catalog,
respectively. Our model (the red contour) prefers a relatively low velocity
dispersion at birth; $\sigma_{v_{\rm opt}}=0.14\;\rm km~{}s^{-1}$, larger
velocity dispersions produce overly diffuse position and proper motion spaces.

In Figure 11, we compare the observed sky positions and proper motions of our
sample today with those produced by our generative model at different values
of $\sigma_{v}$, but still using our fiducial model parameters $R_{\rm opt}$
and $T_{\rm opt}$. Contours contain 95% of the simulated data. Even a modestly
larger velocity dispersion leads to an overly disperse stream, indicating that
Theia 456’s initial velocity dispersion is tightly constrained to be no more
than a few $\sim 0.1\;\rm km~{}s^{-1}$. We note that this velocity dispersion
is lower than what is typically found in young star-forming regions (for
example, NGC 1333 was found to have $\sigma_{v}=0.92\rm~{}km~{}s^{-1}$; Foster
et al., 2015).

Figure 12: Mock phase space of Theia 456 at the present epoch from various
initial conditions, compared with the observed phase space of our sample with
precise radial velocities (black points) and the entire Theia 456 catalog
(gray points). Contours contain 95% of the data and are drawn at different
birth ages, while columns refer to different birth sizes. Red contours in the
middle column most closely represent our best fit model.

Finally, Figure 12 shows the positions and proper motions predicted by our
forward models at different initial birth times and scales, at our derived
velocity dispersion $\sigma_{v_{\rm opt}}$, compared with the observed phase
space of Theia 456 stars today. Each column represents a different
initialization size, while each contour represents a different age. A $10$ pc
ball clearly produces a too-extended stellar distribution in proper motion
space, independent of birth time. For the older age of $300$ Myr, the same is
true, as the position distribution is too extended compared with our observed
sample. In general, sizes and ages larger than our optimal model parameters
are overly influenced by Galactic tides and appear more disperse in phase
space than what is observed in Theia 456 today. Decreasing the age and size
leads to the opposite effect; Galactic tides are not strong enough to
reproduce the size of the distribution observed in phase space today. We
conclude that indeed, the distribution of stars in position and proper motion
space is a sensitive tracer of the cluster’s birth conditions, as our
relatively simple model accurately reproduces the main features of Theia 456
as observed today.

##  5 Discussion

###  5.1 Sensitivity to the Milky Way Potential

Up until this point, we have only considered the effect of observational
uncertainties and how they propagate backwards in time as a function of the
Milky Way potential. To further test the robustness of our model predictions,
we also consider variations in the potential itself based on our imperfect
knowledge of its functional form. Specifically, we modify the characteristic
mass $M$ in the Miyamoto-Nagai (Miyamoto & Nagai, 1975) disk-component of the
MilkyWayPotential, shown below:

$\Phi(R,z)=-\frac{GM}{\sqrt{R^{2}+(a+\sqrt{z^{2}+b^{2}})^{2}}}$ (8)

and compute the differences in model outputs in $\sigma_{v}$, $R$, and $T$
under these modifications. The default value is $M_{\rm def}=6.8\times
10^{10}$ $\rm M_{\odot}$. We consider only variations in the disk component of
the model, and not the bulge, nucleus, or halo, as we find that the disk
component of the potential is the dominant dynamical driver for a thin-disk
structure like Theia 456. We consider variations in $M$ such that the computed
circular velocity of the sun, at an assumed Galactocentric radius of $8.12$
kpc, is reasonably varied. We do not consider variations in $a$ and $b$, as
variations in $a$ produce very similar effects as variations in $M$, and even
large variations in $b$ produce negligible differences in model parameters. At
the fiducial disk parameters in gala, $V_{\rm circ}=231.52\;\rm km~{}s^{-1}$;
we modify $M$ in equation 8 by 20% in either direction, such that $V_{\rm
circ}$ ranges from $[219,244]$ $\rm km~{}s^{-1}$, a reasonable encapsulation
of the current constraints on $V_{\rm circ}$ (Eilers et al., 2019). Table 1
shows the outputted model parameters due to these modifications. Our model
constraints on cluster age are quite robust; reasonable variations to the disk
component of the potential produce only a modest $5\%$ maximum change in $T$.
The velocity dispersion and radius are less robust to variations in the
Galactic potential; $\sigma_{v}$ and $R$ vary by up to $64\%$ and $50\%$,
respectively.

Table 1: Model Outputs for Variations in MilkyWayPotential

$\frac{\displaystyle M}{\displaystyle M_{\rm def}}$ |  $V_{\rm circ}$ (km s-1) |  $\sigma_{v}$ (km s-1) |  $R$ (pc) |  $T$ (Myr)
---|---|---|---|---
0.8 | 219 | $0.23\pm 0.01$ | $3\pm 1$ | $257\pm 2$
1 | $\mathbf{232}$ | $\mathbf{0.14\pm 0.02}$ | $\mathbf{6\pm 1}$ | $\mathbf{245\pm 3}$
1.2 | 244 | $0.11\pm 0.01$ | $5\pm 1$ | $252\pm 2$

††footnotetext: Quoted statistical uncertainties are 1$\sigma$. $V_{\rm circ}$
is the circular velocity of the sun at that given disk mass.

We only consider variations to an axisymmetric and time-independent potential;
more complex considerations, such as perturbations from the spiral arms and
giant molecular clouds (GMCs), require a more focused work on the modeling of
the Milky Way potential. Such baryonic substructure in the disk has been shown
to disrupt passing star clusters (Gieles et al., 2006, 2007; Gieles & Renaud,
2016). These perturbative features, while not necessary to reproduce the bulk
distribution seen today, may be the origin of the particular features of Theia
456, including its diffuse bridge. The bridge is somewhat reminiscent of
observed gaps in the GD-1 and Palomar 5 halo streams which have been
previously explained through the gravitational interaction of stars with Milky
Way substructure (Bonaca et al., 2019; Pearson et al., 2017). We leave the
issue of whether any such encounter with substructure has a significant effect
on the derived dynamical age for future work.

###  5.2 The Stellar Population of Theia 456

In order to obtain an order of magnitude estimate for the total number of
stars and stellar mass of Theia 456, we extrapolate via a Kroupa initial mass
function (Kroupa, 2001). To solve for the normalization constant, we use the
fact that Gaia is essentially complete for $G<17$ (van Leeuwen et al., 2022).
We solve for the normalization constant $N_{0}$ using:

$n_{\rm obs}=N_{0}\int_{M_{\rm min}}^{\infty}m^{\alpha}dm$ (9)

where $n_{\rm obs}$ is the number of Theia 456 stars observed with $G<17$,
$M_{\rm min}$ is the minimum mass in this fully covered region, and $\alpha$
is the Kroupa IMF power law index ($\alpha=-0.3$ for
$0.01\;\text{M}_{\odot}<m<0.08$ M⊙, $\alpha=-1.3$ for
$0.08\;\text{M}_{\odot}<m<0.5$ M⊙, $\alpha=-2.3$ for $m>0.5$ M⊙).

We utilize the FLAME (Final Luminosity and Age Estimator) module of the
astrophysical parameters inference system (Apsis) chain (Bailer-Jones et al.,
2013; Fouesneau et al., 2023) to determine $M_{\rm min}=0.57$ M⊙, simply
taking the smallest mass from the FLAME module for Theia 456 stars with
$G<17$. Having solved for our normalization constant $N_{0}$, we estimate the
total number and mass of Theia 456 by invoking continuity of the IMF at
crossover points to solve for normalization constants in those regimes. Under
these assumptions, we find $n_{\rm total}\simeq 2100$ and $m_{\rm total}\simeq
900\;\rm M_{\odot}$. Note that this mass is roughly twice the mass reported by
Bai et al. (2022), as they only consider stars brighter than $G=18$. At our
derived mass and half-mass radius, our results suggest that Theia 456 was a
low density ($\approx 0.2\;\rm M_{\odot}~{}pc^{-3}$) open cluster at birth. As
a comparison, the Pleiades has a current day mass ($\approx 740$ M⊙) and half-
mass radius ($\approx 3.66$ pc), leading to a characteristic density of
$\approx 1.8$ $\rm M_{\odot}~{}pc^{-3}$ (Moraux et al., 2004). Given the
possibility that the Pleiades is expanding (Converse & Stahler, 2010), its
density at birth may have been even higher. Although we account for the
magnitude limitations of Gaia, incompleteness due to limitations of the
HDBSCAN algorithm used to identify members of Theia 456, particularly near the
edges of the structure, may drive this density to be slightly higher; we
discuss these limitations further in Section 5.5.

###  5.3 Self-Gravity

In our derivation of model parameters to describe Theia 456 at birth, we have
so far ignored any potential effects of self-gravity between stars. The
predominant driver of the early dynamical evolution of stars after their gas
embedded phase comes from the potential of the Milky Way and the cluster’s own
self-gravitational interactions. Even when gas is present, it does not play a
major role in affecting the dynamical state of the young cluster (Sills et
al., 2018), though its expulsion (i.e., the infant mortality paradigm; Lada &
Lada, 2003; Pellerin et al., 2007) is dynamically important, especially in the
formation of tidal tails (Dinnbier & Kroupa, 2020). However, in loosely bound
open clusters and stellar streams, self-gravity may not play a major role in
the bulk dynamical evolution of the system, due to large inter-stellar
separations.

After its initial embedded phase, the acceleration that a star feels is
equivalent to:

$\vec{a}_{\rm star}=-\nabla\Phi_{\rm MilkyWay}-\nabla\Phi_{\rm Cluster}$ (10)

As a simple order-of-magnitude estimate, we take $\Phi_{\rm Cluster}$ to be
that of a Plummer sphere (Plummer, 1911), i.e. $\Phi_{\rm
Cluster}(r)=-GM/\sqrt{r^{2}+a^{2}}$ where $M$ is the total mass of the
cluster, $r$ is the distance to the cluster center, and $a$ is the
characteristic scale radius. The maxima in the potential gradient
$\nabla\Phi_{\rm Cluster}(r)$ occurs at $r=a/\sqrt{2}$, such that
$\nabla\Phi_{\rm Cluster}(a/\sqrt{2})\approx GM/a^{2}$. For our cluster mass
($M\approx 900$ M⊙) and scale ($a\approx 7$ pc, derived using the fact that
$R_{h,\rm opt}\approx 1.3a$), the resulting characteristic acceleration due to
the cluster’s self gravity is $\sim 10^{-10}$ cm s-2. The gradient in our
model Milky Way potential at the predicted birth position of Theia 456 is
$\sim 10^{-8}$ cm s-2.

To further justify our exclusion of self-gravity, we initialize a model run
with our fiducial parameters found in Section 4.2 and run it until present day
without the inclusion of self-gravity. We populate the run with $n_{\rm
total}\simeq 2100$ stars, and assign each star a mass drawn from a Kroupa
IMF222https://github.com/keflavich/imf such that $m_{\rm total}\simeq 900$ M⊙.
At each timestep, we separately compute the acceleration each star feels due
to the Milky Way potential, and due to the sum of gravitational forces from
other stars.333https://github.com/mikegrudic/pytreegrav In Figure 13, we
compare the two accelerations, (taking the median value as a representation of
the gravitational interactions between stars). The typical acceleration due to
the Milky Way potential is orders of magnitude higher than the acceleration
due to self-gravity throughout the forward run, indicating that self-gravity
can safely be ignored when studying the global dynamical properties of Theia
456.

Figure 13: Comparison of the typical acceleration that a model star in Theia
456 feels solely due to the Milky Way potential (black line) and due to the
sum of Nbody interactions (red lines). Shown is the median (for 2166 stars)
magnitude of the acceleration at each timestep. The Milky Way potential is the
dominant driver for the acceleration of stars in Theia 456.

While we can ignore self-gravitational interactions in our orbital
integrations, we again do not rule out that our optimal model parameters may
describe Theia 456 after dynamical relaxation and not strictly at birth; at
these earlier times, dynamical interactions, especially those between massive
binaries, could indeed be important for the cluster’s evolution (Torniamenti
et al., 2021). We also note that Figure 13 shows only the typical
accelerations due to self-gravitational interactions. Close encounters between
individual stars may be dynamically important to accurately reproduce the
position and velocity of individual stars, which is outside of the scope of
this work.

###  5.4 Tracing Birth Conditions using Stellar Positions Today

The idea that stars typically form in bound, dense clusters has been somewhat
challenged as of late (e.g., Gouliermis, 2018; Ward et al., 2020), favoring
the hierarchical collapse of gas across a wide range of surface densities.
This picture implies that the invocation of gas expulsion and/or dynamical
relaxation may not be necessary to explain the wealth of observed, unbound
gas-rid associations within the Milky Way (Wright et al., 2022). Under this
paradigm, for sufficiently low density structures, the Milky Way tides would
play a critical role in the dissolution of the cluster. If one can constrain
the size and initial velocity distribution of a stellar stream, metallicity
gradients in position today may be indicative of metallicity gradients in the
parent molecular cloud. Stellar streams could then provide a probe of the
giant molecular clouds from which stars are born. This motivates high-
resolution spectroscopic observations to obtain precise chemical abundances of
stars to further explore this correlation.

Stellar streams may also pose a powerful tool to study mass segregation;
whether mass segregation is primordial or a dynamical effect is still unclear
(Parker et al., 2015; Dib et al., 2018). If one can constrain the initial size
and velocity distribution of a stellar population, then the presence (or lack
thereof) of mass segregation in streams could help answer this question. In
the case of Theia 456, our model predicts $\sigma_{v_{\rm opt}}=0.14\;\rm
km~{}s^{-1}$, while the virial theorem for a cluster of $m_{\rm total}\approx
900$ M⊙ and $R_{v}\approx 12$ pc (utilizing the fact that $R_{v}\approx 1.7a$
for a Plummer sphere) predicts a velocity dispersion of $\sigma_{v}\approx
1\;\rm km~{}s^{-1}$, indicating that Theia 456 was sub-virial at birth. If one
can quantify that the stream disperses on a timescale shorter than its
crossing time, meaning the typical star does not have time to orbit the
cluster’s center of mass before the cluster begins to desolve, then the
presence of mass segregation points towards the primordial paradigm. Precise
and accurate mass determinations, preferably for a number of stellar streams,
would be necessary to test this hypothesis.

In Figure 14, we explore the potential correlation between birth position and
velocity and the observed phase space today. Using our derived age and birth
scale, we forward-integrate 1,000 model stars in the manner described above
for three different velocity dispersions ($\sigma_{v}=0.14,\;0.25,\;0.50\;\rm
km~{}s^{-1}$), and color-code points by their cylindrical Galactocentric birth
position (top 3 rows) and birth velocities relative to the bulk velocity of
the cluster (bottom 3 rows). For sufficiently low thermal velocities, we see a
clear correlation between initial Galactocentric radius and the position on
the sky today, visible in the top left panel. Somewhere between
$\sigma_{v}=0.25\;\rm km~{}s^{-1}$ and $\sigma_{v}=0.50\;\rm km~{}s^{-1}$,
this gradient in initial position begins to obscure; and a correlation between
the initial angular velocity and position on the sky today emerges, visible in
the bottom middle panel. Qualitatively, this correlation is expected. For low
thermal velocities, the tidal disruption timescale, which we approximate as
$\Delta t_{\rm tidal}=\sqrt{2R_{h,\rm opt}/a_{\rm tidal}}\approx 2$ Myr, is
significantly shorter than the crossing timescale, which at our derived half-
mass radius and velocity dispersion is $\approx 50$ Myr. The opposite occurs
for higher thermal velocities; stars with large relative angular velocities
move quickly to one side of the cluster before tidal forces can disrupt it. At
our derived birth parameters (top row), we see agreement with the claim made
by Bai et al. (2022); it does appear that Theia 456 is dissolving into the
field via differential rotation caused by the Galactic tides.

Figure 14: Mock stars integrated forward, color-coded by their initial
positions (top 3 rows) and velocities (bottom 3 rows) in Galactocentric
cylindrical coordinates, with dynamical members over-plotted as black points.
Rows are sorted by the Maxwell-Boltzmann initialization velocity $\sigma_{v}$.
Positions are measured relative to the Galactic center, and velocities are
measured relative to the bulk cluster velocity.

###  5.5 Future Dispersion of Theia 456

Ever increasing Gaia precision allows us to detect more and more disperse
string-like stellar structures across the sky, which poses the question in
relation to Theia 456: at what point will this structure become so dispersed
in phase space that it is undetectable as a stellar structure of common
origin? Figure 15 shows model streams initialized at different ages (all with
$\sigma_{v}=0.14\;\rm km~{}s^{-1}$, $R=6$ pc) integrated forward in time;
their dispersions in position, proper motion, and parallax increase greatly
with age. Even at large ages, the proper motion dispersion is relatively
small, though an emergent ‘tail’ appears to branch away from the bulk
distribution. Most notably, the dispersion in proper motion for a 350 Myr old
model Theia 456 still falls below the threshold adopted in Cantat-Gaudin &
Anders (2020) of:

$\begin{array}[]{c}\sigma_{pm}=\sqrt{\sigma_{\mu_{\alpha^{\star}}}^{2}+\sigma_{\mu_{\delta}}^{2}}\leq
5\sqrt{2}\frac{\varpi}{4.74},\;\varpi>1\;\rm{mas}\end{array}$ (11)

However, the spread across the sky is so large, $\sim 100^{\circ}$ in right
ascension and $\sim 80^{\circ}$ in declination, that the field of view in
position space for one’s chosen clustering method or algorithm would be
impractical and overwhelmingly filled with contaminants. We conclude that a
stream like Theia 456 will likely be undetectable given our current cluster
finding methods in $\lesssim 100$ Myr. One should be careful when generalizing
these predictions to other stream-like objects, as their diffusion in phase
space due to tidal shredding is greatly affected by their initial velocity
distribution, age, and size, as seen in Figures 11 and 12.

Figure 15: Mock stars (shown as colored points) with different hypothetical
ages, denoted in the top left corner of the top panels, forward-integrated,
and plotted under our catalog in position and proper motion space. The proper
motion dispersion, defined as
$\sigma_{pm}=\sqrt{\sigma_{\mu_{\alpha^{\star}}}^{2}+\sigma_{\mu_{\delta}}^{2}}$,
is shown for each model stream. The inset shows the density distribution of
parallax, and its dispersion in the legend, for each model stream; compared
with the entire Theia 456 catalog in solid black.

The distinct tail emerging in proper motion space shown in Figure 15 points to
limitations within our application of clustering algorithms, like HDBSCAN, to
blindly classify disperse stellar structures. Currently, algorithms
agnostically search for phase space over-densities, and do not consider any
physically motivated phase space morphologies. Application of the dynamical
modelling outlined in this work could be used to refine catalogs generated by
data-driven clustering algorithms. Our results imply that more robust
algorithms, trained on physically generated stellar streams in the process of
tidal disruption, could provide significant improvement over the current
generation of stream-finding algorithms by directly searching for physically
motivated phase space over-densities.

##  6 Conclusions

In this work, we combine Gaia astrometric quantities with MMT-Hectochelle-
derived radial velocities to present a detailed dynamical history of Theia
456, a stellar stream located in the thin disk of the Milky Way. We obtain
observations in six fields from MMT along Theia 456’s extent and derive
precise radial velocities for stars within these fields. We produce a catalog
of 321 stars in Theia 456: 43 with full 6D phase space information and the
remaining 278 with 5D Gaia astrometry. We quantitatively back-track Theia 456
stars with precise radial velocities under a model Milky Potential and find
that its two distinct lobes overlap $\sim 250$ Myr in the past.

To quantify Theia 456’s dynamical evolution, we formulate a Bayesian model,
generate an idealized birth distribution of Theia 456 stars, integrate this
distribution forward under a model Milky Way potential, and compute the
posterior distribution for all 321 stars in Theia 456. We find that Theia
456’s origin is well described as a spherically symmetric ball of stars of age
$245\pm 3$ Myr, a scale of $6\pm 1$ pc (half-mass radius of $9\pm 2$ pc), and
a velocity dispersion of $0.14\pm 0.02\;\rm km~{}s^{-1}$, where all quoted
uncertainties are statistical. We find that parallax uncertainty is the
limiting factor in our statistical model, which should improve with future
Gaia data releases. We also vary model parameters in our assumed Milky Way
potential and find that reasonable variations in the potential contribute an
uncertainty of $5\%$ in age, $50\%$ in radius, and $64\%$ in velocity
dispersion. These uncertainties, especially those on the age, are considerably
smaller than those found in A22 via isochrone modeling and gyrochronology
methods, indicating that dynamical age dating is a very powerful tool,
especially for disperse common-origin structures with high quality kinematics
at both ends.

We then take our generative model and forward integrate those stars to
qualitatively assess our model’s goodness of fit. A simple spherical model for
Theia 456 at birth can reproduce the size and shape of the distribution seen
in phase space today. It does not reproduce some of Theia 456’s unique
features, such as the disperse bridge seen in position space. Reproducing such
features may be an observational tool to probe Milky Way disk substructure.
Theia 456’s appearance as two distinct lobes may also be indicative of
substructure at birth, although other such processes, such as interaction with
a passing giant molecular cloud and the Milky Way’s spiral arms, may be
responsible.

We also derive an order of magnitude estimate of the initial stellar mass and
number of Theia 456. By assuming a Kroupa initial mass function, we find
$m_{\rm total}\approx 900$ M⊙ and $n_{\rm total}\approx 2100$ stars. At our
derived birth scale, Theia 456 formed as a low density open cluster ($\approx
0.2\rm\;M_{\odot}~{}pc^{-3}$). We do not expect, especially at later times in
its evolution, for gravitational interactions between stars in the cluster to
play a significant role in the bulk evolution of the cluster, though they may
be dynamically important for individual stars.

In Figure 14, we show that given a sufficiently low thermal velocity, as found
for Theia 456, the location of stars at birth relative to the Galactic center
may be well correlated with their position at later times. This may present a
way to probe chemical inhomogeneity within molecular clouds, and also mass
segregation. This strategy may be the best method of quantifying inhomogeneity
within molecular clouds.

We also integrate model stars forward in time from the present day, and find
that the Milky Way tides are able to shred these structures quite quickly.
Within another $\sim 100$ Myr, the extent of Theia 456 should roughly double
in right ascension, making this structure unlikely to be identified by the
current generation of clustering algorithms.

Dynamical modelling may be a powerful tool to refine catalogs produced by
data-driven, physics-blind clustering algorithms. It may also be instrumental
in pushing the categorization of stellar streams forward by generating new
algorithms, trained on generated stellar streams, that search for physically
motivated phase space over-densities. Such algorithms may be instrumental in
finding additional disperse structures, like Theia 456, that drive our
knowledge of the Milky Way’s dynamical evolution and star formation history
further forward.

## Acknowledgement

The authors thank Carl Rodriguez and David Hogg for useful conversations. J.C.
acknowledges support from the Agencia Nacional de Investigación y Desarrollo
(ANID), via Proyecto Fondecyt Regular 1231345, and from ANID BASAL project
CATA2-FB210003. Observations reported here were obtained at the MMT
Observatory, a joint facility of the Smithsonian Institution and the
University of Arizona. This paper uses data products produced by the CfA OIR
Telescope Data Center, supported by the Smithsonian Astrophysical Observatory.

Table 2: Theia 456 Catalog

Gaia DR3 ID | $S/N$ | $\frac{V_{\rm rad}}{\rm km~{}s^{-1}}$ | $\frac{V_{\rm rot}}{\rm km~{}s^{-1}}$ | $\frac{\rm Dist.}{\rm pc}$ | $\frac{\varpi}{\rm mas}$ | $\frac{A_{V}}{\rm mag}$ | $\frac{T_{\rm eff}}{\rm K}$ | $\frac{M}{\rm{{M_{\odot}}}}$ | $\frac{\rm[Fe/H]}{\rm dex}$ | $\frac{\rm[\alpha/H]}{\rm dex}$ | $\log\big{(}\frac{g}{\rm cm\;s^{-2}}\big{)}$ | $\log\big{(}\frac{R}{\rm{{R_{\odot}}}}\big{)}$ | $\log\big{(}\frac{\rm Age}{\rm yr}\big{)}$ | $\log\big{(}\frac{L}{\rm{{L_{\odot}}}}\big{)}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
194453957134786176 | 114 |  $-$13.77 | 47.9 | 500 | 1.98 | 0.36 | 6459 | 1.243 | 0.06 |  $-$0.08 | 4.355 | 0.088 | 9.0 | 0.37
194450491096608000 | 18 |  $-$11.90 | 4.5 | 580 | 1.72 | 0.57 | 4774 | 0.746 |  $-$0.27 | 0.16 | 4.660 |  $-$0.144 | 7.8 |  $-$0.62
194445062257874048 | 30 |  $-$13.62 | 0.4 | 530 | 1.88 | 0.45 | 5063 | 0.849 |  $-$0.05 | 0.08 | 4.620 |  $-$0.104 | 8.1 |  $-$0.44
194400944354088960 | 20 |  $-$13.76 | 0.9 | 530 | 1.90 | 0.51 | 4878 | 0.794 |  $-$0.13 | 0.11 | 4.632 |  $-$0.120 | 8.4 |  $-$0.53
194387097379628672 | 32 |  $-$13.97 | 7.6 | 500 | 1.99 | 0.56 | 5489 | 0.933 |  $-$0.05 | 0.07 | 4.572 |  $-$0.074 | 8.7 |  $-$0.24
… | … | … | … | … | … | … | … | … | … | … | … | … | … | …
191391718467413504 | 49 |  $-$13.90 | 7.9 | 510 | 1.95 | 0.33 | 5612 | 0.960 |  $-$0.01 | 0.02 | 4.566 |  $-$0.069 | 8.2 |  $-$0.19
191365566412178048 | 138 |  $-$13.43 | 22.3 | 500 | 2.00 | 0.32 | 6645 | 1.304 | 0.04 |  $-$0.06 | 4.308 | 0.123 | 9.0 | 0.49
191362890644144128 | 10 |  $-$14.25 | 4.9 | 510 | 1.97 | 1.09 | 4635 | 0.774 |  $-$0.16 | 0.22 | 4.635 |  $-$0.121 | 8.7 |  $-$0.62
191357324369704448 | 32 |  $-$13.31 | 0.3 | 530 | 1.90 | 0.35 | 5211 | 0.867 |  $-$0.07 | 0.05 | 4.607 |  $-$0.098 | 8.5 |  $-$0.38
190222319132712832 | 21 |  $-$13.05 | 0.6 | 550 | 1.80 | 0.78 | 4952 | 0.811 |  $-$0.10 | 0.08 | 4.640 |  $-$0.121 | 7.8 |  $-$0.51

††footnotetext: An abbreviated data table for our catalog of stars in Theia
456. Stars listed above come from the subset of stars with MMT-Hectochelle
spectra and derived stellar parameters. A machine readable format is available
for download.

## References

  * Andrews et al. (2022) Andrews, J. J., Curtis, J. L., Chanamé , J., et al. 2022, The Astronomical Journal, 163, 275, doi: 10.3847/1538-3881/ac6952
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Babusiaux et al. (2023) Babusiaux, C., Fabricius, C., Khanna, S., et al. 2023, A&A, 674, A32, doi: 10.1051/0004-6361/202243790
  * Bai et al. (2022) Bai, L., Zhong, J., Chen, L., Li, J., & Hou, J. 2022, Research in Astronomy and Astrophysics, 22, 055022, doi: 10.1088/1674-4527/ac60d2
  * Bailer-Jones et al. (2013) Bailer-Jones, C. A. L., Andrae, R., Arcay, B., et al. 2013, A&A, 559, A74, doi: 10.1051/0004-6361/201322344
  * Barnes (2003) Barnes, S. A. 2003, ApJ, 586, 464, doi: 10.1086/367639
  * Bonaca et al. (2019) Bonaca, A., Hogg, D. W., Price-Whelan, A. M., & Conroy, C. 2019, ApJ, 880, 38, doi: 10.3847/1538-4357/ab2873
  * Bovy (2015) Bovy, J. 2015, ApJS, 216, 29, doi: 10.1088/0067-0049/216/2/29
  * Bressert et al. (2010) Bressert, E., Bastian, N., Gutermuth, R., et al. 2010, MNRAS, 409, L54, doi: 10.1111/j.1745-3933.2010.00946.x
  * Briceño-Morales & Chanamé (2023) Briceño-Morales, G., & Chanamé, J. 2023, MNRAS, 522, 1288, doi: 10.1093/mnras/stad608
  * Campello et al. (2013) Campello, R., Moulavi, D., & Sander, J. 2013, in Advances in Knowledge Discovery and Data Mining, Vol. 7819, 160–172, doi: 10.1007/978-3-642-37456-2_14
  * Cantat-Gaudin & Anders (2020) Cantat-Gaudin, T., & Anders, F. 2020, A&A, 633, A99, doi: 10.1051/0004-6361/201936691
  * Cantat-Gaudin et al. (2018) Cantat-Gaudin, T., Jordi, C., Vallenari, A., et al. 2018, A&A, 618, A93, doi: 10.1051/0004-6361/201833476
  * Cantat-Gaudin et al. (2019) Cantat-Gaudin, T., Krone-Martins, A., Sedaghat, N., et al. 2019, A&A, 624, A126, doi: 10.1051/0004-6361/201834453
  * Cargile et al. (2020) Cargile, P. A., Conroy, C., Johnson, B. D., et al. 2020, ApJ, 900, 28, doi: 10.3847/1538-4357/aba43b
  * Castro-Ginard et al. (2021) Castro-Ginard, A., McMillan, P. J., Luri, X., et al. 2021, A&A, 652, A162, doi: 10.1051/0004-6361/202039751
  * Castro-Ginard et al. (2022) Castro-Ginard, A., Jordi, C., Luri, X., et al. 2022, A&A, 661, A118, doi: 10.1051/0004-6361/202142568
  * Chi et al. (2023) Chi, H., Wang, F., Wang, W., Deng, H., & Li, Z. 2023, ApJS, 266, 36, doi: 10.3847/1538-4365/accb50
  * Conroy et al. (2019) Conroy, C., Naidu, R. P., Zaritsky, D., et al. 2019, ApJ, 887, 237, doi: 10.3847/1538-4357/ab5710
  * Converse & Stahler (2010) Converse, J. M., & Stahler, S. W. 2010, MNRAS, 405, 666, doi: 10.1111/j.1365-2966.2010.16505.x
  * Coronado et al. (2020) Coronado, J., Rix, H.-W., Trick, W. H., et al. 2020, MNRAS, 495, 4098, doi: 10.1093/mnras/staa1358
  * Cui et al. (2012) Cui, X.-Q., Zhao, Y.-H., Chu, Y.-Q., et al. 2012, Research in Astronomy and Astrophysics, 12, 1197, doi: 10.1088/1674-4527/12/9/003
  * Curtis et al. (2020) Curtis, J. L., Agüeros, M. A., Matt, S. P., et al. 2020, ApJ, 904, 140, doi: 10.3847/1538-4357/abbf58
  * Dib et al. (2018) Dib, S., Schmeja, S., & Parker, R. J. 2018, MNRAS, 473, 849, doi: 10.1093/mnras/stx2413
  * Dinnbier & Kroupa (2020) Dinnbier, F., & Kroupa, P. 2020, A&A, 640, A84, doi: 10.1051/0004-6361/201936570
  * Dinnbier & Walch (2020) Dinnbier, F., & Walch, S. 2020, MNRAS, 499, 748, doi: 10.1093/mnras/staa2560
  * Eilers et al. (2019) Eilers, A.-C., Hogg, D. W., Rix, H.-W., & Ness, M. K. 2019, ApJ, 871, 120, doi: 10.3847/1538-4357/aaf648
  * Faherty et al. (2018) Faherty, J. K., Bochanski, J. J., Gagné, J., et al. 2018, ApJ, 863, 91, doi: 10.3847/1538-4357/aac76e
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24, doi: 10.21105/joss.00024
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Foster et al. (2015) Foster, J. B., Cottaar, M., Covey, K. R., et al. 2015, ApJ, 799, 136, doi: 10.1088/0004-637X/799/2/136
  * Fouesneau et al. (2023) Fouesneau, M., Frémat, Y., Andrae, R., et al. 2023, A&A, 674, A28, doi: 10.1051/0004-6361/202243919
  * Gagné et al. (2021) Gagné, J., Faherty, J. K., Moranta, L., & Popinchalk, M. 2021, ApJ, 915, L29, doi: 10.3847/2041-8213/ac0e9a
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Gaia Collaboration et al. (2023) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1, doi: 10.1051/0004-6361/202243940
  * Gieles et al. (2007) Gieles, M., Athanassoula, E., & Portegies Zwart, S. F. 2007, MNRAS, 376, 809, doi: 10.1111/j.1365-2966.2007.11477.x
  * Gieles et al. (2006) Gieles, M., Portegies Zwart, S. F., Baumgardt, H., et al. 2006, MNRAS, 371, 793, doi: 10.1111/j.1365-2966.2006.10711.x
  * Gieles & Renaud (2016) Gieles, M., & Renaud, F. 2016, MNRAS, 463, L103, doi: 10.1093/mnrasl/slw163
  * Goodwin & Bastian (2006) Goodwin, S. P., & Bastian, N. 2006, MNRAS, 373, 752, doi: 10.1111/j.1365-2966.2006.11078.x
  * Gouliermis (2018) Gouliermis, D. A. 2018, PASP, 130, 072001, doi: 10.1088/1538-3873/aac1fd
  * Hao et al. (2022) Hao, C. J., Xu, Y., Wu, Z. Y., et al. 2022, A&A, 660, A4, doi: 10.1051/0004-6361/202243091
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * He et al. (2022) He, Z., Wang, K., Luo, Y., et al. 2022, ApJS, 262, 7, doi: 10.3847/1538-4365/ac7c17
  * Helmi (2020) Helmi, A. 2020, ARA&A, 58, 205, doi: 10.1146/annurev-astro-032620-021917
  * Hunt & Reffert (2023) Hunt, E. L., & Reffert, S. 2023, A&A, 673, A114, doi: 10.1051/0004-6361/202346285
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Jilinski et al. (2008) Jilinski, E., Ortega, V. G., de la Reza, R., Drake, N. A., & Bazzanella, B. 2008, arXiv e-prints, arXiv:0810.1198, doi: 10.48550/arXiv.0810.1198
  * Kalirai & Richer (2010) Kalirai, J. S., & Richer, H. B. 2010, Philosophical Transactions of the Royal Society of London Series A, 368, 755, doi: 10.1098/rsta.2009.0257
  * Kamdar et al. (2019) Kamdar, H., Conroy, C., Ting, Y.-S., et al. 2019, ApJ, 884, 173, doi: 10.3847/1538-4357/ab44be
  * Kerr et al. (2021) Kerr, R. M. P., Rizzuto, A. C., Kraus, A. L., & Offner, S. S. R. 2021, ApJ, 917, 23, doi: 10.3847/1538-4357/ac0251
  * Kounkel & Covey (2019) Kounkel, M., & Covey, K. 2019, AJ, 158, 122, doi: 10.3847/1538-3881/ab339a
  * Kounkel et al. (2020) Kounkel, M., Covey, K., & Stassun, K. G. 2020, AJ, 160, 279, doi: 10.3847/1538-3881/abc0e6
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231, doi: 10.1046/j.1365-8711.2001.04022.x
  * Kruijssen (2012) Kruijssen, J. M. D. 2012, MNRAS, 426, 3008, doi: 10.1111/j.1365-2966.2012.21923.x
  * Kuhn et al. (2019) Kuhn, M. A., Hillenbrand, L. A., Sills, A., Feigelson, E. D., & Getman, K. V. 2019, ApJ, 870, 32, doi: 10.3847/1538-4357/aaef8c
  * Lada & Lada (2003) Lada, C. J., & Lada, E. A. 2003, ARA&A, 41, 57, doi: 10.1146/annurev.astro.41.011802.094844
  * Li et al. (2022) Li, Q.-D., Wang, H.-F., Luo, Y.-P., et al. 2022, ApJS, 262, 20, doi: 10.3847/1538-4365/ac81be
  * Luhman & Esplin (2020) Luhman, K. L., & Esplin, T. L. 2020, AJ, 160, 44, doi: 10.3847/1538-3881/ab9599
  * Luo et al. (2015) Luo, A. L., Zhao, Y.-H., Zhao, G., et al. 2015, Research in Astronomy and Astrophysics, 15, 1095, doi: 10.1088/1674-4527/15/8/002
  * McInnes et al. (2017) McInnes, L., Healy, J., & Astels, S. 2017, Journal of Open Source Software, 2, 205, doi: 10.21105/joss.00205
  * Meingast et al. (2021) Meingast, S., Alves, J., & Rottensteiner, A. 2021, A&A, 645, A84, doi: 10.1051/0004-6361/202038610
  * Miyamoto & Nagai (1975) Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533
  * Moraux et al. (2004) Moraux, E., Kroupa, P., & Bouvier, J. 2004, A&A, 426, 75, doi: 10.1051/0004-6361:20035669
  * Navarro et al. (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563, doi: 10.1086/177173
  * pandas development team (2020) pandas development team, T. 2020, pandas-dev/pandas: Pandas, latest, Zenodo, doi: 10.5281/zenodo.3509134
  * Parker et al. (2015) Parker, R. J., Dale, J. E., & Ercolano, B. 2015, MNRAS, 446, 4278, doi: 10.1093/mnras/stu2393
  * Pearson et al. (2017) Pearson, S., Price-Whelan, A. M., & Johnston, K. V. 2017, Nature Astronomy, 1, 633, doi: 10.1038/s41550-017-0220-3
  * Pellerin et al. (2007) Pellerin, A., Meyer, M., Harris, J., & Calzetti, D. 2007, ApJ, 658, L87, doi: 10.1086/515437
  * Pera et al. (2021) Pera, M. S., Perren, G. I., Moitinho, A., Navone, H. D., & Vazquez, R. A. 2021, pyUPMASK: Unsupervised clustering method for stellar clusters, Astrophysics Source Code Library, record ascl:2101.016
  * Plummer (1911) Plummer, H. C. 1911, MNRAS, 71, 460, doi: 10.1093/mnras/71.5.460
  * Price-Whelan (2017) Price-Whelan, A. M. 2017, The Journal of Open Source Software, 2, doi: 10.21105/joss.00388
  * Sills et al. (2018) Sills, A., Rieder, S., Scora, J., McCloskey, J., & Jaffa, S. 2018, MNRAS, 477, 1903, doi: 10.1093/mnras/sty681
  * Skumanich (1972) Skumanich, A. 1972, ApJ, 171, 565, doi: 10.1086/151310
  * Soderblom (2010) Soderblom, D. R. 2010, ARA&A, 48, 581, doi: 10.1146/annurev-astro-081309-130806
  * Szentgyorgyi et al. (2011) Szentgyorgyi, A., Furesz, G., Cheimets, P., et al. 2011, PASP, 123, 1188, doi: 10.1086/662209
  * Tian (2020) Tian, H.-J. 2020, ApJ, 904, 196, doi: 10.3847/1538-4357/abbf4b
  * Torniamenti et al. (2021) Torniamenti, S., Ballone, A., Mapelli, M., et al. 2021, MNRAS, 507, 2253, doi: 10.1093/mnras/stab2238
  * Ujjwal et al. (2020) Ujjwal, K., Kartha, S. S., Mathew, B., Manoj, P., & Narang, M. 2020, VizieR Online Data Catalog, J/AJ/159/166, doi: 10.26093/cds/vizier.51590166
  * van Leeuwen et al. (2022) van Leeuwen, F., de Bruijne, J., Babusiaux, C., et al. 2022, Gaia DR3 documentation, Gaia DR3 documentation, European Space Agency; Gaia Data Processing and Analysis Consortium. Online at ¡A href=“https://gea.esac.esa.int/archive/documentation/GDR3/index.html”¿https://gea.esac.esa.int/archive/documentation/GDR3/index.html¡/A¿, id. 1
  * Vasiliev & Baumgardt (2021) Vasiliev, E., & Baumgardt, H. 2021, MNRAS, 505, 5978, doi: 10.1093/mnras/stab1475
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Ward et al. (2020) Ward, J. L., Kruijssen, J. M. D., & Rix, H.-W. 2020, MNRAS, 495, 663, doi: 10.1093/mnras/staa1056
  * Waskom (2021) Waskom, M. L. 2021, Journal of Open Source Software, 6, 3021, doi: 10.21105/joss.03021
  * Wright et al. (2022) Wright, N. J., Goodwin, S., Jeffries, R. D., Kounkel, M., & Zari, E. 2022, arXiv e-prints, arXiv:2203.10007, doi: 10.48550/arXiv.2203.10007


# Primitivity and Independent Sets in Direct Products of Vertex-Transitive
Graphs ††thanks: Supported by the National Natural Science Foundation of China
(No.10826084) and Zhejiang Innovation Project (Grant No. T200905).

Huajun Zhang
Department of Mathematics, Shanghai Normal University, Shanghai 200234, China
and
Department of Mathematics, Zhejiang Normal University, Jinhua 321004, P.R.
China
E-mail<EMAIL_ADDRESS>
Abstract. We introduce the concept of the primitivity of independent set in
vertex-transitive graphs, and investigate the relationship between the
primitivity and the structure of maximum independent sets in direct products
of vertex-transitive graphs. As a consequence of our main results, we
positively solve an open problem related to the structure of independent sets
in powers of vertex-transitive graphs.

##  1 Introduction

The direct product $G\times H$ of two graphs $G$ and $H$ is defined by

$V(G\times H)=V(G)\times V(H)$

and

$E(G\times H)=\\{[(u_{1},u_{2}),(v_{1},v_{2})]:[u_{1},v_{1}]\in E(G)\mbox{\
and \ }[u_{2},v_{2}]\in E(H)\\}.$

For a graph $G$, let $G^{n}=G\times\cdots\times G$ denote the $n$-th power of
$G$.

It is clear that if $I$ is an independent set of $G$ (or $H$), then $I\times
H$ (or $G\times I$) is an independent set of $G\times H$. We say that $G\times
H$ is MIS-normal (maximum-independent-set-normal) if each of its maximum
independent sets is of this form. Then the independence number

$\displaystyle\alpha(G\times H)=\max\\{\alpha(G)|H|,\alpha(H)|G|\\}$ (1)

if $G\times H$ is MIS-normal. A product $G_{1}\times G_{2}\times\cdots\times
G_{n}$ is said to be MIS-normal if all of its maximum independent sets are
preimages of projections of maximum independent sets of one of its factors.

This poses two immediate problems: whether (1) holds for all graphs $G$ and
$H$, and whether $G\times H$ is MIS-normal when (1) holds. In general,
however, (1) does not hold for some non-vertex-transitive graphs (see [7]).
So, Tardif [3] asked whether (1) holds for all vertex-transitive graphs $G$
and $H$. Larose and Tardif [2] investigated the relationship between the
projectivity and the structure of maximal independent sets in powers of a
circular graph, Kneser graph, or truncated simplex. Recently, Mario and Vera
[5] proved that (1) holds for some special vertex-transitive graphs, e.g.,
circular graphs and Kneser graphs. In fact, Frankl [6] proved in 1996, one
year before Tardif’s question was posed, that (1) holds for Kneser graphs.
Subsequently, Ahlswede, Aydinian and Khachatrian [8] generalized Frankl’s
result.

In the context of vertex-transitive graphs, the “No-Homomorphism” lemma of
Albertson and Collins [1] is useful to get bounds on the size of independent
sets.

###### Lemma 1.1

([1]) Let $G$ and $H$ be two graphs such that $G$ is vertex-transitive and
there exists a homomorphism $\phi:H\mapsto G$. Then
$\frac{\alpha(G)}{|V(G)|}\leq\frac{\alpha(H)}{|V(H)|},$ and the equality holds
if and only if for any independent set $I$ of cardinality $\alpha(G)$ in $G$,
$\phi^{-1}(I)$ is an independent set of cardinality $\alpha(H)$ in $H$.

By this lemma, it is easy to deduce that $\alpha(G^{n})=\alpha(G)|V(G)|^{n-1}$
for any vertex-transitive graph $G$ and positive integer $n$ (see [2]). So it
is natural to ask whether $G^{n}$ is MIS-normal. Evidently, if $G^{n}$ is MIS-
normal for some $n>2$, so is $G^{2}$. Conversely, Larose and Tardif [2] posed
the following problem.

###### Problem 1.2

(see [2] )Let $G$ be a non-bipartite vertex-transitive graph. If $G^{2}$ is
MIS-normal, is the same for all powers of $G$?

This paper is organized as follows. In the next section, we introduce a
concept of the primitivity of independent sets in a vertex-transitive graph,
and prove that the primitivity can be preserved in direct products under
certain conditions. Based on these results we establish in section $3$ a
direct product theorem on the MIS-normality. As a consequence, Problem 1.2 is
solved.

##  2  Primitivity of independent sets

In the sequel of this paper, let $G$ and $H$ be vertex-transitive graphs. By
$I(G)$ we denote the set of all maximum independent sets of $G$. For any
subset $A$ of $V(G)$, let $\alpha(A)$ denote the independence number of the
induced subgraph of $G$ by $A$, and we define

$N_{G}(A)=\\{b\in G:\mbox{$(a,b)\in E(G)$ for some $a\in A$}\\},$
$N_{G}[A]=N_{G}(A)\cup A\ \mbox{and}\ \overline{N}_{G}[A]=G-N_{G}[A].$

In Lemma 1.1, by taking $H$ as an induced subgraph of $G$ and $\phi$ as the
embedding mapping, we obtain the following lemma (cf. [4]).

###### Lemma 2.1

$\frac{\alpha(G)}{|V(G)|}\leq\frac{\alpha(B)}{|B|}$ holds for all $B\subseteq
V(G)$. Equality implies that $|S\cap B|=\alpha(B)$ for every $S\in I(G)$.

A graph $G$ is said to be non-empty if $E(G)\neq\emptyset$. Lemma 2.1 implies
that $\alpha(G)\leq|V(G)|/2$ for all non-empty vertex-transitive graphs.
Equality holds if and only if $G$ is bipartite, which we state as a corollary
for reference.

###### Corollary 2.2

Let $G$ be a non-empty vertex-transitive graph. Then
$\frac{\alpha(G)}{|G|}\leq\frac{1}{2}$, and equality holds if and only if $G$
is bipartite.

###### Proposition 2.3

Let $A$ be an independent set of $G$. Then
$\frac{|A|}{|N_{G}[A]|}\leq\frac{\alpha(G)}{|V(G)|}$. Equality implies that
$|S\cap N_{G}[A]|=|A|$ for every $S\in I(G)$, and in particularly $A\subseteq
S$ for some $S\in I(G)$.

Proof. Since $A$ is an independent set, clearly

$\frac{|A|+\alpha(\overline{N}_{G}[A])}{|N_{G}[A]|+|\overline{N}_{G}[A]|}\leq\frac{\alpha(G)}{|V(G)|}.$

By Lemma 2.1 we see that
$\frac{\alpha(\overline{N}_{G}[A])}{|\overline{N}_{G}[A]|}\geq\frac{\alpha(G)}{|V(G)|}$,
so $\frac{|A|}{|N_{G}[A]|}\leq\frac{\alpha(G)}{|V(G)|}$. Equality in the
latter implies equality in the former. In this case any $S\in I(G)$ must be
the union of a maximum independent set in $\overline{N}_{G}[A]$ and an
independent set of size $|A|$ in $N_{G}[A]$, and thus $|S\cap N_{G}[A]|=|A|$.
$\Box$

An independent set $A$ in $G$ is said to be imprimitive if $|A|<\alpha(G)$ and
$\frac{|A|}{|N_{G}[A]|}=\frac{\alpha(G)}{|V(G)|}$. We say that $G$ is IS-
imprimitive if $G$ has an imprimitive independent set. In the other case, $G$
is _IS-primitive_.

###### Proposition 2.4

Let $A$ be a maximum imprimitive independent set of $G$. Set
$B=\overline{N}_{G}[A]$. Then $\frac{\alpha(B)}{|B|}=\frac{\alpha(G)}{|V(G)|}$
and $\\{\sigma(B)|\sigma\in\mbox{Aut}(G)\\}$ forms a nontrivial partition of
$V(G)$, i.e., $\sigma(B)\cap B=\emptyset$ or $B$ for each
$\sigma\in\mbox{Aut}(G)$.

Proof. Clearly
$\frac{|A|+\alpha(B)}{|N_{G}[A]|+|B|}\leq\frac{\alpha(G)}{|V(G)|}$. Combining
the condition of $A$ and Lemma 2.1, we have
$\frac{\alpha(B)}{|B|}=\frac{\alpha(G)}{|V(G)|}$. By definition,
$N_{G}[\sigma(A)]=\sigma(N_{G}[A])$ for all $\sigma\in\mbox{Aut}(G)$. Suppose
that there exists a $\sigma\in\mbox{Aut}(G)$ such that $\sigma(B)\neq B$ and
$\sigma(B)\cap B\neq\emptyset$. Then $\sigma(N_{G}[A])\neq N_{G}[A]$ and

$\displaystyle|V(G)|>|N_{G}[A]\cup\sigma\big{(}N_{G}[A]\big{)}|>|N_{G}[A]|.$
(2)

Let $C=\sigma(A)\cup(A-N_{G}[\sigma(A)])$. Then $C$ is also an independent set
and

$N_{G}[C]\subseteq N_{G}[A]\cup\sigma(N_{G}[A]).$

By Proposition 2.3, $|S\cap N_{G}[A]|=|A|$ for all $S\in I(G)$, which implies
that $(S-N_{G}[A])\cup A\in I(G)$ for all $S\in I(G)$. Similarly,

$\displaystyle((S-N_{G}[A])\cup A)-N_{G}[\sigma(A)])\cup\sigma(A)$
$\displaystyle=$ $\displaystyle(S-N_{G}[A]\cup
N_{G}[\sigma(A)])\cup(A-N_{G}[\sigma(A)])\cup\sigma(A)$ $\displaystyle=$
$\displaystyle(S-N_{G}[A]\cup N_{G}[\sigma(A)])\cup C$

is also a maximum independent set of $G$, which implies $|S\cap(N_{G}[A]\cup
N_{G}[\sigma(A)])|=|C|$ for all $S\in I(G)$.

Given a $u\in V(G)$, suppose that there are $r$ $S$’s in $I(G)$ such that
$u\in S$. Since $G$ is vertex-transitive, the number $r$ is independent of the
choice of $u$. Thus $r|V(G)|=\alpha(G)|I(G)|$. On the other hand, since
$|S\cap(N_{G}[A]\cup N_{G}[\sigma(A)])|=|C|$ for all $S\in I(G)$,
$|C||I(G)|=r|N_{G}[A]\cup N_{G}[\sigma(A)]|$. Combining the above two
equalities, we have $\frac{|C|}{|N_{G}[A]\cup
N_{G}[\sigma(A)]|}=\frac{\alpha(G)}{|V(G)|}$. Thus, by Proposition 2.3 we have

$\frac{\alpha(G)}{|V(G)|}\geq\frac{|C|}{N_{G}[C]}\geq\frac{|C|}{|N_{G}[A]\cup
N_{G}[\sigma(A)]|}=\frac{\alpha(G)}{|V(G)|},$

which implies $N_{G}[C]=N_{G}[A]\cup N_{G}[\sigma(A)]$ and
$\frac{|C|}{|N_{G}[C]|}=\frac{\alpha(G)}{|V(G)|}$. By (2), we have
$|A|<|C|<\alpha(G)$, contradicting the maximality of $|A|$. This completes the
proof. $\Box$

The concept of primitivity comes from permutation groups: A permutation group
$\Gamma$ acting on a set $X$ is called primitive if $\Gamma$ preserves no
nontrivial partition of $X$. In the other case, $\Gamma$ is imprimitive. As
usual (see e.g. [2]), a vertex-transitive graph $G$ is called primitive if its
automorphism group, as a permutation group on $V(G)$, is primitive. By
Proposition 2.4 we see that if $G$ is primitive, then $G$ is IS-primitive. But
the converse is not true.

For any $S\subseteq V(G)\times V(H)$, $a\in G$ and $u\in H$, define

$\partial_{G}(u,S)=\\{b\in G:(b,u)\in S\\},\ \ \partial_{H}(a,S)=\\{v\in
H:(a,v)\in S\\},$

and

$\partial_{G}(S)=\\{b\in G:\partial_{H}(b,S)\neq\emptyset\\},\ \
\partial_{H}(S)=\\{v\in H:\partial_{G}(v,S)\neq\emptyset\\}.$

By definition we see that $\partial_{G}(S)$ and $\partial_{H}(S)$ are in fact
the projections of $S$ on $G$ and $H$, respectively.

###### Lemma 2.5

Suppose $G\times H$ is MIS-normal and
$\frac{\alpha(H)}{|H|}\leq\frac{\alpha(G)}{|G|}$. If $G\times H$ is IS-
imprimitive, then one of the following two possible cases holds:

  * (i)

$\frac{\alpha(H)}{|H|}=\frac{\alpha(G)}{|G|}$, and one of them is IS-
imprimitive or both $G$ and $H$ are bipartite;

  * (ii)

$\frac{\alpha(H)}{|H|}<\frac{\alpha(G)}{|G|}$, and $G$ is IS-imprimitive or
$H$ is disconnected.

Proof. Throughout this proof, we denote $N_{G\times H}[A]$ by $N[A]$ for
brevity. Suppose that $G\times H$ is IS-imprimitive and let $A$ be a maximum
imprimitive independent set of $G\times H$. Clearly, $\alpha(G\times
H)=\alpha(G)|V(H)|$, and thus $\frac{|A|}{|N[A]|}=\frac{\alpha(G\times
H)}{|V(G\times H)|}=\frac{\alpha(G)}{|V(G)|}$. If $E(G)=\emptyset$, the result
is trivial, so we suppose $E(G)\neq\emptyset$, then Corollary 2.2 implies that
$\frac{\alpha(H)}{|V(H)|}\leq\frac{\alpha(G)}{|V(G)|}\leq\frac{1}{2}$. By
Proposition 2.3, there exists some $S\in I(G\times H)$ such that $A=S\cap
N[A]$. Since $G\times H$ is MIS-normal, we may assume that $S=S^{\prime}\times
H$ for some $S^{\prime}\in I(G)$. Thus $A=(S^{\prime}\times H)\cap N[A]$. Set
$B=\overline{N}[A]$. Then, by Proposition 2.4, $\sigma(B)\cap B=\emptyset$ or
$B$ for every $\sigma\in\mbox{Aut}(G\times H)$.

Set $C=\partial_{G}(B)$. For every pair $a$ and $b$ of $C$, select
$u\in\partial_{H}(a,B)$ and $v\in\partial_{H}(b,B)$. Since $G$ and $H$ are
vertex-transitive, there exist $\gamma\in\mbox{Aut}(G)$ and
$\tau\in\mbox{Aut}(H)$ such that $a=\gamma(b)$ and $u=\tau(v)$. It is clear
that $\sigma=(\gamma,\tau)\in\mbox{Aut}(G\times H)$ and
$(a,u)=\sigma(b,v)\in\sigma(B)\cap B$. By Proposition 2.4, we conclude that
$\sigma(B)=B$. Thus, we have $\partial_{H}(a,B)=\tau(\partial_{H}(b,B))$.
Therefore, $|\partial_{H}(a,B)|=|\partial_{H}(b,B))|$ for any $a,b\in C$. In
the following, we will complete the proof by two cases.

Case $1$: $C\neq V(G)$. Set $\overline{C}=(V(G)-C)$. Then $(\overline{C}\times
H)\cap B=\emptyset$, and thus $\overline{C}\times H\subseteq N[A]$. For every
$S^{\prime\prime}\in I(G)$, it is clear that $S^{\prime\prime}\times H$ is a
maximum independent set of $G\times H$. Since
$\frac{\alpha(B)}{|B|}=\frac{\alpha(G\times H)}{|G\times
H|}=\frac{\alpha(G)}{|V(G)|}$ and $|\partial_{H}(a,B)|=|\partial_{H}(b,B)|$
for all $a,b\in\partial_{G}(B)$, from Lemma 2.1 and the MIS-normality of
$G\times H$ it follows that

$\frac{|(S^{\prime\prime}\times H)\cap B|}{|B|}=\frac{|S^{\prime\prime}\cap
C|}{|C|}=\frac{\alpha(G)}{|V(G)|}.$

Thus for every $S^{\prime\prime}\in I(G)$,

$\displaystyle\frac{\alpha(G)}{|V(G)|}=\frac{|S^{\prime\prime}|}{|V(G)|}=\frac{|S^{\prime\prime}\cap
C|+|S^{\prime\prime}\cap\overline{C}|}{|C|+|\overline{C}|}=\frac{|S^{\prime\prime}\cap\overline{C}|}{|\overline{C}|}=\frac{|S^{\prime\prime}\cap
C|}{|C|}.$ (3)

Recall that $\overline{C}\times H\subseteq N[A]$ and $A\subseteq
S^{\prime}\times H$, it is easy to see that $A=N[A]\cap(S^{\prime}\times H)$
and $\partial_{G}(A\cap(\overline{C}\times H))=S^{\prime}\cap\overline{C}$.
Setting $F=S^{\prime}\cap\overline{C}$, we have that $a\times H\subseteq A$
for every $a\in F$. If $N_{G}[F]\cap C\neq\emptyset$, then there exist $a\in
F$ and $b\in{C}$ such that $(a,b)\in E(G)$. Since $B=\overline{N}[A]$ and
$a\times H\subseteq A$, by definition, $(b,u)\subseteq\overline{N}[a\times H]$
for every $u\in\partial_{H}(b,B)$. Hence $N_{H}[H]\neq\emptyset$ and
$E(H)=\emptyset$, which contradicts that
$\frac{\alpha(H)}{|H|}\leq\frac{1}{2}$. Thus $N_{G}[F]\cap C=\emptyset$, i.e.,
$N_{G}[F]\subseteq\overline{C}$. By Proposition 2.3 and (3),

$\frac{\alpha(G)}{|V(G)|}\geq\frac{|F|}{|N_{G}[F]|}=\frac{|{S^{\prime}}\cap\overline{C}|}{|N_{G}[F]|}\geq\frac{|{S^{\prime}}\cap\overline{C}|}{|\overline{C}|}=\frac{\alpha(G)}{|V(G)|}.$

Therefore $\frac{|F|}{|N_{G}[F]|}=\frac{\alpha(G)}{|V(G)|}$, so $G$ is IS-
imprimitive and (i) holds.

Case $2$: $C=V(G)$. Since $|\partial_{H}(a,B)|=|\partial_{H}(b,B))|$ for all
$a,b\in V(G)$, we have $\partial_{G}(N[A])=V(G)$ and
$|\partial_{H}(a,N[A])|=|\partial_{H}(b,N[A])|<|H|$ for all $a,b\in V(G)$.
Since $A=(S^{\prime}\times H)\cap N[A]$,
$\partial_{H}(a,N[A])\subseteq\partial_{H}(a,S^{\prime}\times H)$ for all
$a\in\partial_{G}(A)$. Thus $\partial_{H}(a,A)=\partial_{H}(a,N[A])$ for all
$a\in\partial_{G}(A)$. Select two vertices $a$ and $b$ of $V(G)$ such that
$a\in\partial_{G}(A)$ and $(a,b)\in E(G)$. Then, for every
$u\in[V(H)-\partial_{H}(b,N[A])]$ and $v\in\partial_{H}(a,N[A])$, it is clear
that $[(b,u),(a,v)]\not\in E(G\times H)$, so $(u,v)\not\in E(H)$. This means
$u\not\in N_{H}(\partial_{H}(a,N[A]))$, that is,

$\displaystyle V(H)-\partial_{H}(b,N[A])\subseteq
V(H)-N_{H}(\partial_{H}(a,N[A])).$ (4)

If $\partial_{H}(b,N[A])=\partial_{H}(a,N[A])$, it follows from (4) that $H$
is disconnected, and so either (i) or (ii) holds.

Suppose that $\partial_{H}(b,N[A])\neq\partial_{H}(a,N[A])$ and set
$D=\partial_{H}(a,N[A])-\partial_{H}(b,N[A])$. It is easy to check that

$\displaystyle
2|D|=|\partial_{H}(a,N[A])\cup\partial_{H}(b,N[A])-\partial_{H}(a,N[A])\cap\partial_{H}(b,N[A])|.$

Since $D\subseteq H-\partial_{H}(b,N[A])$ and
$D\subseteq\partial_{G}(a,N[A])$, by (4), we have

$D\subseteq V(H)-\partial_{H}(b,N[A])\subseteq
V(H)-{N}_{H}(\partial_{H}(a,N[A]))\subseteq V(H)-{N}_{H}(D).$

So $D$ is an independent set of $H$ and

$\displaystyle N_{H}[D]$ $\displaystyle\subseteq$ $\displaystyle
D\cup[\partial_{H}(b,N[A])-\partial_{H}(a,N[A])]$ $\displaystyle=$
$\displaystyle\partial_{H}(a,N[A])\cup\partial_{H}(b,N[A])-\partial_{H}(a,N[A])\cap\partial_{H}(b,N[A]),$

which implies that
$\frac{1}{2}\geq\frac{\alpha(H)}{|V(H)|}\geq\frac{|D|}{|N_{H}[D]|}\geq\frac{1}{2}$.
Thus $\frac{\alpha(G)}{|V(G)|}=\frac{\alpha(H)}{|V(H)|}=\frac{1}{2}$. By
Corollary 2.2, $G$ and $H$ are both bipartite, so (i) holds and the proof
completed. $\Box$

###### Theorem 2.6

Let $G$ and $H$ be two non-bipartite vertex-transitive graph such that
$\frac{\alpha(H)}{|V(H)|}=\frac{\alpha(G)}{|V(G)|}$. If $G\times H$ is MIS-
normal, then $G$, $H$ and $G\times H$ are all IS-primitive.

Proof. First, suppose that $G$ is IS-imprimitive and let $A$ be an imprimitive
independent set in $G$. For any $S\in I(H)$, let
$S^{\prime}=(\overline{N}_{G}[A]\times S)\cup(A\times H)$. It is clear that
$S^{\prime}$ is an independent set of $G\times H$ and

$\displaystyle|S^{\prime}|$ $\displaystyle=$
$\displaystyle|\overline{N}_{G}[A]\alpha(H)|+|A||V(H)|=(|\overline{N}_{G}[A]|+|N_{G}[A]|)\alpha(H)$
$\displaystyle=$ $\displaystyle|V(G)|\alpha(H)=\alpha(G\times H),$

i.e., $S^{\prime}$ is a maximum independent set of $G\times H$, contradicting
the MIS-normality of $G$. Therefore, $G$ is IS-primitive. Similarly, $H$ is
also IS-primitive. By Lemma 2.5, $G\times H$ is IS-primitive. $\Box$

##  3 MIS-normality of the Products of Graphs

The following theorem is the main result on the MIS-normality of products of
vertex-transitive graphs in this paper.

###### Theorem 3.1

Let $G$ and $H$ be two vertex-transitive graphs. Suppose that there exists an
induced subgraph $G^{\prime}$ of $G$ such that $G^{\prime}\times H$ is MIS-
normal and
$\frac{\alpha(G^{\prime})}{|V(G^{\prime})|}=\frac{\alpha(G)}{|V(G)|}$. Then
either: (i) $G\times H$ is MIS-normal, or (ii)
$\frac{\alpha(G)}{|V(G)|}=\frac{\alpha(H)}{|V(H)|}$ and $G$ is IS-imprimitive,
or (iii) $\frac{\alpha(G)}{|V(G)|}<\frac{\alpha(H)}{|V(H)|}$ and $G$ is
disconnected.

Proof. If $E(H)=\emptyset$, the result is obvious, so we assume that
$E(H)\neq\emptyset$. By Lemma 2.1 and the MIS-normality of $G^{\prime}\times
H$, we have the following inequality

$\frac{\alpha(G\times H)}{|V(G)||V(H)|}\leq\frac{\alpha(G^{\prime}\times
H)}{|V(G^{\prime})||V(H)|}=\max\left\\{\frac{\alpha(G)}{|V(G)|},\frac{\alpha(H)}{|V(H)|}\right\\}\leq\frac{\alpha(G\times
H)}{|V(G)||V(H)|},$

yielding

$\displaystyle\frac{\alpha(G\times
H)}{|V(G)||V(H)|}=\frac{\alpha(G^{\prime}\times
H)}{|V(G^{\prime})||V(H)|}=\max\left\\{\frac{\alpha(G)}{|V(G)|},\frac{\alpha(H)}{|V(H)|}\right\\}.$
(5)

For every $\sigma\in\mbox{Aut}(G)$, it is clear that $\sigma(G^{\prime})\times
H$ is MIS-normal. Let $S$ be a maximum independent set of $G\times H$. By
Lemma 2.1 and (5), $S\cap(\sigma(G^{\prime})\times H)$ is a maximum
independent set of $\sigma(G^{\prime})\times H$. Clearly, for each
$a\in\partial_{G}(S)$, there is a $\sigma\in\mbox{Aut}(G)$ such that
$a\in\sigma(G^{\prime})$. We therefore have that $|\partial_{H}(a,S)|=|H|$ or
$\alpha(H)$ for each $a\in\partial_{G}(S)$. In the following we distinguish
three cases to complete the proof.

Case $1$: $|\partial_{H}(a,S)|=|V(H)|$ for every $a\in\partial_{G}(S)$. By
(5), we obtain that $|\partial_{G}(S)|=\alpha(G)$. Since $E(H)\neq\emptyset$,
$\partial_{G}(S)$ is an independent set of $G$. This implies that
$S=\partial_{G}(S)\times H$.

Case $2$: $|\partial_{H}(a,S)|=\alpha(H)$ for every $a\in\partial_{G}(S)$. By
(5), we have that $\partial_{G}(S)=G$,
$\frac{\alpha(H)}{|V(H)|}\geq\frac{\alpha(G)}{|V(G)|}$ and $\partial_{H}(a,S)$
is a maximum independent set of $H$ for every $a\in G$. Let $a$ be a fixed
vertex of $G$, and set

$C=\\{c\in G:\partial_{H}(c,S)=\partial_{H}(a,S)\\}.$

If $C=G$, then $S=G\times\partial_{H}(a,S)$. If $C\neq G$, then choose $d\in
G-C$ and $c\in C$. Since $\partial_{H}(c,S)\neq\partial_{H}(d,S)$, there
exists $u\in\partial_{H}(c,S)$ and $v\in\partial_{H}(d,S)$ such that $(u,v)\in
E(H)$ and $[(c,u),(d,v)]\not\in E(G\times H)$. This implies that $(c,d)\in
E(G)$ and thus G is disconnected.

Case $3$: $|\partial_{H}(a,S)|=|V(H)|$ and $|\partial_{H}(b,S)|=\alpha(H)$ for
some $a,b\in\partial_{G}(S)$. By (5),
$\frac{\alpha(H)}{|V(H)|}=\frac{\alpha(G)}{|V(G)|}$ and $\alpha(G\times
H)=\alpha(G)|V(H)|=\alpha(H)|V(G)|$. Set

$C=\\{c\in G:|\partial_{H}(c,S)|=|V(H)|\\}\mbox{\ and \ }D=\\{d\in
G:|\partial_{H}(d,S)|=\alpha(H)\\}.$

Since $E(H)\neq\emptyset$, it is clear that $C$ is an independent set of $G$
and $(c,d)\not\in E(G)$ for every $c\in C$ and $d\in D$. So $N_{G}[C]\subseteq
V(G)-D$. Moreover,

$|S|=\alpha(H)|V(G)|=|C||V(H)|+|D|\alpha(H).$

Thus
$\frac{|C|}{|N_{G}[C]|}\geq\frac{|C|}{|V(G)|-|D|}=\frac{\alpha(H)}{|V(H)|}=\frac{\alpha(G)}{|V(G)|}$.
By Proposition 2.3, $\frac{|C|}{|N_{G}[C]|}=\frac{\alpha(G)}{|V(G)|}$, that
is, $G$ is IS-imprimitive.

This completes the proof. $\Box$

The following Corollary solves Problem 1.2 in a bit more general setting.

###### Corollary 3.2

Let $G$ be a vertex-transitive, non-bipartite graph. If $G^{2}$ is MIS-normal,
then $G^{n}$ is also MIS-normal and IS-primitive for all $n\geq 3$.

Proof. We prove by induction on $n$. Since $G^{2}$ is MIS-normal, by Theorem
2.6, $G$ and $G^{2}$ are both IS-primitive. Assume that $G^{d}$ is MIS-normal
and IS-primitive for all $d=2,\ldots,n-1$. We now prove that $G^{n}$ is MIS-
normal and IS-primitive. Note that $G^{n}=G^{2}\times G^{n-2}$. Let
$G^{\prime}$ be some subgraph of $G^{2}$ that is isomorphic to $G$, for
instance, the subgraph induced by the set of vertices $\\{(u,u):u\in V(G)\\}$.
It is clear that

$\frac{\alpha(G^{\prime})}{|V(G^{\prime})|}=\frac{\alpha(G)}{|V(G)|}=\frac{\alpha(G^{2})}{|V(G^{2})|}$

and $G^{\prime}\times G^{n-2}$ is isomorphic to $G^{n-1}$. Thus by assumption,
$G^{\prime}\times G^{n-2}$ is MIS-normal. By Theorem 3.1 and Theorem 2.6, it
is easy to see that $G^{n}$ is MIS-normal and IS-primitive. This completes the
proof. $\Box$

Acknowledgement The author is greatly indebted to the anonymous referees for
giving useful comments and suggestions that have considerably improved the
manuscript. He is grateful also for many valuable discussions with Professor
J. Wang and Professor C.J. Zhou.

## References

  * [1] M.O. Albertson and K.L. Collins, Homomorphisms of $3$-chromatic graphs, Discrete Math., 54 (1985) 127-132.
  * [2] B. Larose and C. Tardif, Projectivity and independent sets in powers of graphs, J. Graph Theory, 40 (2002) 162-171.
  * [3] C. Tardif, Graph products and the chromatic difference sequence of vertex-transitive graphs, Discrete Math., 185 (1998) 193-200.
  * [4] P.J. Cameron and C.Y. Ku, Intersecting families of permutations, European J. Comb., 24 (2003) 881-890.
  * [5] V.P. Mario and J. Vera, Independent and coloring properties of direct products of some vertex-transitive graphs, Discrete Math., 306 (2006) 2275-2281.
  * [6] P. Frankl, An Erdős-Ko-Rado Theorem for direct products, European J. Combin., 17 (1996) 727-730.
  * [7] P.K. Jha and S. Klavz̆ar, Independence in direct-product graphs, Ars Combin., 50 (1998) 53-60.
  * [8] R. Ahlswede, H. Aydinian and L.H. Khachatrian, The intersection theorem for direct products, European J. Combin., 19 (1998) 649-661.