Individual chaotic behaviour of the S-stars in the Galactic centre

Located at the core of the Galactic Centre, the S-star cluster serves as a remarkable illustration of chaos in dynamical systems. The long-term chaotic behaviour of this system can be studied with gravitational $N$-body simulations. By applying a small perturbation to the initial position of star S5, we can compare the evolution of this system to its unperturbed evolution. This results in the two solutions diverging exponentially, defined by the separation in position space $\delta_{r}$, with an average Lyapunov timescale of $\sim$420 yr, corresponding to the largest positive Lyapunov exponent. Even though the general trend of the chaotic evolution is governed in part by the supermassive black hole Sagittarius $\rm A^{*}$ (Sgr $\rm A^{*}$), individual differences between the stars can be noted in the behaviour of their phase-space curves. We present an analysis of the individual behaviour of the stars in this Newtonian chaotic dynamical system. The individuality of their behaviour is evident from offsets in the position space separation curves of the S-stars and the black hole. We propose that the offsets originate from the initial orbital elements of the S-stars, where Sgr $\rm A^{*}$ is considered in one of the focal points of the Keplerian orbits. Methods were considered to find a relationship between these elements and the separation in position space. Symbolic regression turns out to provide the clearest diagnostics for finding an interpretable expression for the problem. Our symbolic regression model indicates that $\left\langle\delta_r\right\rangle \propto e^{2.3}$, implying that the time-averaged individual separation in position space increases rapidly with the initial eccentricity of the S-stars.


Introduction
In the pursuit of understanding the nature of chaos and its manifestation in the universe, the S-star cluster is used as a laboratory for experimenting with the underlying astrophysics that give rise to chaos.The orbital data of 27 S-stars orbiting the supermassive black hole Sagittarius A * from (Gillessen et al. 2009) are employed as initial conditions for a Newtonian gravitational body simulation of the S-star cluster.The orbital evolution of the system was tracked over 10 4 yr by (Portegies Zwart et al. 2023;Boekholt et al. 2023).When perturbing the initial conditions of the star S5 by displacing its initial  coordinate by  = 15m, the evolution of the perturbed solution diverges exponentially from the canonical solution, implying that the system is chaotic.
In this paper, we present the methods used to quantify chaos in the simulation of the S-star cluster and delve into the individual differences in the chaotic behaviour of the S-stars and Sgr A*.For more details on the simulation, refer to (Portegies Zwart et al. 2023).

Measuring chaos
In the gravitational -body problem, chaos can be measured with the Lyapunov timescale.The Lyapunov timescale represents the timescale on which the system becomes unpredictable.We describe the chaotic S-star orbital evolution with phase-space distance.Phase-space distance as a function of time is defined as ★ E-mail: beckers@strw.leidenuniv.nl★★ E-mail: cpoppelaars@strw.leidenuniv.nl where ì  and ì  are the position and velocity vectors of an Sstar, and  and  denote the canonical and perturbed solutions, respectively.We define an (initial) perturbation at  = 0, (0).The evolution of () is then approximately described by an exponential function with time dependency, where  is the maximum positive Lyapunov exponent.The growth factor,   (), is the value of  at some time  as a fraction of the initial perturbation, (0), i.e.
From this equation, one can find  = ln(G  )/t, which is the reciprocal of the Lyapunov timescale, . (4)

Separation in position space
We calculated the time evolution of the separation in position space between the canonical and perturbed solution,   , for each S-star and Sgr A*.We show this time evolution in the left panel of  subtracted from the time evolution of log 10 (  ) for each S-star.The grey shading shows the region between the non-reduced maxima and minima of   over the whole system at each timestep, indicating the magnitude of reduction in the spread of the curves, though its lower section is obscured by the coloured curves.
Figure 1.The separation in position space grows approximately exponentially, as expected from Equation 2. From a least-squares fit to the curve of each star and the black hole, we infer an ensemble mean Lyapunov timescale of   ≃ 420 yr.The black hole is less sensitive to the perturbation than most stars, as seen from its   values, which are generally lower in magnitude than those of the S-stars.We attribute this to the mass difference of order 10 5 M ⊙ between the central body and the stars, leading to much higher inertia.This stabilizes the black hole's position, even in the presence of perturbations.Perturbations caused by close encounters between stars are propagated through the entire system, driven by feedback from the black hole (Portegies Zwart et al. 2023).Hence, big events with multiple close encounters, such as the one at  = 2876 yr, influence the general trend of all curves.However, it is not immediately evident why there are vertical offsets in the S-star curves, as the S-stars have been given identical masses in this simulation.Furthermore, while differences in the shapes of the individual curves could account for an offset, the curves as a whole exhibit a noticeable shift relative to each other.

