11institutetext:  European Space Agency (ESA), European Space Research and
Technology Centre (ESTEC), Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands,
Emails<EMAIL_ADDRESS>22institutetext: European Southern
Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching bei München
Emails<EMAIL_ADDRESS><EMAIL_ADDRESS>
# The 800pc long tidal tails of the Hyades star cluster

Possible discovery of candidate epicyclic overdensities from an open star
cluster

Tereza Jerabkova 11 Henri M.J. Boffin 22 Giacomo Beccari 22 Guido de Marchi 11
Jos H. J. de Bruijne 11 Timo Prusti 11

(Received ; accepted )

The tidal tails of stellar clusters provide an important tool for studying the
birth conditions of the clusters and their evolution, coupling, and
interaction with the Galactic potential. The Gaia satellite, with its high-
quality astrometric data, opened this field of study, allowing us to observe
large-scale tidal tails. Theoretical models of tidal-tail formation and
evolution are available. However, the exact appearance of tidal features as
seen in the Gaia catalogue has not yet been studied. Here we present the
$N-$body evolution of a Hyades-like stellar cluster with backward-integrated
initial conditions on a realistic 3D orbit in the Milky Way (MW) galaxy
computed within the AMUSE framework. For the first time, we explore the effect
of the initial cluster rotation and the presence of lumps in the Galactic
potential on the formation and evolution of tidal tails. For all of our
simulations we present Gaia observables and derived parameters in the
convergent point (CP) diagram. We show that the tidal tails are not naturally
clustered in any coordinate system and that they can span up to 40 km/s
relative to the cluster centre in proper motions for a cluster age of 600-700
Myr. Models with initial rotation result in significant differences in the
cluster mass loss and follow different angular momentum time evolution. Thus
the orientation of the tidal tails relative to the motion vector of the
cluster and the current cluster angular momentum constrain the initial
rotation of the cluster. We highlight the use of the standard CP method in
searches for co-moving groups and introduce a new compact CP (CCP) method that
accounts for internal kinematics based on an assumed model. Using the CCP
method, we are able to recover candidate members of the Hyades tidal tails in
the Gaia Data Release 2 (DR2) and early Data Release 3 (eDR3) reaching a total
extent of almost 1kpc. We confirm the previously noted asymmetry in the
detected tidal tails. In the eDR3 data we recovered spatial overdensities in
the leading and trailing tails that are kinematically consistent with being
epicyclic overdensities and thus would present candidates for the first such
detection in an open star cluster. We show that the epicyclic overdensities
are able to provide constraints not only on the cluster properties, but also
on the Galactic potential. Finally, based on $N-$body simulations, a close
encounter with a massive Galactic lump can explain the observed asymmetry in
the tidal tails of the Hyades.

###### Key Words.:

parallaxes – proper motions – stars: kinematics and dynamics – open clusters
and associations: Hyades cluster – Galaxy: stellar content – Galaxy: structure

##  1 Introduction

Stars and stellar clusters form in sub-parsec (pc) dense filamentary regions
of molecular clouds (Kroupa, 1995a, c; Lada & Lada, 2003; Marks & Kroupa,
2012; André et al., 2014). The detailed small-scale physics of star formation
and its initial conditions coupled to the galaxy-wide scales through
gravitational forces and stellar populations still remains to be fully
understood. Tidal tails of star cluster and large-scale relics of star
formation in filaments (Jerabkova et al., 2019b; Beccari et al., 2020) are
unique structures providing a physical link between the initial conditions of
cluster formation and the Galactic potential, its tides, and shears. These
structures are in transition, becoming part of the galactic field (Baumgardt &
Makino, 2003; Kroupa, 2005; Dinnbier & Kroupa, 2020a, b; Meingast et al.,
2020; Dinnbier & Walch, 2020) and creating an additional link to galactic
stellar population signatures (Marks & Kroupa, 2011; Jeřábková et al., 2018).

Embedded star clusters that survive the process of the expulsion of their
residual gas re-virialise to open clusters (OCs) with half-mass radii of about
$3\,$pc (Kroupa et al., 2001; Banerjee & Kroupa, 2017) and as OCs evolve
through mass loss from their dying stars and through loss of members driven by
the energy equipartition or two-body relaxation processes. This process causes
a small increase in binding energy per dynamical time for a small fraction of
stars, such that the stars leave the cluster and orbit the Galaxy (Baumgardt &
Makino, 2003). If the cluster is on a circular orbit, stars that leave with a
lower velocity will fall towards the Galactic centre and will thus speed ahead
as the angular velocity increases on their eccentric orbit with apogalacticon
at the cluster orbit. A star that leaves with a higher velocity will fall
behind it as its orbital angular velocity decreases on its eccentric orbit
with perigalacticon at the cluster orbit. The stars drifting from the cluster
form a symmetrical S-shaped stellar distribution in the immediate surroundings
of the cluster.

Because the large majority of stars leave an OC through the energy
equipartition process (rather than being ejected through three-body or binary-
binary encounters, Oh et al., 2015; Oh & Kroupa, 2016), their velocities
relative to the cluster are within a few km/s. The lost stars slowly drift
from the cluster following very similar orbits. The cluster thus develops a
thin leading tidal tail closer to the Galactic centre and a thin trailing
tidal tail farther away. The tails are expected to compose a symmetric
structure around the cluster (Thomas et al., 2018). While the stars initially
drift away from the cluster potential, they later partially and repeatedly
return as tidal-tail stars, and their cluster of origin moves on orbits with
slightly different eccentricities and consequently with different orbital
frequencies. This drifting-away and returning motion leads to overdensities
along the tidal tails that are known as epicyclic overdensities (Küpper et
al., 2008; Just et al., 2009; Küpper et al., 2012). The tidal tails of
clusters on eccentric orbits have more complicated structures (Küpper et al.,
2010), but they still display regularly spaced epicyclic overdensities, the
spacing of which depends on the Galactic potential. The distance between the
cluster and the Küpper overdensities decreases as the cluster looses mass
(Küpper et al., 2010), allowing a very precise measurement of the Galactic
potential (Küpper et al., 2015).

Much progress has been made on the detection and study of tidal tails around
globular clusters (GC) in the halo of the MW (e.g. Carballo-Bello et al. 2020;
Piatti & Carballo-Bello 2020 or the identification of the long stellar stream
of the massive globular cluster $\omega$ Centauri by Ibata et al. 2019). The
reason is mostly that because the typical orbit of a GC is very high above the
Galactic plane, contamination by field stars is removed. Moreover, GCs are
massive ($>10^{4}\,M_{\odot}$) and have long lifetimes ($\approx 13$ Gyr),
allowing the tidal tails to be well populated by the cluster stars and to
display large extensions on the sky. A landmark discovery, for example, is
given by the long tidal tails of the globular cluster Palomar 5 (Odenkirchen
et al., 2003).

Detecting the tidal tails around OCs is more challenging because of their
young ages (from a few hundred million years up to a few billion years at
most) and their low mass ($\approx 10^{2}-10^{3}$ M⊙). Moreover, OCs are
mostly confined to the Galactic disc, which makes identifying their tidal
tails challenging because they contain few stars and are heavily contaminated
by field stars. However, in the same way as for GC, it is expected that tidal
tails around OCs display an S-shape and epicyclic overdensities, the existence
of which is related to the orbit of the cluster, the initial conditions at
formation, and the Galactic potential. Complications arise because OCs may
have orbits that take them somewhat out of the Galactic mid-plane, and the bar
and spiral patterns of the Galaxy may lead to non-axisymmetric perturbations,
which may affect the properties of the tails.

The advent of the ESA Gaia satellite for the first time allows the study of
young stellar populations on large scales (¿100pc) in six-dimensional (6D)
position and velocity space. For example, using Gaia DR2 (Gaia Collaboration
et al., 2018b), Kounkel et al. (2018) and Zari et al. (2019) demonstrated that
the Orion star-forming region is composed of a variety of populations with
different spatial and kinematic properties each that are all likely generated
in multiple events instead of in a progressive star formation history (see
also Beccari et al., 2017; Jerabkova et al., 2019a; Kroupa et al., 2018,
discussing three bursts of star formation in the Orion Nebula cluster). Zari
et al. (2018) described the 6D properties and age structures in which young
stars are found within a region of 500 pc around the Sun, and a number of
other studies focusing on individual star-forming regions have been reported.
One example is the large-scale picture of the Gamma Velorum region (Beccari et
al., 2018; Cantat-Gaudin et al., 2019). The Gaia DR2 catalogue also led to the
discovery of stellar relic filaments. These are spatial structures of a few pc
that are more than 90 pc long and consist of stars of equal ages that are
younger than a few hundred million years (Jerabkova et al., 2019b; Beccari et
al., 2020). A Galaxy-wide survey using the Gaia DR2 catalogue reveals many
elongated structures made of coeval stars, many of which are tidal tails of
dissolving star clusters (Kounkel & Covey, 2019a). Meingast et al. (2020)
described the nearby extended spatial distribution of stars that co-move with
their clusters but are not bound to them. After the release of the Gaia DR2
catalogue, tidal tails have been found around four nearby ($<300\,$pc) OCs,
namely Blanco 1 ($\approx 100\,$ Myr, Zhang et al. 2020), the Hyades ($\approx
600-700\,$Myr; Meingast & Alves, 2019; Röser et al., 2019; Douglas et al.,
2019; Gossage et al., 2018; Reino et al., 2018; Lodieu, 2020; Gaia
Collaboration et al., 2018a), Coma Berenices ($\approx 750\,$Myr; Tang et al.,
2019; Fürnkranz et al., 2019), and Praesepe ($\approx 800\,$Myr; Röser &
Schilbach, 2019).

The Hyades cluster has attracted much attention because it is only about
$45\,$pc distant from the Sun (e.g. see Chumak et al. 2005 for a
visualisation), such that its stellar content and its tidal tails can most
likely be mapped out with high precision. Using 6D phase-space constraints
(Roeser et al., 2010), Röser et al. (2011) found the half-mass radius of the
Hyades to be $4.1\,$pc and its tidal radius to be $9\,$pc with a bound mass of
$275\,M_{\odot}$. The total mass of Hyades stars found so far is about
$435\,M_{\odot}$. Interestingly, the measured bulk (3D) velocity dispersion
within $9\,$pc is $\approx 0.8\pm 0.1\,$km/s, but the theoretical value for
Hyades stellar mass is expected to be $0.36\,$km/s (table 3 in Röser et al.
2011). The Hyades therefore show a possible mass discrepancy of a factor of
four; part of this may be due to unresolved multiple stars (Dabringhausen et
al., 2016) and/or to the fact that the cluster is not in dynamical equilibrium
(Oh & Evans, 2020). The cluster contains stars with masses between about
$0.1\,M_{\odot}$ and $2.6\,M_{\odot}$, distributed according to a present-day
mass function (PDMF) that is broadly consistent with evaporation from the
cluster, a canonical initial mass function (IMF), and a majority of stars that
have been born in binaries (Kroupa, 1995b, 2011; Ernst et al., 2011).

The astrophysical evolution of this cluster has been studied in depth using
$N-$body methods linked with stellar evolution codes (Portegies Zwart et al.,
2001; Madsen, 2003), confirming the expected theoretical velocity dispersion.
Chumak et al. (2005) in addition provided the first detailed $N-$body study of
the expected tidal-tail structure of the Hyades. They showed that each of the
two symmetrical tails should be about $500\,$pc long and also documented the
degeneracy between the initial mass and the initial radius of the cluster for
the models to fit its present-day bulk properties. The symmetrical tidal tails
contain very pronounced Küpper overdensities about $200\,$pc from the cluster.
Ernst et al. (2011) based their models on the Röser et al. (2011) Hyades
survey constraints, thereby verifying and significantly improving the details
of the cluster evolution. The resultant best-fitting initial mass of a Plummer
model is $1230\,M_{\odot}$ with an initial half-mass radius of $2.62\,$pc,
which leads to an average mass-loss rate of $1.4\,M_{\odot}$/Myr. The present-
day tidal tails reach out to a length of $\approx 800\,$pc according to this.
These derived quantities depend on the choice of the Galactic potential,
however, whereby the current position and velocity vectors of the cluster are
well constrained.

In this paper we perform $N-$body simulations of the Hyades cluster in a
realistic Galactic potential with the aim to investigate the expected extent
of the tidal tails. With this in mind, we develop a new recipe to find these
stellar structures using the Gaia DR2 and eDR3 data with the final goal of
reporting evidence for or against the expected $\approx 800\,$ pc long tidal
tails. This work thus concentrates on the problem of developing a new method
(called the compact convergent point (CCP) method) for finding the leading and
trailing tails in their full expected length.

This work is structured as follows: We introduce Hyades-like cluster evolution
on a realistic orbit and reveal the structure of its tidal tails in various
Gaia spaces in Sec. 2. In the same section we also introduce the new CCP
method in order to project the tidal tails that are otherwise spread in
spatial and velocity coordinates into a compact configuration. In Sec. 3 we
use the information gained in Sec. 2 and identify candidate Hyades tidal-tail
members in the Gaia data spanning a length of 800 pc. In Sec. 4 we study the
effect of an initial angular moment and of a lumpy Galactic potential on the
formation and evolution of tidal tails for the first time. The last two
sections, Sec. 5 and Sec. 6, present a discussion and some conclusions.

