Modelling the M68 stellar stream with realistic mass loss and frequency distributions in angle-action coordinates

TL;DR

This paper introduces a new angle-action coordinate method for simulating stellar streams from globular clusters, capturing variable mass loss and frequency, and applies it to the M68 stream to estimate its age and mass-loss rate.

Contribution

The paper presents a novel, efficient method for simulating stellar streams using angle-action coordinates that accounts for realistic mass loss and frequency distributions, improving upon traditional N-body simulations.

Findings

Stream width exceeds what is expected from cluster stars alone.

Estimated stream age is approximately 3 Gyr.

Cluster mass-loss rate is about 0.496 M_sun per Myr per arm.

Abstract

We develop a new method for simulating stellar streams generated by globular clusters using angle-action coordinates. This method reproduces the variable mass-loss and variable frequency of the stripped stars caused by the changing tidal forces acting on the cluster as it moves along an eccentric orbit. The model incorporates realistic distributions for the stripping angle and frequency of the stream stars both along and perpendicular to the stream. The stream is simulated by generating random samples of stripped stars and integrating them forward in time in angle-frequency space. Once the free parameters are calibrated, this method can be used to simulate the internal structure of stellar streams more quickly than N-body simulations, while achieving a similar level of accuracy. We use this model to study the surface density of the stellar stream produced by the globular cluster M68 (NGC 4590). We select $291$ stars from the Gaia-DR3 catalogue along the observable section that are likely to be members of the stream. We find that the width of the stream is too large to be explained by stars stripped from the cluster alone. We simulate the stream using the present method and include the Gaia selection function and observational errors, and the process of separating the stream stars from the foreground. By comparing these results with the observed data, we estimate the age of the stream, or equivalently the cluster accretion time, to be $3.04_{-0.29}^{+5.63}$ Gyr, and the mass-loss of the cluster to be $0.496 \pm 0.030$ M$_{\odot}$ Myr$^{-1}$ arm$^{-1}$.

Tables

Table 1. Table 1 : Parameters of the M68 stellar stream model.

	​​​​​​ $A$	​​​​​​ $C$	​​​ ${\tau }_1$	​​ ${\tau }_2$
	​​​​​​ ($\ast$)	​​​​​​ ($\ast$)	​​​ (rad)	​​ (rad)
$N_s$	​​​​​​ $128.736$	​​​​​​ $10.957$	​​​ $0.660$	​​ $0.248$
	​​​​​​ (mrad Myr-1)	​​​​​​ (mrad Myr-1)	​​​ (rad)	​​ (rad)
$\mathrm{mean}\left(\Delta {\overline{\Omega }}_1\right)$	​​​​​​ $0.295$	​​​​​​ $0.261$	​​​ $0.812$	​​ $0.456$
$\mathrm{std}\left(\Delta {\overline{\Omega }}_1\right)$	​​​​​​ $0.099$	​​​​​​ $0.021$	​​​ $0.524$	​​ $0.451$

Table 2. (a) Best-fitting parameters of the double exponential model (Eq.  3 ) for the number of stripped stars N s N_{\rm s} , where the units of A A and C C are: ∗ ≡ stars × ( 4 π / 139 rad ) − 1 \ast\equiv{\rm stars}\times(4\pi/139\,{\rm rad})^{-1} , and mean and standard deviation of the frequency along the principal axis of the stream of the stream stars Δ ​ Ω ¯ 1 \Delta\bar{\varOmega}_{1} .

	​​​​​​ $A$	​​​​​​ $C$	​​​ ${\tau }_1$	​​ ${\tau }_2$
	​​​​​​ ($\ast$)	​​​​​​ ($\ast$)	​​​ (rad)	​​ (rad)
$N_s$	​​​​​​ $128.736$	​​​​​​ $10.957$	​​​ $0.660$	​​ $0.248$
	​​​​​​ (mrad Myr-1)	​​​​​​ (mrad Myr-1)	​​​ (rad)	​​ (rad)
$\mathrm{mean}\left(\Delta {\overline{\Omega }}_1\right)$	​​​​​​ $0.295$	​​​​​​ $0.261$	​​​ $0.812$	​​ $0.456$
$\mathrm{std}\left(\Delta {\overline{\Omega }}_1\right)$	​​​​​​ $0.099$	​​​​​​ $0.021$	​​​ $0.524$	​​ $0.451$

Table 3. (b) Mean and covariance matrix of the best-fitting multivariate Gaussian distribution to the stripping points in angle space.

	$\mu$ (mrad)	$\Sigma$ (mrad2)
$\Delta \overline{\alpha }$	$\left(\begin{array}{c}-0.103\\ -1.602\\ -2.255\end{array}\right)$

Table 4. (c) Mean and standard deviation of the best-fitting univariate Gaussian distribution to the frequencies perpendicular to the principal axis of the stream.

	$\mu$	$\sigma$
	(mrad Gyr-1)	(mrad Gyr-1)
$\Delta {\overline{\Omega }}_2$	$-2.465$	$4.150$
$\Delta {\overline{\Omega }}_3$	$-6.883$	$4.239$

Table 5. Table 2 : Scale factors of the Gaussian distributions used to compute the bundle of orbits for M68.

	$r_h$	$\delta$	$\alpha$	$v_r$	${\mu }_{\delta }$	${\mu }_{\alpha \ast }$
	(kpc)	(deg)	(deg)	(km s-1)	(mas yr-1)	(mas yr-1)
${\epsilon }_{\mu }$	0.25	2.5	2.5	0.2	0.05	0.05
	${\rho }_h$	$a_h$	$q_h$	${\beta }_h$
	(M⊙ kpc-3)	(kpc)	-	-
${\epsilon }_{\varphi }$	${10}^6$	$2$	$0.25$	$0.1$

Equations

$$\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{int_orb}}{i}nt_{o}rb}\Delta\bar{\theta}_{j}\mbox{\footnotesize$(t)$}=\kappa\Bigl[\Delta\bar{\varOmega}_{j}\;\!(t-t_{\rm s})+\Delta\bar{\alpha}_{j}\Bigr]\hfill t\geqslant t_{\rm s},$$

$$\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{orbit_gc}}{o}rbit_{g}c}\theta_{i}^{\scalebox{0.5}{\>\!GC}}=\varOmega_{i}^{\scalebox{0.5}{\>\!GC}}t+\theta_{i}^{\scalebox{0.5}{\>\!GC}}\mbox{\footnotesize$(0)$},$$

$$\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{double_exp}}{d}ouble_{e}xp}f\mbox{\footnotesize$(\theta_{r}^{\scalebox{0.5}{\>\!M}})$}=A\left[\exp\!\left(-\frac{\theta_{r}^{\scalebox{0.5}{\>\!M}}}{\tau_{1}}\right)+\exp\!\left(\frac{\theta_{r}^{\scalebox{0.5}{\>\!M}}-2\pi}{\tau_{2}}\right)\right]+C,$$

$$\text{addcontentsline}llasection\text{numberline}\text{string}\text{crtrefnumber}equa{l}_{e}xpequa{l}_{e}xp{\theta }_r^{\text{scalebox}0.5L}\equiv \frac{2\pi {\tau }_1}{{\tau }_1+{\tau }_2}.$$

$$\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{stripped_stars_pdf}}{s}tripped_{s}tars_{p}df}p\mbox{\footnotesize$(\theta_{r}^{\scalebox{0.5}{\>\!M}}\>\!|\>\!A,C,\tau_{1},\tau_{2})$}=N_{\rm c}\,f\mbox{\footnotesize$(\theta_{r}^{\scalebox{0.5}{\>\!M}})$},$$

$$t_{\rm s}=\dfrac{\theta_{r}^{\scalebox{0.5}{\>\!S}}-\theta_{r}^{\scalebox{0.5}{\>\!GC}}\mbox{\footnotesize$(0)$}}{\varOmega_{r}^{\scalebox{0.5}{\>\!GC}}}-t_{\scalebox{0.5}{\>\!D}}.$$

$$N_{\text{scalebox}0.5S}=2T{N}_T.$$

$$\text{addcontentsline}llasection\text{numberline}\text{string}\text{crtrefnumber}linea{r}_{r}elationlinea{r}_{r}elation\Delta {\bar{\Omega }}_i\approx {\lambda }_i\Delta {\bar{J}}_i,$$

$${\theta }_r^{\text{scalebox}0.5RS}\equiv {\theta }_r^{\text{scalebox}0.5S}\mathrm{rem}2\pi .$$

$$\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{bayes}}{b}ayes}p\mbox{\footnotesize$(T|\>\!d)$}=\frac{\mathcal{L}\mbox{\footnotesize$(d\>\!|T)$}\,p\mbox{\footnotesize$(T)$}}{p\mbox{\footnotesize$(d\>\!)$}},$$

$$\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{likelihood}}{l}ikelihood}\mathcal{L}\mbox{\footnotesize$(d\>\!|T)$}\equiv\prod_{\Delta\bar{\theta}_{1}\in d}\,p_{\scalebox{0.5}{\>\!SD}}\mbox{\footnotesize$(\Delta\bar{\theta}_{1}|T)$},$$

$$\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{max_dkl}}{m}ax_{d}kl}\hat{T}\equiv\underset{T}{\mathrm{argmax}}\>\!\big(p\mbox{\footnotesize$(T|\>\!d)$}\big).$$

$$R\simeq \left(\begin{array}{ccc}0.689& -0.483& 0.540\\ -0.336& 0.448& 0.829\\ -0.642& -0.752& 0.147\end{array}\right).$$

$$R_{\text{scalebox}0.65TM}\simeq \left(\begin{array}{ccc}0.693& -0.468& 0.549\\ -0.367& 0.426& 0.827\\ -0.620& -0.775& 0.124\end{array}\right).$$

$$\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{sky_R}}{s}ky_{R}}\begin{pmatrix}\cos\!\mbox{\footnotesize$(\phi_{1})$}\cos\!\mbox{\footnotesize$(\phi_{2})$}\\ \sin\!\mbox{\footnotesize$(\phi_{1})$}\cos\!\mbox{\footnotesize$(\phi_{2})$}\\ \sin\!\mbox{\footnotesize$(\phi_{2})$}\\ \end{pmatrix}=R_{1}\times\begin{pmatrix}\cos\!\mbox{\footnotesize$(\alpha)$}\cos\!\mbox{\footnotesize$(\delta)$}\\ \sin\!\mbox{\footnotesize$(\alpha)$}\cos\!\mbox{\footnotesize$(\delta)$}\\ \sin\!\mbox{\footnotesize$(\delta)$}\\ \end{pmatrix}.$$

$$L_{1}\mbox{\footnotesize$(\alpha_{x},\alpha_{y})$}\equiv\sum\big(\phi_{2}-\operatorname{mean}\!\mbox{\footnotesize$(\phi_{2})$}\big)^{2},$$

$$R_1\simeq \left(\begin{array}{ccc}-0.997& -0.073& 0.025\\ -0.007& -0.237& -0.971\\ 0.077& -0.969& 0.236\end{array}\right).$$

$$\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{lk_2}}{l}k_{2}}L_{2}\mbox{\footnotesize$(\alpha_{x},\alpha_{y})$}\equiv\sum\phi_{2}^{2}.$$

$$R_2\simeq \left(\begin{array}{ccc}-0.978& -0.198& 0.062\\ 0.014& -0.364& -0.931\\ 0.207& -0.910& 0.359\end{array}\right).$$

$$\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{int_period}}{i}nt_{p}eriod}\int_{\theta_{r}^{\scalebox{0.5}{\>\!Inf}}}^{\theta_{r}^{\scalebox{0.5}{\>\!Sup}}}\!f\mbox{\footnotesize$(\theta_{r}^{\scalebox{0.5}{\>\!M}})$}\,d\theta_{r}^{\scalebox{0.5}{\>\!M}}=Ah\tau_{1}+Ag\tau_{2}\exp\!\left(-\frac{2\pi}{\tau_{2}}\right)+C(\theta_{r}^{\scalebox{0.5}{\>\!Sup}}-\theta_{r}^{\scalebox{0.5}{\>\!Inf}}),$$

$$\begin{array}{rl}h& \equiv \exp \left(-\frac{{\theta }_r^{\text{scalebox}0.5Inf}}{{\tau }_1}\right)-\exp \left(-\frac{{\theta }_r^{\text{scalebox}0.5Sup}}{{\tau }_1}\right)\\ g& \equiv \exp \left(\frac{{\theta }_r^{\text{scalebox}0.5Sup}}{{\tau }_2}\right)-\exp \left(\frac{{\theta }_r^{\text{scalebox}0.5Inf}}{{\tau }_2}\right).\end{array}$$

$$N_{\rm c}^{-1}=\int_{\theta_{r}^{\>\!\scalebox{0.6}{$-T$}}}^{\theta_{r}^{\>\!\scalebox{0.61}{$0$}}}f\mbox{\footnotesize$(\theta_{r}^{\scalebox{0.5}{\>\!M}})$}\,d\theta_{r}^{\scalebox{0.5}{\>\!M}}.$$

$$N_{\rm c}^{-1}=\int_{2\pi-\theta_{r}^{\scalebox{0.5}{\>\!R}}}^{2\pi}f\mbox{\footnotesize$(\theta_{r}^{\scalebox{0.5}{\>\!M}})$}\,d\theta_{r}^{\scalebox{0.5}{\>\!M}}+N_{\rm p}\!\int_{0}^{2\pi}\!f\mbox{\footnotesize$(\theta_{r}^{\scalebox{0.5}{\>\!M}})$}\,d\theta_{r}^{\scalebox{0.5}{\>\!M}}+\int_{0}^{\theta_{r}^{\>\!\scalebox{0.61}{$0$}}}\!f\mbox{\footnotesize$(\theta_{r}^{\scalebox{0.5}{\>\!M}})$}\,d\theta_{r}^{\scalebox{0.5}{\>\!M}},$$

