© The Authors 2023

1. Introduction

In general, stars form spatially confined in the densest regions of molecular clouds (Lada & Lada 2003; Allen et al. 2007). After their formation, different processes lead to the loss of stellar members. We briefly describe these processes below.

Early gas expulsion. The gas in the central part of the newly formed star cluster has not been completely converted into stars. When the most massive stars have ignited, the ionising radiation heats up the remaining gas, leading to its removal from the star cluster. The initially virialised mixture of stars and gas turns into a dynamically hot and expanding star cluster. After loosing a substantial amount of members, the remaining star cluster revirialises, now having a larger diameter than at its birth. This process only takes a few crossing times, which are typically of the order of a few million years for open star clusters (e.g., Baumgardt & Kroupa 2007).

Stellar ejections. In the Galactic field, OB stars are observed to move with much higher velocities than the velocity dispersion of the young stellar component in the Galactic field. They are assumed to be ejected form young star clusters either by the disintegration of massive binaries, where one component explodes in a supernova (supernova ejection), leaving the second component with its high orbital velocity, or by close dynamical interactions of multiple stellar systems (e.g., in binary-binary encounters) with energy transfer between the components. By orbital shrinkage of one binary, potential energy is transferred to the other binary, leading to its disintegration and leaving the two components with high kinetic energy behind (e.g., Poveda et al. 1967; Pflamm-Altenburg & Kroupa 2006; Oh & Kroupa 2016). In both cases, the velocity of the star is higher than the escape velocity of the star cluster, allowing the particular star to escape from the star cluster into the Galactic field.

Stellar evolution. Due to stellar evolution, the total mass of the star cluster decreases continuously, reducing the binding energy and the tidal radius of the star cluster. The negative binding energy of stars that are only slightly bound may turn into a positive value.

Evaporation and tidal loss of stars. Gravitationally bound stellar systems embedded in an external gravitational field lose members due to the tidal forces. In the frame of a star cluster, its potential well is lowered by the tidal forces at two opposite points, known as Lagrange points. The velocities of stars in the Maxwell tail of the velocity distribution are sufficiently high to escape through these two Lagrange points. Each of both streams of escaping stars forms an elongated structure pointing in opposite directions; these are the tidal tails. If the size of the star cluster is small enough compared to the spatial change of the external force field, then the two depressions of the star cluster potential at the Lagrange points are equal. This is typically the case for star clusters in the solar vicinity orbiting the Galactic centre. Thus, the two streams of stars through the Lagrange points are equal, and the tidal pattern is expected to be symmetric with respect to the star cluster (e.g., Küpper et al. 2010). Therefore, the two tidal arms should be equal.

Recently, the detailed analysis of Gaia data revealed asymmetries in the tidal tails of the four open star clusters Hyades (Jerabkova et al. 2021), Preasepe, Coma Berenices (Jerabkova et al., in prep.) and NGC 752 (Boffin et al. 2022), challenging the common assumption of symmetric tidal tails. However, as it cannot be predicted through which Lagrange point a star will escape from the star cluster into one of the tidal arms, evaporation can be interpreted as a stochastic process. Thus, different numbers of members in the two tidal tails of an individual star cluster are expected. It is the aim of this work to quantify the degree of asymmetry in tidal tails that is due to the stochastic population of tidal tails and to compare with the observations.

In Sect. 2 we describe how (a)symmetry in tidal tails of star clusters is quantified throughout this work. The data of the observed star clusters and the derived asymmetries are presented in Sect. 3. The evaporation process of the numerical Monte Carlo model is outlined in Sect. 4, including the determination of statistical asymmetries due to stochastic evaporation. Section 5 compares the numerical asymmetries with a theoretical model that interprets the evaporation of stars through the two Lagrange points into the tidal tails as a Bernoulli experiment.

2. Measuring the (a)symmetry of tidal tails

In order to explore the degree of (a)symmetry in the tidal tails, the distance distribution of tail members was quantified. Küpper et al. (2010) performed N-body simulations and calculated the stellar number density as a function of the distance to the star cluster along the star cluster orbit. In their simulations, 21 models covering a range of different initial conditions such as the galactocentic orbital radius or the inclination of the orbit were computed with 65536 particles initially. The authors used the NNODY4 code (Aarseth 1999, 2003), which performs a full force summation over all particles. This high number of particles led to a smooth population of the tidal arms. Thus, only small statistical fluctuations were expected. In order to explore the statistical fluctuations of the member number of tidal tails of open star clusters with only a few hundred stars in the tidal tails, multiple simulations of star clusters resembling the observed ones are required.