Position space separation offsets
The systematic offsets found in the position space separation curves must be an effect of some underlying astrophysics.The (initial) S-star orbits can introduce variations in the evolution of   of each star, so it is worthwhile to investigate the Keplerian elements.

A relationship between the offsets in position space separation and Keplerian elements
To map the separations of S-star curves and Sgr A*, we subtract the black hole curve from each S-star curve, and take the temporal mean for each S-star; we define this as Δ S-star BH log 10 (  ) .To find a relationship between this and the Keplerian elements of the stars, we use the initial orbital parameters from (Gillessen et al. 2009).We expect to see that a low semi-major axis () corresponds to a higher magnitude in   , as stars that are close to the black hole should be more sensitive to the changes inside the system.Moreover, eccentricity () and   should be directly correlated, since a highly elliptical orbit has a closer pericenter to the black hole.In Figure 2, we show the relationship between ,  and Δ S-star BH log 10 (  ) .The pattern that emerges in the colour gradient in Figure 2 hints that eccentricity and semi-major axis play a major role in the curve separations.The influence of eccentricity and semi-major axis on the mean separation are both important.The other orbital parameters can be compared with each other in a similar fashion, however, no apparent pattern is found.This observation suggests that their influence on Δ S-star BH log 10 (  ) is negligible, compared to that of the eccentricity and semi-major axis.Furthermore, the stars with the highest colour values in Figure 2 are S21 and S29. (Portegies Zwart et al. 2023) demonstrate that S21 is one of the stars that were part of the big event at 2876 yr.Moreover, S29 was shown to have 6 close encounters during the simulation.In addition, S67 has 10 close encounters, the most of any star.Therefore, despite its relatively low eccentricity and large semi-major axis, it is not green in Figure 2 as opposed to its nearest neighbour, S87.

Symbolic regression
We adopt PySR (Cranmer 2023), a symbolic regression Python package for discovering interpretable analytical equations that describe an underlying pattern in a dataset.Only semi-major axis and eccentricity are used in the model, since symbolic regression may not provide accurate results when the dimensionality of the data is high (Matchev et al. 2022).Moreover, PySR does not correctly interpret the evolution of   ,  and  for each star, suggesting it does not support time-series data1.Therefore, we provide the model with Δ S-star BH log 10 (  ) and the initial semi- major axis and eccentricity.The results can be found in Table 1.
We use this to estimate   approximately as Subtracting Equation 5 from log 10 (  ) reduces the offsets significantly, as shown in the right panel of Figure 1.

Conclusions
We aimed to characterize the individual behaviour of stars during the evolution of a Newtonian chaotic dynamical system.Using a gravitational -body simulation of the S-star cluster, we derive the Lyapunov timescale of the system over 10 4 yr to be ∼420 yr.
In the position space of the S-stars, we find a vertical offset between their curves.A relationship with the Keplerian orbits of the S-star was proposed to explain this offset in   .We adopt symbolic regression to find the relation ⟨  ⟩ = 2 • 10 4 10 −0.00018446/au  2.3 , where  and  are taken from the initial orbit.We conclude that ⟨  ⟩ ∝  2.3 ; the time-averaged individual phase-space distance with respect to the black hole increases rapidly with the star's orbital eccentricity.

Fig. 1 .
Fig. 1. (left) Time evolution of separation in position space for each S-star and the central black hole.(right)  = −0.00018446•/au+log()+4.3103subtracted from the time evolution of log 10 (  ) for each S-star.The grey shading shows the region between the non-reduced maxima and minima of   over the whole system at each timestep, indicating the magnitude of reduction in the spread of the curves, though its lower section is obscured by the coloured curves.

Fig. 2 .
Fig. 2. Initial semi-major axis and eccentricity of each S-star, with Δ S-starBH log 10 (  ) indicated by colour.Each S-star is labelled, according to the Gillssen catalogue(Gillessen et al. 2009).A colour gradient can be seen, ranging from orange in the bottom right to pink in the top left.