$$\begin{split}N_{\rm p}&\equiv\big|\theta_{r}^{\>\!\scalebox{0.6}{$-T$}}\big|\,\operatorname{mod}\,2\pi\\[2.5pt] \theta_{r}^{\scalebox{0.5}{\>\!R}}&\equiv\big|\theta_{r}^{\>\!\scalebox{0.6}{$-T$}}\big|\;\!\operatorname{rem}\;\!2\pi.\end{split}$$

$$G_{\rm L}\mbox{\footnotesize$(x\>\!|\>\!\mu,\sigma)$}\equiv\sigma^{-1}\exp\!\big(\>\!z-e^{z}\big),$$

$$\mathrm{norm}(\Delta {\bar{\Omega }}_1)\equiv \frac{\Delta {\bar{\Omega }}_1-\mathrm{mean}(\Delta {\bar{\Omega }}_1)}{\mathrm{std}(\Delta {\bar{\Omega }}_1)}.$$

$$\int_{\hat{\mu}}^{\lambda_{\rm sup}}p\mbox{\footnotesize$(x)$}\,dx=\sigma\int_{\hat{\mu}}^{\infty}p\mbox{\footnotesize$(x)$}\,dx,$$

$$\int_{\lambda_{\rm inf}}^{\hat{\mu}}p\mbox{\footnotesize$(x)$}\,dx=\sigma\int^{\hat{\mu}}_{-\infty}p\mbox{\footnotesize$(x)$}\,dx.$$

$$\operatorname{logN}\!\mbox{\footnotesize$(x|\mu,\sigma,s)$}\equiv\frac{1}{\sqrt{2\pi}\!\>(x-\mu)\!\>s}\exp\!\left[\frac{-\log\!\left(\frac{x-\mu}{\sigma}\right)^{2}}{2s^{2}}\right],$$

Full text

††pubyear: 2025††pagerange: Modelling the M68 stellar stream with realistic mass loss and frequency distributions in angle-action coordinates–J

1 Introduction

Stellar streams are formed by stars that have been stripped from a progenitor cluster by tidal forces caused by the host galaxy. Recently, thanks to the data provided by the Gaia mission (Gaia Collaboration et al., 2016), many streams have been discovered in the Solar neighborhood (Bonaca & Price-Whelan, 2025). These streams can provide valuable information about the history and evolution of the Milky Way (e.g. Deason & Belokurov, 2024) and its characteristics, especially its potential (e.g. Hunt & Vasiliev, 2025). In addition, stellar streams can be used to constrain the properties of the Galactic satellites such as the Large Magellanic Cloud (LMC) (e.g. Koposov et al., 2023) or the spatial and mass distribution of the sub-halos predicted by cosmological simulations (e.g. Bullock & Boylan-Kolchin, 2017).

Stellar streams can also be used to study the formation and evolution of their progenitors, as well as to estimate their current properties. For this purpose, the streams generated by globular clusters offer a unique opportunity because they are dynamically cold, the gas within the stream is expected to be negligible, and in general, they present a simpler structure than the large streams generated by dwarf galaxies such as Sagittarius (e.g. Ibata et al., 1994). The most promising streams are long, close to the Sun, and have no significant perturbations. The best examples are the streams of M68 (Palau & Miralda-Escudé, 2019), NGC 3201 (Palau & Miralda-Escudé, 2021) and Palomar 5 (e.g. Odenkirchen et al., 2001). In such cases, we can use the length and density of the stream to estimate the accretion time and to describe the process of mass loss that the cluster has undergone during its evolution within the potential of the Milky Way.

For a cluster following a general eccentric orbit in the halo of its host galaxy, the density of its stream depends on several factors. Firstly, the shape of the stream changes along the orbit. Its curvature is maximum at the apocentres and minimum at the pericentres. Similarly, the volume of the stream is compressed and expanded at these locations, respectively (e.g. Hozumi & Burkert, 2015). In addition, the density of the stream depends on the mass loss of the progenitor cluster. This is maximal at the pericentres, where the tidal perturbations are strongest. This implies that the stripped stars form clumps that are generated at regular intervals along the cluster orbit (e.g. Dehnen et al., 2004). The mass loss is also a function of the mass of the progenitor and its density distribution, which evolve as stars are lost along the orbit. Another significant factor is the relative motion of the stream stars with respect to the cluster orbit. Due to small variations in angular phase between the stream stars, their relative motion create over- and under-densities within the stellar stream. This phenomenon is known as epicyclic density variations (e.g. Just et al., 2009; Küpper et al., 2010, 2012). Finally, the relative velocity of the stream stars with respect to the progenitor makes the stream longer and less dense over time.

Boundary cases allow us to minimise or eliminate some of these factors, making it easier to study the remaining ones. For example, in a circular orbit, the curvature and compression of the stream along the orbit as well as the tidal forces acting on the cluster are constant. In this case, we can isolate the variations in mass loss caused by the evolution of the mass and density profile of the cluster, the epicyclic density variations, and the growth of the stream length over time. A more general approach is to describe the stream in angle-action coordinates (Binney & Tremaine, 2008, hereafter BT08). In this case, the stream has a rectilinear shape (Eyre & Binney, 2011). This greatly simplifies its analysis by eliminating the geometry variations along the orbit that are present in other coordinate systems. In addition, by rotating the angle-action coordinates to align them with the principal axes of the stream, the projection distortions are minimised and each coordinate can be studied independently. Furthermore, since the stars follow uniform rectilinear motions, the stream does not present epicyclic density variations. Consequently, the density of the stream depends on the mass loss of the cluster, the frequency distribution of the stream stars, and the time since its inception.

Previous studies have modelled stellar streams using angle-action coordinates. For example, Sanders (2014) and Bovy (2014) modelled a stellar stream under the assumption of uniform mass loss from the cluster. Similarly, they assumed that the frequency distribution of the stripped stars is independent of the stripping time. These assumptions are reasonable, given that these models were created to constrain the Galactic potential. In such cases, it is unnecessary to accurately model the details of the stream surface density, since only its large-scale structure needs to be considered. Conversely, Fardal et al. (2015) included a variable mass loss from the progenitor cluster. This variation is caused by the changing tidal forces acting on the cluster as it moves along an eccentric orbit. However, the stream was not fully modelled in angle-action space. Building upon these studies, we present a new methodology for modelling stellar streams in angle-action coordinates. This takes into account both the variable mass loss of the progenitor cluster and the variable frequency distributions of the stripped stars. Such improvements are generally necessary to enable the model to be compared with real stellar streams.

In this paper, we use this methodology to study the stream density of the globular cluster M68 (NGC 4590). We use the methodology, notation, and the N​​ -body simulation introduced in (Palau et al., 2025, hereafter P25). In this simulation, the globular cluster is modelled using a realistic stellar population that follows an isochrone generated by PARSEC/COLIBRI (§2.1, P25). This population is then distributed in phase-space according to a generalised King model (§2.2, P25). The orbit of the cluster is integrated backwards during $T=1.5$ Gyr from the current position of the cluster within an axisymmetric mass model of the Milky Way (§2.4, P25). This model consists of a spherical bulge, a Miyamoto-Nagai disc, and a dark halo modelled with a Navarro, Frenk & White (NFW) density profile. The stars are located at the final position of the inverse integration and then integrated forward using the N​​ -body code PeTar (§2.3, P25). The cluster follows an eccentric regular S-Tube orbit and passes through three pericentre passages (§2.5, P25).

In this paper we denote the action-angle coordinates by $(\theta_{i},J_{i})$ and the frequencies by $\varOmega_{i}$, where $i=\{r,\phi,z\}$. We compute their value using the Averaged Generating Function (AvGF) method implemented in the Python package galpy (§3, P25). When the coordinates and frequencies are given relative to the centre of the globular cluster, they are denoted by $(\Delta\theta_{i},\Delta J_{i})$ and $\Delta\varOmega_{i}$. We apply a rotation in phase-space to align angles, actions and frequencies with the principal axes of the stream (Tremaine, 1999). The rotation is described in Appendix A. The rotated coordinates are denoted by $(\Delta\bar{\theta}_{j},\Delta\bar{J}_{j})$ and the frequencies by $\Delta\bar{\varOmega}_{j}$, where $j=1$ corresponds to the direction aligned with the principal axis of the stream and $j=\{2,3\}$ correspond to the perpendicular directions.

In Section 2, we describe the mass loss of the globular cluster along its orbit and present the internal structure of the stellar stream based on the N​​ -body simulation. In Section 3, we model the stellar stream in angle-action coordinates. The first application of this model is to determine the accretion time and the mass loss of the M68 globular cluster. To achieve this, we compare the surface density of the stream predicted by the model to the observed stream star data. In Section 4 we present our selection of stream star candidates from the Gaia catalogue. Section 5 describes how we generate mock observations from the stream model. These mocks include the Gaia selection function, observational uncertainties, and the process of separation of stream stars from the foreground that we use to identify the real stream stars. In Section 6 we compare the observed surface density with the mock distribution and estimate the accretion time of the cluster and the mass loss with a Likelihood technique. Finally, we present our conclusion in Section 7.

2 Mass loss along an eccentric orbit

The mass loss of the progenitor cluster is proportional to the tidal forces exerted by the host potential (§7.5.6, BT08). For a cluster with an eccentric orbit, the tidal forces increase to a relative maximum at the pericentres. Thus, the cluster undergoes a tidal shock at each pericentre passage and the number of stars stripped peaks. The tidal forces generated by the disc crossings are stronger, but they act for a short time and only result in an increase in the velocity dispersion of the cluster (§8.2.2.f, BT08). For a spherical Galactocentric radius r\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\gtrsim}}}}}R_{\scalebox{0.7}{\odot}}, we consider this effect negligible (Vesperini & Heggie, 1997).

In the top left panel of Figure 1 we show the mass of M68 $M_{\rm gc}$ as a function of the simulation time $t$, where $t=-1.5$ Gyr corresponds to the start of the simulation, and $t=0$ corresponds to the current time. In the lower left panel we show the mass loss of the progenitor cluster $M_{\rm loss}$ in bins of $7.5$ Myr. The mass loss as a function of $t$ is estimated from the stripping time $t_{\rm s}$ of each stream star (Eq. 8, P25). The red, blue, and green vertical dashed lines mark the times of the three pericentre passages that the cluster has undergone since the start of the simulation $t_{\rm peri}\equiv(-1353.59,\,-896.13,\,-439.64)$ Myr. The mass loss begins when the cluster approaches the first pericentre and is heated by the tidal shock. Previously, the number of stripped stars is zero. Between the pericentres, the mass loss is approximately constant at \mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\sim}}}}}0.5 M{}_{\scalebox{0.7}{\odot}} Myr*-1*, and it increases to \mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\sim}}}}}5.5 M{}_{\scalebox{0.7}{\odot}} Myr*-1* at the peaks. The mass loss over the last $50$ Myr is not accounted for, as the stripped stars are close to the cluster and are not classified as stream stars (App. D, P25). During the $1.5$ Gyr of the simulated orbit, the cluster loses about $0.34$ per cent of its initial mass ($M_{\rm ini}=1.25\times\!10^{5}$ M{}_{\scalebox{0.7}{\odot}}) at each pericentre passage, and about $1.27$ per cent in total.

In general, for a cluster with an eccentric orbit, the stripped stars at each pericentre generate a pair of symmetric stellar streams (e.g. Dehnen et al., 2004). Each pair expands along the principal axis of the stream, one at a higher frequency (leading arm) and the other at a lower frequency (trailing arm) with respect to the progenitor cluster. In this way, the tidal stream is formed by a series of small internal streams. In cartesian Galactocentric coordinates, the internal streams appear separated from each other when the progenitor is close to the apocentre, and they overlap when the cluster is close to the pericentre (Hozumi & Burkert, 2015). In angle-action coordinates aligned with the stream, the internal components are clearly visible in the subspace $\{\Delta\bar{\theta}_{1},\Delta\bar{\varOmega}_{1}\}$, which corresponds to the angle and frequency aligned with the principal axis of the stream. In contrast, the perpendicular frequencies $\Delta\bar{\varOmega}_{2}$ and $\Delta\bar{\varOmega}_{3}$ have a weaker or no internal structure. Their distribution along the stream is shown in Appendix B.

In the right panel of Figure 1, we show the angle and the frequency along the principal axis of the simulated stream. The centre of the cluster is located at the origin, and it is marked as a large grey dot. We also show the area containing 68 per cent of the cluster star distribution with a grey contour line. The stream stars are shown as small black dots for the leading ($\Delta\bar{\theta}_{1}>0$) and trailing ($\Delta\bar{\theta}_{1}<0$) arms. These stars are approximately stripped from the centre of the cluster with a frequency in the range |\Delta\bar{\varOmega}_{1}|\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\sim}}}}}[0.25\,\text{,}\,0.75] rad Gyr*-1*. Once the gravitational potential of the progenitor is negligible, the stripped stars follow a linear motion with constant frequency along the principal axis of the stream (Eq. 5, P25). Since most of the stars are stripped at regular intervals near the pericentres, the stars appear concentrated in peaks. The centres of the three peaks are marked by three coloured dashed lines. Each line marks the frequency of test particles ejected from the cluster at the time of each pericentre passage, with the frequencies distributed uniformly as a function of the angle.