In order to avoid effects by uncertainties of the Galactic tidal field, here the actual velocity vector of the star cluster was used as the reference (Fig. 1) and not the distance to the cluster along the orbit of the star cluster. The membership and distance of each star was determined as follows: For each star, the distance was given by its position vector, r_i = |r_i|, in the star cluster reference frame. If its distance was larger than the tidal radius, r_t, it was a tidal tail star. The orientation angle, φ_i, is the angle enclosed by the position vector, r_i, of the star and the actual velocity vector of the star cluster, v_cl. A star was considered to be a member of the leading arm if 0 ≤ φ_i ≤ 90°, and to be a member of the trailing arm if φ_i >  90°. Here, the membership of a star of belonging to either the leading or the trailing arm is a sharp criterion without any probability. For stars that are located very close to the border separating the leading from the trailing arm, the errors of the observed data of the positions of the stars were not considered.

Fig. 1. Sketch of the method for quantifying the (a)symmetry in tidal tails. The angle of the orientation, φi, of each star is the angle enclosed by the position vector, ri, of the particular star with respect to the centre of the star cluster and the actual velocity vector, vcl, of the star cluster.

The normalised asymmetry, ϵ, is given by the difference of the number of members in the leading tail, n_l, and the number of members in the trailing tail, n_t, divided by the total number of tidal tail stars, n,

$$\begin{array}{c}\hfill ϵ=\frac{n_{\mathrm{l}}-n_{\mathrm{t}}}{n_{\mathrm{l}}+n_{\mathrm{t}}}=\frac{n_{\mathrm{l}}-n_{\mathrm{t}}}{n}.\end{array}$$(1)

3. Observed open star clusters

In this section, the structure of the tidal tails of the Hyades, Praesepe, Coma Berenices, and NGC 752 are analysed using the method described in Sect. 2. The present-day stellar distributions have been obtained by the Jerabkova et al. (2021) compact convergent point (CCP) method, which allows the tidal tails to be mapped to their tips. Table 1 summarises the observed values of the star clusters that are of interest for this work (see Kroupa et al. 2022 for details).

The results of the analysis of the observational data are shown in Figs. 2–8. The data for the Hyades are published in Jerabkova et al. (2021), and the detailed analysis of its tidal tails can be seen in Fig. 2. The upper left panel shows the spatial distribution of all stars associated with the Hyades. These are all stars that are identified to be cluster members, that is, have a distance to the cluster centre smaller than the tidal radius (Table 1), and those stars that are identified to be tidal tail members, that is, have a distance to the cluster centre larger than the tidal radius. The leading arm is clearly more populated than the trailing arm. Furthermore, the leading tidal arm also has a higher surface density of tidal tail members, as shown in Fig. 3.

Fig. 2. Hyades. (upper left) Spatial distribution of the members of the Hyades (black), the leading tidal arm (green), and the trailing tidal arm (red) from Jerabkova et al. (2021). The coordinates xMW and yMW refer to a non-rotating rest frame of the Milky Way with a shifted origin at which the position of the Sun lies (0, 0). The positive x-axis points towards the Galactic centre, and the positive y-axis points towards Galactic rotation. The arrow indicates the motion of the star cluster centre in the Galactic rest frame. Upper right: positions of cluster and tidal tidal members in polar coordinates calculated as described in Sect. 2. lower left: radial cumulative distribution of the number of stars in the leading and trailing tidal arm of the Hyades (lower right:) Radial cumulative evolution of the normalised asymmetry (Eq. (1)) considering tidal tail stars with a distance to the cluster centre smaller than or equal to r.

Fig. 3. Surface number density of stars of the Hyades (cluster plus the two tidal arms). Each dot represents one star. The local number density is calculated with the sixth-nearest neighbour method (Casertano & Hut 1985) and is colour-coded. The black arrow indicates the velocity vector of the Hyades star cluster in the Galactic rest frame. The centre of the star cluster is located at the intersection of the two thin dashed lines.

Fig. 4. Praesepe from Jerabkova et al., (in prep.): Similar to Fig. 2.

Fig. 5. Praesepe: Similar to Fig. 3.

Fig. 6. Coma Berenices from Jerabkova et al. (in prep.): Similar to Fig. 2.

Fig. 7. Coma Berenices: Similar to Fig. 3.

Fig. 8. NGC 752 from Boffin et al. (2022): Similar to Fig. 2.

In the upper right panel of Fig. 2, the orientation angle and the distance of all stars are shown using the method described in Sect. 2. Stars with an orientation angle smaller than 90° formally belong to the leading arm, and stars with an angle larger than 90° belong to the trailing arm. Within the tidal radius, the cluster is well represented by a Plummer model (Röser et al. 2011). Thus, in the φ-r diagram, the region between 0 and 180° is fully populated. The region with distances between the tidal radius of 9 pc and ≈50 pc is still uniformly populated, but much sparser than within the bound cluster. This corresponds to a relatively spherically shaped region around the Hyades cluster, which is sometimes called the stellar corona of a star cluster. At distances larger than ≈50 pc, two distinct, sharply confined tidal arms are visible.