##  2 Simulated Gaia view of tidal tails

The main aim of this section is to present for the first time the results of a
simulation in the observable (and derived from observable) parameters of the
Gaia catalogue of the tidal tails of this star cluster. We then determine
whether these simulations can be used to detect the full or at least a larger
extent of the tidal tails of the Hyades star cluster.

We used the code called astronomical multipurpose software environment, AMUSE
(Portegies Zwart & McMillan, 2018; Portegies Zwart et al., 2013; Pelupessy et
al., 2013; Portegies Zwart et al., 2009). It allowed us to couple the
evolution of a star cluster and the description of a Galactic potential to
simulate the evolution of a Hyades-like star cluster on a realistic orbit.

###  2.1 Present-day parameters of the Hyades and coordinate setup

We summarise in Tab. 1 some of the present-day properties of the Hyades as
derived in Röser et al. (2011), Goldman et al. (2013), and Ernst et al.
(2011). In particular, Röser et al. (2011) and Goldman et al. (2013) provided
a catalogue of Hyades members down to 0.116 $M_{\odot}$ and 0.1 $M_{\odot}$,
respectively. The reported parameters are still discrepant because the method
that are used in the different studies differ. As an example, the age of the
Hyades is estimated to be $625\pm 50$ Myr in Perryman et al. (1998), $648\pm
50$ Myr in De Gennaro et al. (2009), 500-650 Myr in Lebreton et al. (2001),
675-700 Myr (Douglas et al., 2019), 680 Myr (Gossage et al., 2018), 675-700
Myr (Reino et al., 2018), 650 $\pm$ 70 Myr (Lodieu, 2020), and 795 Myr (Gaia
Collaboration et al., 2018a). For the purpose of this study, we adopted a
nominal age value of 655 Myr and considered an age spread of 620 Myr to 695
Myr (see Tab. 1) for the models we present. The velocity dispersion of the
Hyades star cluster has been estimated in a number of studies. Madsen (2003)
and Makarov et al. (2000) reported a velocity dispersion of about $0.3$ km/s
based on the cluster virial mass. Douglas et al. (2019) concluded that the
upper limit on the velocity dispersion is about $0.5$ km/s. Röser et al.
(2011) and Oh et al. (2015) derived an even higher value for the velocity
dispersion of about 0.8 km/s. When we discuss our result, we consider all
these values.

Table 1: Present-day properties of the Hyades star cluster based on Röser et al. (2011), Ernst et al. (2011), Goldman et al. (2013), Gaia Collaboration et al. (2018a), Douglas et al. (2019), Gossage et al. (2018), Reino et al. (2018), and Lodieu (2020). ** See the text for a discussion of the velocity dispersion.  parameter | value
---|---
age |  $\approx$ 600-700 Myr
number of stars ($r<30$ pc) | 724
stellar mass ($r<30$ pc) | 469 $M_{\odot}$
3D velocity dispersion ($r<9$ pc) |  $\approx 0.3-0.8$ km/s **

We used the Astropy SkyCoord package to compute all coordinate transformations
in this work (Astropy Collaboration et al., 2018). Galactic and equatorial
coordinates are centred on the barycenter of the Solar System. For the
galactocentric coordinates we used the ’pre-v4.0’ Astropy frame default setup,
which adopts the coordinates of the Sun in the galactocentric Cartesian system
as

$\displaystyle[X_{\odot},Y_{\odot},Z_{\odot}]=[8300.0\,\mathrm{pc},0.0\,\mathrm{pc},27.0\,\mathrm{pc}],$

and the velocities as

$\displaystyle[V_{X,\odot},V_{Y,\odot},V_{Z,\odot}]=[11.1,\,232.24,\,7.25]\,\mathrm{km/s}.$

This framework was chosen because it is one of the default settings in the
widely used tool AstroPy (Astropy Collaboration et al., 2013, 2018). The
circular velocity at the distance of the Sun agrees well with observations
(e.g. Reid et al., 2014). The adopted distance of the Sun from the Galactic
mid-plane is 10 pc larger than the distance found in recent studies (Yao et
al., 2017; Karim & Mamajek, 2017). This difference does not significantly
affect the interpretation of the computed star cluster, in addition, Bland-
Hawthorn & Gerhard (2016) estimated the distance of the Sun from the Galactic
plane to be $25\pm 5$pc. Thus the default setting from AstroPy was used here.

When we refer to the Hyades star cluster, we use the cluster parameters
published by the Gaia Collaboration et al. (2018a, their Table A.3). That is,
the velocity of the cluster in equatorial barycentric spherical coordinates is

$\displaystyle[\mu^{*}_{\alpha},\mu_{\delta},V_{r}]=[101.005\,\mathrm{mas/yr},-28.490\,\mathrm{mas/yr},39.96\,\mathrm{km/s}].$

In their Table A.3, Gaia Collaboration et al. (2018a) published the values of
the equatorial barycentric Cartesian velocities ($U^{\prime}$, $V^{\prime}$,
$W^{\prime}$) and parallax, $\varpi$, which allowed us to derive the cluster
R.A. and Dec coordinates by transforming between the spherical and Cartesian
coordinates,

$\displaystyle[\mathrm{R.A.},\mathrm{Dec},\varpi]=[67.985\,\mathrm{deg},17.012\,\mathrm{deg},21.052\,\mathrm{mas}].$

###  2.2 $N$-body setup

The formation and evolution of tidal tails mostly depend on the interplay
between the dynamical and stellar evolution of the host star cluster and the
interaction with the Galactic potential. The complexity of the study is
therefore mostly driven by the complexity of the physical processes at play,
the uncertainties of the initial conditions of a star cluster, and the
parameters describing the Galaxy. The computations we present here were all
performed with AMUSE (Portegies Zwart & McMillan, 2018; Portegies Zwart et
al., 2013; Pelupessy et al., 2013; Portegies Zwart et al., 2009), which
allowed us to couple the star cluster evolution and the description of a
Galactic potential in a convenient way.

####  2.2.1 Star cluster setup and evolution

To integrate the internal dynamical evolution of star clusters, we used the
Huayno code (Pelupessy et al., 2012; Jänes et al., 2014), which is a class of
$N$-body integrators that implements a variety of kick-drift-kick algorithms
through the Hamiltonian splitting strategy of adjustable order (see Portegies
Zwart & McMillan, 2018; Pelupessy et al., 2012; Jänes et al., 2014, for
detailed descriptions). To integrate the evolution of a Hyades-type star
cluster, the Huayno code was implemented within the AMUSE environment, which
allowed us to couple the direct N-body integration with the stellar evolution
code SeBa (Portegies Zwart & Verbunt, 1996; Toonen et al., 2012) and the
integration along the orbit (see Sec. 2.2.2 for the details). We used the
optimal softening length,

$\varepsilon=2GM_{\mathrm{cl}}^{2}/(N|U|),$ (1)

where $M_{\mathrm{cl}}$ is the stellar mass of the cluster, $G$ is the
gravitational constant, $N$ is the number of stars, and $U$ is the potential
energy of the cluster. For our simulations the value of the softening length
was $\varepsilon=0.004\,\mathrm{pc}$ for the initial cluster parameters. The
mean stellar separation was $\approx 0.2$ pc, which is significantly larger
than the softening length and ensured a sufficiently correct evaluation of
two-body relaxation. The physical origin of the optimal softening length
described by Eq. (1) is $\varepsilon=4R_{\mathrm{virial}}/N$. For the other
parameters, we used the default Huayno values such as the type 8 integrator
HOLD_DKD (for details, see Pelupessy et al., 2012; Jänes et al., 2014;
Portegies Zwart & McMillan, 2018). The parameters we adopted to define the
initial setup for the cluster are given in Tab. 2. In particular, the stellar
initial mass function (IMF) $\xi_{\star}(m)$ was implemented as a two-part
power-law function based on Kroupa et al. (2001),

$\xi_{\star}(m)=\left\\{\begin{array}[]{ll}\,\,\,\,\,\,\,\,\,m^{-1.3}\hskip
46.94687pt0.08\leq m/M_{\odot}<0.50\,,\\\ 0.5\,m^{-2.3}\hskip
46.94687pt0.50\leq m/M_{\odot}<120\,,\\\ \end{array}\right.$ (2)

where $\mathrm{d}N_{\star}=\xi_{\star}(m)\,\mathrm{d}m$ is the number of stars
in the mass interval $m$ to $m+\mathrm{d}m$. We did not include primordial
binaries in the initial cluster setup.

Essentially, the initial model can be seen as constituting a post-gas-
expulsion re-virialised model (Kroupa et al., 2001). We note that it is
interesting and important to study the effect of the exact initial conditions
in detail (e.g. relation of initial mass to radius that matches the observed
present-day properties of the Hyades, the gas expulsion phase, the stellar
IMF, and the initial binary fraction). In this respect, Dinnbier & Kroupa
(2020a, b) have recently studied the effects of gas expulsion on the early
evolution stage ($100\,$Myr) of a tidal tail. With our setup we are limited to
study the overall morphology of tidal tails; see Sec. 5.2, where the
limitations of our method are discussed at length.

To ensure that our setup provides a reliable output comparable with the Hyades
and the morphology of its tidal tails, we proceeded as follows: 1) For a
number of simulations we ran the same setup with the PhiGRAPE code (Harfst et
al., 2007; Gaburov et al., 2009) that is available in AMUSE and was also used
by Kharchenko et al. (2009) to study open clusters and their tidal tails. We
confirmed that these simulations give comparable results. 2) We used the same
initial setup as Ernst et al. (2011), who used the Aarseth Nbody6 code to
simulate the tidal tails of the Hyades, to allow a comparison.

Table 2: Initial parameters of the Hyades star cluster based on Röser et al. (2011) and Ernst et al. (2011). M1, M2, and M3 are the names of our models parameter | value
---|---
Cluster profile | Plummer mode
initial mass $M_{cl,0}$ | $1230\,M_{\odot}$
half-mass radius, $r_{h}$ | 2.62 pc
IMF | canonical IMF, Eq. (2)
angular momentum, $L_{z}$ |  M2 $L_{z}$¡0, M1 $L_{z}$=0, M3 $L_{z}$¿0
initial cluster coordinates | (galactocentric Cartesian f.)
$X,Y,Z$ | 5.605, 7.382, -0.045 kpc
$v_{X},v_{Y},v_{Z}$ | 163.8, -114.0, 3.4 km/s

####  2.2.2 Galaxy parameters and computation setup

To follow the orbital evolution of stellar clusters in the Galaxy, we used a
detailed model of the Galactic potential from the Galaxia module (Martínez-
Barbosa et al., 2016, 2017) available in AMUSE. The Galaxia module consists of
a semi-analytical axisymmetric potential given by a Galactic bulge, disc, and
halo (Allen & Santillan, 1991) and allows additional non-axisymetric
components, such as a bar and spiral arms. We did not implement any non-
axisymetric components. The adopted physical values of the MW are summarised
in Tab. 3. The parameters that are not listed are the same as the Galaxia
module default values. For the MW potential we considered two sets of
parameters from Allen & Santillan (1991) and Irrgang et al. (2013). Allen &
Santillan (1991) has been used in recent studies (e.g. Dinnbier & Kroupa,
2020b) and therefore serves as a good comparison benchmark. Irrgang et al.
(2013) presented updated values of the MW potential, see Tab. 3 for
comparison. We used the sixth-order rotating BRIDGE (to bridge integration
within the cluster to orbit integration of stars in the Galaxy) by setting the
option SPLIT_6TH_SS_M13 with 13 symmetric evaluations of the force per time
step, $\Delta t=0.25$ Myr. We ensured that the relative energy error per time
step was always below $10^{-7}$ (Martínez-Barbosa et al., 2016, 2017;
Portegies Zwart & McMillan, 2018).

Table 3: MW parameters. Values are taken from Allen & Santillan (1991) for M1-M4 and from Irrgang et al. (2013) for M5.  parameter | value M1-4 | value M5
---|---|---
Mass of the bulge ($M_{b}$) | $1.41\cdot 10^{10}\,M_{\odot}$ | $0.95\cdot 10^{10}\,M_{\odot}$
Bulge scale-length ($b_{1}$) | 0.03873 kpc | 0.23 kpc
Disc mass ($M_{d}$) | $8.56\cdot 10^{10}\,M_{\odot}$ | $6.6\cdot 10^{10}\,M_{\odot}$
Disc scale-length 1 ($a_{1}$) | 5.31 kpc | 4.22 kpc
Disc scale-length 2 ($a_{2}$) | 0.25 kpc | 0.292 kpc
Halo mass ($M_{h}$) | $1.07\cdot 10^{11}\,M_{\odot}$ | $0.24\cdot 10^{11}\,M_{\odot}$
Halo scale-length ($a_{3}$) | 12 kpc | 2.565 kpc
lumpy galaxy | NO (M1,2,3) | NO
| YES (M4) |
– lump mass function | power-law (slope -1.6), |
– minimum lump mass | $M_{min}\in(10,10^{3})\,M_{\odot}$ |
– maximum lump mass | $M_{max}\in(10^{4},10^{8})\,M_{\odot}$ |
– lump size | 10 –100 pc |