Each peak can be separated from the others by grouping the stripped stars into three intervals defined by the time boundaries $[-1500,\,-1055,\,-610,\,-200]$ Myr (App. D, P25). Each group forms an internal stream. In the top panel of Figure 2, we show the number of stars in each internal component along the principal axis of the stream. Close to the cluster, between |\Delta\bar{\theta}_{1}|\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\sim}}}}}[0.1\,\text{,}\,0.4] deg, we can observe two peaks in the total star distribution (black histogram), corresponding to the third (blue) and second (green) internal streams. The internal streams expand in time in angle space as the stars move along the principal axis. As a result, the density of these streams decreases and they become indistinguishable from the others at large angles. This is the case of the internal stream generated during the first pericentre passage (red). In the bottom panel of Figure 2 we show a projection in angle space of each component. The vertical separation between the internal streams is arbitrary and has been chosen to show each component separately. In the next section, we present a methodology for modelling the mass loss of a globular cluster that reproduces this internal structure. We then use this model to describe the internal dynamics of the stream.

3 Model of the stream in angle-frequency space

In this section, we develop a model of a stellar stream in angle-action coordinates. The main advantage in comparison to an N​​ -body simulation is that this method generates a stream model much faster, especially when the integration time is long and the cluster is composed of a large number of stars. In addition, it allow us to generate large star samples. Using these samples, we can estimate the density along the stream with negligible statistical error and convolve it with observational uncertainties and selection functions. Large samples are in general impractical to obtain with N​​ -body simulations, since they require repeating the simulation many times with different initial conditions. The main drawback of modelling in angle-action coordinates is that the model depends on several free parameters that have to be estimated from a theory or calibrated from a simulation. Furthermore, computing angle-action coordinates requires static potentials with a high degree of symmetry (§3.5, BT08).

We model the stream in the reference frame aligned with the principal axes of the stream with the globular cluster centre located at the origin ($\Delta\bar{\theta},\Delta\bar{\varOmega}$). In these coordinates, the orbit of a stripped star is derived from Eq. 1 of P25:

[TABLE]

where $j=\{1,2,3\}$, $t$ is the simulation time, $t_{\rm s}$ is the stripping time of the stream star (Eq. 8, P25), and the stripping angle is defined as \Delta\bar{\alpha}\equiv\Delta\bar{\theta}\mbox{\footnotesize(t_{\rm s})}. In this model, we assume that the leading and trailing arms are symmetric and introduce the parameter $\kappa=(1,-1)$, which gives the sign of the integration for the leading and trailing arms respectively. The orbit of a stream star, as described by Eq. 1, is equivalent to a rectilinear uniform motion from an initial position $\Delta\bar{\alpha}$ at a velocity $\Delta\bar{\varOmega}$. For this reason, it is convenient to model the stream in angle-frequency space.

In the following subsections, we model the distribution of stripped stars along the cluster orbit in order to estimate the stripping time of each star. We take the cluster orbit from an initial time $T$ to the present. This initial time $T$ corresponds to the age of the stream and is considered equivalent to the accretion time of the globular cluster. We also model the distribution of stripping angles, which are the initial conditions for the motion, as well as the distribution of frequencies of the stripped stars along the principal axes of the stream. To increase the sample size of the N​​ -body simulation, and therefore to minimise the statistical error when determining the probability distributions, we multiply the angles and frequencies of the stream stars by the parameter $\kappa$. This reverses the sign of the angles and frequencies in the trailing arm, merging the two arms together.

The model of the stellar stream is generated using the following procedure:

Determine the distribution of the number of stripped stars along the cluster orbit from an initial time $t=-T$ to the present time $t=0$. 2. 2.

From this distribution and the total number of stripped stars, generate a random sample of stripping times $t_{\rm s}$. 3. 3.

For each star, randomly generate the initial stripping angle $\Delta\bar{\alpha}$ and frequency $\Delta\bar{\varOmega}$ according to the estimated distributions. 4. 4.

Integrate each star forward in time using Eq. 1.

The actions of the stream stars are estimated from their frequencies, which are then used to compute the inverse transformation from angle-actions to Galactocentric cartesian coordinates. This transformation is necessary to evaluate the Gaia selection function and simulate observational errors. We then validate the generated model by comparing it with the N​​ -body simulation. Finally, we use the model to study the surface density of the stream as a function of $T$.

3.1 Distribution of the number of stripped stars along the cluster orbit

In general, the number of stripped stars depends on the tidal forces acting on the cluster. These forces, in turn, depend on the position of the cluster along its orbit. In angle-action coordinates, the orbit of the cluster is given by Eq. 1 of P25:

[TABLE]

where the values of the frequency $\varOmega_{i}^{\scalebox{0.5}{\>\!GC}}$ and angle at the current time \theta_{i}^{\scalebox{0.5}{\>\!GC}}\mbox{\footnotesize(0)} of the M68 cluster are listed in Eq. 3 of P25. Figure 3 shows the number of stripped stars $N_{\rm s}$ as a function of the radial angle of the globular cluster $\theta_{r}^{\scalebox{0.5}{\>\!GC}}$, within bins of $4\pi/139\simeq 0.09$ rad. The coordinate $\theta_{r}$ is defined such that the pericentres are located at $\theta_{r}=2\pi n$ rad, where $n$ is an integer. Similarly, the apocentres are located at $\theta_{r}=(2n+1)\pi$ rad. In this simulation, the pericentres correspond to $\theta_{r}^{\scalebox{0.5}{\>\!GC}}=\{-4\pi,-2\pi,\,0\}$ rad, and the apocentres to $\theta_{r}^{\scalebox{0.5}{\>\!GC}}=\{-3\pi,-\pi,\,\pi\}$ rad. In Figure 3, the pericentres are marked with vertical red, green and blue dashed lines. The stripping time $t_{\rm s}$ of the stream stars is also indicated on the upper ordinate axis, and the start and end of the simulation are marked with vertical grey dashed lines.

There is a small delay between the pericentres and the position of the peaks of an angle $\theta_{r}^{\scalebox{0.5}{\>\!D}}\simeq 82.5$ mrad, or equivalently $t_{\scalebox{0.5}{\>\!D}}\simeq 6$ Myr. This delay may be because the stars require a period of time to leave the cluster after the tidal shock, or because the stripping time estimation does not account for the influence of the cluster potential on the stripped stars during the stripping process. We introduce the angle $\theta_{r}^{\scalebox{0.5}{\>\!M}}\equiv\theta_{r}^{\scalebox{0.5}{\>\!GC}}+\theta_{r}^{\scalebox{0.5}{\>\!D}}$ to correct the radial angle of the cluster for the delay. When expressed as a function of $\theta_{r}^{\scalebox{0.5}{\>\!M}}$, the peaks coincide with the beginning and end of the angular radial period. The N​​ -body simulation completed two periods corresponding to $\theta_{r}^{\scalebox{0.5}{\>\!M}}\in[-4\pi\,\text{,}\,0)$. Due to the periodicity of this evolution, the stars stripped within these two periods can be grouped, and $N_{\rm s}$ can be expressed as a single period of $\theta_{r}^{\scalebox{0.5}{\>\!M}}\in[0\,\text{,}\,2\pi)$. In the bottom panel of Figure 3, we plot the $N_{\rm s}$ as black dots as a function of $\theta_{r}^{\scalebox{0.5}{\>\!M}}$, where the angle corresponds to the centre of the bin.

We model the evolution of $N_{\rm s}$ using a double exponential function:

[TABLE]

where $\theta_{r}^{\scalebox{0.5}{\>\!M}}\in[0\,\text{,}\,2\pi)$. The first exponential approximates the peak after the tidal shock at the pericentre. Its scale length in radians is modelled by the parameter $\tau_{1}$. The second exponential models the approach to the pericentre, and its scale length is determined by $\tau_{2}$. The parameter $A$ determines the relative amplitude of the peak with respect to the baseline $C$. We plot the best-fitting f\mbox{\footnotesize(\theta_{r}^{\scalebox{0.5}{>!M}})} as a solid red line in Figure 3, where the parameters are estimated following the procedure described in Appendix C.1. The values of the parameters are listed in Table 1(a). In the bottom panel of Figure 3, we indicate the limits of the radial period with vertical dashed grey lines and mark with a vertical grey solid line the angle where the two exponentials are equal:

[TABLE]

We also show in grey the second exponential at the beginning of the period so the asymmetry of the slopes is clearly visible.

The probability of a star being stripped as a function of $\theta_{r}^{\scalebox{0.5}{\>\!M}}$ is obtained by normalising Eq. 3:

[TABLE]

where $N_{\rm c}$ is the normalisation constant. This probability density function (PDF) is periodic, and can therefore be generalised for an orbit of arbitrary length by wrapping the angle $\theta_{r}^{\scalebox{0.5}{\>\!M}}$ within the interval $[0\,\text{,}\,2\pi)$. In Appendix C.2, we provide a formal definition of the normalisation constant for a general orbit. Using Eq. 5, we can generate random samples of stripping angles of the stream stars $\theta_{r}^{\scalebox{0.5}{\>\!S}}$. In Appendix C.3, we introduce the methodology for generating random samples from this distribution. The stripping time of the stream stars $t_{\rm s}$ is obtained by subtracting the delay $\theta_{r}^{\scalebox{0.5}{\>\!D}}$ from $\theta_{r}^{\scalebox{0.5}{\>\!S}}$ and isolating the time from Eq. 2:

[TABLE]

Given an orbit of duration $T$, we define the number of stripped stars per unit of time per arm as $N_{T}$. The total number of stripped stars along the orbit is therefore given by:

[TABLE]

For the N​​ -body simulation of M68, the time between the two complete periods is approximately $913.95$ Myr, and the total number of stripped stars is $4071$. Assuming these stars are equally distributed between the two arms, we fix the parameter $N_{T}\simeq 2.23$ stars Myr*-1* arm*-1* to reproduce the results of the simulation. In Section 6, we determine the value of this parameter from the observational data of M68.

The total mass loss is estimated by multiplying $N_{\scalebox{0.5}{\>\!S}}$ by the mean of the star mass distribution of the cluster. The synthetic population of the cluster follows a Kroupa initial mass function with a mean of approximately $0.289$ M{}_{\scalebox{0.7}{\odot}} (§2.1, P25). For the stream simulated with this model, the estimated mass lost is about $1930.94$ M{}_{\scalebox{0.7}{\odot}} for $T=1.5$ Gyr, which corresponds to \mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\sim}}}}}1.55 per cent of the initial mass of the cluster (Table 1, P25). In this model, we assume that this mass loss is small and neglect variations in mass and density distribution of the cluster along the orbit. We also neglect the variations in the stellar mass function along the stream caused by star mass segregation (Webb & Bovy, 2022). We assume that this simplification is valid for orbits of a similar duration to that of the simulation.

3.2 Stripping points distribution

We determine the distribution of stripping points $\Delta\bar{\alpha}$ by integrating the orbit of the stream stars backwards in time until they return to the globular cluster. The integration time is estimated using the ratio of the relative angle and frequency of each stream star to the globular cluster (Eq. 7 of P25). This procedure is analogous to the Inverse Time Integration (invi) method introduced in §4 of P25. The resulting distribution is roughly flat along the principal axis of the stream and approximately a correlated Gaussian in the perpendicular plane (Fig. 4, P25).

We model this distribution using a multivariate Gaussian distribution. Due to the presence of outliers in the coordinate $\Delta\bar{\alpha}_{1}$ (Left panel Fig. 4, P25), we estimate the mean using the median. The values of the parameters of the best-fitting distribution are listed in Table 1(b). In the left panel of Figure 4, we plot in black the distribution of stripping points along the principal axis of the stream $\Delta\bar{\alpha}_{1}$, estimated with an histogram. We also plot the first component of the best-fitting multivariate Gaussian distribution as a solid red line. Due to the asymmetric distribution of $\Delta\bar{\alpha}_{1}$, the Gaussian model overestimates the positive values and underestimates the peak. However, these differences are negligible because variations in \Delta\bar{\alpha}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\sim}}}}}0.5 mrad are much smaller than the length of an arm of a stream. In the case of the N​​ -body simulation of M68, the length of an arm is \mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\sim}}}}}1 rad (Fig 3c, P25). Nevertheless, in order to provide a complete theoretical study, we present a list of distributions that better fit $\Delta\bar{\alpha}_{1}$ in Appendix D.

In the right panel of Figure 4, we plot the stripping points in the perpendicular plane of the stream ($\Delta\bar{\alpha}_{2},\Delta\bar{\alpha}_{3}$). The $1$, $2$, and $3$-$\sigma$ levels of the best-fitting multivariate Gaussian distribution are shown as solid red lines. The median of the distribution is marked with a large red dot and the centre of the globular cluster is marked with a large grey dot. As can be seen, this model provides a good fit to the distribution of stripping points in this plane.

3.3 Frequency distribution along the stream

In the reference frame aligned with the principal axes of the stream, the frequencies are approximately proportional to the actions, according to the following relation:

[TABLE]