The lower left panel of Fig. 2 shows the cumulative number of stars in the two tidal arms separately as a function of the distance to the cluster centre. Up to ≈20–30 pc, the two distributions are nearly equal. Beyond ≈30 pc, the distributions start to deviate from each other. Up to ≈100 pc, the cumulative distributions increase nearly constantly. The total number of stars in the leading arm increases faster than the number of stars in the trailing arm. At ≈100 pc, the two distributions flatten and again have a constant but shallower slope. The lower right panel in Fig. 2 shows the cumulative normalised asymmetry, ϵ(≤r), calculated by Eq. (1).

Figure 4 shows the same data analysis for the Praesepe star cluster (data are taken from Jerabkova et al., in prep.). The lower left panel shows that the radial cumulative distance distribution at small distances to the star cluster centre shows the same functional behaviour as in the case of the Hyades. Within the tidal radius, the cumulative distributions rise rapidly and abruptly become shallower in slope at the tidal radius. However, the cumulative distributions of the two arms are identical up to a distance of ≈50–70 pc. Beyond the stellar corona, the distributions diverge continuously from each other, whereas the trailing arm contains more members than the leading arm, which is different from the Hyades. The surface density is shown in Fig. 5.

Figure 6 shows the analysis of the Coma Berenices star cluster (data were taken from Jerabkova et al., in prep.). The radial cumulative distributions have the same qualitative trend as those of the Hyades. Within the star cluster and the stellar corona, the two distribution functions are identical. Again, beyond a distance of ≈50–70 pc to the cluster centre, the distributions diverge continuously from each other. In this case, the leading arm contains more members. The surface density is shown in Fig. 7.

Figure 8 shows the analysis of the NGC 752 star cluster (data were taken from Boffin et al. 2022). Because the errors of the parallaxes increase with increasing distance of the star cluster, the analysis of the tidal tails in three dimensions is affected by a distortion of the sample in the x–y plane (Boffin et al. 2022). Therefore, the stellar sample is analysed in projection on the sky.

4. Monte Carlo simulations

In order to compare the observed asymmetries of tidal tails with the degree of asymmetry that is due to the stochastic evaporation of stars from star clusters embedded in a Galactic tidal field, a large number of test particle integrations were performed in the Galactic gravitational potential.

4.1. Numerical model

In the simulations, a star cluster was set up as a Plummer phase-space distribution (Plummer 1911; Aarseth et al. 1974), with initial parameters b_Pl being the Plummer radius and the total mass, M_Pl. The centre of mass of the model, r_Pl, moved in a Galactic potential as given in Allen & Santillan (1991). The orbit of each stellar test particle was integrated in the two gravitational fields, the Plummer potential as the gravitational proxy of the star cluster and the full Galactic gravitational potential. Thus, the equations of motion of the star i and of the cluster model are

$$\begin{array}{c}\hfill {\mathit{a}}_i=-\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{Pl}}-\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{MW}},\end{array}$$(2)

$$\begin{array}{c}\hfill {\mathit{a}}_{\mathrm{Pl}}=-\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{MW}}.\end{array}$$(3)

The test particles were distributed in the model cluster potential according to the Plummer phase-space distribution function and were integrated in time in the static Galactic potential and a static Plummer potential, whose origin was integrated as a test particle in the Galactic potential. In a self-gravitating system, these particles evaporate from the cluster, which gained sufficient energy by energy redistribution between the gravitationally interacting particles to exceed the binding energy to the cluster. In the model cluster, all particles were treated as test particles and were integrated in a static potential without gravitational interaction between the particles. As energy redistribution does not occur in this model, the evaporation process was simulated as follows: If the Plummer sphere were set up in isolation, all particles would have negative energy and would be bound to the cluster. As the model cluster was positioned in the Galactic external potential, the combined potential was lowered in two opposite points on the intersecting line defined by the Galactic origin and the cluster centre; these are the Lagrange points. As a consequence, a fraction of the initial set of stars have positive energy with respect to the Lagrange points and lead to a continuous stream of escaping stars across the tidal threshold of the clusters (or práh, according to Kroupa et al. 2022), with individual escape timescales up to a Hubble time. This method was successfully tested in Fukushige & Heggie (2000).