We estimated the initial position of the Hyades in the Galaxy using its
current position as described in Sec. 2.1 and by integrating backwards (using
the integration method described above, i.e. the same setup as for the forward
cluster integration the orbit of a point mass for the nominal age of the
Hyades of 655 Myr; see Tab. 1). This was done for the two potential MW
parametrisations we considered. The derived trajectory is shown in
galactocentric Cartesian coordinates in Fig. 1 (see Tab. 2 for the initial
coordinates for the Hyaes).

With Galactic parameters we computed the tidal radius of the star cluster. The
tidal radius, $R_{t}$ is here defined as (see Binney & Tremaine 1987)

$R_{t}=\frac{M(R<R_{t})}{3M_{\mathrm{Gal}}}\cdot R_{\mathrm{Gal}}\,,$ (3)

where $M(R<R_{t})$ is the stellar mass within the cluster tidal radius,
$M_{\mathrm{Gal}}$ is the Galactic mass within the cluster orbit, and
$R_{\mathrm{Gal}}$ is the distance of the star cluster from the Galactic
centre. We used an iterative procedure to estimate $R_{t}$ and $M(R<R_{t})$
for a given position of the cluster in the Galactic potential.

Figure 1: Backward-integrated orbit in galactocentric Cartesian coordinates
from the present-day coordinates for the Hyades (red circle) to its likely
birth location (green circle). The backward-integration was performed using a
point mass in the same Galactic potential as the simulated Galactic cluster.

To facilitate comparison of our simulation results with the observed Hyades
star cluster, we followed several steps: 1) The backward-integrated initial
conditions ensured that the simulated star cluster was almost exactly at the
position of the present-day Hyades. However, because the cluster is losing
stars, its orbit is slightly different from the orbit of a single particle. To
correct for this difference, we shifted each simulation snapshot to the
present-day position of the Hyades star cluster and aligned their velocity
vectors, which necessitated a rotation of the system. This operation thus
makes the necessary corrections in projection on the sky without changing
internal kinematics. This is standard practice, but previous papers have
usually shifted the star cluster to the centre (zero) position and aligned the
velocity vector with one of the Cartesian axes such that no rotation was
applied. The advantage of our method is that our results are directly
comparable with the Gaia data in all parameter spaces. 2) The age of the
cluster is uncertain. This presents a clear obstacle because, at a given age,
the cluster might have a certain special feature in one of the projections
that shifts slightly with age. To account for this uncertainty and ensure that
we did not discard valuable stars from the Gaia catalogue, we stacked several
snapshots spanning over 75 Myr from the simulation. 3) We computed the cluster
and its tails in the Galactic (spherical and Cartesian) and galactocentric
(spherical, cylindrical, and Cartesian) coordinate systems and representations
(as CP projected velocities), which are directly comparable with Gaia data. We
note that to the best of our knowledge, this is the first time that
simulations are reported that are fully comparable with Gaia measurements. We
emphasize that the internal kinematics of a star cluster and of its tail in
general depends on the acceleration that is generated along the orbit. For an
orbit with a large variation of the acceleration, that is, for highly
eccentric orbits, the simulation snapshots would therefore differ greatly for
different spatial positions. In the case of the Hyades star cluster, the orbit
is close to circular. The stacking of snapshots along the orbit is therefore
not significantly affected by changes in acceleration along the orbit.

We now discuss the resulting star cluster and its tidal tails for combined
snapshots from 620 Myr to 695 Myr with a time step of 5 Myr time-step for the
non-rotating cluster in the smooth Galactic potential (model M1). This is
plotted for several parameters measured by the Gaia satellite. In Fig. 2 the
star cluster and its tidal tails are plotted. The time-stacking of several
snapshots shows the epicyclic motions of individual stars. Otherwise, the
extent and structure of the tails are comparable to previous simulations
(Ernst et al., 2011; Chumak et al., 2005). We note that the realistic star
cluster trajectory with excursions out of the Galactic plane does not have a
large effect on the physical appearance of the tails. However, it is very
important to consider the projection effects caused by the present-day
position of the cluster and its tail on the sky. That is, whether the tails
are in or out of the Galactic plane changes their appearance on the sky and in
other parameter spaces significantly (proper motions and their derived
projections).

Figure 2: Hyades-like simulated star cluster at the series of snapshots at
620, 625, 630, 635, 640, 645, 650, 655, 660, 665, 670, 675, 680, 685, 690, and
695 Myr in Galactic Cartesian coordinates. The Sun is marked as a yellow
point. The grey arrow points to the Galactic centre, and the black arrow is
the cluster velocity vector in the corresponding coordinates. The time
stacking of snapshots shows the movement of individual stars to and from the
epicyclic overdensities.

Meingast & Alves (2019) used galactocentric cylindrical coordinates centred on
the Hyades cluster with a search radius of 1.5 km/s in order to search for its
tidal tails. Fig. 3 shows that in galactocentric cylindrical coordinates, the
two tidal tails in their entirety indeed extend for only a few km/s. This is
because galactocentric cylindrical coordinates eliminate projection effects
introduced by the extended nature of the tails. This coordinate system
presents the most compact representation of the members of the tail and is
independent of the model.

Figure 3: Simulated star cluster and its tidal tails (M1) plotted in velocity
space in galactocentric cylindrical coordinates. The red cross marks the
central velocity value of the Hyades.

However, Meingast & Alves (2019) recovered only a small part of the tidal
tails in comparison to Röser et al. (2019), who used the CP method (discussed
in more detail below). One reason for this is that radial velocity
measurements are required to compute velocities in galactocentric coordinates
from the Gaia astrometric data. This presents a significant limitation because
only bright stars have radial velocity measurements in the Gaia catalogue
($G<15$ mag), and thus they constitute a very incomplete data sample to work
with. The other aspect is the extension of the tail in velocity space. As
shown in Fig. 3, the central region is dense and thus can be recovered by a
standard clustering algorithm such as DBSCAN (e.g. Jerabkova et al., 2019b),
but the outer parts drop rapidly in density while the contribution from
contamination grows (because the extended tails are farther from the Sun),
making it difficult to distinguish the tails from the field contamination.

Fig. 4 shows the projection of the model tidal tails into space defined by
some selected combinations of Gaia parameters. The core of the cluster is
represented by a dense clump in proper motions (on the sky $\mu^{*}_{\alpha}$,
$\mu_{\delta}$, or distance-corrected $v^{*}_{\alpha}$, $v_{\delta}$) and in
radial velocities. However, the tidal tail is then projected in various shapes
with several clumps or overdensities. Without the knowledge given to us by the
model, it would be very difficult if not impossible to recover these
structures by searching for overdensities in proper motions and radial
velocities in the real dataset. The shape of a cluster and its tails in proper
motions and radial velocities depends on the physical properties of the tails
and on the position in the galaxy. This means that the very same tidal tails
would be represented by different structures in proper motions in the sky
projections if they were located elsewhere.

Figure 4: Modeled (M1) Hyades star cluster and its tail at time-stacked
snapshots from 620 Myr to 695 Myr. Top panel: Proper motion distribution in
physical km/s units. The central clump represents the star cluster, and the
extended distribution (up to 80 km/s) consists of members of the tidal tail.
Middle panel: Proper motion in R.A. and radial velocity distribution showing
the same 80 km/s span for the tidal tails members. Bottom panel: On-the-sky
distribution of proper motions of the modeled cluster and its tails. Because
of the proximity of Hyades to the Sun, the spread in proper motion is large
(see the loose lump representing the cluster members), and the two
overdensities in the shape of tentacles caution against using a blind
clustering algorithm in order to search for tidal tails in the data.

Models are expected to be able to easily test for the epicyclic overdensities,
which may also be used to identify tail members in the Gaia data, as has been
pointed out in a recent study of Oh & Evans (2020). In order to understand the
position of the epicyclic overdensities in different spatial and velocity
space, we used the simulated M1 model and identified the overdensities in
$X-Y$ Galactic Cartesian coordinates, the see upper left panel of Fig. 5. In
the other panels of the same figure, the epicyclic overdensities are plotted
with the same colour-coding.

Figure 5: Epicyclic overdensities depicted in colour in different plots
showing various phase-space coordinates for model M1. We used the $X-Y$ phase-
space to identify epicyclic overdensities. Stars belonging to the epyciclic
overdensities are plotted in colour and appear in various Gaia parameter
spaces.

###  2.3 Compact convergent point method: Recovery of the full extent of tidal
tails

The conclusion so far is that tidal tails are not in general represented by a
simple overdensity in velocity space. While plots showing modelled proper
motions and radial velocities of an open cluster and its tail have been
documented here for the first time, the difficulties of detecting extended co-
moving structures have been known and explored in detail for quite a while
(e.g. van Leeuwen, 2009). Two effects contribute to the appearance of tidal
tails in velocity space. We describe them in detail below.

The first effect is a projection effect. In an imagined set of perfectly co-
moving stars on a circular orbit in the Galaxy, these points will have the
same Cartesian velocities in the Galactic coordinates. However, the on-the-sky
projection effects mean that stars with different positions on the sky will
have different proper motion values. A good solution to this approach is the
well-established CP method (Strömberg, 1939; de Bruijne, 1999; van Leeuwen,
2009, and references therein). The CP method is based on correcting proper
motions for projection effects of velocities on the celestial sphere that are
caused by the spatial extent of a stellar structure. The proper motions (in
km/s taking individual parallaxes into account) of the constituent stars in a
spatially co-moving group point to a so-called CP on the sky. It is then
possible to compute the expected or predicted parallel and perpendicular
velocity components that point to the CP,
$v_{||\mathrm{pred}},v_{\bot\mathrm{pred}}$, for a star with a given position
in space or on the sky by defining the CP, $v_{\bot\mathrm{pred}}=0$. These
values are then compared with the projected values of measured proper motions
(in km/s), $v_{||\mathrm{obs}},v_{\bot\mathrm{obs}}$. For stars that are co-
moving, the differences $v_{||\mathrm{pred}}-v_{||\mathrm{obs}}$ and
$v_{\bot\mathrm{obs}}$ will be close to zero.

The second effect on the appearance of the tidal tails in velocity space
consist of the actual velocity differences between the members of the tail
that are in general caused by their drifting apart because they are on
slightly different orbits in the Galaxy. Individual stars in tidal tails are
indeed not on exactly the same orbits and have non-zero relative velocities
with respect to the cluster. These effects are amplified with time as the
members drift apart from each other. These differences are not corrected for
by the CP method and complicate the detection of the outer parts of tidal
tails and other extended structures. This is shown in the top panel in Fig. 6,
which plots the projected modelled star cluster and its tail (M1) using the
the CP method (in black and for different ages in colour).

Figure 6: Hyades cluster and its tidal tails used to explain the transition
from the classical CP metod to newly introduced CCP method. Top panel:
Classical CP method diagram showing the velocity distribution of M1-simulated
Hyades and its tidal tails projected onto the parallel (horizontal axes) and
perpendicular direction (vertical axes) to the CP. The predicted velocity
component is removed based on the on-the-sky position of individual objects
(see text for more details). The snapshots beginning at 620 Myr to 695Myr are
plotted, each placed at the Hyades position in spatial and velocity space. In
red we show snapshots at 620, 625, 630, 635, 640, 645, and 650 Myr, in black
those at 655 Myr, and in blue the snapshots at 660, 665, 670, 675, 680, 685,
and 690 Myr. This shows the scatter in the CP diagram that is caused by the
time evolution of the tidal tail. The extent of the points on the horizontal
axes reaches 40 km/s. Middle panel: Dependence of
$v_{||\mathrm{pred}}-v_{||\mathrm{obs}}$ as used in the CP method (horizontal
axes in the top panel) as a function of the distance from the cluster centre
($R_{cl}$). The objects farther from the cluster centre move relatively
faster. The red line shows the approximately linear relation denoted as the
function $\Lambda(R-R_{cl})$. Bottom panel: The modelled points in the CP
diagram (top panel) without the velocity trend from the middle panel described
by the function $\Lambda(R-R_{cl}).$ This operation compacts the whole
structure (star cluster plus the tidal tail) within a region spanning a few
km/s. We call this method the CCP method and use it in Sec. 3.1 to search for
the full extent of the Hyades tidal tails in the Gaia data.

The power of the standard CP method has been demonstrated by Röser et al.
(2019). These authors used the Gaia DR2 and found that the tidal tails of the
Hyades cluster extend up to 200 pc from the cluster. The top panel of Fig. 6
can be compared with Fig. 1 of Röser et al. (2019). We used the age stacking
method again and plot several time-snapshots from our simulation in one plot.
The CP plot shows the central dense part (the stellar cluster), and that the
tidal tails extend up to a distance of 20 km/s in velocity space. Röser et al.
(2019) were only able to detect the overdensity above the background to an
extent of up to 5 km/s, and their detected structure is not symmetrical. This
is again a reminder that it is difficult to find tidal tails because they are
most likely not dense enough to be distinguished from the background
contamination. The structure of the tidal tails is a weak function of age in
this diagram.