where $\lambda_{i}$ are the eigenvalues of a linear approximation of the Hamiltonian (Eq. 4, P25). Note that this is only valid for a progenitor of size and mass of a globular cluster (Amorisco, 2015). For the globular cluster M68 and the potential model of the Milky Way used in the N​​ -body simulation, the absolute values of the eigenvalues satisfy $\lambda_{1}\gg\lambda_{2},\lambda_{3}$ (Eq. 6, P25). Since the three actions of the stream stars are similar in magnitude (Fig. 3a, P25), it follows that $\Delta\bar{\varOmega}_{1}\gg\Delta\bar{\varOmega}_{2},\Delta\bar{\varOmega}_{3}$. Consequently, the length and the surface density of the stream are determined by $\Delta\bar{\varOmega}_{1}$. We therefore model this frequency in more detail, on the assumption that it is uncorrelated with the frequencies perpendicular to the stream.

In Figure 5, we plot the dependency of the frequency along the principal axis of the stream of the stripped stars $\Delta\bar{\varOmega}_{1}$ on the radial angle of the globular cluster $\theta_{r}^{\scalebox{0.5}{\>\!GC}}$. We observe that, on average, the frequencies of the stars stripped at pericentres are larger than those stripped at apocentres. Similarly, the dispersion of the distribution also changes along the cluster orbit, being maximum at the pericentres and minimum at the apocentres. This is because the tidal forces are stronger at the pericentres, resulting in stars being stripped with greater relative velocity with respect to the cluster. The frequency peaks appear with a delay relative to the pericentres, which are indicated by vertical coloured dashed lines. We therefore express the frequency as a function of the corrected angle $\theta_{r}^{\scalebox{0.5}{\>\!M}}$, similarly as we do for the number of stripped stars (Section 3.1).

In the top left panel of Figure 6, we plot the frequency $\Delta\bar{\varOmega}_{1}$ of the stars stripped during the angular interval $\theta_{r}^{\scalebox{0.5}{\>\!M}}\in[-4\pi\,\text{,}\,0]$. This interval corresponds to the two complete periods present in the simulation. The angle is wrapped within $\theta_{r}^{\scalebox{0.5}{\>\!M}}\in[-\pi\,\text{,}\,\pi)$, with the pericentre at $\theta_{r}^{\scalebox{0.5}{\>\!GC}}=0$ rad. To visualise the variations in the mean and standard deviation (std) of the frequency along the cluster orbit, we divide $\theta_{r}^{\scalebox{0.5}{\>\!M}}$ in 41 bins of approximately $0.153$ rad and compute the mean and standard deviation of $\Delta\bar{\varOmega}_{1}$ of the stars within each bin. We plot the mean and standard deviation as a large red dot with error-bars in the bottom left panel of Figure 6. We also show the frequency values as small grey dots for reference. We observe that the mean and the standard deviation are approximately constant close to the apocentre at $\theta_{r}^{\scalebox{0.5}{\>\!M}}=\{-\pi,\,\pi\}$ rad, and increase at the pericentre. We also note that the slopes of the peak are asymmetric, with the incident slope being steeper than the exiting.

In order to model the evolution of the distribution of $\Delta\bar{\varOmega}_{1}$, we divide the interval $\theta_{r}^{\scalebox{0.5}{\>\!M}}\in[-\pi\,\text{,}\,\pi)$ into 70 bins of $2\pi/70\simeq 0.09$ rad, and compute the mean and standard deviation of the frequency for the stars within each bin. In the top and bottom panels of Figure 7, we plot the results for the mean and standard deviation respectively, where the angle $\theta_{r}^{\scalebox{0.5}{\>\!M}}$ is taken as the centre of the bin. In the case of the mean, the data shows a slightly asymmetric peak with low dispersion along the radial angle. On the other hand, the standard deviation peak is approximately symmetric and the dispersion is larger. We model the evolution of the mean and standard deviation using the double exponential function introduced in Eq. 3. We calculate the best-fitting parameters using the method described in Appendix C.1, and list the results in Table 1(a). The best-fitting models are shown as a solid red lines in Figure 7. We observe that the models are in very good agreement with the data obtained from the N​​ -body simulation.

We model the frequency distribution of $\Delta\bar{\varOmega}_{1}$ of the stripped stars within each bin as a Gaussian distribution, with the mean and standard deviation given by the best-fitting double exponential models. We verify this assumption in Appendix D and propose a more accurate model that depends on additional free parameters. Thus, the frequency distribution of $\Delta\bar{\varOmega}_{1}$ depends on the stripping angle of the stream stars along the cluster orbit. For a given unwrapped stripping angle $\theta_{r}^{\scalebox{0.5}{\>\!S}}$, computed using the method described in Section 3.1, we evaluate the double exponential models for the remainder, which is defined as follows:

[TABLE]

3.4 Frequency distributions perpendicular to the stream

In the top centre and top right panels of Figure 6, we plot the frequencies of the stripped stars along the axes perpendicular to the stream $\Delta\bar{\varOmega}_{2}$ and $\Delta\bar{\varOmega}_{3}$. These plots are analogous to the plot for $\Delta\bar{\varOmega}_{1}$ in the top left panel of the same figure. The prevalence of negative frequencies is a consequence of the choice of a right-handed basis for a clockwise rotating cluster, or equivalently for a cluster of $J_{\phi}\equiv L_{z}<0$. The two frequencies present a peak at $\theta_{r}^{\scalebox{0.5}{\>\!M}}=0$, corresponding to the pericentre. These peaks are less prominent than the peak of $\Delta\bar{\varOmega}_{1}$. Similarly, they are more populated than the rest of the period due to the greater mass loss near the pericentres (Section 3.1).

The perpendicular frequencies are about $50$ times smaller than $\Delta\bar{\varOmega}_{1}$ (Section 3.3), so their impact on the length and surface density of the stream is less significant. Therefore, we neglect the variations in $\theta_{r}^{\scalebox{0.5}{\>\!M}}$ and model the frequencies integrated within the interval $\theta_{r}^{\scalebox{0.5}{\>\!M}}\in[-\pi\,\text{,}\,\pi)$ with Gaussian distributions. The distributions of $\Delta\bar{\varOmega}_{2}$ and $\Delta\bar{\varOmega}_{3}$, estimated with a histogram, are shown as grey lines in the bottom centre and bottom right panels of Figure 6. We also plot the best-fitting univariate Gaussian distributions for each frequency in red and list the values of the best-fitting parameters in Table 1(c). We verify the compatibility of the Gaussian model with the simulated data using a Kolmogorov-Smirnov test. We obtain p ​-values of approximately $0.44$ and $0.12$ for $\Delta\bar{\varOmega}_{2}$ and $\Delta\bar{\varOmega}_{3}$, respectively. This indicates that the hypothesis that the data are generated by Gaussian distributions is not rejected at the 95 per cent confidence level ($\alpha=0.05$).

3.5 Actions and inverse coordinate transformation

In the reference frame aligned with the principal axes of the stream, the actions are proportional to the frequencies, as given by Eq. 8. The proportionality constants are the eigenvalues. These can be evaluated at the position of the globular cluster using the Torus Mapper (TM) algorithm111Galpy (Bovy, 2015) provides a Python front end for this functionality (Binney & McMillan, 2016). Alternatively, the eigenvalues can be estimated from the N​​ -body simulation, as described in Appendix F of P25. The actions are required for the inverse transformation from angle-action coordinates to Galactocentric cartesian coordinates. This transformation is also computed numerically using the Torus Mapper code. With a tolerance $\texttt{tol}=2\times\!10^{-4}$, the precision achieved in sky coordinates is \mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\sim}}}}}0.2 deg from the true value.

In Figure 8, we plot the actions of the stream stars of both arms as a small black dots, and the position of the cluster centre as a large blue dot. The $1$, $2$, and $3\text{-}\sigma$ levels of the distribution obtained from the model of the frequencies are indicated as red lines. The left panels show the actions along the principal axis of the stream $\Delta\bar{J}_{1}$, versus the perpendicular actions $\Delta\bar{J}_{2}$ and $\Delta\bar{J}_{3}$. We note that these variables are correlated. These correlations have been neglected when assuming that the frequency distributions are uncorrelated (Section 3.3). Nevertheless, the model accurately describes the central tendency and variance of the simulated stars. On the other hand, the perpendicular frequencies to the stream are uncorrelated, and the model accurately describes their distribution.

3.6 Comparison with the N​​ -body simulation

We validate the model by comparing it with the N​​ -body simulation of the M68 globular cluster. We simulate $10^{3}$ streams of $T=1.5$ Gyr and select the stars within $\theta_{r}^{\scalebox{0.5}{\>\!GC}}\in[-4\pi\,\text{,}\,7\pi/4]$ rad or equivalently $t_{\rm s}\in[-1353.28\,\text{,}\,-39.62]$ Myr. This approximately corresponds to the range in which the stars are stripped in the N​​ -body simulation (Section 2). We estimate the distribution of angles along the principal axis of the stream $\Delta\bar{\theta}_{1}$ of each sample using a histogram with bins of $6/625$ rad. In the top panel of Figure 9, we plot the mean of the distributions in red and mark the $1$ and $2\text{-}\sigma$ levels with a shaded red area. We also plot the histogram for the stars of the N​​ -body simulation in black for reference. The bottom panel shows the distributions of the internal components of the stream, as defined in Section 2. Here, we mark the mean of the model with coloured lines and the N​​ -body with light coloured lines. The components follow the same colour code as in Figure 2.

We observe a slight overestimation of the first peak, which corresponds to the last pericentre passage (blue). We note that we have not used data within $\theta_{r}^{\scalebox{0.5}{\>\!GC}}\in[0\,\text{,}\,2\pi)$ rad to calibrate the model, since the last period is incomplete. As shown in Figure 5, the stars stripped during this period have a slightly lower mean and a slightly larger standard deviation of $\Delta\bar{\varOmega}_{1}$ than the stripped stars in previous periods. We also observe a slight underestimation of the second peak, which corresponds to the second pericentre passage (green) plus the slowest stars from the first pericentre passage (red). The N​​ -body simulation is within the $2\text{-}\sigma$ level of the model, so this discrepancy could be explained by a statistical bias in the N​​ -body simulation sample. Alternatively, a slight overestimation of the standard deviation of $\Delta\bar{\varOmega}_{1}$ makes the stars expand faster and therefore appear less compact near the peak. We consider these differences to be negligible and conclude that the model is in very good agreement with the simulation.

3.7 Dynamical evolution

This model allow us to describe the evolution of a stream from an arbitrary accretion time $T$. In Figure 10, we show an stream of nine complete angular periods, with limits $\theta_{r}^{\scalebox{0.5}{\>\!GC}}=-(2n+1)\pi$ rad, where $n$ is an integer such that $-1\leqslant n\leqslant 8$. These periods are indicated by a colour gradient, with red representing the first period and dark purple representing the last. This angular length corresponds to an accretion time $T\simeq 4323.29$ Myr. In the top panel, we show the distribution of stripped stars along the principal axis of the stream $\Delta\bar{\theta}_{1}$ up to the period indicated by the colour of the line. We also include the stars stripped during \theta_{r}^{\scalebox{0.5}{\>\!GC}}\in[\pi\,\text{,}\,\theta_{r}^{\>\!\scalebox{0.61}{0}}) rad in black, located at \Delta\bar{\theta}_{1}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\sim}}}}}[0\,\text{,}\,0.15] rad. In the bottom panel we show each period separately. In both panels, the dashed vertical lines indicate the angle corresponding to the maximum value of the distribution of each period.

In general, the stream stars are stripped in bursts during the tidal shocks at the pericentres (Section 2 and 3.1). These stars form peaks in angular space $\Delta\bar{\theta}_{1}$ that expand along the principal axis of the stream. The stream stars stripped close to the pericentre have large frequencies $\Delta\bar{\varOmega}_{1}$ with high dispersion (Section 3.3), so they move and expand quickly. On the other hand, the stars stripped closer to the apocentre have lower frequencies with lower dispersion and remain close together for a longer period of time. In the top panel of Figure 10, we can see the peak corresponding to the last pericentre passage (dark purple). In this case, the faster stars have not had time to expand, so these stars together with the slowest form a clear overdensity. The section between the peak and the cluster, which is located at the origin of coordinates, has a density of about one third that of the peak. This region is therefore visible as an underdensity.

As can be seen in the bottom panel, the stars stripped one period earlier (dark blue) show a much clearer separation between the fast and slow stars. In this case, the slower stars form an overdensity that includes the slowest stars from the previous period (light blue). Once the peaks have expanded significantly (green and orange periods), the slow and fast stars are no longer clearly separable. The peaks overlap in such a way that the density saturates, resulting in a constant value along the stream. The grey horizontal line in the top panel shows the value at which the density saturates. For an older stream, additional periods will increase the length of the constant density region but not its saturation value. The total length of the stream is determined by the fastest stars stripped during the first pericentre passage (red). In this case, each arm of the stream is about $4.62$ rad long.

In general, streams generated by globular clusters following an eccentric orbit similar to M68 will present symmetric arms with the following density structure:

Large underdensity close to the cluster. 2. 2.

A concentrated overdensity corresponding to the last pericentre passage. 3. 3.

Density oscillations of smaller amplitude, depending on the frequency distribution of the stripped stars. 4. 4.

Constant density of length proportional to the number of pericentre passages undergone by the progenitor. 5. 5.

A long density decline of a length determined by the fastest stars stripped during the first pericentre passage.

In practice, we only observe sections of streams, particularly when they are close to the Sun. Therefore, the section of a stream that is observed determines the research methodology and the possible conclusions that can be obtained. In the following section, we use Gaia data to determine the visible part of the M68 stream and its surface density.

4 GDR3 M68 stream star selection