The aim of this work is to quantify the expected distribution of the asymmetry between the two tidal tails for the observed total number of tidal tail members that is due to stochastic evaporation through the Lagrange points. As not all test particles escape from the star cluster into the tidal arms, a few test runs were required to calibrate the total number of test particles, such that the number of escaped stars was nearly equal to the observed number of tidal tail members.

As the stellar test particles perform many revolutions around the star cluster centre before evaporating into the Galactic tidal field, the equations of motion were integrated with a time-symmetric Hermite method (Kokubo et al. 1998). The expressions of the accelerations and the corresponding time derivatives are listed in Appendix A for completeness.

As the total mass of the Plummer model is constant during the simulation, models with minimum and maximum star cluster mass were calculated in order to enclose the mass loss of real star clusters. Röser et al. (2011) estimated an initial stellar mass of the Hyades of 1100 M_⊙ and determined a current stellar mass that is gravitationally bound within the tidal radius of 275 M_⊙. Therefore, in the minimum model, the mass of the Plummer sphere was given by the current stellar mass of the star cluster (Table 1), and in the maximum model, the mass of the Plummer sphere was set to four times the current stellar mass. The Plummer parameter, b_Pl, was identical in both models. Furthermore, these two models also take the possible range of tidal radii into account: a minimum model (current mass) with the smallest tidal radius, and a maximum model (here four times the current mass) with a plausible maximum tidal radius.

The orbits of the four star clusters were integrated backwards in time over their assumed lifetime from their current position (Table 1) to their location of formation. At this position, 1000 randomly created Plummer models were set up for each maximum and minimum cluster. Each configuration was integrated forwards in time over the assumed age of the respective star cluster. Finally, all particles outside the tidal radius were assigned to their corresponding tidal arm using the method described in Sect. 2, and the normalised asymmetry, ϵ, was calculated.

4.2. Statistical asymmetry

The numerically obtained distribution of the normalised asymmetry of all four star clusters is shown by the histograms in Fig. 9 for the minimum models and in Fig. 10 for the maximum models. The resulting mean asymmetry, μ, and the dispersion, σ, of the Gaussian fits are listed in Table 2. In the last column in Table 2, the observed asymmetry is tabulated in units of the respective model dispersion. For example, in the case of the maximum model, the observed asymmetry of the Praesepe cluster is a 1.7σ event.

Fig. 9. Monte Carlo results for the minimum model. Normalised distribution of the asymmetry in the Monte Carlo simulations for each star cluster model. The solid curve is the best-fitting Gaussian distribution function. The vertical solid line marks the observed value (Table 1).

Fig. 10. Similar to analysis as in Fig. 9, but for the maximum model.

Dispersion and mean value of the minimum and maximum agree well. We conclude that the mass of the Plummer model has only a small influence. In the case of the Praesepe, the Coma Berenices, and the NGC 752 star cluster, the observed asymmetries have a probability lower than 3σ and can be interpreted as pure statistical events. If the observed asymmetric tidal tails of the Hyades were solely the result of stochastic evaporation, however, then the asymmetry would at least be a 6.7σ event.

5. Theoretical considerations

In this section, a theoretical stochastic evaporation model is developed and compared with the results of the numerical models of Sect. 4.

5.1. Theoretical distribution function, f(ϵ), of the asymmetry

The evaporation of a star into one of the tidal tails can be treated as a Bernoulli experiment. Let p_l be the probability for a star to end up in the leading tail, and (1 − p_l) the probability to end up in the trailing tail. For a total number n of tidal tail members, the probability that n_l stars are located in the leading arm is given by a binomial distribution,

$$\begin{array}{c}\hfill b_n(n_{\mathrm{l}})=\left(\begin{array}{c}n\\ n_{\mathrm{l}}\end{array}\right)p_{\mathrm{l}}^{n_{\mathrm{l}}}{(1-p_{\mathrm{l}})}^{n-n_{\mathrm{l}}},\end{array}$$(4)

with an expectation value of E(n_l)=np_l and a variance Var(n_l)=np_l(1 − p_l).

According to the theorem by de Moivre-Laplace, the binomial distribution converges against the Gaussian distribution,

$$\begin{array}{c}\hfill g(n_{\mathrm{l}})=\frac{1}{\sqrt{2\pi {\sigma }_{\mathrm{l}}^2}}{e}^{-\frac{{(n_{\mathrm{l}}-{\mu }_{\mathrm{l}})}^2}{2{\sigma }_{\mathrm{l}}^2}},\end{array}$$(5)

for increasing n_l with an expectation value E(n_l)=μ_l = np_l and a variance $\mathrm{V}\mathrm{a}\mathrm{r}(n_{\mathrm{l}})={\sigma }_{\mathrm{l}}^2=n{p}_{\mathrm{l}}(1-p_{\mathrm{l}})$.