The middle panel of Fig. 6 shows the dependence of the velocity difference,
$v_{||\mathrm{pred}}-v_{||\mathrm{obs}}$, in the CP method (i.e. the
horizontal axis of the top panel in Fig. 6) on the spatial distance from the
cluster centre. We demonstrate that it is possible to approximate this
relation by the linear function $\Lambda(R-R_{cl})$ for distances in
Galactocentric coordinates, $R_{cl}-R$ denotes distance difference of
individual stars from clusters centre. For the M1 model at a simulation time
of $625$ Myr, the relation has the form

$\Lambda(R-R_{cl})=-0.026\,(\pm 0.005)\,(R-R_{cl})-0.06\,(\pm 0.01)\,.$ (4)

While there is some variation for different time snapshots (see the
uncertainty spread in the brackets above in eq.4), the variations are small
enough for this method to be robust when the tail is developed (i.e. the age
of the cluster exceeds 100 Myr). This means that the farther away the location
of the two tidal tails from the star cluster centre ($R_{cl}$), the larger the
velocity difference. We verified that this correlation is robust and not
sensitive to the type of orbit and age (after detectable tidal tails have
developed) for a cluster at a fixed position and orientation. This is shown in
the bottom panel of Fig. 6: Even when single relation is used for a number of
snapshots that span almost 100 Myr, the final distribution transformed with
eq.4 remains compact. This trend notably remains robust, and Eq. 4 is also
valid for model M5 when the updated MW potential from Irrgang et al. (2013) is
used.

In order to make the most compact representation of the tidal tails, we
removed the $\Lambda(R-R_{cl})$ correlation from the CP method velocity
difference $v_{||\mathrm{pred}}-v_{||\mathrm{obs}}$. The bottom panel of Fig.
6 thus shows the newly introduced compact CP method (CCP). Because Röser et
al. (2019) and other studies have reported that the models (comparable with
our M1 model) agree well with the data, we use the CCP method in order to try
to recover the full extent of the Hyades tidal tails in Sec. 3.1.

We emphasize that the results above (mainly the appearance of tidal tails in
velocity space) demonstrate that tidal tails of open clusters are not in
general represented by any model-independent overdensity in proper motions or
derived quantities (e.g. using the CCP method). This needs to be kept in mind
when searching for tidal tails in Gaia data. A numerical model might even be
necessary in order to be able to find the extended parts and interpret them
correctly.

##  3 Searching for the tidal tails in the Gaia data

With the modelling described above, the newly gained knowledge, and the
developed CCP method, we now attempt to recover the full extent of the
theoretically expected tidal tails around the Hyades in the Gaia data. To do
so, we downloaded Gaia DR2 (Gaia Collaboration et al., 2018b) and the eDR3
data (Gaia Collaboration et al., 2020) from the ESA
archive111https://gea.esac.esa.int/archive/; see the ADQL set of commands in
Appendix A.. The chosen range of parallaxes has the potential that the Hyades
tidal tail members can be identified up to a distance of $500\,$pc from the
Sun, that is, the full length of the tail up to 1000 pc. We queried only
targets with high-quality parallax measurement, paralax_over_error > 10,
allowing us to compute the distance $d$ as

$d/[\mathrm{pc}]=\frac{1000}{\varpi/[\mathrm{mas}]}\,.$

Moreover, we used the RUWE quality criterion to filter out objects with
spurious astrometric solutions (Lindegren, document GAIA-C3-TN-LU-LL-124-01).
We note that previous studies such as Röser et al. (2019) applied a smaller
range of parallaxes, limiting the search volume in the Gaia catalogue and thus
the detectable extent of the tail. The smaller volume, however, also means
lower contamination. The models computed above show that the tail does not
produce a simple overdensity in velocity space, but that the extended parts of
the tails are much harder to recover because the contamination increases with
cubic distance from the Sun. Our search for the full extent of the Hyades
tidal tails, which are comprised of a few hundred stars in the catalogue of
about $10^{8}$ targets is thus very challenging.

The new approach used here therefore is to employ the $N-$body simulations
described above in order to separate the expected members of the Hyades star
cluster and its tails from the Galactic field. The standard CP method already
uses model expectations to define the tidal tails. By using the $N-$body
models to formulate the CCP method, we therefore merely generalise this
approach. For the purpose of this initial study, we used the benchmark M1
model. See also Röser & Schilbach (2020) who extended CP method based on
empirical correlation in a similar way as the CCP method applied here.

Röser et al. (2019) reported that the detected parts of the Hyades tidal tails
agree very well with previously published $N-$body models that are basically
identical to our M1 model. This means that we should at least be able to
detect the parts of the tail reported by Röser et al. (2019) with the M1
model, and that any deviations from the model in the data will indicate the
need for more realistic models.

###  3.1 Hyades and tidal tail selection method

In order to recover the Hyades cluster and its tidal tails from our initial
catalogue, we followed a multi-step procedure, comprising steps S1, S2, and
S3. We describe the steps in detail below.

S1 We used the CCP method as described in Sec. 2.3. In the CCP method, the
stellar cluster and its extended tidal tails are mapped to a clump with a
radius of $1.5$ km/s; see the bottom panel in Fig. 6. The S1 criterion selects
objects from the initial catalogue that have
$-1.5<(v_{\mathrm{||obs}}-v_{\mathrm{||pred}})-\Lambda(R-R_{cl}),v_{\perp\mathrm{obs}}<1.5$
km/s. This procedure reduces our initial DR2 catalogue of 11,223,898 sources
to 22,909 objects and our initial eDR3 catalogue of 13,896,271 sources to
28,218 objects. In spatial coordinates, the leading tail of the Hyades cluster
can be distinguished from the contamination, see Fig. 7. This shows that while
contamination is still present after the S1 step, the method is successful
because it allows the cluster and the tail to emerge from it. In the initial
catalogue, the overdensity visible in the standard CP diagram is not evident,
as is the case for the smaller data set shown in Fig. 1 in Röser et al.
(2019).

Figure 7: Gaia data sample after the S1 CCP criterion was applied, which
filters out $\approx$ 22 000 stars from the initial $\approx$ 11 000 000. Left
panel: Galactic Cartesian coordinates. The yellow point shows the (0,0)
position of the Sun, and the arrows show the velocity of the Hyades in the
$X-Y$ plane and the direction to the GC. The black Gaia data points are
highlighted with the grey model to guide the eye. This suggests an overdensity
to the eye. Right panel: Classical CP method plot. The overdensity at (0,0) is
already visible after the S1 selection.

S2 We used the model-predicted proper motions (top panel of Fig. 4) and
selected only objects with a distance from any modelled proper motion within
1.5 km/s. We used tangential velocities in km/s. The coordinate system can be
either equatorial or galactic. For the purpose of this study, we chose to make
the cut at 1.5 km/s so that targets with an uncertainty of 0.5 km/s on their
proper motion were still accounted for within $3\sigma$. The inclusiveness or
completeness of the selected sample changes when this value is changed, and it
will be interesting to explore this in more detail in a specific future study.
This selection further reduces the number of objects from 22,909 to 2102. We
plot these targets as open orange circles in Fig. 8. They perfectly overlap
(by construction) with the model predictions in Fig. 4. However, some of these
points are still contamination, which is shown in the proper motion plot (the
lump of points near (0,0) mas/yr in Fig. 9). For the eDR3 we used the
benchmark model M5 and the corresponding tangential velocities. The number of
objects from the S1 selection value of 28,218 is further reduced to 1,774.
They are plotted as open blue circles in Fig. 10 and Fig. 14.

Figure 8: Gaia data showing selection criteria S1 and S2. The black background
points show the M1 model, the open yellow circles are Gaia objects that are
consistent with the CCP modelled distribution (i.e. after the S1 selection was
applied) and with the modelled proper motion values (S2 selection). The full
orange points show the Gaia objects that in addition to this, are also
consistent with the spatial distribution of the M1 model (S2 selection). Top
panel: Diagram of the CCP method. Bottom panel: Proper motion distribution in
equatorial spherical coordinates, R.A. and Dec, in physical km/s units.

S3 To clean the sample from the contaminants that were still present after S2
was applied, we used the modelled spatial coordinates and chose only the
targets that overlap with the 3D tidal tail shape within their 3$\sigma$
uncertainty. The finally selected DR2 data sample is plotted as full orange
circles (Fig. 8,9,10) and contains 1109 objects, 411 of which belong to a
leading tail and 331 to a trailing tail (assuming the tidal radius is 9 pc).
For the eDR3 data sample, for which we used the M5 model, the final catalogue
contains 862 objects, 293 of which belong to the leading tail and 166 to the
trailing tail. The data are plotted as filled blue circles in Fig. 10 and Fig.
14. The final data sets can be acquired through the CDS222Centre de Données
astronomiques de Strasbourg, https://cds.u-strasbg.fr/ service and the tables
showing the data format are available in the appendix, see Tab. 5 and Tab. 6.

Figure 9: On-the-sky proper motions in equatorial spherical coordinates R.A.
and Dec in mas/yr. See Fig. 8 and text for the detailed description. The black
points are the M1 model points, the open orange circles show the remaining
objects after the S1+S2 selection was applied, and the filled orange circles
show the objects after the S1+S2+S3 selection.

Figure 10: Tidal tails in Galactic Cartesian coordinates. The large yellow
point at (0,0) in all panels marks the Sun. The black arrow shows the velocity
vector in the given coordinate system, and the grey arrow points to the GC.
See Fig. 8 and text for the detailed description. The black points are the M1
model points, the open orange circles show the remaining points after the
S1+S2 selection was applied, and the filled orange circles show the objects
after the S1+S2+S3 selection.

###  3.2 Members of the tidal tail: RV and CaMD check.

We present a catalogue of candidate members of the extended tidal tails of the
Hyades. The selection after the CCP method was applied is based on model M1,
as described above. We did not use radial velocities (RV) in the process of
tail selection, however. Because we used the underlying $N-$body model, the RV
values of the selected candidate members are expected to follow the
measurements. For the few stars that have Gaia RV measurements, this is indeed
the case, as shown in Fig. 11.

Figure 11: Tangential velocity in R.A. in km/s on the vertical axis and the
radial velocities, RV, on the horizontal axis. See Fig. 8 and text for the
detailed description. The black points are the M1 model points, the open
orange circles are obtained after the S1+S2 selection was applied, and the
filled orange circles remain after the S1+S2+S3 selection criterion was
applied. The radial velocity information is only available for a small number
of objects.

The extended tidal tails can also be detected using the colour-absolute
magnitude diagram (CaMD). This procedure is completely independent from the
method used above. Fig. 12 shows the CaMD for the S1+S2-selected,
S1+S2+S3-selected and the Gaia Hyades CaMD from Gaia Collaboration et al.
(2018a). The agreement of the Hyades Gaia DR2 sample (Gaia Collaboration et
al., 2018a) and our selected objects is encouraging.

Figure 12: Colour-absolute magnitude diagram for S1+S2 -selected objects (open
orange circles), S1+S2+S3-selected objects (full orange circle), and the Gaia
Hyades CaMD from Gaia Collaboration et al. (2018a) as open grey circles. For
the selection criteria S1, S2, and S3, see Sec.3.1. The PARSEC isochrone is
plotted as a black line. The missing binary sequence in our selected points is
most likely due to the RUWE filtering, which partially filters out unresolved
binary stars. Mean error bars computed in $M_{G}$ bins computed based on Gaia
photometric fluxes (see Gaia DR2 primer, issue 1.5) are shown in red at the
right.

In the upper part of the CaMD (for targets brighter than $M_{G}=10$ mag) that
shows the overall more precise photometric values, the scatter in the S3-final
selected Hyades and the corresponding data sample for the tidal tails is
clearly smaller. The S3-final data sample presents candidate members whose
membership of the Hyades cluster should be explored through detailed chemical
tagging. However, the decrease in scatter from selections S2 to S3 provides
supportive evidence that they are associated with the Hyades star cluster.

###  3.3 Detection of epyciclic overdensities

Fig. 5 shows the positions of epicyclic overdensities in various spaces.
Because the tails contain a relatively small number of stars, the
overdensities, if present, might be difficult to spot in real data. In order
to quantify the density structure along the tail in a robust way, we used the
spatial distribution in $X-Y$ space. The histograms along the tail in $X-Y$
space are shown in Fig. 13, accompanied by models M1 and M5. We used the
AstroPy histogram function with the implemented Freedman-Diaconis rule to
choose the optimal bin widths.

The simulation results of M1 and M5 follow qualitatively similar evolution
paths. This results in three clear, symmetric, epicyclic overdensities in the
leading and trailing tail. The detailed tail structure and exact position of
the epicyclic overdensities differ in the models.

In order to recover the Hyades cluster and its tail in the Gaia data, we
followed a multi-step procedure. We discussed above that the CCP method that
implements the velocity gradient along the tail is robust and does not change
from model M1 to M5. The subsequent steps that employ the detailed shape of
the distribution of proper motions then differ among the models, however. Thus
we used models M1 and M5 to search for Hyades tidal tails in the data. The
results are compared in the bottom four panels in Fig. 13. For the case of
DR2, some substructures in the tails are suggested by the data. These become
much clearer for the eDR3 data. Because the position of the overdensities is
better matched by the M5 model using the MW potential parameters from Irrgang
et al. (2013), the signal becomes clearer for this selection.