In this section, we improve the M68 star selection by Palau & Miralda-Escudé (2019) using data from the Gaia Data Release 3 (GDR3) catalogue (Gaia Collaboration et al., 2023). We apply the following cuts to reduce the amount of data to be analysed and remove low-quality measurements:

parallax $<1/0.3$ mas 2. 2.

$-20\leqslant$ dec $\leqslant 80$ deg 3. 3.

$170\leqslant$ ra $\leqslant 302$ deg 4. 4.

b $>15$ deg 5. 5.

bp_rp $\neq$ Null 6. 6.

ruwe $<1.2$ 7. 7.

visibility_periods_used $\geqslant 10$ 8. 8.

duplicate_source $=$ False

For each star that passes these cuts, we determine its intersection with a phase-space volume surrounding the orbit of the cluster. This volume is defined by a bundle of orbits computed using different phase-space coordinates of the cluster and different potentials. This method is described in detail in Appendix C of Palau & Miralda-Escudé (2019), and Appendix E provides the values of the free parameters used in this case. We pre-select the $6614$ stars with the highest intersection and plot their sky coordinates as black dots in the top panel of Figure 11. We observe that the stars are distributed along a path on the sky pointing towards the cluster, which is marked by a large green dot. The star density is significantly higher in the regions where the cluster orbit goes beyond \mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\sim}}}}}5 kpc from the Sun (Fig. 9, P25 and Section 5) and when its projection is close to the disc. We mark the limits of the disc with two dashed grey lines at Galactocentric latitude of $b=\pm 15$ deg.

To determine the absolute magnitude of the stars, we estimate their distance from the closest point of the orbit of the cluster, as described in §5.3 of P25. We also apply the dust reddening correction introduced in App. B of P25. In the top panel of Figure 12, we show the Hertzsprung–Russell or Colour Magnitude Diagram (CMD) of the pre-selected stars. We also include in red a synthetic population of M68 generated by PARSEC/COLIBRI (§2.1, P25). We observe that most of the stars are distributed in two clumps, both of which have an extension for brighter stars. The smaller clump, located at 0.2\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\lesssim}}}}}{G_{\rm BP}\!-\!G_{\rm RP}}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\lesssim}}}}}1.2 mag, has an extension at M_{\rm G}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\sim}}}}}4 mag that fits with the main sequence of the M68 globular cluster. Most of these stars belong of the stream. In order to select the stars that are compatible with the CMD of M68, we define three polygons that correspond to the main sequence (red), the red giant branch (green), and the horizontal branch (blue). We have excluded a section of the turn off because it is highly contaminated by stars from the extension of the second clump. The polygons are shown in the bottom panel of Figure 12. The stars located within one of these polygons are selected.

In the middle panel of Figure 11, we plot as small dots the $601$ pre-selected stars compatible with the CMD of M68. The stream clearly emerges as an elongated overdensity compared to the background. We eliminate the $58$ isolated stars from the selection, defined as the stars with fewer than $3$ neighbours within an angular radius of $2.5$ deg. These stars are shown in blue in Figure 11. The sections of the stream that are close to the cluster and close to the disc are highly contaminated by foreground stars. We eliminate these sections by defining the limits of the observable part of the stream:

$\delta\,\geqslant-3.5$ deg 2. 2.

$\alpha\leqslant 290\,$ for $\,\delta\geqslant 64$ deg 3. 3.

$\alpha\leqslant 270\,$ for $\,\delta<64$ deg.

These limits are marked as blue dashed lines and define the boundaries of the blue shaded areas in the bottom panel of Figure 11. The sections of the stream beneath these areas are difficult to observe using data from the GDR3 catalogue alone. Conversely, the observable section of the stream is projected onto the white area of the sky. The main peak of the stream, which corresponds to the last pericentre passage, lies below the $\delta=-3.5$ deg limit. Its location and the possibility of observing it are discussed in Appendix F. After removing the isolated stars and those within the blue areas, we obtain $291$ stream stars in the final GDR3 selection. These stars are plotted as small coloured dots in the bottom panel of Figure 11, where black indicates that the star is consistent with being originated from the globular cluster, and red indicates otherwise (Section 4.1).

The GDR3 catalogue provides radial velocities for $2$ stars in the final selection. We checked the Survey of Surveys (SoS) DR1 compilation (Tsantaki et al., 2022), which includes radial velocities from several sources, and found $2$ star from the LAMOST DR5 catalogue (Zhao et al., 2012) and $6$ stars from the SEGUE catalogue (Yanny et al., 2009). Additionally, the STREAMFINDER catalogue (Ibata et al., 2021) includes $11$ stars corresponding to the Fjörm stream222The M68 stream is named Fjörm and labelled 22 in the catalogue.. However, we decided not to include radial velocities in the analysis, as all the sources together only provide radial velocities for about $7$ per cent of the final selection. Furthermore, the DESI Milky Way Survey (MWS) (Cooper et al., 2023) is expected to provide spectroscopic measurements for \mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\sim}}}}}100 stream stars in future releases of the DESI catalogue. This will significantly improve the separation between stream and foreground stars, as well as enabling a more detailed spectroscopic analysis of the M68 stream.

4.1 Stream width

In order to determine the width of the stream, we rotate the sky coordinates to minimise the distortion caused by the projection. In Appendix A we introduce the rotation matrix from $\{\delta,\alpha\}$ to the coordinates $\{\phi_{1},\phi_{2}\}$, which minimise the curvature of the stream. In the top panel of Figure 13, we plot the rotated sky coordinates of the GDR3 star selection. The width of the stream is approximately $6$ deg and remains constant along the observable section of the stream. Only the part farthest from the globular cluster (\phi_{1}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\gtrsim}}}}}80 deg) has a width smaller than $3$ deg, most likely because this section is the most distant from the Sun and only the brightest stars are visible.

In the bottom panel of Figure 13, we compare the width of the stream with the N​​ -body simulation. The simulated stars are plotted as grey dots and the GDR3 star selection is plotted as coloured dots. We observe that the simulation is approximately contained within the limits of -11\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\lesssim}}}}}\phi_{2}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\lesssim}}}}}-8 deg, while the observed stream is about twice as wide spanning about -12\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\lesssim}}}}}\phi_{2}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\lesssim}}}}}-6 deg. Additionally, the simulated leading arm of the stream extends beneath the cluster orbit, which is shown as a solid green line in the bottom panel of Figure 13. Approximately half of the width of the stream is located above the orbit. We note that the orbit depends on the Milky Way potential and that its exact location is especially sensitive for \phi_{1}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\gtrsim}}}}}80 deg. However, variations in the potential do not significantly affect the position of the orbit close to the cluster, for \phi_{1}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\sim}}}}}0 deg.

We identify the stars that can be traced back to the globular cluster, and define them as the main component of the stream. The remaining stars are defined as members of the envelope. The stars in the main component are those that fulfil one of the following conditions:

$0\leqslant\phi_{1}\leqslant 85$ deg and $-11.5\leqslant\phi_{2}\leqslant-8.2$ deg 2. 2.

$\phi_{1}>85$ deg

We do not apply the constraint on $\phi_{2}$ in the second condition because, in general, a long stream presents curvature in the rotated coordinate frame, as can be seen in the N​​ -body simulation (grey dots) in the bottom panel of Figure 13. In the same panel, the $199$ stars of the GDR3 selection that fulfil the above conditions are shown as black dots, and the $92$ stars of the envelope are shown as red dots. The sky coordinates of these stars in the ICRS reference frame are also shown in the bottom panel of Figure 11 with the same colour code. Most of the envelope stars are located above the stream ($\phi_{2}>-8.2$ deg) and are distributed uniformly along it. A small fraction is located below the stream, close to the cluster, but mainly within the range 60\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\lesssim}}}}}\phi_{1}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\lesssim}}}}}90.

In Figure 14, we plot the proper motions of the GDR3 star selection. The stars of the main component (black dots) fit with the N​​ -body simulation (grey dots). The larger dispersion is explained by the observational uncertainties. In the top right corner, we illustrate the size of the uncertainties using a star with a Gaia magnitude of $G\simeq 19.17$ mag, which is equal to the mean magnitude of the selected stars. There is no a clear distinction between the main component and the envelope (red dots). Nevertheless, the envelope has systematically larger proper motions in the range -2\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\lesssim}}}}}\mu_{\alpha\ast}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\lesssim}}}}}2 mas yr*-1* than the main component of the stream. This suggests that the stars in the envelope may be slightly closer to the Sun than the main component.

Other streams present a similar structure consisting of a main component and a fainter envelope. For example Jhelum (Woudenberg et al., 2023; Awad et al., 2024) or GD-1 (Malhan et al., 2019; Valluri et al., 2025; Tavangar & Price-Whelan, 2025). In the case of M68, the size of the globular cluster is well constrained (§2.2, P25), and variations within the observational uncertainties do not significantly increase the width of the simulated stream. Therefore, we present a list of possible causes that could explain this discrepancy in width, as well as similar features of other streams:

Prior to accretion, the progenitor cluster orbited a Milky Way satellite galaxy, generating an stellar stream due to the tidal forces of the satellite. After accretion, this stream is mixed with a new one generated by the cluster due to the tidal forces of its new host galaxy (Carlberg, 2020; Malhan et al., 2021; Qian et al., 2022). 2. 2.

A cluster close to a chaotic phase-space region or resonant trapped orbit family can develop a stream with a wide range of star frequencies. This causes the width of the stream to expand, forming a fan-shaped feature. (Yavetz et al., 2021, 2023). 3. 3.

Perturbations from a tilting disc can produce widening of stellar streams (Nibauer et al., 2024). 4. 4.

Interactions with Milky Way satellites, such as the Sagittarius dwarf galaxy, the LMC (e.g. El-Falou & Webb, 2022; Dillamore et al., 2022) or the Galactic bar (e.g. Price-Whelan et al., 2016; Thomas et al., 2023) could alter the structure of the stream and explain the observed widening. 5. 5.

According to simulations, interacting with hypothetical dark matter subhaloes can increase the width of the stream (e.g. Carlberg & Agler, 2023; Carlberg et al., 2024).

Another explanation is that the M68 stream comprises several star populations of different origins. In fact, the metallicity of the stream is broader than that of similar globular cluster streams, such as Palomar 5 (Martin et al., 2022). Therefore, it is possible that some stream stars, especially those in the envelope, are remnants of an accreted system that included M68 and was formed from a dwarf galaxy and other globular clusters (Section 6.1). This question remains open and will require more precise and complete spectroscopic surveys in future studies.

4.2 Stream surface density

As the observed width of the stream cannot be explained by an axisymmetric potential, such as that used in the N​​ -body simulation, and given the possibility that the stream consists of multiple components, we restrict our analysis to the stars from the main component of the stream. These stars are shown as black dots in the bottom panels of Figure 11, 13, and 14. To eliminate projection distortions, especially prevalent for long streams, we estimate the angle along the principal axis of the stream relative to the globular cluster $\Delta\bar{\theta}_{1}$ for each stream star. This angle is estimated from the closest point of the cluster orbit, using the same method as that used to estimate the Heliocentric distances (§5.3, P25). There is a systematic offset of $5.47$ mrad between the exact and the estimated values of $\Delta\bar{\theta}_{1}$ in the observable section of the stream. After the correction, the standard deviation of the error is $5.12$ mrad, which we consider to be negligible.

In the top panel of Figure 15, we plot the observed surface density of the main component of the stream, estimated using a histogram with bins of approximately $10$ mrad. The limits of the observable section $\Delta\bar{\theta}_{1}\simeq[0.284\,\text{,}\,0.967]$ rad are shown as blue dashed vertical lines. The surface density presents a double-peak structure close to the cluster. It is the highest at \Delta\bar{\theta}_{1}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\sim}}}}}0.3 rad, and peaks at \Delta\bar{\theta}_{1}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\sim}}}}}0.45 rad, which corresponds to the point of closest approach to the Sun. For larger angles, the observed surface density declines as the stream moves away from the Sun.

5 Mock star selection

In this section, we introduce the process of generating mock final selections of stream stars. These mocks are used to determine the PDF of the surface density of the stream, taking into account the effect of the selection function, observational uncertainties and the process of separating the stream stars from the foreground. Additionally, we use a mock random sample of observed stream stars to test the methods for determining the age and mass loss of the progenitor cluster. This process consists of ten steps:

1. (1)

Simulation of the stellar stream of age $T$ in angle-frequency space (Section 3). 2. (2)

Conversion from the angle-action coordinates to Galactocentric cartesian coordinates using the Torus Mapper code (Subsection 3.5). 3. (3)

Estimation of the angle $\Delta\bar{\theta}_{1}$ from the cluster orbit (Subsection 4.2). 4. (4)

Generation of a synthetic stellar population compatible with the globular cluster (§2.1, P25). 5. (5)

Application of the dust reddening extinction (App. B, P25). 6. (6)

Application of the GDR3 selection function (§5.1, P25). 7. (7)

Estimation of the GDR3 astrometric, photometric, and spectroscopic observational uncertainties (§5.2, P25). 8. (8)

Generation of random samples following Gaussian distributions given by the simulated values and the observational uncertainties (§5.4, P25). 9. (9)

Selection process (Section 4). 10. (10)

Determination of the stars of the main component of the stream (Subsection 4.1).

The following summary outlines the process of selecting stream stars from the GDR3 catalogue (9) introduced in Section 4:

1. (9.1)

Application of the cuts to reduce the total amount of data. The quality cuts in ruwe, visibility_periods_used, and duplicate_source are not included in this simulation. 2. (9.2)