The normalised asymmetry, ϵ, is related to the number of members of the leading tail, n_l, by

$$\begin{array}{c}\hfill ϵ=\frac{n_{\mathrm{l}}-n_{\mathrm{t}}}{n}=\frac{2n_{\mathrm{l}}}{n}-1.\end{array}$$(6)

The relation between the distribution function of the normalised asymmetry, f(ϵ), and the distribution of the member number of stars in the leading tail, g(n_l), is given by

$$\begin{array}{c}\hfill f(ϵ)\mathrm{d}ϵ=g(n_{\mathrm{l}})\mathrm{d}{n}_{\mathrm{l}}.\end{array}$$(7)

The distribution function of the asymmetry can then be calculated by

$$\begin{array}{c}\hfill f(ϵ)=g(n_{\mathrm{l}}(ϵ))\left|\frac{\mathrm{d}{n}_{\mathrm{l}}}{\mathrm{d}ϵ}(ϵ)\right|=\frac{1}{\sqrt{2\pi {\sigma }_{ϵ}^2}}{e}^{-\frac{{(ϵ-{\mu }_{ϵ})}^2}{2{\sigma }_{ϵ}^2}},\end{array}$$(8)

with dn_l/dϵ following from Eq. (6) and

$$\begin{array}{c}\hfill {\sigma }_{ϵ}^2=\frac{4{\sigma }_{\mathrm{l}}^2}{n^2}=\frac{4p_{\mathrm{l}}(1-p_{\mathrm{l}})}{n},\end{array}$$(9)

and

$$\begin{array}{c}\hfill {\mu }_{ϵ}=\frac{2{\mu }_{\mathrm{l}}}{n}-1=2p_{\mathrm{l}}-1.\end{array}$$(10)

5.2. Symmetric evaporation

In the case of a symmetric population of the tidal tails, the evaporation probabilities into the two arms are identical, $p_{\mathrm{l}}=p_{\mathrm{t}}=\frac{1}{2}$. The expectation value and the variance are

$$\begin{array}{c}\hfill {\mu }_{ϵ}=0,{\sigma }_{ϵ}=\frac{1}{\sqrt{n}}.\end{array}$$(11)

For all four clusters, the theoretically expected asymmetry due to stochastic evaporation, if the tidal tails are symmetrically populated, is listed in Table 2. In all four cases, the theoretical values are close to the numerically obtained ones.

5.3. Asymmetric evaporation

We considered the case that the evaporation and the distribution processes within the vicinity of the star cluster into the tidal arms are asymmetric. According to Eq. (10), the expectation value of the asymmetry, μ_ϵ, increases, that is, more stars evaporate into the leading arm than into the trailing arm, if the population probability, p_l, of the leading arm increases.

For a given observed asymmetry, ϵ_obs, the probability of this event can be calculated as multiples, k, of the dispersion for the assumed evaporation probability p_l,

$$\begin{array}{c}\hfill k=\frac{{ϵ}_{\mathrm{obs}}-{\mu }_{ϵ}(p_{\mathrm{l}})}{{\sigma }_{ϵ}(p_{\mathrm{l}})}=\frac{\sqrt{n}({ϵ}_{\mathrm{obs}}-2p_{\mathrm{l}}+1)}{2\sqrt{p_{\mathrm{l}}(1-p_{\mathrm{l}})}}.\end{array}$$(12)

Figure 11 shows this function for all four star clusters. The vertical dashed line marks the case of symmetrically populated tidal tails. For example, the Hyades have an observed asymmetry of ϵ_obs = 0.298 (Table 1). In the case of a symmetric evaporation, $p_{\mathrm{l}}=\frac{1}{2}$, the theoretical dispersion is σ = 0.043. Thus, the probabillity of this event is 0.298/0.043 = 6.9σ. This point is marked in Fig. 11 by the intersection of the solid line labelled ’Hyades’ and the vertical dashed line.

Fig. 11. Probability of the observed asymmetry of all four star clusters as a multiple of the dispersion in dependence of the assumed evaporation probability into the leading arm, pl. See Sect. 5.3 for details.

On the other hand, in order to increase the event probability (decreasing k), the probability of evaporation into the leading arm needs to be increased. For a given k, Eq. (12) can be solved for the required evaporation probability, p_l, into the leading arm. The emerging quadratic equation leads to two solutions of p_l,

$$\begin{array}{c}\hfill p_{1,2}=-\frac{a}{2}\pm \sqrt{{\left(\frac{a}{2}\right)}^2-b},\end{array}$$(13)

where