We used the M5 model to select tidal-tail candidates in the eDR3 data that
belong to the detected overdensities. Their stellar kinematics are shown in
Fig. 14 as open red circles. The position of the stars that belong to the
overdensities is as expected from the epicyclic overdensities. Further
analysis achieved with a large modelling effort would provide more insights,
but the detected substructures in the tidal tails of the Hyades appear to
represent the first observed epicyclic overdensities in an open star cluster
based on the spatial and kinematical properties.

Figure 13: Each pair of panels (from top to bottom) shows the cluster and its
tidal tails in $X-Y$ Galactic coordinates rotated so that the $V_{Y}$
component is horizontal. This allows us to plot the histogram along the tail
and to reduce projection effects caused by the spatial alignment of the tail.
The upper two double panels show models $M1$ and $M5$ at the age of 655 Myr.
The lower two double panels show the data selected based on the $M1$ and $M5$
models.

Figure 14: Candidate epicyclic overdensities plotted in proper motions and
tangential velocities together with the eDR3 final data selection. The top
panel shows tangential velocities, and the bottom panel shows proper motions.
The underlying black points belong to model M5 at an age of 655 Myr. Open blue
circles show the data after the S2 selection was applied based on model M5.
The filled blue points are the final data after the S3 selection. In Fig. 13
two spatial overdensities are clearly visible (highlighted by semi-transparent
vertical blue bars), and stars that belong to them are plotted as filled red
circles (for the overdensity in the leading part of the tail) and as dark red
circles (for the overdensity in the trailing part of the tail).

##  4 Effect of the initial angular momentum and lumpy Galactic potential

We have discussed the importance of using detailed kinematic modelling in
order to collect and interpret the most useful Gaia data. With our novel
model-based method, we are able to identify candidate members of the Hyades
tidal tails up to the largest distances of the full extent of 800 pc and can
confirm the detected asymmetry already reported by Röser et al. (2019). The
$N$-body model and the observations differ in several instances, however. We
appear to be unable to detect the epicyclic overdensities that are a firm
prediction in the evolution of a tail in the Galactic potential. Furthermore,
Röser et al. (2011) pointed out that the cluster has an elevated (super-
virial) velocity dispersion. This has been confirmed by Oh & Evans (2020), but
see the discussion of the value of the velocity dispersion in Sec. 2.1. The
aim of this section is to explore whether the initial angular momentum and/or
lumpy Galactic potential could remove some of the discrepancies between model
and observation.

###  4.1 Initial angular momentum: Cluster rotation

In addition to the initial setup, we also considered simulations with a non-
zero $Z$-component of angular momentum, $L_{z}$, such that the star cluster
initially rotates around its $Z$-axis perpendicular to the Galactic plane. To
initialise the rotating cluster, we followed several steps: 1) We set up a
Plummer model identical to the description in Tab. 2 with average zero angular
momentum. 2) The sign of $L_{z}$ was chosen: when $L_{z}<0,$ the cluster
rotation sense is counter-orbit, and when $L_{z}>0,$ it is co-rotating with
the cluster orbit. 3) We verified the value of the angular momentum for each
star, $L_{z}^{\star}$. When the sign of $L_{z}^{\star}$ is opposite to the
chosen direction of rotation in step 2, the velocity signs are changed to the
opposite values. This procedure ensures that the total energy of the star
cluster remains unchanged while allowing it to have a non-zero total angular
momentum. The aim here is to assess for the first time how and whether an
initially rotating cluster may affect the properties of its tidal tails.

###  4.2 Lumpy Galactic disc

In addition to the smooth semi-analytical Galactic potential, we also
considered local potential fluctuations that might be caused by molecular
clouds or other star clusters, for example. The question we would like to
address with this setup has been illustrated by Röser et al. (2019). These
authors pointed out that the detected tails of the Hyades cluster are
asymmetric: the trailing tail is shorter and less populated. This might be
caused by the interactions with potential perturbations, such as those
generated by giant molecular clouds or star clusters.

Introducing lumps to the smooth Galactic potential described above adds a
number of additional parameters to be considered, such as their number,
density distribution, their orbits, mass, and size distribution. This in-depth
study is beyond the scope of this paper, in which we primarily seek a novel
method to extract likely members of non-compact 6D structures. The effect of a
clumpy Galaxy on these structures deserves an individual consideration based
on the exact setup and the interactions the cluster encounters during its
evolution. For the purpose of assessing the possible effect of lumps on the
tidal tails of the Hyades, we considered the following model. The smooth
Galactic potential and initial stellar cluster setup was identical to the
simulations above, see Tab. 2 and Tab 3. The basic model of the lumpy Galaxy
presented in the AMUSE book (Portegies Zwart & McMillan, 2018) was used for
guidance. We initialised lumps with a power-law number distribution slope of
$-1.6$ (the Salpeter value would be $-2.35$), sampling between 4,000 to 10,000
lumps that are randomly distributed in a region in the Galaxy from 3.5 kpc to
8.5 kpc distance from the GC on corresponding circular orbits. We varied the
minimum, $M_{min}\in(10,10^{3})\,M_{\odot}$, and maximum,
$M_{max}\in(10^{4},10^{7})\,M_{\odot}$, masses of the distribution. Because we
set up the lumps as point masses, their physical size is represented by a
softening length, which we varied from 10 pc up to 100 pc.

These parameters allowed us to study various interactions of tidal tails with
lumps, from close encounters of massive lumps to a large number of small
perturbations. We emphasize that the aim of this study is not to simulate a
realistic distribution of GMCs in the MW, but to understand whether
gravitational interactions with lumps can produce a non-symmetric tidal tail
reminiscent of the observed one.

###  4.3 Initial spin angular momentum of the cluster: Simulation results

The mass loss and thus also the mass and length of the tidal tails are higher
for the cluster model that has an initial spin angular momentum vector in the
same direction as its orbital angular momentum vector in comparison to the
initially non-rotating cluster. On the other hand, the mass loss and the mass
and length of the tidal tails are lower for the cluster that initially rotates
in the opposite direction compared to its orbital angular momentum. The mass
loss for the different models is shown in Fig. 15. The initial spin angular
momentum in our models is small because we initialised the Plummer models,
therefore the effects on the length and orientation of the leading and
trailing tidal tails relative to the cluster centre are evident but barely
statistically significant.

Figure 15: Angular momentum and stellar mass of the tidal tails as a function
of time for simulated star clusters. The left panel shows the evolution of the
$z$ component of the spin angular momentum, $L_{z}$, of the star cluster with
time. We plot models M1, M2, and M3, which were evolved from the same initial
conditions (see Tab. 2) but with different initial spin angular momenta (see
Sec. 2.2). Model M4 considers the same MW potential as M1-M3 with the addition
of lumps, and model M5 uses the MW potential from Irrgang et al. (2013). The
$L_{z}$ values are always computed within the tidal radius of the cluster. The
right panel shows the mass outside of the tidal radius ($R_{t}$) of the
cluster at a give time for the same models (M1, M2, and M3). The black lines
are linear fits to the individual tracks and are used to compute the average
mass-loss rate. The angular momentum of all our modelled clusters given by the
motion in the Galaxy has a positive value.

The initial loss of stars from the cluster carries away most of the initial
spin angular momentum because the cluster immediately begins to de-spin, and
thus the stars farthest from the cluster, those that were lost first, carry
most of the information on the initial cluster spin. Steps towards
understanding the depopulation of the phase-space distribution function of
stars in a rotating satellite have been performed by Piatek & Pryor (1995) and
Kroupa (1997), but much more work on this physical process is needed to fully
understand the configuration of initial spin and subsequent tidal tail. A
large number of simulations needs to be analysed in order to quantify the
correlation between the initial cluster spin and the orientation and extent on
the tidal-tail tips. The results obtained here are thus clearly suggestive of
a potentially powerful method for constraining the initial cluster spin by
measuring the orientation of the tail tips relative to the cluster velocity
vector and their distance from the cluster position, but this particular
problem is beyond the aims of this work, which are to introduce the new CCP
method. We emphasise that while there is an effect caused by cluster rotation
on the formation of tidal tails, there is no indication of asymmetric tidal
tails.

###  4.4 Lumpy Galactic potential: Results

In order to assess whether encounters with lumps might in principle be able to
account for the observed Hyades tidal-tail asymmetry and for its potentially
high velocity dispersion, we performed ten simulations with varied lump
parameters. The aim here was not to perform a detailed analysis of all
possible cases, but to broadly assess the possibilities. In the range of
parameters we explored (see Tab. 3), the main results based on the evolution
of the tidal tails can be divided into three statistical categories that we
list below.

1) No discernible effect on the tidal tails is detected. This is the case for
large effective sizes of the lumps, with only a few or no encounters with
massive ($>10^{5}\,M_{\odot}$) lumps.

2) A discernible effect on the tidal tails becomes evident, that is, the tidal
tails show a wider and more dispersed distribution of their members than in
the smooth Galactic potential model, with an effect on the cluster itself.
This is the result in simulations with many encounters, for which the
effective sizes is not very important, without any close interaction with a
massive ($>10^{7}\,M_{\odot}$) lump. No significant asymmetries between the
leading and trailing tails are produced, and the velocity dispersion of the
cluster is not significantly affected.

3) An essentially complete destruction of the star cluster created by a close
encounter with a massive lump, which leads to a temporal tail asymmetry and an
increased velocity dispersion of the cluster. In one of the simulations for
model M4, a close encounter with a massive ($10^{7}\,M_{\odot}$) lump has been
detected around the age of the cluster at 640 Myr. The star cluster winds up
around the lump and is destroyed. The nominal age of the Hyades is 625 Myr,
therefore this case represents a good comparison, noting that it would mean
that the cluster is in the process of disruption (see as well Oh & Evans,
2020). In the following paragraph we discuss this event in more detail and use
model M4 for comparison with the simulations without lumps discussed above.

Figure 16: Encounter of the cluster and its tidal tails with a massive
Galactic lump. This asymmetrically stretches the tail (note the displaced
positions of the epicyclic overdensities; compare with Fig. 2) and changes the
velocity distribution in the star cluster and in its tails. As the star
cluster is disrupted by the lump, its velocity distribution increases. The
snapshot here combines lump distance and mass so that the velocity dispersion
of the model cluster is increased to the potentially higher observed value in
the Hyades star cluster (see text for more details). This affects the
appearance of the star cluster and of its tail in different parameter spaces.
The change when the CP method is applied in the bottom panel is clearly
visible. The disruption process thus can not only produce a temporary
asymmetry in the tidal tails (prior to complete disruption), but can affect
the detection using the CP method, which can lead to an incomplete recovery of
the tail.

The simulated cluster in the lumpy Galaxy model M4 is destroyed at an age of
640 Myr by being gravitationally torn and sheared by a massive lump. However,
the most interesting part of the simulation from the point of view of the
Hyades stellar cluster is not the actual act of cluster destruction, but the
moments before it, because this part can be compared with the observed Hyades
regarding the reported potentially higher velocity dispersion and tidal-tail
asymmetry. Because the simulated stellar cluster is close to being disrupted,
its velocity dispersion increases with time as the lump approaches, as shown
in the last two columns in Tab. 4. For the purpose of comparison with the
Hyades, we therefore chose the snapshot at an age of 641 Myr for which the
velocity dispersion within 9 pc of the cluster centre is 0.8 km/s (see Tab.
1). This is closest to the observed potentially higher value. Oh & Evans
(2020) have concluded that based on the present-day velocity dispersion, the
Hyades cluster is in a stage of disruption according to kinematic forward-
modelling.

This snapshot is plotted in Fig. 16 showing the $X-Y$ positions in the Galaxy
(Sun centred), tangential velocities in R.A. and Dec in km/s, and in the CP
diagram. These plots can be compared with the same coordinate representation
of simulation M1 (see Fig. 10 for the $X-Y$ positions of the cluster and tail
stars in the Galaxy, Fig. 4 for the proper motion plot, and Fig. 6 for the CP
diagram). In addition to the elevated velocity dispersion, the close presence
of the lump already started stretching the tidal tail. Fig. 16 clearly shows
an asymmetry such that i) the trailing tail is lagging behind more strongly,
with the first epicyclic overdensity being stretched out almost completely.
The proper motion distribution is also affected, and in combination, an
asymmetric CP velocity projection becomes evident. The bottom panel in Fig. 16
shows the effect on the CP velocity distribution. The contamination from the
background increases with the relative velocity difference from the cluster
centre because this corresponds to regions that are more distant from the Sun
because of the large extent of the tail. It is therefore straightforward to
see that because of the asymmetry, the observed distribution would be
reflected in the detected tidal tails when the CP method is used. Fig. 5 shows
that the leading and trailing tails are reversed in the CP distribution. This
means that while the close-lump encounter in the simulation in model M4 is
consistent with the observation of a potentially higher velocity dispersion,
it does not appear to explain the observed asymmetry. For reference, in Fig. 1
of Röser et al. (2019) the first 200 pc of the Hyades leading tail is seen in
the CP diagram between 0 and $-$5 km/s. Our diagram instead shows an under-
density in the same region. Nevertheless, it is possible that if a wider
parameter space or more cluster–lump encounters are explored in the future,
the CP diagram of the Hyades might be matched. More precise astrometric
parameters with future Gaia data releases will provide firmer data
constraints.