Pre-selection of stars compatible with a phase-space volume following the orbit of the cluster. 3. (9.3)

Estimation of the distance of the stars from the cluster orbit. 4. (9.4)

Selection of stars compatible with the CMD of the progenitor cluster. 5. (9.5)

Elimination of the sections of the stream obscured by foreground contamination.

Figure 16 illustrates several steps of this process. In the top panel, we plot in black the surface density of both arms of a simulated stream of age $T=3039.7$ Myr (1). This stream corresponds to the configuration that best fits the GDR3 star selection, as discussed in Section 6. The surface density is estimated from a random sample of approximately $1.3\times\!10^{7}$ stars using a histogram of $6.25$ mrad bins, and is expressed as a function of the angle along the principal axis of the stream $\Delta\bar{\theta}_{1}$. The position of the globular cluster is indicated by a green dashed line and the boundaries of the observable section are marked by blue dashed lines. This stream has symmetric arms and four well-defined peaks, with the saturation density covering approximately two-thirds of the observable section. For reference, the surface density of a simulated stream that has reached saturation for the entire depicted angular range is shown as a grey histogram.

In the middle-top panel of Figure 16, we plot the stars that pass the GDR3 selection function (6). This function depends on the sky coordinates and the $G$ magnitude of the stars. The determining factor is the distance. In the bottom panel of Figure 16, we plot the Heliocentric distance of the stream stars as a function of $\Delta\bar{\theta}_{1}$, and mark the cluster position at $r_{\rm h}=10.404$ kpc with a large green dot. The entire trailing arm of the stream lies beyond $10$ kpc, so only the brightest stars are observable. Conversely, the leading arm is closer to the Sun, at up to \mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\sim}}}}}5 kpc. The surface density appears as a wide peak within the interval \Delta\bar{\theta}_{1}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\sim}}}}}[0.4\,\text{,}\,0.65] rad. For \Delta\bar{\theta}_{1}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\sim}}}}}0.35,\ {\rm and}\ 0.5 rad, the underdensities between the real peaks are visible. For values \Delta\bar{\theta}_{1}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\gtrsim}}}}}0.65 rad, the density declines due to the increasing distance of the stream stars.

In the middle-bottom panel of Figure 16, we plot the stars that pass the selection process (9) and belong to the main component of the stream (10). Approximately $2.6\times\!10^{5}$ stars pass through these cuts. This number is large enough to account for the full distribution of observational uncertainties of the stream stars, and therefore ensures that the surface density has converged to its true value. The mock observed distribution peaks at \Delta\bar{\theta}_{1}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\gtrsim}}}}}0.4 rad. This peak does not correspond to the angle of closest approach of the stream, which is $\Delta\bar{\theta}_{1}\simeq 0.56$ rad. It is possible that some stars around \Delta\bar{\theta}_{1}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\sim}}}}}0.56 rad are lost during the selection process. In particular, stars with large proper motions uncertainties have a smaller intersection with the volume following the cluster orbit. For angles \Delta\bar{\theta}_{1}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\gtrsim}}}}}0.55 rad, the density declines almost linearly up to the limit of the observable section.

6 Age and mass loss estimation of the M68 stellar stream

We estimate the age $T$ of the M68 stellar stream, or equivalently the accretion time of the cluster, using the posterior distribution p\mbox{\footnotesize(T|d)} given by the Bayes’ theorem:

[TABLE]

where $d$ is the GDR3 selection of observed stars obtained in Section 4 and plotted in the top panel of Figure 15. The likelihood function \mathcal{L}\mbox{\footnotesize(d>!|T)} is defined as:

[TABLE]

where the product is over all elements of the dataset $d$ and p_{\scalebox{0.5}{\>\!SD}}\mbox{\footnotesize(\Delta\bar{\theta}_{1}|T)} is the PDF of the surface density within the observable limits. This distribution is estimated by combining $10^{3}$ samples of mock observed stream stars, generated following the procedure described in Section 5. In total, the samples contain a number of stars sufficiently large to cover the distributions of observational uncertainties completely, and therefore ensure that $p_{\scalebox{0.5}{\>\!SD}}$ has converged to the true probability distribution. We estimate the PDF from the samples using a histogram with $150$ bins within the observed range $\Delta\bar{\theta}_{1}\in[0.284\,\text{,}\,0.967]$ rad, which corresponds to bins of approximately $4.55$ mrad. We verified that this number of bins yields the optimal fit by taking a quarter of the star sample to generate the histogram and using the remaining three quarters to determine the best fit.

The prior distribution p\mbox{\footnotesize(T)} is taken as a uniform distribution defined from the present time $T=0$ to $T=12$ Gyr, which is the mean of the estimated age of the cluster (Table 1, P25). The best-fitting configuration $\hat{T}$ is obtained by maximising the posterior distribution:

[TABLE]

This methodology is tested in Appendix G using the N​​ -body simulation. In this case, we recover the age of the stream $T$ with an error of approximately $2$ Myr, and estimate the measurement uncertainty to be approximately $10$ per cent of the best-fitting value.

Figure 17 shows the posterior distribution as a function of the stream age $T$. The radial angle of the globular cluster $\theta_{r}^{\scalebox{0.5}{\>\!GC}}$ is also indicated on the upper ordinate axis. Each pericentre passage is marked with a blue dashed line. An interval of $\pm 70$ Myr centred at the pericentres is shaded in blue to represent the approximate size of the peak. In the top panel, we plot the posterior distribution using a logarithmic scale. Moving from younger to older ages, we can see how the posterior distribution increases with each pericentre passage, reaching a maximum at $\hat{T}=3039.7$ Myr. This value corresponds to the best-fitting configuration, and is marked in the figure with a red dashed line. This behaviour is explained by the fact that the stream grows in bursts caused by tidal shocks at each pericentre passage, rather than linearly over time. For $T<\hat{T}$, the value of the posterior distribution decreases between the peaks because stars released after the pericentre passage move slower and only populate the central part of the observable section of the stream. Consequently, the central density of the stream increases, but its length remains unchanged. This worsens the fit because the observed stream is longer than the simulated one.

We also show the posterior distribution using a linear scale in the bottom panel of Figure 17. After the maximum $\hat{T}$, the posterior distribution decreases for the subsequent two pericentre passages. These two additional tidal shocks populate the region \Delta\bar{\theta}_{1}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\gtrsim}}}}}0.8 rad, which has low surface density in the observed data, as shown in the top panel of Figure 15. For T\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\gtrsim}}}}}4 Gyr, the posterior distribution remains constant. This is explained because the stream reaches saturation within the observable interval. For older streams, the density within the observable section does not change, and therefore the posterior distribution has a constant value.

We plot the PDF of the surface density for the best-fitting configuration p_{\scalebox{0.5}{\>\!SD}}\mbox{\footnotesize(\Delta\bar{\theta}_{1}|\hat{T})} in red in the bottom panel of Figure 15, and compare it to the distribution of GDR3 selected stream stars, shown in grey. This distribution accurately describes several features of the data, including the presence of two peaks and the small underdensity between them, the length of the central peak, and the linear decay at larger angles. The main discrepancy is the excess of observed stars in the range \Delta\bar{\theta}_{1}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\sim}}}}}[0.4\,\text{,}\,0.5] rad. This can be explained by a number of envelope stars being projected within the limits of the main component, by the inclusion of foreground stars, or by statistical biases present in the observed sample. Accurate modelling of these density fluctuations is necessary in order to estimate the age of the stream correctly. In Appendix J, we demonstrate that simpler models assuming constant mass loss and a constant frequency distribution for the stripped stars cannot reproduce these density variations, and consequently bias the estimated age of the stream. Nevertheless, we consider the discrepancy between our model and the observed sample to be small, and we conclude that the best-fitting configuration is a good representation of the data. The PDF of the best-fitting configuration is also shown in Figure 16 as a black histogram.

We estimate the age of the simulated stream to be $T=$3.0397_{-0.287,049,79}^{+5.632,032,24} Gyr, where the uncertainties are calculated using the methodology introduced in Appendix H. The lower uncertainty is small because the stream cannot be shorter than the observed GDR3 star sample. Conversely, the upper uncertainty is large. This is due to the small difference between a stream of $T=3039.7$ Myr and a stream that has reached saturation. This difference can be seen in the middle-bottom panel of Figure 16 for angles \Delta\bar{\theta}_{1}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\gtrsim}}}}}0.8 rad, where the saturated stream is shown as a grey histogram. In this region, it is difficult to determine the surface density with high precision in practice. This is because the stream is located at a distance of \mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\sim}}}}}8 kpc from the Sun. At this distance, the selection procedure can erroneously exclude stars with large proper motion errors. Furthermore, once it has reached saturation within the observable section, the surface density is insensitive to the age of the stream. These two phenomena explain the large upper uncertainty of our estimate.

In order to determine the mass loss of the cluster, we adjust the parameter $N_{T}$ (Section 3.1) so that the mean of the mock observed star samples (Section 5) equals the number of stars in the observed sample. We note that the model used to generate the mock star samples assumes that the cluster mass and density distribution remain constant throughout its evolution (Section 3.1). Therefore, we determine $N_{T}$ assuming that it is independent of the properties of the cluster. We tested this methodology in Appendix G using the N​​ -body simulation, obtaining uncertainties of about $10$ per cent with respect to the true value. In the GDR3 star selection, there are $199$ stars in the main component of the stream (Section 4). Adjusting for this number, we obtain $N_{T}\simeq 1.718\pm 0.104$ stars Myr*-1* arm*-1*. Multiplying this by the mean of the Kroupa distribution (Section 3.1), we obtain a mass loss of approximately $0.496\pm 0.030$ M{}_{\scalebox{0.7}{\odot}} Myr*-1* arm*-1*. This implies a total mass lost of $3016.5$ M{}_{\scalebox{0.7}{\odot}} since accretion, corresponding to approximately $2.4$ per cent of the cluster initial mass. These values are sufficiently small to be compatible with the model assumption that the mass of the cluster is constant.

6.1 Comparison to previous studies

The age and mass loss of M68 have been studied by Chen et al. (2025), using the Fjörm stream333The M68 stream is also named Fjörm. from the STREAMFINDER catalogue (Ibata et al., 2024). They estimate $T=2.34\pm 0.82$ Gyr, $N_{T}=1.7\pm 0.5$ stars Myr*-1* arm*-1*, and mass loss equal to $0.5\pm 0.13$ M{}_{\scalebox{0.7}{\odot}} Myr*-1* arm*-1*. These results are consistent with ours, with the estimated age being at the lower end of our uncertainty range. In Appendix I, we discuss whether this age discrepancy could be explained by excluding the stars in the envelope of the stream. We find that including all stars in the selection reduces the estimated age of the stream, which partially explains the discrepancy in the results.

Our estimate of the age indicates that the globular cluster M68 was accreted about T\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\gtrapprox}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\gtrapprox}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\gtrapprox}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\gtrapprox}}}}}3 Gyr ago. This is broadly consistent with several major accretion events. Some studies have associated M68 with the Canis Major dwarf (Kruijssen et al., 2019), Helmi streams (e.g. Massari et al., 2019; Kruijssen et al., 2020; Callingham et al., 2022), Sagittarius dwarf Bonaca et al. (2021), and Cetus Polar stream (Malhan et al., 2022). None of these associations are conclusive. However, according to our best-fitting result, it is possible that M68 was accreted more recently. As suggested by Sun et al. (2023), it is possible that M68 is associated with a minor accretion event. Furthermore, integrating the cluster orbit backwards in time over a period $\hat{T}$ reveals that the cluster was accreted at a radius $R\simeq 27.6$ kpc, which is close to the apocentre $R_{\rm apo}\simeq 31.8$ kpc. This is consistent with the cluster having being accreted. However, the potential has not been optimised for this particular stream. Therefore, we postpone drawing any conclusions until a more realistic potential is available.

7 Conclusion

The dynamics of a stream is greatly simplified in angle-action coordinates. This allows us to simulate a stream with a high degree of accuracy. The developed model focuses on reproducing the effects of variable tidal forces acting on the cluster as it moves along an eccentric orbit. We find that the mass loss of the cluster, as well as the mean and standard deviation of the frequencies of the stripped stars, can be modelled using a simple double exponential model. This model has four free parameters for each variable and depends on the cluster radial angle. Due to its simplicity, the double exponential opens up the possibility of explaining the perturbations observed in the N​​ -body simulation using a semi-analytical approach. Once the free parameters are calibrated, the model can reproduce the internal structure of a stream obtained with the N​​ -body simulation. This model is an improvement on previous methodologies for simulating streams in angle-action coordinates, which neglected the time dependency of mass loss and stripping frequencies (Sanders, 2014; Bovy, 2014).

This model incorporates several simplifications. The most significant of these, especially for long streams, is the assumption that the cluster mass and its density distribution are constant. Modelling this aspect is crucial for old streams, such as Palomar 5, and for predicting the dissolution time of the progenitor clusters. Furthermore, the mass lost by ejection and evaporation have been neglected. Although this is expected to be less than the mass lost by tidal shocks, it is not negligible throughout the life of a globular cluster. Additionally, the frequencies perpendicular to the principal axis of the stream of the stripped stars are considered to be independent of the radial angle of the cluster. Their correlations with the frequency along the stream have also been neglected. Eliminating these assumptions could marginally improve the accuracy of the simulated surface density.