$$\begin{array}{c}\hfill a=-\frac{k^2+n({ϵ}_{\mathrm{obs}}+1)}{k^2+n}\mathrm{and}b=\frac{n{({ϵ}_{\mathrm{obs}}+1)}^2}{4k^2+4n}.\end{array}$$(14)

If the observed asymmetry of the Hyades were a 3σ event, then an evaporation probability into the leading arm of approximately p_l = 0.585 would be required.

6. Discussion and conclusions

The common assumption that the two tidal tails of star clusters that move on nearly circular obits around the Galactic centre evolve equally has recently faced a challenge because the analysis of Gaia data reveals asymmetries in the tidal tails of four nearby open star clusters.

Because the evaporation of stars can be treated as a stochastic process, the normalised difference of the number of member stars of the two tails should follow a distribution function. This distribution has been quantified here by use of Monte Carlo simulations of test particle configurations integrated in the full Galactic potential and compared with a theoretical approach. The theoretical and numerical results are found to agree with each other.

Comparing the individual distribution functions of the asymmetry with the observed ones, we conclude that the observed asymmetry of Praesepe, Coma Berenices, and NGC 752 might be the result of the stochastic nature of the evaporation of stars through the two Lagrange points. On the other hand, the asymmetry of the Hyades is a 6.7σ event. In order to interpret the asymmetry as a 3σ event, asymmetric evaporation probabilities into the leading arm of 58.5% and into the trailing arm of 41.5% are required.

It might be speculated that external effects might lead to an additional broadening of the distribution function of the asymmetry. Assuming a different value for the Galactic rotational velocity or a different position of the Sun than used in this work might not lead to a larger scatter as the 50/50% evaporation probabilities at the two Lagrange points are not expected to vary in an almost flat rotation curve.

A stronger effect on the asymmetry of the tidal tails might be due to local deviations from a logarithmic potential (as required in the case of a flat rotation curve). These local variations can be a result of an interaction with a Galactic bar or spiral arms (cf. Bonaca et al. 2020; Pearson et al. 2017). How strong this effect will be can hardly be estimated and will be explored in further numerical studies. However, the main result here remains solid: the pure evaporation of stars through the Lagrange points basically is a simple Bernoulli process.

If the asymmetric evaporation probabilities have an internal origin, then larger asymmetries in the star cluster potential and the kinematics of the evaporation process are required than Newtonian dynamics can provide (Kroupa et al. 2022). On the other hand, the asymmetry might be due to an external perturbation, for example through the encounter with a molecular cloud (Jerabkova et al. 2021). However, the detailed analysis of the asymmetry in the tidal tails of the Hyades revealed that the process leading to the asymmetry must affect both arms equally in terms of the qualitative structure. The bending of the radial number distributions occurs at the same distance to the star cluster, but with different strength (Fig. 2, lower left panel). If an encounter with an external object had occurred on one side of the star cluster, then it is expected that more than only the total number of tidal tail members is reduced. Instead, the functional form of the cumulative number distribution in both arms would be completely different, and rapid changes of the slopes of the cumulative number distribution would not be expected to occur at the same distance from the star cluster centre in the two tidal tails on opposite sites of the star cluster, as is observed in the tidal tails of the Hyades (Fig. 2, lower left panel).

Acknowledgments

PK acknowledges support through the Grant Agency of the Czech Republic under grant number 20-21855S and the DAAD-East-European-Exchange programme at the University of Bonn.

References

Appendix A: Hermite formulae for the Allen-Santillan potential

The total Galactic gravitational field receives contributions from three components: the Galactic bulge, disk, and halo. The rotationally symmetric gravitational potentials are taken from Allen & Santillan (1991). The time-symmetric Hermite integrator from Kokubo et al. (1998) requires the accelerations a, and jerks j = ȧ of the three components of the Galactic gravitational field, which are listed below as expressions of the position vector, r = (x, y, z), the decomposition vectors ρ = (x, y, 0) and z = (0, 0, z), and the corresponding velocity vectors ṙ = v = (ẋ,ẏ,ż), $\dot{\mathit{\rho }}=(\dot{x},\dot{y},0)$, and ż = (0,0,ż).

The total Galactic potential,

$$\begin{array}{c}\hfill \mathrm{\Phi }(\rho ,z)={\mathrm{\Phi }}_{\mathrm{bulge}}(\rho ,z)+{\mathrm{\Phi }}_{\mathrm{disk}}(\rho ,z)+{\mathrm{\Phi }}_{\mathrm{halo}}(\rho ,z),\end{array}$$(A.1)

the total acceleration,

$$\begin{array}{c}\hfill \mathit{a}={\mathit{a}}_{\mathrm{bulge}}+{\mathit{a}}_{\mathrm{disk}}+{\mathit{a}}_{\mathrm{halo}}=-\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{bulge}}-\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{disk}}-\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{halo}},\end{array}$$(A.2)