More research is clearly required to reach firm conclusions. Based on our
results, however, we suggest that an encounter between the Hyades and a lump
with a mass of about $10^{7}\,M_{\odot}$ appears to be a possible explanation
not only for the measured potentially higher velocity dispersion of the Hyades
(Röser et al., 2011; Oh & Evans, 2020), but possibly also for the asymmetry of
the detected tidal tails (Röser et al., 2019, and also this work). This
explanation may turn out to be implausible, however, because it would require
a lump of this mass to be about 117 pc away from the real Hyades cluster. The
Hyades are located close to the Sun and in an inter-arm region of the Galaxy
that largely lacks massive molecular clouds. While the closest star-forming
molecular clouds are about 100 pc away, their masses are about four orders of
magnitude lower.

###  4.5 Quantitative comparison to our simulations

To make a quantitative comparison of the models M1, M2, M3, M4 described
above, we list in Tab. 4 time-averaged quantities in a Hyades-compatible age
span (620-690 Myr for M1, M2, and M3 and prior disruption times 615-635 Myr
for M4). We show the tidal radius ($R_{t}$), the 3D velocity dispersion
($\sigma$), the stellar mass of the cluster ($M$), and the number of stars and
stellar remnants ($N$) within a specified radius. For the lumpy Galaxy model
we also include the time that is closer to the disruption event caused by the
close interaction with a massive lump. The increase in the velocity dispersion
of the cluster is clearly visible.

parameter | value – M1 | value – M2 | value – M3 | value – M4⋆ | value – M4† | value – M5
---|---|---|---|---|---|---
| $L_{z}^{0}=0$ | $L_{z}^{0}<0$ | $L_{z}^{0}>0$ |  $L_{z}^{0}=0$, lumps |  $L_{z}^{0}=0$, lumps | $L_{z}^{0}=0$
|  |  |  |  | (during disruption) |
$R_{t}$ [pc] | 9.2 $\pm$ 0.1 | 9.5 $\pm$ 0.1 | 7.9 $\pm$ 0.1 | 5.2 $\pm$ 0.2 | 4.9 $\pm$ 0.7 | 9.1 $\pm$ 0.1
$\sigma(r<R_{t})$ [km/s] | 0.55 $\pm$ 0.01 | 0.57 $\pm$ 0.02 | 0.53 $\pm$ 0.01 | 0.85 $\pm$ 0.07 | 0.95 $\pm$ 0.09 | 0.52 $\pm$ 0.01
$\sigma(r<3\,\mathrm{pc})$ [km/s] | 0.62 $\pm$ 0.01 | 0.68 $\pm$ 0.02 | 0.66 $\pm$ 0.05 | 1.2 $\pm$ 0.1 | 1.17 $\pm$ 0.06 | 0.61 $\pm$ 0.01
$\sigma(r<9\,\mathrm{pc})$ [km/s] | 0.55 $\pm$ 0.01 | 0.58 0.01 | 0.52 $\pm$ 0.01 | 0.68 $\pm$ 0.02 | 0.93 $\pm$ 0.27 | 0.53 $\pm$ 0.01
$M(r<R_{t})]$ | 394 $\pm$ 11 | 436 $\pm$ 6 | 249 $\pm$ 2 | 60 $\pm$ 5 | 56 $\pm$ 2 | 380 $\pm$ 7
$M(r<30\,\mathrm{pc})]$ | 536 $\pm$ 7 | 568 $\pm$ 4 | 423 $\pm$ 8 | 191 $\pm$ 10 | 175 $\pm$ 4 | 660 $\pm$ 5
$N(r<R_{t})$ | 704 $\pm$ 26 | 845 $\pm$ 14 | 347 $\pm$ 8 | 70 $\pm$ 7 | 63 $\pm$ 3 | 690 $\pm$ 13
$N(r<30\,\mathrm{pc})$ | 1092 $\pm$ 17 | 1220 $\pm$ 14 | 770 $\pm$ 15 | 363 $\pm$ 20 | 317 $\pm$ 12 | 980 $\pm$ 9
average mass loss [$M_{\odot}$/yr] ⋆⋆ | 0.9 $\pm$ 0.1 | 0.8 $\pm$ 0.1 | 1.1 $\pm$ 0.1 | 1.3 $\pm$ 0.2 | ($=$) 1.3 $\pm$ 0.2 | 0.9 $\pm$ 0.1
Table 4: Simulated present-day properties (mean values between ages of 620 and
695 Myr for M1, M2 and M3 and 615 to 635 Myr for M4; time before cluster
disruption) of the Hyades star cluster. ⋆ The mass loss in model M4 is more
complex than in the other models, see the right panel of Fig. 15, thus the
linear mass-loss average only serves as an orientation value. ⋆⋆ The average
mass loss of each of the models has been calculated by a linear fit to the
cluster mass loss over the simulation time (see Fig. 15, left panel).
$\dagger$ As described in Sec. 4, the M4 model is disrupted by a close
encounter with a $10^{7}\,M_{\odot}$ lump. The disruption event occurs
approximately 645 Myr after the start of the simulation, and we present the
last column of values, which are averages of the values between 638 Myr and
643 Myr. This column shows the increase in velocity dispersion before the
cluster is disrupted.

The left panel of Fig. 15 shows the evolution in the angular momentum of the
star cluster models within different setups (M1, M2, M3, and M4). The two
initially non-rotating setups develop rotation in the direction of orbital
motion (the cluster spin angular momentum points in the same direction as the
orbital angular momentum). For the lumpy Galaxy model, the cluster has a
higher mass loss and thus a lower cluster mass at the nominal age of the
Hyades. This model also develops a non-zero spin angular momentum that faster
returns to the zero value than for the M1 model in conjunction with its
significantly higher mass loss.

The right panel shows the time evolution of the stellar mass that is outside
of the tidal radius, $R_{t}$. As expected, the model in a lumpy Galaxy loses
mass at the highest rate and with the highest dispersion. To compute the
average mass loss, which is also provided in Tab. 4, we approximated each
relation using a linear fit, and the average mass-loss rate is then
represented by its slope.

##  5 Discussion

###  5.1 Comparison with the study of Röser et al. (2019)

Röser et al. (2019) used the CP method to recover a co-moving group of objects
based on correcting their proper motions for on-the-sky projection. They
considered the region of 200 pc around the Sun. In their sample, an
overdensity was found. They used a clustering algorithm to recover objects
belonging to this overdensity. These authors published their data set, which
comprises the leading tail extending up to 170 pc from the centre of the
Hyades with 292 stars, and a trailing tail up to 70 pc from the cluster,
containing 237 stars.

Based on the Gaia IDs, we cross-matched our S1-selected data sample with the
published data set of Röser et al. (2019) and recovered 891 of 1316 targets
(with our full initial catalogue, we recovered 1182 targets). The RUWE quality
filtering removes about 100 targets from the data set of Röser et al. (2019).
The additional missing targets are removed by the S1 cut because they are the
most kinematically incompatible with being tidal tails members based on our
model. One simple example of this behaviour are targets with an opposite
kinematic pattern based on their location in the leading or trailing tail.

In Fig. 17 we show our baseline model (as used in all plots so far) over-
plotted with the cross-matched data set (Röser et al., 2019) and our
S1-selected data sample. It is immediately visible that while there is general
agreement, a number of points from Röser et al. (2019) do not overlap with the
model. As an example, we chose the clustering of points extending to (0,0) and
plotted them in a different colour. This clustering is scattered around the
Hyades and the tidal tail.

Figure 17: Comparison of models and data selection done by Röser et al.
(2019). Top panel: Proper motion distribution in mas/yr plotted for the M1
model as black background points. In red we show targets identified as members
of the Hyades and its tail by Röser et al. (2019) in common with our
S1-selected sample, and yellow points represent a selected subset of these
that creates a clump in proper motion distribution that is offset from the
model prediction. Bottom panel: Spatial distribution of the same targets as in
the top panel, shown in Galactic Cartesian coordinates with its central sun
plotted as a yellow point with a black line. The selected clump in proper
motion (yellow points) is scattered over the tidal tails and partially
overlaps it. The points are selected based on their position in measured
proper motion without the need to introduce any model.

The CaMD does not provide a difference between the clump and the other data
points from Röser et al. (2019). More precise age and chemical composition
estimates are required to investigate the association of these points with the
Hyades. However, the spatial distribution, the proper motions (on-the-sky and
in km/s, and also the radial velocity measurements available for a data sub-
sample) in Fig.17 and Fig.18 suggest that these are more likely contamination
objects and not members of the Hyades cluster and its tail.

Figure 18: Same set of points to complement Fig. 17 (black for model M1 , red
for Röser et al. (2019) cross-matched with our S1-selected data sample, and
yellow for the proper motion in mas/yr selected clump) plotted in tangential
velocities in km/s (top panel) and proper motion in R.A. vs. radial velocities
(bottom panel).

###  5.2 Limitation of the $N$-body integrator

The quantification of the tidal-tail models we performed relies on the
$N-$body integrator Huayno (Pelupessy et al., 2012; Jänes et al., 2014), which
uses an optimal softening length (as described in Sec. 2.2.1, see eq. 1) to
allow integrating close encounters between stars. This is a reasonable
approximation to the overall energy-equipartition process coming from the many
weak stellar encounters, but it fails to take the strong close gravitational
interactions between stars into account. Because these are very rare and tend
to cause stars to be ejected, the overall results as quantified in terms of
the stellar mass loss from the cluster through evaporation are expected to be
reasonably accurate over the time of the simulations. This evaporation, rather
than the strong three-body encounters that eject stars, is the relevant
physical process that populates the tidal tails. We tested whether our models
are able to reproduce the tidal tails reasonably accurately by comparing the
mass-loss rate from our cluster models with models that were computed with the
very advanced and precise Aarseth $N-$body codes (e.g. Ernst et al., 2011).
The models presented here and those of Ernst et al. (2011) have a comparable
mass loss, but the general trend of our models with softening is that the
mass-loss values are slightly lower. This is expected because softening makes
the energy equipartition process less efficient (in the limit of a completely
relaxation-free code, such as is used in galactic dynamics simulations, there
would be no energy equipartition through two-body encounters and thus no mass
loss through this process). The overall morphology of the tidal tails computed
here with the Küpper epicyclic overdensities and length indicates, however,
that the energetics of the evaporating stars are calculated sufficiently
accurately for our purpose in comparison with the Aarseth code.

###  5.3 Effect of the chosen MW potential parameters

We have considered two different MW-like potentials, those by Allen &
Santillan (1991) and Irrgang et al. (2013). In general, a direct comparison of
cluster evolution in different potentials is difficult because the clusters
will have different trajectories and will be exposed to different tidal field.
We have started with the present-day Hyades position and integrated back in a
given potential for the same amount of time to obtain the initial conditions.
We used the same cluster initial mass and radius. Fig. 15 shows that the
models M1 and M5 follow the same angular momentum evolution and have the same
mass loss.

Interestingly, there is a difference in the structure of the tidal tails
between model M1 and M5 for the same time snapshot, as shown in Fig. 13. This
demonstrates the great potential of epicyclic overdensities to constrain the
Galactic potential, as described by Küpper et al. (2008, 2015). The updated
model of the MW potential by Irrgang et al. (2013) reproduces the position of
the detected overdensities in the eDR3 data better than the Allen & Santillan
(1991) potential. This opens the possibility for future work to significantly
constrain the Galactic potential.

###  5.4 Limitations of the model-dependent selection

The standard way of searching for star clusters is to use a clustering
algorithm that identifies overdensities in velocity (often proper motion)
space (Cantat-Gaudin et al., 2020). This method not only allowed the recovery
of stellar clusters, but also of some extended (100 pc long) stellar
structures (Jerabkova et al., 2019b; Kounkel & Covey, 2019b; Beccari et al.,
2020). Removing projection effects in the standard CP method (van Leeuwen,
2009) showed to be very promising in recovering up to 200 pc long (in length
extending from the cluster) tidal tails of open clusters (e.g. Hyades, Röser
et al., 2019). The $N-$body modelling of Hyades-like star clusters on a
realistic orbit in the Galactic potential here has demonstrated that the full
extent of the tidal tails cannot be recovered with any of the methods
mentioned above because it is represented by an extended (overall up to $\pm
40\,$km/s) structure in different velocity spaces, and because the extent
depends on the initial conditions such as the age and position of the cluster
and the tail in the Galaxy relative to the Sun. We therefore developed the new
CCP method that takes the modelled parameters into account (see Sec. 2.3).

The model-dependent CCP search method comes with the same caveat as the
standard CP method for the velocities to recover star clusters. When the
observed system deviates from the model, for example, because an environmental
effect is not taken into account in the model, important parts of the stellar
population may not be evident. This has been demonstrated in the lumpy Galaxy
model M4, which deviates from the smooth Galactic model M1 in all phase
spaces. While this can present a certain obstacle in the search for tidal
tails, the comparison with models becomes necessary in order to understand the
reason for this when the search is incomplete (as it appears to be in the case
of Hyades).