This model depends on several free parameters. The total number can be reduced by making additional assumptions. For example, the distribution of stripping angles of the stream stars has a minimal impact on the morphology of the stream and can therefore be simplified. In addition, the peak in frequency dispersion of the stripped stars can assumed to be symmetric. In principle, all of these free parameters can be expressed as a function of the properties of the cluster, such as its mass and density distribution, as well as the potential of the Galaxy. Further developments of this model will include such relations.

This model can be used to study the surface density of a stellar stream along its principal axis and its dynamical evolution. The main difficulty in applying this method to study real streams lies in calculating the transformation between cartesian and angle-action coordinates. This is only possible for static potentials, and it is computationally expensive, which limits its range of applicability. The transformation is also less accurate for orbits within strong resonant trapped families or close to these resonances. This model can be used to study stellar streams with a known progenitor similar to globular clusters. Examples include the streams of Palomar 5, NGC 3201, M2, and M5. The first application is the study of the stream generated by the globular cluster M68.

Using data from the GDR3 catalogue, we improved upon previous star selections of the M68 stream, obtaining $291$ stars that are likely to be members of the stream. The width of the stream is too large to be explained only by stars stripped from the cluster. We present several explanations for this phenomenon, including the effects of complex and evolving potentials, or the presence of remnants from an accretion event. This feature can potentially introduce biases in the constraints on the Milky Way potential, particularly for methods that use a simulation of the stream based on the properties of the cluster, or based on integrating stars backwards to bring them back to their progenitor. In addition, the observed surface density of the stream presents several density variations in the region closest to the cluster. These variations may be due to the small size of the star sample. In general, they cannot be reproduced by a model that assumes a constant mass loss and a constant frequency distribution for the stripped stars. If these density variations are real, modelling the internal structure of the M68 stream is necessary to correctly estimate its age and the mass loss of the cluster.

We use an static and axisymmetric potential model for the Milky Way. This potential has not been optimised to fit the M68 stream. Nevertheless, the discrepancies between the simulation and the observations are of about one degree in declination within the observable section of the stream. These discrepancies do not affect the estimation of the surface density. By selecting stars consistent with having been stripped from the cluster, we determine the age of the stream, obtaining $T=$3.0397_{-0.287,049,79}^{+5.631,979,31} Gyr. The uncertainties are dominated by the fact that we only observe a section of the stream. As a consequence that the surface density reaches a saturation level within the observable section and therefore the stream only grows in length, we cannot accurately estimate the upper limit of its age. Conversely, we can precisely estimate the lower limit because the stream cannot be shorter than the observed section. The exact value of the best-fitting configuration depends on the details of the observed surface density. In particular, it depends on the section of the stream that is farthest from the globular cluster, as well as on the separation between the main component and the envelope. Improving the observed star selection could potentially reduce the uncertainty in this estimate. Our current best-fitting age of $T=3039.7$ Myr suggests that M68 was accreted more recently than the major accretion events that dominated the formation of the Milky Way.

We also estimate the number of stripped stars $N_{T}\simeq 1.717\pm 0.104$ stars Myr*-1* arm*-1* and the mass loss to be $0.496\pm 0.030$ M{}_{\scalebox{0.7}{\odot}} Myr*-1* arm*-1*. These values are consistent with the modelling assumption that the mass of the cluster is constant. While these results are compatible with previous research, they are approximately $20$ per cent lower than those obtained using the N​​ -body simulation. This suggests an inconsistency between the observed properties of the cluster and those inferred from the stream. Alternatively, this could indicate that a significant number of stars have been lost during the selection process and the determination of the main component of the stream. Further research is required to address this discrepancy.

Acknowledgements

This work is supported by NSFC (12573022, 12273021), National Key R&D Program of China (2023YFA1607800, 2023YFA1607801, 2023YFA1605600, 2023YFA1605601), 111 project (No. B20019), and the science research grants from the China Manned Space Project (No. CMS-CSST-2021-A03). We also acknowledge the Yangyang Development Fund. The computation of this work is done on the Gravity supercomputer at the Department of Astronomy, Shanghai Jiao Tong University.

Data Availability

A presentation summarising the content of this paper is available on Zenodo https://zenodo.org/records/17330448. The Python package invi containing the codes used in this paper is available on GitHub: https://github.com/cgpalau-astro/invi. The following files related to the M68 stellar stream are available on Zenodo https://zenodo.org/records/17020518:

N​​ -body simulation $T=1500$ Myr 2. 2.

N​​ -body simulation $T=3040$ Myr 3. 3.

Mock Gaia-DR3 star catalogue 4. 4.

Gaia-DR3 star catalogue 5. 5.

Additional data for the Gaia-DR3 star catalogue

Software

The following Tools and Python packages are used in this research:

[TABLE]

Appendix A Rotation matrices

A.1 Angle-Action coordinates

The rotation matrix from the angle-action coordinates relative to the globular cluster ($\Delta\theta,\Delta J$) to the angle-action coordinates aligned with the principal axes of the stellar stream ($\Delta\bar{\theta},\Delta\bar{J}$) is estimated from the N​​ -body simulation of the stream as described in Appendix F of P25. The obtained rotation matrix, which is used in this paper, is:

[TABLE]

We verify that this matrix is approximately equal to the rotation matrix obtained by diagonalising a linear approximation of the Hamiltonian function at the position of the globular cluster (§3, P25). We use the Torus Mapper algorithm (Section 3.5) with tolerance $\texttt{tol}=2\times\!10^{-4}$ to numerically compute the Hessian matrix. By inverting the matrix whose columns are the eigenvectors of the Hessian matrix ordered in descending order of absolute eigenvalue magnitude, we obtain the following rotation matrix:

[TABLE]

We consider these two methods to be approximately equivalent, since the error in all values in $R_{\scalebox{0.65}{TM}}$ is less than $20$ per cent of $R$.

A.2 Sky coordinates

The rotation from the ICRS sky coordinates ($\delta,\alpha$) to the sky coordinates aligned with the stream ($\phi_{1},\phi_{2}$) is given by the following expression:

[TABLE]

The rotation matrix $R_{1}$ is expressed in terms of the Euler angles $\alpha_{x},\alpha_{y}\in[-\pi/2\,\text{,}\,\pi/2)$ and $\alpha_{z}\in[0\,\text{,}\,2\pi)$. It is estimated by minimising the following loss function:

[TABLE]

where the summation is defined for all stream stars in the GDR3 selection (Section 4) and $\phi_{2}$ is computed using Eq. 15 fixing the value of $\alpha_{z}=0$. The optimal values obtained for the angles are $(\alpha_{x},\alpha_{y})\simeq(-1.324,-0.316)$ rad. Then, a horizontal rotation of $\alpha_{z}=-3.440$ rad is included to position the closest extreme of the stream to the cluster at the origin of the coordinates. The resulting matrix is:

[TABLE]

Typically, the rotation is defined so that the average $\phi_{2}$ is approximately zero. This is achieved using the following loss function:

[TABLE]

For the sake of completeness, we provide the rotation matrix corresponding to this definition, even though we only use $R_{1}$ in this paper. The optimal values obtained for the angles from Eq. 18 and the horizontal rotation are $(\alpha_{x},\alpha_{y},\alpha_{z})\simeq(-1.144,-0.524,-3.640)$ rad. The resulting matrix is:

[TABLE]

Appendix B Frequency distributions along the stream

Figure 18 shows the frequencies of the simulated stream stars perpendicular to the stream as a function of the angle along the principal axis of the stream $\Delta\bar{\theta}_{1}$. The frequency space $\Delta\bar{\varOmega}_{2}$ (top panel) do not present a significant structure that allows us to identify the three internal components of the stream. Conversely, the internal components generated during the last two pericentre passages can be identified in the frequency space $\Delta\bar{\varOmega}_{3}$ (bottom panel). However, the component generated during the first pericentre passage is much more extended and populates the extremes of each arm. This makes it difficult to separate, as it is merged with the previous internal components.

Appendix C Double exponential function

The double exponential function is defined in Eq. 3 as a function of the radial angle of the progenitor cluster, corrected for the delay between the mass-loss peaks and the pericentre positions over the periodic domain $\theta_{r}^{\scalebox{0.5}{\>\!M}}\in[0\,\text{,}\,2\pi)$ (Section 3.1). This function depends on the parameters $A$ and $C$, as well as the scale lengths $\tau_{1}$ and $\tau_{2}$. Its integral over a general interval such that $0\leqslant\theta_{r}^{\scalebox{0.5}{\>\!Inf}}\leqslant\theta_{r}^{\scalebox{0.5}{\>\!Sup}}<2\pi$ has an analytical solution:

[TABLE]

where $h$ and $g$ are defined as:

[TABLE]

C.1 Parameter estimation

We use the curve_fit function provided by the Python package scipy (Virtanen et al., 2020) to compute the best-fitting parameters for a given two-dimensional data set. This function implements the Trust Region Reflective algorithm or trf method.

C.2 Normalisation constant

For a general cluster orbit, defined from $t=-T$ to $t=0$, the angle \theta_{r}^{\scalebox{0.5}{\>\!M}}\in[\theta_{r}^{\>\!\scalebox{0.6}{-T}}\!\!\>,\>\!\theta_{r}^{\>\!\scalebox{0.61}{0}}], where 0\leqslant\theta_{r}^{\>\!\scalebox{0.61}{0}}<2\pi. For this orbit, the inverse of the normalisation constant in Eq. 5 is:

[TABLE]

When \theta_{r}^{\>\!\scalebox{0.6}{-T}}\geqslant 0, this integral can be solved analytically using Eq. 20. When \theta_{r}^{\>\!\scalebox{0.6}{-T}}<0, the inverse of the normalisation constant can be expressed in terms of integrals within a single angular period. When expressing \theta_{r}^{\>\!\scalebox{0.6}{-T}} in radians and unwrapping by a period of $2\pi$:

[TABLE]

where $N_{\rm p}$ is the number of periods within \theta_{r}^{\scalebox{0.5}{\>\!M}}\in[\theta_{r}^{\>\!\scalebox{0.6}{-T}}\!\!\>,\>\!0] and $\theta_{r}^{\scalebox{0.5}{\>\!R}}$ is the remainder:

[TABLE]

C.3 Random sample

We generate random samples from the probability density function defined in Eq. 5 using an inverse transform sampling algorithm. First, we unwrap the angle $\theta_{r}^{\scalebox{0.5}{\>\!M}}$ by a period of $2\pi$ and evaluate p\mbox{\footnotesize(\theta_{r}^{\scalebox{0.5}{>!M}})} at regular spacing within \theta_{r}^{\scalebox{0.5}{\>\!M}}\in[\theta_{r}^{\>\!\scalebox{0.6}{-T}}\!\!\>,\;\!\theta_{r}^{\>\!\scalebox{0.61}{0}}]. We then compute the cumulative sum of these evaluations and normalise the result by subtracting the minimum value and dividing by the maximum. This ensures that the cumulative sum is contained within the interval $[0\,\text{,}\,1]$. The result obtained is a numerical evaluation of the cumulative distribution function (CDF). We compute the inverse of the CDF using the cubic spline method implemented by the scipy function make_interp_spline. The inverse is obtained when the parameter x is the CDF and y is the angle $\theta_{r}^{\scalebox{0.5}{\>\!M}}$. Finally, we generate a random sample by evaluating the spline with a random sample from a uniform distribution defined on the interval $[0\,\text{,}\,1]$.

Appendix D Improved models

We provide an improved model of the stripping points and frequency distribution along the principal axis of the stream using the left-skewed Gumbel distribution:

[TABLE]

where $z\equiv(x-\mu)/\sigma$. This distribution provides better fits than a Gaussian distribution with the same number of free parameters.

D.1 Angle $\boldsymbol{\Delta\bar{\alpha}_{1}}$

As shown by the grey histogram in the left panel of Figure 19, the distribution of stripping points along the principal axis of the stream $\Delta\bar{\alpha}_{1}$ is skewed towards negative values. We also show the first component of the multivariate Gaussian distribution introduced in Section 3.2 as a red line. The distribution G_{\rm L}\mbox{\footnotesize(x>!|-0.072,>!0.174)}, plot in blue, provides a better fit than the Gaussian distribution, especially for positive values. Using distributions with more free parameters improves the fit to the peak region. For example, the Johnson SU and Generalised Logistic distributions provide better fits with four and three parameters, respectively.

D.2 Frequency $\boldsymbol{\Delta\bar{\varOmega}_{1}}$

The frequency distribution of the stripped stars along the principal axis of the stream $\Delta\bar{\varOmega}_{1}$ is modelled using the mean and standard deviation of the simulated data within bins defined by the radial angle $\theta_{r}^{\scalebox{0.5}{\>\!M}}$ (Section 3.3). We normalise the values of $\Delta\bar{\varOmega}_{1}$ by subtracting the mean and dividing by the standard deviation within each bin as follows:

[TABLE]

We then integrate all the data within the limits $\theta_{r}^{\scalebox{0.5}{\>\!M}}\in[-\pi\,\text{,}\,\pi]$. The resulting distribution, estimated using a histogram, is plotted in grey in the right panel of Figure 19. The best-fitting Gumbel distribution G_{\rm L}\mbox{\footnotesize(x>!|>!0.480,>!0.873)} is plotted in blue, and a Gaussian distribution centred at zero with a standard deviation of one is plotted in red. The Gumbel distribution provides a slightly better fit because the data is left-skewed. Nevertheless, the Gaussian distribution is a good approximation for modelling the surface density of the stream.

Appendix E Pre-selection parameters for M68