and the total jerk,

$$\begin{array}{c}\hfill \mathit{j}={\mathit{j}}_{\mathrm{bulge}}+{\mathit{j}}_{\mathrm{disk}}+{\mathit{j}}_{\mathrm{halo}},\end{array}$$(A.3)

are calculated as the sum of all three Galactic components.

A.1. Bulge

The gravitational potential of the Galactic bulge is described by

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{bulge}}(\rho ,z)=-G{M}_1\frac{1}{\sqrt{{\rho }^2+z^2+b_1^2}}.\end{array}$$(A.4)

The corresponding acceleration is given by

$$\begin{array}{c}\hfill {\mathit{a}}_{\mathrm{bulge}}=-G{M}_1{({\rho }^2+z^2+b_1^2)}^{-\frac{3}{2}}\mathit{r}\end{array}$$(A.5)

and the jerk is given by

$$\begin{array}{c}\hfill {\mathit{j}}_{\mathrm{bulge}}=-G{M}_1{({\rho }^2+z^2+b_1^2)}^{-\frac{3}{2}}(\mathit{v}-3\frac{\dot{\mathit{r}}\bullet \mathit{r}}{{\rho }^2+z^2+b_1^2}\mathit{r}).\end{array}$$(A.6)

A.2. Disk

The gravitational potential of the Galactic disk is modeled by a Myamoto-Nagai disk,

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{disk}}(\rho ,z)=-G{M}_2\frac{1}{\sqrt{{\rho }^2+{(a_2+\sqrt{z^2+b_2^2})}^2}}.\end{array}$$(A.7)

This potential leads to an acceleration of

$$\begin{array}{c}\hfill {\mathit{a}}_{\mathrm{disk}}=-G{M}_2{({\rho }^2+{(a_2+{(z^2+b_2^2)}^{\frac{1}{2}})}^2)}^{-\frac{3}{2}}\end{array}$$(A.8)

$$\begin{array}{c}\hfill (\mathit{\rho }+(1+\frac{a_2}{\sqrt{z^2+b_2^2}})\mathit{z})\end{array}$$(A.9)

and a jerk of

$$\begin{array}{cc}& {\mathit{j}}_{\mathrm{disk}}=3M_{2}G{({\rho }^2+{(a_2+{(z^2+b_2^2)}^{\frac{1}{2}})}^2)}^{-\frac{5}{2}}·\hfill \end{array}$$(A.10)

$$\begin{array}{cc}\hfill & (\mathit{\rho }\bullet \dot{\mathit{\rho }}+(1+\frac{a_2}{\sqrt{z^2+b_2^2}})\mathit{z}\bullet \dot{\mathit{z}})(\mathit{\rho }+(1+\frac{a_2}{\sqrt{z^2+b_2^2}})\mathit{z})\hfill \end{array}$$(A.11)

$$\begin{array}{cc}\hfill & -M_{2}G{({\rho }^2+{(a_2+{(z^2+b_2^2)}^{\frac{1}{2}})}^2)}^{-\frac{3}{2}}·\hfill \end{array}$$(A.12)

$$\begin{array}{cc}\hfill & (\dot{\mathit{\rho }}-\frac{a_2\mathit{z}\bullet \dot{\mathit{z}}}{{(z^2+b_2^2)}^{\frac{3}{2}}}\mathit{z}+(1+\frac{a_2}{\sqrt{z^2+b_2^2}})\dot{\mathit{z}}).\hfill \end{array}$$(A.13)

A.3. Halo

The gravitational potential of the Galactic halo is described by a modified logarithmic potential,

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{halo}}=-\frac{G{M}_{\mathrm{halo}}(\le r)}{r}-\frac{G{M}_3}{(\gamma -1)a_3}{[-\frac{\gamma -1}{1+{\left(\frac{r}{a_3}\right)}^{\gamma -1}}+ln(1+{\left(\frac{r}{a_3}\right)}^{\gamma -1})]}_r^{100},\end{array}$$(A.14)

where

$$\begin{array}{c}\hfill M(\le r)=\frac{M_3{\left(\frac{r}{a_3}\right)}^{\gamma }}{1+{\left(\frac{r}{a_3}\right)}^{\gamma -1}}\end{array}$$(A.15)

is the total halo mass within the Galactocentric radius r. The associated acceleration is

$$\begin{array}{c}\hfill {\mathit{a}}_{\mathrm{halo}}=-\frac{GM(\le r)}{r^3}\mathit{r}\end{array}$$(A.16)

and the associated jerk is