In addition, the identified candidate members of the Hyades and its tail need
to be further verified in detailed chemical studies. This is especially
important for tidal tails with ages of about 500 Myr to 1 Gyr because the
identification based on isochrones becomes very difficult in these cases; see
for example the study by Hawkins et al. (2020) using the LAMOST333LAMOST
survey contains two main parts: the LAMOST ExtraGAlactic Survey (LEGAS), and
the LAMOST Experiment for Galactic Understanding and Exploration (LEGUE)
survey of MW stellar structure. spectroscopic survey to establish the age for
the nearby Pisces-Eridanus stellar stream as 120 Myr and not 1 Gyr, as
estimated before through isochrone fitting.

###  5.5 Relation of period distribution to age using Gaia-derived periods

One way to test whether candidate stars are likely past Hyades members is to
determine whether their spin periods are similar to those of the current
Hyades members, that is, determine whether the spin speeds of these stars are
consistent with the age of the Hyades stars. Curtis et al. (2019) used data
from the Transiting Exoplanet Survey Satellite (TESS) to estimate period
distributions for the Pisces-Eridanus stellar stream. They compared the period
distribution of the stream members with period distributions of stars
belonging to nearby open clusters of different ages. They acquired spin period
distributions for the 120 Myr old Pleiades, for the 670 Myr old Praesepe
cluster, and for the 1000 Myr old NGC 6811. Based on this comparison, they
concluded that the age of the Pisces-Eridanus stellar stream is comparable to
that of the Pleiades, that is, about 120 Myr. This method thus significantly
increases the precision of age estimates for stellar groups with ages
$>100\,$Myr for which low-mass stars (which lie below the turnoff in the
colour-magnitude diagram) can no longer be distinguished from the main
sequence.

Together with chemical tagging (see e.g. Hawkins et al., 2020), this method is
very important for the purpose of testing the membership of the identified
stars in the Hyades tidal tail. There is a difference between estimating the
age of a star cluster and the age of a large-scale co-moving group of stars:
For the cluster, the membership probability is given by the spatial and
kinematic signature, and because of the compactness in phase-space, these
constraints usually allow reducing contamination from the field. For large-
scale structures such as relic filaments and tidal tails, which are close to
becoming the field population, this is not the case. Data in addition to CaMD-
age estimates, such as spin-period–age data, are therefore extremely useful.

A separate dedicated study would be needed to derive the period distribution
from the TESS data, as was done by Curtis et al. (2019). However, rotational
periods have been estimated by Lanzafame et al. (2018) for Gaia DR2 data. We
acquired all the periods from the Lanzafame et al. (2018)
gaiadr2.vari_rotatio_modulation catalogue in our initial data sample, see the
black points in Fig. 19. The orange points in this figure show the available
periods in the catalogue for our final Hyades and tail candidate members: only
a handfull of points are available at present. All these points belong to the
cluster area and are not from the tidal tail. For comparison we plot the
periods derived from TESS data by Curtis et al. (2019).

Figure 19: Stellar spin-rotation periods (for tail candidates after selection
S3), either estimated from TESS data by Curtis et al. (2019) or taken from
Lanzafame et al. (2018), as a function of the Gaia-estimated effective
temperature.

Unfortunately, the incompleteness in the data does not allow us to reach any
conclusions on the tidal tails candidates. However, it is interesting to note
that the Hyades cluster members that have Gaia-estimated spin periods appear
to be more consistent with the period distribution of the Pleiades than with
Praesepe, which is closer in age to the Hyades. The clear conclusion from this
exercise is the need for further analysis and more spin-period and chemical
composition measurements of Hyades and its tail member candidates.

###  5.6 Future prospects and general context

The Gaia early-data release 3 (EDR3), available since December 2020, promises
an improved precision in the determination of the astrometric parameters. This
will significantly improve the confidence in detecting structures in proper
motion values and will help us to better constrain distances. One key question
that might be answered in a follow-up study with more precise astrometric
parameters is whether the epicyclic overdensities are truly absent from the
data.

On the other hand, there is no expectation that radial velocity information
will be available for faint targets ($G>15$ mag) in the future Gaia data
releases. This will always present a limitation of the detection method and
thus enforces the use of the CP method and the newly introduced CCP method,
neither of which relies on the radial velocity information.

Future surveys such as the ESO 4MOST (4-metre Multi-Object Spectroscopic
Telescope) and the ING WEAVE (William Herschel Telescope Enhanced Area
Velocity Explorer), with expected first data approximately in 2025, will
change the situation regarding radial velocities. While these surveys will
provide radial velocities for more targets than Gaia, neither of these
instruments will deliver a full-sky RV survey complete down to $G\approx 20$
mag.

##  6 Conclusions

We performed $N-$body simulations of Hyades-like star clusters on a realistic
orbit in the Galactic potential within the AMUSE computational environment. We
explored the effect of initial cluster rotation and presence of Galactic lumps
in addition to the smooth Galactic potential on the evolution of the tidal
tail. The main results are summarized below.

  1. 1.

The $N-$body simulations show that tidal tails are extended low-density
structures. The contamination by Galactic-field stars means that they are
difficult to recover in their full extent using standard methods that rely on
finding clusters and overdensities in the data. The most naturally clustered
tidal tails appear to be in galactocentric cylindrical velocities, but the
required RV measurements for computing the stellar velocities in
galactocentric coordinates are (and will remain) highly incomplete in the Gaia
catalogue.

  2. 2.

The most powerful method for detecting co-moving structures is the so-called
CP method (van Leeuwen, 2009; Röser et al., 2019). It does not rely on the RV
measurements a priori and removes the projection effect arising from the large
extent on the sky. However, as demonstrated in the top panel of Fig. 6, the
full extent of the tail (out to a distance of 600 pc from the cluster for the
Hyades) is extended (by $\pm 20\,$km/s) in the CP plot. Because the
contamination increases with the third power of the distance from the Sun in
the case of the Hyades, the CP method will therefore not reliably recover the
overdensity for the full extent of the tidal tail.

  3. 3.

We suggest using the model-dependent CCP method. The CCP method uses the
velocity versus distance properties of tidal tail models to formulate a
transformation such that the extended 5D structure becomes a compact
collection of points (see Fig. 6). Using this method, we were able to recover
candidate members of the Hyades tidal tails at its full expected extent of 800
pc. The CCP method is only marginally model dependent because the
transformation remains robust and largely invariant in all models we used.
Moreover, the CP method assumes that the structure being sought consists of a
group of stars with the same velocity vector direction, which introduces less
realistic model assumptions.

  4. 4.

We confirmed the previously noted asymmetry in the detected tidal tail: The
leading tail is significantly more pronounced than the trailing tail, which is
almost absent. We showed that a close encounter with a massive Galactic lump
(or dark matter substructure, as suggested for some halo streams, Bonaca et
al. 2019) could lead to an asymmetric tidal tail.

  5. 5.

In the eDR3 data we recovered spatial overdensities in leading and trailing
tails that are kinematically consistent with being epicyclic overdensities and
thus may present the first such detection for an open star cluster. By using
two different MW potentials, Allen & Santillan (1991) and Irrgang et al.
(2013) in our N-body simulations, we showed that the epicyclic overdensities
of open star clusters are able to provide constraints not only on the cluster
properties, but also on the Galactic potential.

Future Gaia data releases will deliver more precise astrometric measurements,
and with upcoming spectral surveys (4MOST and WEAVE) a significantly improved
6D phase-space map of the solar neighbourhood is expected. With this
information in hand, we will have the opportunity to explore the full
potential of extended Galactic structures. These structures provide not only
valuable clues to the detailed physics of star-formation, but are also a link
to the large scale of the Galactic potential, its tides, and shears. We showed
that detailed modelling is fruitful, if not necessary, for interpreting and
finding extended tidal tails in the Gaia data and for distinguishing them from
the Galactic field star populations that the members of tidal tails are
becoming part of.

###### Acknowledgements.