For the case of M68, the volume is constructed with $N=510$ Gaussian distributions computed using a bundle of $M=2\times\!10^{4}$ orbits of length $l=60$ Myr. The means of the Gaussian distributions used to compute the bundle of orbits are taken from Table 3 of P25 for the position of the globular cluster, and from Table 2 of P25 for the potential of the Milky Way. The scale factors of the Gaussian distributions are listed in Table 2. Stars with an intersection P_{\mathchoice{\raisebox{0.0pt}{\resizebox{12.70836pt}{3.5pt}{\hbox{\raisebox{0.0pt}{\displaystyle\rm REG}}}}}{\raisebox{0.0pt}{\resizebox{12.70836pt}{3.5pt}{\hbox{\raisebox{0.0pt}{\textstyle\rm REG}}}}}{\raisebox{0.0pt}{\resizebox{11.27542pt}{3.5pt}{\hbox{\raisebox{0.0pt}{\scriptstyle\rm REG}}}}}{\raisebox{0.0pt}{\resizebox{11.2754pt}{3.5pt}{\hbox{\raisebox{0.0pt}{\scriptscriptstyle\rm REG}}}}}}\geqslant 6\times\!10^{-6} ${\rm yr}^{3}\,{\rm deg}^{-2}\,{\rm pc}^{-1}\,{\rm mas}^{-3}$ are pre-selected.

Appendix F ICRS surface density

It is interesting to observe the peaks of the M68 stellar stream because we can use their location to estimate the time since the pericentre passages. This allows us to constrain the potential of the Milky Way by requiring the cluster to return to the pericentres after integrating backwards over the estimated time periods. To determine the location of the peaks and the expected number of stars, we simulate a stream with an age $T=3039.7$ Myr and a mass loss of $N_{T}\simeq 1.718$ stars Myr*-1* arm*-1* using the method introduced in Section 3.

In the left panel of Figure 20, we show the projected surface density in sky coordinates of the sample of $10\,442$ simulated stream stars. The declination has been scaled to divide the surface into equal square bins of area $16$ µsr. We mark the position of the globular cluster M68 as a large red dot at $(\delta,\alpha)=(-26.744,189.867)$ deg. The observable limit at $\delta=-3.5$ deg is marked by a solid grey line, and the area below this limit is shaded in blue. The stream is not visible in the blue area using only GDR3 data because it is obscured by foreground star contamination (Section 4). In the right panel, we show a sample of $645$ stars that pass the Gaia selection function (§5.1, P25).

The coloured dashed lines indicate the positions along the stream of the peaks, following the colour code in Figure 1 and 2. The peaks corresponding to the last pericentre passage (blue dashed lines) are clearly visible in the simulated data (left panel). Their location can be precisely determined because they mark the end of the underdensity between the cluster and the peaks in both arms. Due to projection distortion, the underdensity of each arm appears uneven in length. The exact location of the main peaks depends on the potential of the Milky Way, but generally, they always appear below the observational threshold.

Conversely, the second (green dashed line) and first (red dashed line) peaks of the leading arm lie within the observable area. However, as can be seen in the right panel of Figure 20, the stream has a low surface density after the GDR3 selection function cut is applied. Additionally, the selection process and the presence of foreground contamination, neither of which is depicted in this figure, will complicate the identification of the peaks. Therefore, we conclude that identifying the peaks of the stream from its surface density using GDR3 data alone is problematic. However, it may be possible to identify the peaks if the signal-to-noise ratio of the star selection is improved. This could be achieved, for example, by incorporating radial velocity measurements into the selection methodology.

Appendix G Test fitting methodology

We test the methodology for constraining the age of the stellar stream and the mass loss of the globular cluster using a mock GDR3 star selection obtained from the N​​ -body simulation (§5.4, P25). This sample consists of $116$ stars with magnitudes $G<20$ mag located in the main component of the stream. We compute the posterior distribution defined in Eq. 10 as described in Section 6 and plot the results in Figure 21.

We observe that the posterior distribution peaks between approximately $T\in[1.3\,\text{,}\,1.8]$ Gyr, with a maximum at $\hat{T}=1351.3$ Myr. The value of the posterior distribution is approximately zero for streams with fewer than three pericentre passages. This is because the streams are too short. The maximum value approximately coincides with the first tidal shock experienced by the cluster at the pericentre at $t=-1353.28$ Myr, with a difference of $1.98$ Myr. For older streams, the posterior distribution decreases until the fourth pericentre passage. The relative maximums occur because the slower stars, released before and after the pericentre passages, populate the range \Delta\bar{\theta}_{1}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\sim}}}}}[0.35\,\text{,}\,0.45] rad of the observable section of the stream. This increases the surface density of the stream but not its length, thereby improving the fit. For streams with more than four pericentre passages, the posterior distribution is approximately zero, indicating that the streams are too long.

We estimate the age of the simulated stream to be $T=$1.3513_{-0.008,896,91}^{+0.171,221,85} Gyr using the methodology introduced in Appendix H to determine asymmetric uncertainties. The relative uncertainties are of approximately $3$ and $12$ per cent of the estimated value, which we consider to be small. The number of stripped stars per unit of time per arm $N_{T}$ (Section 3.1) is determined as in Section 6. We obtain $N_{T}\simeq 2.002\pm 0.173$ stars Myr*-1* arm*-1*, which is approximately $10$ per cent smaller than the value of $2.23$ stars Myr*-1* arm*-1* that best reproduces the N​​ -body simulation. We conclude that the discrepancies in determining the age and mass loss of the cluster are small, and that the methodology used is valid.

Appendix H Uncertainty determination

Given a PDF p\mbox{\footnotesize(x)} where $x\in(-\infty\,\text{,}\,\infty)$, and its maximum ${\hat{\mu}}$, we determine the upper limit of the $1\sigma$-level by solving the following equation for $\lambda_{\rm sup}\geqslant{\hat{\mu}}$:

[TABLE]

where $\sigma\simeq 0.683$. Similarly, we determine the lower limit by solving for $\lambda_{\rm inf}\leqslant{\hat{\mu}}$:

[TABLE]

The superior $1\sigma$-uncertainty of the variable $x$ is defined as $\epsilon_{\rm sup}\equiv\lambda_{\rm sup}-{\hat{\mu}}$, and the inferior as $\epsilon_{\rm inf}\equiv{\hat{\mu}}-\lambda_{\rm inf}$. The same procedure applies when $x$ is limited, except that the infinities are replaced by the limits.

Appendix I Age estimate using the entire GDR3 star selection

The best-fitting age of the M68 stream estimated using only the stars in the main component is $\hat{T}=3039.7$ Myr (Section 6).This measurement is compared to the estimate using all the stars in the GDR3 selection, which includes both the main and the envelope components. In this case, we obtain $\hat{T}=2637.4$ Myr. This difference is due to the fact that the stars in the envelope populate the region \Delta\bar{\theta}_{1}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\lesssim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\lesssim}}}}}0.8 rad of the observable section of the stream. Therefore, a stream with a larger surface density close to the cluster is favoured. This is achieved by a younger stream that has undergone one fewer pericentre passage. Using the method of Appendix H, we estimate the errors, obtaining $T=$2.6374_{-0.148,382,6}^{+0.257,618,37} Gyr. This result is closer to the $T=2.34\pm 0.82$ Myr obtained by Chen et al. (2025), indicating that the difference is partially caused by the exclusion of the stars in the envelope of the stream.

Appendix J Constant mass loss and frequency distribution

We eliminate the angular dependency of the double exponential model (Eq. 3) by setting the amplitude $A=0$. In this case, the model is reduced to f\mbox{\footnotesize(\theta_{r}^{\scalebox{0.5}{>!M}})}=C. This implies that the mass loss of the cluster is constant, and that the frequency distribution of the stripped stars along the principal axis of the stream $\Delta\bar{\varOmega}_{1}$ is Gaussian (Section 3.3). The values that best reproduce the results of the N​​ -body simulation are given by: $N_{\rm s}\simeq 29.371$ ${\rm stars}\times(4\pi/139\,{\rm rad})^{-1}$, $\operatorname{mean}(\Delta\bar{\varOmega}_{1})=0.391$ mrad Myr*-1*, and $\operatorname{std}(\Delta\bar{\varOmega}_{1})=0.115$ mrad Myr*-1*. Hereafter, we refer to this as the ‘Gaussian Model’, and to the model introduced in Section 3 as the ‘Variable Model’.

The top panel of Figure 22 shows the distribution of $\Delta\bar{\varOmega}_{1}$ of the stream stars from the N​​ -body simulation as a grey histogram. We estimate this distribution using only the stars stripped during the two complete periods of the simulation, which correspond to $\theta_{r}^{\scalebox{0.5}{\>\!M}}\in[-4\pi\,\text{,}\,0)$. This distribution has a sharp peak at $\Delta\bar{\varOmega}_{1}\in[0.25\,\text{,}\,0.3]$ rad Gyr*-1*, corresponding to the slow stars stripped around the apocentres. The stars stripped at the pericentres move faster and populate the tail of the distribution. The best-fitting Gaussian distribution is shown in red. This model cannot reproduce the skewness of the distribution or its peak. However, the Gaussian distribution accurately reproduces the far upper tail of the distribution, which ensures that a stream simulated with the Gaussian Model will be of a similar length to that obtained with the N​​ -body simulation.

In the middle panel of Figure 22, we show the PDF of the surface density of a simulated stream as a function of the angle $\Delta\bar{\theta}_{1}$. As in Figure 16, the black histogram corresponds to a stream of age $T=3039.7$ Myr, and the grey histogram corresponds to a saturated stream over the entire angular range depicted. Both streams are simulated using the Variable Model. The red histogram shows the surface density of a simulated stream of age $T=3039.7$ Myr computed using the Gaussian Model. This model cannot reproduce the underdensity next to the globular cluster ($\Delta\bar{\theta}_{1}=0$) or the peaks. In addition, it does not produce a constant saturation density, but a curvilinear surface density that underestimates \Delta\bar{\theta}_{1}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\sim}}}}}[0.3\,\text{,}\,0.9] rad and overestimates \Delta\bar{\theta}_{1}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\gtrsim}}}}}0.9 rad with respect to the Variable Model (black histogram).

In the bottom panel of Figure 22, we plot the distribution of stars from the main component of the simulated stream that pass the GDR3 selection function and the selection process (Section 5). The Gaussian Model (red histogram) presents a flat PDF within the range \Delta\bar{\theta}_{1}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\sim}}}}}[0.3\,\text{,}\,0.5] rad, which appears incompatible with the observational data shown in Figure 15. We test whether this model can accurately determine the accretion time of the simulated globular cluster by using a mock star sample obtained from the N​​ -body simulation, as described in Appendix G. We obtain an age of $T=1.51\pm 0.08$ Gyr, which overestimates the time since the first pericentre passage by approximately $159$ Myr.

The model with constant mass loss can be improved assuming that the frequency distribution of the stream stars along the principal axis of the stream $\Delta\bar{\varOmega}_{1}$ is a Log-Normal distribution of equation:

[TABLE]

with parameters $\mu\simeq 0.125$ mrad Myr*-1*, $\sigma\simeq 0.242$ mrad Myr*-1*, and $s\simeq 0.436$. We plot this distribution as a green line in the top panel of Figure 22. Although the Log-Normal distribution accurately describes the skewness of the distribution, it does not reproduce the peak corresponding to the slow stars stripped at the apocentre. Furthermore, the Log-Normal distribution overestimates the fraction of stars with \Delta\bar{\varOmega}_{1}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\gtrsim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\gtrsim}}}}}0.7 rad Gyr*-1*, resulting in simulated streams that are longer than those obtained with the N​​ -body simulation.

The middle panel of Figure 22 shows as a green histogram the surface density of a simulated stream of age $T=3039.7$ Myr computed assuming a constant mass loss and assuming that $\Delta\bar{\varOmega}_{1}$ follows the Log-Normal distribution. We refer to this model as the ‘Log-Normal Model’. As in the Gaussian case, this model does not reproduce the density variations close to the cluster or the density distribution slope around \Delta\bar{\theta}_{1}\mathrel{\mathchoice{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\displaystyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\textstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptstyle\sim}}}}{\raise 0.0pt\hbox{\scalebox{0.8}{\raise 0.0pt\hbox{\scriptscriptstyle\sim}}}}}1 rad obtained with the Variable Model (black histogram). However, it reproduces the constant density saturation level predicted by the Variable Model. In the bottom panel of Figure 22, we plot the distribution of the final star selection obtained using the Log-Normal Model as a green histogram. This distribution is similar to that obtained using the Gaussian Model. Therefore, the estimated age of the stream is also similar, being $T=1.60^{+0.08}_{-0.11}$ Myr, which overestimates the age of the N​​ -body stream by approximately $246$ Myr.

We conclude that the models assuming constant mass loss provide a good approximation of the stream surface density far from the progenitor cluster, where the density peaks are negligible. While the Gaussian Model can reproduce the length of the stream, it cannot reproduce the constant saturation density. Furthermore, it overestimates the age of the N​​ -body simulation of the M68 stream due to the inaccurate modelling of the density peaks. On the other hand, the Log-Normal Model requires an additional free parameter compared to the Gaussian Model, but it can accurately reproduce the stream saturation level. However, as with the Gaussian Model, it overestimates the age of the stream. In addition, it overestimates the length of the stream by a factor of approximately two. Therefore, in order to accurately estimate the age of the M68 stream, we consider it necessary to model the density variations caused by variable mass loss and variable frequency distribution of the stream stars.