$$\begin{array}{c}\hfill {\mathit{j}}_{\mathrm{halo}}=-\frac{GM(r)}{r^3}(\mathit{v}-\frac{\mathit{r}\bullet \mathit{v}}{r^2}(\gamma -3-\frac{(\gamma -1){\left(\frac{r}{a_3}\right)}^{\gamma -1}}{1+{\left(\frac{r}{a_3}\right)}^{\gamma -1}})\mathit{r}).\end{array}$$(A.17)

A.4. Units and parameters

In Allen & Santillan (1991), the spatial variables of the potential have the unit kiloparsec (kpc), the velocity has the unit 10 km s⁻¹, and the potential has the unit 100 km² s⁻². Thus, the accelerations need to be multiplied by 0.104572 in order to have the unit pc Myr⁻² and their time derivatives, j, by 1.06936 × 10⁻³ to have the unit pc Myr⁻³. The mass has the unit 2.3262 × 10⁷ M_⊙ and is scaled such that the gravitational constant is G = 1. The parameters of the potentials are summarised in Table A.1.

A.5. Hermite formulae for the Plummer sphere

The gravitational potential of the star cluster is described by a Plummer sphere,

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{Pl}}(r)=-G{M}_{\mathrm{Pl}}\frac{1}{\sqrt{r^2+b_{\mathrm{Pl}}^2}},\end{array}$$(A.18)

with an acceleration of

$$\begin{array}{c}\hfill {\mathit{a}}_{\mathrm{Pl}}=-G{M}_{\mathrm{Pl}}{(r^2+b_{\mathrm{Pl}}^2)}^{-\frac{3}{2}}\mathit{r}\end{array}$$(A.19)

and a jerk of

$$\begin{array}{c}\hfill {\mathit{j}}_{\mathrm{Pl}}=-G{M}_{\mathrm{Pl}}{(r^2+b_{\mathrm{Pl}}^2)}^{-\frac{3}{2}}(\mathit{v}-3\frac{\dot{\mathit{r}}\bullet \mathit{r}}{r^2+b_{\mathrm{Pl}}^2}\mathit{r}).\end{array}$$(A.20)

All Tables

All Figures

	Fig. 1. Sketch of the method for quantifying the (a)symmetry in tidal tails. The angle of the orientation, φi, of each star is the angle enclosed by the position vector, ri, of the particular star with respect to the centre of the star cluster and the actual velocity vector, vcl, of the star cluster.
In the text

Fig. 2. Hyades. (upper left) Spatial distribution of the members of the Hyades (black), the leading tidal arm (green), and the trailing tidal arm (red) from Jerabkova et al. (2021). The coordinates xMW and yMW refer to a non-rotating rest frame of the Milky Way with a shifted origin at which the position of the Sun lies (0, 0). The positive x-axis points towards the Galactic centre, and the positive y-axis points towards Galactic rotation. The arrow indicates the motion of the star cluster centre in the Galactic rest frame. Upper right: positions of cluster and tidal tidal members in polar coordinates calculated as described in Sect. 2. lower left: radial cumulative distribution of the number of stars in the leading and trailing tidal arm of the Hyades (lower right:) Radial cumulative evolution of the normalised asymmetry (Eq. (1)) considering tidal tail stars with a distance to the cluster centre smaller than or equal to r.

In the text

Fig. 3. Surface number density of stars of the Hyades (cluster plus the two tidal arms). Each dot represents one star. The local number density is calculated with the sixth-nearest neighbour method (Casertano & Hut 1985) and is colour-coded. The black arrow indicates the velocity vector of the Hyades star cluster in the Galactic rest frame. The centre of the star cluster is located at the intersection of the two thin dashed lines.

In the text

	Fig. 4. Praesepe from Jerabkova et al., (in prep.): Similar to Fig. 2.
In the text

	Fig. 5. Praesepe: Similar to Fig. 3.
In the text

	Fig. 6. Coma Berenices from Jerabkova et al. (in prep.): Similar to Fig. 2.
In the text

	Fig. 7. Coma Berenices: Similar to Fig. 3.
In the text

	Fig. 8. NGC 752 from Boffin et al. (2022): Similar to Fig. 2.
In the text

	Fig. 9. Monte Carlo results for the minimum model. Normalised distribution of the asymmetry in the Monte Carlo simulations for each star cluster model. The solid curve is the best-fitting Gaussian distribution function. The vertical solid line marks the observed value (Table 1).
In the text

	Fig. 10. Similar to analysis as in Fig. 9, but for the maximum model.
In the text

	Fig. 11. Probability of the observed asymmetry of all four star clusters as a multiple of the dispersion in dependence of the assumed evaporation probability into the leading arm, pl. See Sect. 5.3 for details.
In the text