TJ would like to wholeheartedly acknowledge the encouragement and positiveness
of astronomer, observer and friend Rebeca Galera Rosillo ($\star
1988$-$\dagger 2020$) especially during the pandemic lock-down at La Palma:
Gracias por todo, tu recuerdo permanecerá en el cielo nocturno para siempre.
The authors are thankful to the anonymous referee for useful comments and
suggestions that improved the manuscript. TJ is grateful to Siegfried Röser
for useful discussions and clarification of the CP method implementation. TJ
acknowledges support through the European Space Agency fellowship programme,
the Erasmus+ programme of the European Union under grant number 2017-1-CZ01-
KA203-035562 on La Palma and the support from the HISPK in Bonn. This research
has made use of the SIMBAD database and of the VizieR catalogue access tool,
operated at CDS, Strasbourg, France. This work has made use of data from the
European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia),
processed by the Gaia Data Processing and Analysis Consortium (DPAC, https:
//www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been
provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement.

## References

  * Allen & Santillan (1991) Allen, C. & Santillan, A. 1991, Rev. Mexicana Astron. Astrofis., 22, 255
  * André et al. (2014) André, P., Di Francesco, J., Ward-Thompson, D., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 27
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., SipHocz, B. M., et al. 2018, aj, 156, 123
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Banerjee & Kroupa (2017) Banerjee, S. & Kroupa, P. 2017, A&A, 597, A28
  * Baumgardt & Makino (2003) Baumgardt, H. & Makino, J. 2003, MNRAS, 340, 227
  * Beccari et al. (2020) Beccari, G., Boffin, H. M. J., & Jerabkova, T. 2020, MNRAS, 491, 2205
  * Beccari et al. (2018) Beccari, G., Boffin, H. M. J., Jerabkova, T., et al. 2018, MNRAS, 481, L11
  * Beccari et al. (2017) Beccari, G., Petr-Gotzens, M. G., Boffin, H. M. J., et al. 2017, A&A, 604, A22
  * Binney & Tremaine (1987) Binney, J. & Tremaine, S. 1987, Galactic dynamics (Princeton University Press)
  * Bland-Hawthorn & Gerhard (2016) Bland-Hawthorn, J. & Gerhard, O. 2016, ARA&A, 54, 529
  * Bonaca et al. (2019) Bonaca, A., Hogg, D. W., Price-Whelan, A. M., & Conroy, C. 2019, ApJ, 880, 38
  * Cantat-Gaudin et al. (2020) Cantat-Gaudin, T., Anders, F., Castro-Ginard, A., et al. 2020, A&A, 640, A1
  * Cantat-Gaudin et al. (2019) Cantat-Gaudin, T., Mapelli, M., Balaguer-Núñez, L., et al. 2019, A&A, 621, A115
  * Carballo-Bello et al. (2020) Carballo-Bello, J. A., Salinas, R., & Piatti, A. E. 2020, MNRAS[arXiv:2009.11320]
  * Chumak et al. (2005) Chumak, Y. O., Rastorguev, A. S., & Aarseth, S. J. 2005, Astronomy Letters, 31, 308
  * Curtis et al. (2019) Curtis, J. L., Agüeros, M. A., Mamajek, E. E., Wright, J. T., & Cummings, J. D. 2019, AJ, 158, 77
  * Dabringhausen et al. (2016) Dabringhausen, J., Kroupa, P., Famaey, B., & Fellhauer, M. 2016, MNRAS, 463, 1865
  * de Bruijne (1999) de Bruijne, J. H. J. 1999, MNRAS, 306, 381
  * De Gennaro et al. (2009) De Gennaro, S., von Hippel, T., Jefferys, W. H., et al. 2009, ApJ, 696, 12
  * Dinnbier & Kroupa (2020a) Dinnbier, F. & Kroupa, P. 2020a, A&A, 640, A84
  * Dinnbier & Kroupa (2020b) Dinnbier, F. & Kroupa, P. 2020b, A&A, 640, A85
  * Dinnbier & Walch (2020) Dinnbier, F. & Walch, S. 2020, MNRAS[arXiv:2008.08602]
  * Douglas et al. (2019) Douglas, S. T., Curtis, J. L., Agüeros, M. A., et al. 2019, ApJ, 879, 100
  * Ernst et al. (2011) Ernst, A., Just, A., Berczik, P., & Olczak, C. 2011, A&A, 536, A64
  * Fürnkranz et al. (2019) Fürnkranz, V., Meingast, S., & Alves, J. 2019, A&A, 624, L11
  * Gaburov et al. (2009) Gaburov, E., Harfst, S., & Portegies Zwart, S. 2009, New A, 14, 630
  * Gaia Collaboration et al. (2018a) Gaia Collaboration, Babusiaux, C., van Leeuwen, F., et al. 2018a, A&A, 616, A10
  * Gaia Collaboration et al. (2018b) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018b, A&A, 616, A1
  * Gaia Collaboration et al. (2020) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2020, arXiv e-prints, arXiv:2012.01533
  * Goldman et al. (2013) Goldman, B., Röser, S., Schilbach, E., et al. 2013, A&A, 559, A43
  * Gossage et al. (2018) Gossage, S., Conroy, C., Dotter, A., et al. 2018, ApJ, 863, 67
  * Harfst et al. (2007) Harfst, S., Gualandris, A., Merritt, D., et al. 2007, New A, 12, 357
  * Hawkins et al. (2020) Hawkins, K., Lucey, M., & Curtis, J. 2020, MNRAS, 496, 2422
  * Ibata et al. (2019) Ibata, R. A., Bellazzini, M., Malhan, K., Martin, N., & Bianchini, P. 2019, Nature Astronomy, 3, 667
  * Irrgang et al. (2013) Irrgang, A., Wilcox, B., Tucker, E., & Schiefelbein, L. 2013, A&A, 549, A137
  * Jänes et al. (2014) Jänes, J., Pelupessy, I., & Portegies Zwart, S. 2014, A&A, 570, A20
  * Jerabkova et al. (2019a) Jerabkova, T., Beccari, G., Boffin, H. M. J., et al. 2019a, A&A, 627, A57
  * Jerabkova et al. (2019b) Jerabkova, T., Boffin, H. M. J., Beccari, G., & Anderson, R. I. 2019b, MNRAS, 489, 4418
  * Jeřábková et al. (2018) Jeřábková, T., Hasani Zonoozi, A., Kroupa, P., et al. 2018, A&A, 620, A39
  * Just et al. (2009) Just, A., Berczik, P., Petrov, M. I., & Ernst, A. 2009, MNRAS, 392, 969
  * Karim & Mamajek (2017) Karim, M. T. & Mamajek, E. E. 2017, MNRAS, 465, 472
  * Kharchenko et al. (2009) Kharchenko, N. V., Berczik, P., Petrov, M. I., et al. 2009, A&A, 495, 807
  * Kounkel & Covey (2019a) Kounkel, M. & Covey, K. 2019a, AJ, 158, 122
  * Kounkel & Covey (2019b) Kounkel, M. & Covey, K. 2019b, AJ, 158, 122
  * Kounkel et al. (2018) Kounkel, M., Covey, K., Suárez, G., et al. 2018, AJ, 156, 84
  * Kroupa (1995a) Kroupa, P. 1995a, MNRAS, 277, 1491
  * Kroupa (1995b) Kroupa, P. 1995b, MNRAS, 277, 1522
  * Kroupa (1995c) Kroupa, P. 1995c, MNRAS, 277, 1507
  * Kroupa (1997) Kroupa, P. 1997, New A, 2, 139
  * Kroupa (2005) Kroupa, P. 2005, in ESA Special Publication, Vol. 576, The Three-Dimensional Universe with Gaia, ed. C. Turon, K. S. O’Flaherty, & M. A. C. Perryman, 629
  * Kroupa (2011) Kroupa, P. 2011, in Computational Star Formation, ed. J. Alves, B. G. Elmegreen, J. M. Girart, & V. Trimble, Vol. 270, 141–149
  * Kroupa et al. (2001) Kroupa, P., Aarseth, S., & Hurley, J. 2001, MNRAS, 321, 699
  * Kroupa et al. (2018) Kroupa, P., Jeřábková, T., Dinnbier, F., Beccari, G., & Yan, Z. 2018, A&A, 612, A74
  * Küpper et al. (2015) Küpper, A. H. W., Balbinot, E., Bonaca, A., et al. 2015, ApJ, 803, 80
  * Küpper et al. (2010) Küpper, A. H. W., Kroupa, P., Baumgardt, H., & Heggie, D. C. 2010, MNRAS, 401, 105
  * Küpper et al. (2012) Küpper, A. H. W., Lane, R. R., & Heggie, D. C. 2012, MNRAS, 420, 2700
  * Küpper et al. (2008) Küpper, A. H. W., MacLeod, A., & Heggie, D. C. 2008, MNRAS, 387, 1248
  * Lada & Lada (2003) Lada, C. J. & Lada, E. A. 2003, ARA&A, 41, 57
  * Lanzafame et al. (2018) Lanzafame, A. C., Distefano, E., Messina, S., et al. 2018, A&A, 616, A16
  * Lebreton et al. (2001) Lebreton, Y., Fernandes, J., & Lejeune, T. 2001, A&A, 374, 540
  * Lodieu (2020) Lodieu, N. 2020, Mem. Soc. Astron. Italiana, 91, 84
  * Madsen (2003) Madsen, S. 2003, A&A, 401, 565
  * Makarov et al. (2000) Makarov, V. V., Odenkirchen, M., & Urban, S. 2000, A&A, 358, 923
  * Marks & Kroupa (2011) Marks, M. & Kroupa, P. 2011, MNRAS, 417, 1702
  * Marks & Kroupa (2012) Marks, M. & Kroupa, P. 2012, A&A, 543, A8
  * Martínez-Barbosa et al. (2016) Martínez-Barbosa, C. A., Brown, A. G. A., Boekholt, T., et al. 2016, MNRAS, 457, 1062
  * Martínez-Barbosa et al. (2017) Martínez-Barbosa, C. A., Jílková, L., Portegies Zwart, S., & Brown, A. G. A. 2017, MNRAS, 464, 2290
  * Meingast & Alves (2019) Meingast, S. & Alves, J. 2019, A&A, 621, L3
  * Meingast et al. (2020) Meingast, S., Alves, J., & Rottensteiner, A. 2020, arXiv e-prints, arXiv:2010.06591
  * Odenkirchen et al. (2003) Odenkirchen, M., Grebel, E. K., Dehnen, W., et al. 2003, AJ, 126, 2385
  * Oh & Evans (2020) Oh, S. & Evans, N. W. 2020, MNRAS, 498, 1920
  * Oh & Kroupa (2016) Oh, S. & Kroupa, P. 2016, A&A, 590, A107
  * Oh et al. (2015) Oh, S., Kroupa, P., & Pflamm-Altenburg, J. 2015, ApJ, 805, 92
  * Pelupessy et al. (2012) Pelupessy, F. I., Jänes, J., & Portegies Zwart, S. 2012, New A, 17, 711
  * Pelupessy et al. (2013) Pelupessy, F. I., van Elteren, A., de Vries, N., et al. 2013, A&A, 557, A84
  * Perryman et al. (1998) Perryman, M. A. C., Brown, A. G. A., Lebreton, Y., et al. 1998, A&A, 331, 81
  * Piatek & Pryor (1995) Piatek, S. & Pryor, C. 1995, AJ, 109, 1071
  * Piatti & Carballo-Bello (2020) Piatti, A. E. & Carballo-Bello, J. A. 2020, A&A, 637, L2
  * Portegies Zwart & McMillan (2018) Portegies Zwart, S. & McMillan, S. 2018, Astrophysical Recipes; The art of AMUSE (Iop Publishing Ltd)
  * Portegies Zwart et al. (2009) Portegies Zwart, S., McMillan, S., Harfst, S., et al. 2009, New A, 14, 369
  * Portegies Zwart et al. (2013) Portegies Zwart, S., McMillan, S. L. W., van Elteren, E., Pelupessy, I., & de Vries, N. 2013, Computer Physics Communications, 184, 456
  * Portegies Zwart et al. (2001) Portegies Zwart, S. F., McMillan, S. L. W., Hut, P., & Makino, J. 2001, MNRAS, 321, 199
  * Portegies Zwart & Verbunt (1996) Portegies Zwart, S. F. & Verbunt, F. 1996, A&A, 309, 179
  * Reid et al. (2014) Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2014, ApJ, 783, 130
  * Reino et al. (2018) Reino, S., de Bruijne, J., Zari, E., d’Antona, F., & Ventura, P. 2018, MNRAS, 477, 3197
  * Roeser et al. (2010) Roeser, S., Demleitner, M., & Schilbach, E. 2010, AJ, 139, 2440
  * Röser & Schilbach (2019) Röser, S. & Schilbach, E. 2019, A&A, 627, A4
  * Röser & Schilbach (2020) Röser, S. & Schilbach, E. 2020, A&A, 638, A9
  * Röser et al. (2019) Röser, S., Schilbach, E., & Goldman, B. 2019, A&A, 621, L2
  * Röser et al. (2011) Röser, S., Schilbach, E., Piskunov, A. E., Kharchenko, N. V., & Scholz, R. D. 2011, A&A, 531, A92
  * Strömberg (1939) Strömberg, G. 1939, Popular Astronomy, 47, 172
  * Tang et al. (2019) Tang, S.-Y., Pang, X., Yuan, Z., et al. 2019, ApJ, 877, 12
  * Thomas et al. (2018) Thomas, G. F., Famaey, B., Ibata, R., et al. 2018, A&A, 609, A44
  * Toonen et al. (2012) Toonen, S., Nelemans, G., & Portegies Zwart, S. 2012, A&A, 546, A70
  * van Leeuwen (2009) van Leeuwen, F. 2009, A&A, 497, 209
  * Yao et al. (2017) Yao, J. M., Manchester, R. N., & Wang, N. 2017, MNRAS, 468, 3289
  * Zari et al. (2019) Zari, E., Brown, A. G. A., & de Zeeuw, P. T. 2019, A&A, 628, A123
  * Zari et al. (2018) Zari, E., Hashemi, H., Brown, A. G. A., Jardine, K., & de Zeeuw, P. T. 2018, A&A, 620, A172
  * Zhang et al. (2020) Zhang, Y., Tang, S.-Y., Chen, W. P., Pang, X., & Liu, J. Z. 2020, ApJ, 889, 99

##  Appendix A Gaia queries

To query Gaia data, we used the ESA
archive444https://gea.esac.esa.int/archive/ to acquire Gaia DR2 data in order
to search for the Hyades tidal tails using the following ADQL set of commands:

⬇

SELECT gs.source_id,gs.ra,gs.ra_error,

gs.dec,gs.dec_error,gs.parallax,

gs.parallax_error,gs.pmra,gs.pmra_error,

gs.pmdec,gs.pmdec_error,gs.phot_g_mean_mag,

gs.bp_rp,gs.radial_velocity,

gs.radial_velocity_error, r.ruwe

FROM gaiadr2.gaia_source AS gs,

gaiadr2.ruwe AS r

WHERE (gs.source_id = r.source_id

AND gs.parallax>=2.0

AND gs.parallax_over_error>=10

AND r.ruwe < 1.4)

RUWE in the download script above notes the renormalised unit weight error
that is used to filter out objects with spurious astrometric solutions
(Lindegren, document GAIA-C3-TN-LU-LL-124-01). We used several parallax bins
to query the full range of parallaxes ($\varpi>2\,\mathrm{mas}$), resulting in
a total number of queried objects of 11,223,898.

To query Gaia eDR3 data, we used following set of ADQL commands:

⬇

SELECT gs.source_id,gs.ra,gs.ra_error,

gs.dec,gs.dec_error,gs.parallax,

gs.parallax_error,gs.pmra,gs.pmra_error,

gs.pmdec,gs.pmdec_error,gs.phot_g_mean_mag,

gs.bp_rp,gs.dr2_radial_velocity,

gs.dr2_radial_velocity_error, gs.ruwe

FROM gaiaedr3.gaia_source AS gs,

WHERE (gs.parallax>=2.0

AND gs.parallax_over_error>=10

AND gs.ruwe < 1.4)

The changed catalogue name for the radial velocity parameter and the RUWE
parameter is part of the official eDR3 table. This query results in 13,896,271
objects.

###  A.1 Query the Gaia rotation periods

⬇

SELECT gs.source_id,gs.ra,gs.ra_error,

gs.dec,gs.dec_error,gs.parallax,

gs.parallax_error,gs.pmra,gs.pmra_error,

gs.pmdec,gs.pmdec_error,gs.phot_g_mean_mag,

gs.bp_rp,gs.radial_velocity,

gs.radial_velocity_error, gs.teff_val, r.ruwe,

gv.source_id, gv.best_rotation_period, gv.max_activity_index

FROM gaiadr2.gaia_source AS gs,

gaiadr2.ruwe AS r,

gaiadr2.vari_rotation_modulation as gv

WHERE (gs.source_id = r.source_id

AND gs.source_id = gv.source_id

AND gs.parallax>=2.0

AND gs.parallax_over_error>=10

AND r.ruwe < 1.4)

##  Appendix B Catalogue tables

source ID | R.A. [deg] | Dec [deg] |  $\mu_{\alpha^{*}}$ [mas/yr] |  $\mu_{\delta}$ [mas/yr]
---|---|---|---|---
2171661054009007616 | 322.706 | 51.214 | 13.535 | 11.171
2171661054003013504 | 322.705 | 51.214 | 12.889 | 11.188
1978513488498261888 | 322.009 | 47.876 | 14.500 | 13.828
⋮ | ⋮ | ⋮ | ⋮ | ⋮
Table 5: Final selection of tidal tail members in Gaia DR2 data based on model M1. The full table is available at the CDS. source ID | R.A. [deg] | Dec [deg] |  $\mu_{\alpha^{*}}$ [mas/yr] |  $\mu_{\delta}$ [mas/yr]
---|---|---|---|---
5305130139316222080 | 135.090 | -55.267 | -16.274 | 19.302
2174289161681454720 | 325.689 | 54.0383 | 10.528 | 9.471
5311109828992983168 | 137.422 | -54.213 | -16.763 | 17.373
⋮ | ⋮ | ⋮ | ⋮ | ⋮
Table 6: Final selection of tidal tail members in Gaia eDR3 data based of
model M5. The full table is available at the CDS.