Punctuated chaos and the unpredictability of the Galactic Centre S-star orbital evolution

ABSTRACT

We investigate the chaotic behaviour of the S-star cluster in the Galactic Centre using precise N-body calculations, free from round-off or discretization errors. Our findings reveal that chaos among the Galactic Centre S-stars arises from close encounters, particularly among pairs and near the massive central body. These encounters induce perturbations, causing sudden changes in the orbital energies of the interacting stars. Consequently, neighbouring solutions experience roughly exponential growth in separation. We propose a theory of ‘punctuated chaos’ that describes the S-star cluster’s chaotic behaviour. This phenomenon results from nearly linear growth in the separation between neighbouring orbits after repeated finite perturbations. Each participating star’s orbit experiences discrete, abrupt changes in energy due to the perturbations. The cumulative effect of these events is further amplified by the steady drift in orbital phase. In the Galactic Centre, perturbations originate from coincidental encounters occurring within a distance of ≲100 au between at least two stars (in some cases, three stars). Our model satisfactorily explains the observed exponential growth in the 27 S-stars cluster. We determine that the S-star system has a Lyapunov time-scale of approximately 462 ± 74 yr. For the coming millennium, chaos in the S-star cluster will be driven mainly by a few of the closest orbiting stars: S2, S5, S6, S8, S9, S14, S18, S31, S21, S24, S27, S29, and S38.

1 INTRODUCTION

Newton’s laws of motion lead to chaos. This chaotic behaviour is often quantified by measuring the growth in time of the separation, δ, between two neighbouring solutions, i.e. by solving the equations of motion of the multibody system twice, once with and once without a small change in the initial conditions. If the evolution of δ is roughly exponential, but δ is still small by the end of the calculations, we call the evolution chaotic.

Chaotic behaviour can be quantified using the Lyapunov time-scale, which (in a more rigorous treatment) is the reciprocal of the maximum positive Lyapunov exponent. The Lyapunov time-scale has been measured for only a limited number of multibody systems, because these measurements are expensive in terms of computer time. It is not even clear that a specific system, within the uncertainty of its initial conditions, leads to a unique Lyapunov time-scale, because the available parameter space may have an irregular structure with stable and chaotic regions (Hayes 2008).

In the Solar system, chaos is mainly driven by resonant overlap (Chirikov 1979; Wisdom 1980; Tamayo et al. 2021; Mogavero & Laskar 2022; Rath, Hadden & Lithwick 2022), which has some resemblance with resonant relaxation in the galactic centre (Sridhar & Touma 2016), or even with violent relaxation (Kandrup, Vass & Sideris 2003). Previous studies, however, found that repeated weak scatterings among minor bodies also form a major driver for changes in the orbital parameters, for example in studies on the orbit of mercury (Laskar & Gastineau 2009), but also in Halley’s comet (Boekholt et al. 2016) to be chaotic on a time-scale of less than 3000 yr, as it interacts with Venus and Jupiter. In numerical studies of planet–planet scattering experiments in hypothetical multiplanet systems chaotic behaviour is also recognized to be driven by encounters rather than resonant overlap (Chatterjee et al. 2008; Jurić & Tremaine 2008). This event-driven chaotic behaviour does not comply with Chirikov’s resonant overlap paradigm.

This finding hints towards an event-driven process, where resonances are not required, but where chaos is driven by close encounters. We informally refer to such behaviour as ‘punctuated chaos’. A similar chaotic behaviour was found in stellar clusters (Goodman, Heggie & Hut 1993), where mutual interactions between stars seemed to drive chaos in the system. A system consisting of many eccentric and inclined orbits behaves differently from a flat system with non-crossing orbits and a dominant central body. The S-star cluster in the Galactic Centre (Ghez et al. 2008; Gillessen et al. 2009) may be a good example where resonant overlap does not form the major driver for chaos.

Understanding chaos in self-gravitating systems is important for understanding a wide variety of astronomical phenomena, including the sources of gravitational waves, extreme mass ratio inspirals, and the probability that a minor body hits the planet Earth. Our picture of punctuated chaos is also important for other Hamiltonian systems, such as the multibody pendulum (Hoover & Griswold Hoover 2017), and it may be responsible for a slow down of chaotic behaviour in critical fluctuations (Das & Green 2019).

In our analysis of the S-stars, we will neglect low-mass stars, unobserved compact objects, and other interstellar material (planets, asteroids, dust, and gas). The low-mass components are expected to be copious and important, but we limit ourselves to study chaos in an idealized system of black hole with S-stars. Chaos in the studied idealized system then probably only provides an upper limit to the actual Lyapunov time-scale.

2 CHAOS AS AN EVENT-DRIVEN PROCESS

2.1 The Lyapunov time-scale

Consider a Keplerian two-body system, e.g. a star with a planet or supermassive black hole with a star. The difference between two neighbouring solutions, started with slightly different initial conditions, grows approximately linearly with time. The orbital phases of the two solutions gradually diverge, until their separation δ has grown to the size of the orbits.

In a three-body system consisting of a binary and a third body, the binary orbits are repeatedly perturbed by the third body. The effect of these events on two neighbouring binary orbits depends on their separation δ, and can cause the two solutions to diverge more quickly than in a two-body system. A sequence of perturbations, each with a subsequent period of linear divergence until the next perturbation, can then drive exponential divergence between two initially almost identical systems (Section 2.2). The degree of chaos, measured in terms of the Lyapunov time-scale, can then be derived by the frequency and strength of perturbations in the system (Section 4.1).

In direct N-body calculations we can quantify chaos by measuring the Lyapunov time-scale. Instead of a measure over the system’s entire lifetime, which would be the setting for a rigorous treatment, a measure over a finite time interval is more practical in our case. We therefore define the growth factor G_δ(t) following from an initial separation δ(0) after some time t. The evolution of the separation is then described as an exponential function of time δ(t) = δ(0) e^λt, where λ is the Lyapunov exponent. The growth factor can then be written as G_δ ≡ δ(t)/δ(0) = e^λt, and λ = log (G_δ)/t, or in terms of the Lyapunov time-scale, t_λ = 1/λ = t/log (G_δ).

2.2 Punctuated chaos

Based on the event-driven process described above, we derive approximate expressions for the Lyapunov time-scale. We start again by considering a particle in a Kepler orbit around a much more massive body. With an initial semimajor axis a₀ and total mass m, the initial orbital frequency |$\omega _0 = \sqrt{Gm/a^3_0}$|⁠. Let a neighbouring solution be separated by an infinitesimal displacement δx₀ at time t = 0. This displacement leads to a small difference in the semimajor axis of the same order, i.e. δa₀ ∼ δx₀. The resulting difference in frequency is

The separation of the two motions along the orbit grows with time t > 0 according to

The growth of the initial displacement is linear with time from t₀ to t, but such that δa remains constant as the growth is along the orbit, i.e. the growth is in the orbital phase rather than in energy (or angular momentum).

Now suppose that an instantaneous perturbation acts on the motion at time t₁, causing the velocity of the Kepler motion to receive a slight kick.

There is a resulting change in energy, and therefore a change in semimajor axis Δa₁. But since the perturbation will depend on the position of the body, the difference in semimajor axis between the two motions, δa, will also change, by an amount that we call δa₁. We suppose for the sake of argument that δa₁ ∼ δx₁, i.e. the difference in position at the time of the event. (In Section 4.1, we give a specific example in which it is easy to see that δa₁ ∝ δx₁.) The new difference in semimajor axis implies a difference in orbital frequency δω₁ ∼ δa₁ω₁/a₁ at time t₁. Thus for t > t₁ the displacement varies as

If a second perturbation occurs at time t₂ > t₁, we can see from equations (2) and (3) that the displacement is

If these perturbations recur at roughly comparable intervals Δt, and if ω does not change by a large factor, it can be seen that the displacement at some large time t will be

Thus the linear growth of equation (2) transforms into exponential growth. The corresponding Lyapunov exponent is |${\cal O}(\omega)$| if ωΔt ≲ 1; it is of order the reciprocal of the crossing time. The case ωΔt ≳ 1 is also of interest and leads to a smaller estimate of order ln (ωΔt)/Δt.

Though we here considered perturbed Keplerian motion with one dominant body, which is relevant to our discussion of the S-stars in Sections 3 and 4, the conclusions are valid more widely. As an example, we consider resonant three-body scattering events, which can be viewed as a prolonged sequence of perturbations of the Kepler motion. So long as the three particles remain democratic (i.e. at comparable distances) and are of comparable mass, the perturbations in any of the three two-body motions will be of order unity and occur at intervals of order the crossing time, t_cr. Therefore, as in the above discussion, the Lyapunov exponent will be of order 1/t_cr. Numerical examples of Dejonghe & Hut (1986) indeed show that the separation of neighbouring solutions grows roughly exponentially until dissolution of the resonance (in the sense of a democratic three-body behaviour). On the other hand, if the evolution of a triple system is dominated by protracted excursions of the third body, of order T ≫ t_cr, then the estimate will decrease to one of order 1/T (following the result for the case ωΔt ≳ 1, and neglecting a logarithm). Usually, the evolution is a mix of prolonged excursions interspersed with periods of frequent interplay (Szebehely 1971), and the Lyapunov exponent λ will be 1/T ≲ λ < 1/t_cr, where T is the duration of the longest excursion. The lifetime of the Pythagorean problem, for example, is about 16t_cr (Aarseth 2003, their p. 238), and the growth of the separation of neighbouring solutions in this time is about 8.5 dex (Dejonghe & Hut 1986). Thus T ≲ 16t_cr and λ ≃ 1.2/t_cr, and the finite-time Lyapunov exponent lies in the range 1/T ≲ λ ≲ 1/t_cr.

The result of the model (that the Lyapunov exponent λ is of order 1/t_cr for comparable masses without lengthy excursions) is consistent with the results in Goodman et al. (1993), who considered the general N-body problem. This is a rather independent confirmation, as their model was also based on assuming that the separation between neighbouring solutions grows as a result of two-body encounters (see Table 1), but otherwise did not follow the same approach as ours.

3 APPLICATION OF PUNCTUATED CHAOS

We now study punctuated chaos in an actual astrophysical system, namely the S-star cluster in the Galactic Centre. All stars affect each other and we therefore stretch the assumptions of extreme mass ratio, the approximation of only two orbiting particles, and the three-dimensional nature of the problem.

Our analysis starts by acquiring converged solutions for the dynamical evolution of the S-stars. Acquiring a converged solution is important because round-off and discretization errors grow exponentially, rendering a non-converged solution inappropriate for studying chaos. We acquire converged solutions to the N-body problem by integrating Newton’s equations of motion. The dynamics near supermassive black holes is best described with general relativity rather than by Newtonian dynamics. Still we perform our simulations using a Newtonian approximation and ignore the theory of general relativity. We consider this appropriate because the highest speed in our simulations (the star S2 at pericentre) does not exceed about ∼0.017 times the speed of light. At these relatively low speeds, the system behaves like a Newtonian system, at least within our brief simulation time of 10⁴ yr (Portegies Zwart et al. 2022). In addition, our model for punctuated chaos was derived for Newtonian systems. In the presence of the black hole, in future it will be worth exploring the effect of general relativity on the chaotic behaviour of the S-star cluster.

The results presented here are acquired using converged solutions. These result from a procedure in which the length of the mantissa, controlling precision, and the accuracy of the integrator are improved at each iteration step, as is explained in Portegies Zwart & Boekholt (2018). Such converged solutions can be achieved using brutus. In brutus, we control round-off by extending the numerical precision, and accuracy by reducing the tolerance in the Bulirsch–Stoer integrator (Bulirsch & Stoer 1964). By repeating the same calculation with higher precision and better accuracy, we eventually reach a solution for which the results remain identical to a pre-determined number of decimal places; we call this the converged solution.

3.1 Measuring chaos in the S-star cluster

We consider a realization of the S-star cluster, consisting of 27 early-type stars that orbit the supermassive black hole in the Galactic Centre. We adopted the orbital parameters of the 27 S-stars reported in Gillessen et al. (2009, their table 7). The numbering of these stars is not simply S1–S27, because we only use those stars for which an orbit is determined.

The initial conditions are generated using the amuse (Pelupessy et al. 2013; Portegies Zwart et al. 2013) routine generate_Sstar_cluster.py from Portegies Zwart & McMillan (2018). This routine adopts the cited orbital elements for the S-stars, and solves Kepler’s equation to acquire Cartesian coordinates and velocities of the S-stars. We assume each S-star to have a mass of 20 M_⊙, and adopt the central massive black hole mass of 4.154 × 10⁶ M_⊙ (Gravity Collaboration 2019). We integrate this system for 10 000 yr until the solution converges using brutus (Boekholt & Portegies Zwart 2015).

In Fig. 1, we present the converged solution of the projected orbits of the S-stars integrated over 10 000 yr. This solution was obtained using a word length L_{|$\mathrm{ w}$|} = 128 bits, and a tolerance in the Bulirsch–Stoer integrator of ϵ = 10⁻²⁴ leading to a total final phase-space error less than 1/10⁷, and a relative energy error in the integration of dE/E ≃ 10⁻¹⁵. Each calculation took about 20 h on the single core of a Xeon W-11855M processor. For L_{|$\mathrm{ w}$|} = 116 bits precision and a tolerance of ϵ = 10⁻²¹ the solution is converged to four decimal places, which would have sufficed for this study. Regular double precision, however, would have been insufficiently precise and would not have led to the required accuracy.

Although the orbits in Fig. 1 appear to follow a Keplerian orbit nicely, in Fig. 2 we show that this is not the case. Here we show a detail of the orbits of several stars from Fig. 1. Each orbital revolution is indicated with a thin line. The lines do not fully overlap, indicating that their orbits change with time.

After having established a converged simulation for the S-stars (see Fig. 1), we introduce a relative shift by translating the Cartesian x coordinate of star S5 by 10⁻¹⁰. This amounts to moving the initial position of S5 by ∼15 m in the x direction.

Instead of the separation in position space for each individual S-star, as illustrated in Fig. 1, we can also present the more usual total phase-space distance as a function of time, which we define as (Dejonghe & Hut 1986)

Here 0 and 1 refer to the original and shifted solution, respectively. In Fig. 3, we present the evolution of the phase-space distance between the two solutions. The overplotted dotted line (in red) corresponds to a Lyapunov time-scale of t_λ ≃ 450 yr.

In Fig. 4, we present the evolution of the root-mean-square difference in the semimajor axes between the two solutions for the S-stars. Here we define

The events are visible at discrete times followed by a constant difference δ_a. In Fig. 3, a particularly strong event is visible in the phase-space distance evolution around t = 2876 yr. According to punctuated chaos the orbital separation changes abruptly during events. We identify several with the colours and the vertical dashed lines (see Fig. 4). We automatically recognize these events in the simulation data using the pelt algorithm (Killick, Fearnhead & Eckley 2012), using the piecewise constant model. Both free parameters in the pelt algorithms (the controls the minimum distance between change points and the grid of possible change points) were set to unity.

To further investigate which stars are involved, we apply the pelt algorithm to individual stars. Three examples are presented in Fig. 5. Here the colours are meant to guide the eye and indicate the events already identified on the full set of 27 stars in Fig. 4. The vertical red-dashed lines are the events as identified using the pelt algorithm for each individual star.

In Table 1, we present the times at which events were detected, which stars are associated with these events, and at what distance from the black hole (in terms of au and in distance with respect to the sum of the stars’ Hill radii). The stars S2 and S18 are frequent participants, and they form the main drivers of chaos. We also present in Table 1 the values for r_ij, and r_ij/r_H. The former is defined as the distance between the two stars i and j nearest to the black hole at the moment of closest approach to the black hole. The latter gives the relative distance in terms of the radius of the Hill sphere of the two participating stars (second column). The last column gives the relative velocity between the two S-stars at the moment of closest mutual approach.

In Fig. 6, we present the results of applying the theory of punctuated chaos to the evolution of the S-star cluster over 10 000 yr. Here we overplotted the actual data and the filtered curve, as in Fig. 3. The red curve results from measuring the time and magnitude of the events in semimajor axis (see Fig. 4 and Table 1), which are then applied directly to our model for punctuated chaos (see Section 2.2).

In Appendix A, we present the algorithm based on equation (6) that is used to draw the red curve in Fig. 6. The events (moments of close encounter and the mutual distance between the two S-stars at peribothron) are measured in the N-body simulations, and presented in Table 1. Note that several punctuated events listed in Table 1 are not directly related to the jumps in the semimajor axis presented in Fig. 4 (see also equation 7). For individual stars these jumps are visible, but they tend to be obscured by the ensemble averaging. Each of the events listed in Table 1, however, do relate to a jump in the semimajor axis in the S-stars associated with the close encounter. These jumps are visible in the evolution of the semimajor axis for individual S-stars; as examples, we present this evolution for S6, S27, and S66 in Fig. 5, which also show the early (13, 408, and 807 yr) and late (7007 yr) jumps listed in Table 1 (red dashed lines).

The comparison between the theory and the simulations is striking. Although the red curve in Fig. 6 does not perfectly match the filtered curve (black), it is astonishingly similar. A slight deviation where the red curve overshoots the black curve, around t = 7007 yr, seem to be caused by several close encounters between S2 and S18. In our model, we anticipate on standalone encounters, and the effect of multiple encounters in short succession between the same stars can also cause the evolution of δ to decay.

3.2 The big event at 2876 years

The largest event happens around t = 2876 yr after the start of the simulations.

In Fig. 7, we present the distance of the nearest star to the black hole (blue), and the distance between this star and its next-nearest stellar neighbour (orange), around the moment of the big event at t = 2876 yr. Along the top we identify the stars involved in their close proximity. It starts at t = 2870 yr with S2 being the closest to the black hole, and S14 being closest to S2. At t ≃ 2872 yr star S6 takes over the closest position to S2 from S14. At t ≃ 2873 yr star S6 becomes the closest star to the black hole, and at t ≃ 2874 yr star S21 becomes the closest neighbour of star S6, until a little before t = 2876 yr star S6 and S21 have their mutual close encounter with each other and with the black hole.

In Fig. 8, we present for the stars S6 and S21 the evolution of δ_r (left-hand panels), i.e. the spatial separation between the two orbits, the semimajor axis (middle panels), and the difference in orbital frequency (right-hand panels) around the moment of the big event at t = 2876 yr. In the evolution of the spatial separation (left-hand panel) the consequences of the event are visible. The main driver of this evolution is the almost instantaneous change in orbital energy (middle panel, expressed here in semimajor axis) and the separation in orbital frequency (right-hand panel).

4 DYNAMICS OF PUNCTUATED CHAOS IN THE S-STARS

The picture of punctuated chaos in Section 2.2 lacks any dynamical characterization of the events. They are simply circumstances that change the frequency of the orbital motion. Even at the close of Section 3.1, in the construction of Fig. 6, the data on the changes in frequency are simply read from a full numerical simulation. In this section, we attempt to explain, in order of magnitude, how these changes arise, on the basis of few-body interactions.

For the purposes of this discussion, let us reduce the problem to a three-body system consisting of a black hole (mass M, particle 0) and two S-stars (mass m, particles 1, 2). Let r_i be the position vector of star i relative to the black hole, and focus on star 1. Its equation of motion is

with obvious notation r_i, r₁₂. The first term on the right is the acceleration of the motion of star 1 relative to the black hole, the second term is the direct perturbation by particle 2, and the third term is the indirect perturbation, caused by the acceleration of the black hole by particle 2. (If we were to write down the equation for the case of N S-stars, the second and third terms on the right would be replaced by sums from j = 2 to N, with the index 2 replaced by j.)

4.1 Direct perturbations

From the previous section (especially Fig. 7 and Table 1), we know empirically that the major events are associated with a close approach (to a distance of order 10–50 au) between two stars (see Table 1), at times where the distances to the black hole are larger by an order of magnitude. If we label the stars in the encounter as 1 and 2 (ignoring all other stars for the time being), then it is clear that the perturbation on the motion of particle 1 is dominated by the second term on the right of equation (8), i.e. the direct perturbation, due to the gravitational attraction of star 2. For the moment we also ignore the third term, but return to it briefly in Section 4.2.

We now consider the growth of chaos in the single S-star that we have hitherto labelled as star 1, and suppose that it has a succession of close encounters with other stars (in ‘events’) at some interval Δt. Let δx denote the variation in some quantity x, such as energy or position, between two neighbouring solutions. Specifically, let δE₀ be the variation in specific energy E ≡ |v²/2 − GM/r| at a time just before an event, where v is the speed, and let δr₀ be the variation in position at the same time. (Henceforth the subscripts denote an index of events in a sequence of events.)

Then we have

where δE₁ is the variation just before the next event, which will be the same as the value just after the first event, and ΔE is the change in E at that first event. Now we estimate this as

where d is the distance of closest approach in the encounter, and V is the relative speed of the two stars in the encounter. Both these parameters are listed for the punctuations measured in the S-star cluster in Table 1. (This estimate is made by taking Δv to be the maximum perturbation Gm/d² multiplied by the time-scale of the encounter, d/V.) We use the same estimate |$\sqrt{GM/r}$| for V as for v, so that

Its variation can be estimated as

and so equation (9) becomes

After this event the variation in the orbital frequency will be δω₁, and so

where a is the semimajor axis, just as in the simplified model of punctuated chaos described in Section 2.2. Since ω²a³ = GM and E ∼ GM/a (though strictly here, M should be replaced by M + m, but we ignore the stellar mass compared with the black hole mass), we can re-express this as

By equation (13), this in turn becomes

Equations (13) and (16) are explicit estimates that allow us to map the effect of an encounter on δE and δr. It has matrix

The larger eigenvalue of this matrix gives the factor by which the variation, in either E or r, is multiplied as a result of an event.

To get a feel for the magnitude of this evolution, consider the six encounters in Table 1 that involve star S2. The mean distance of closest approach in these is about 40 au, which we take for the value of d. Since the duration of the simulation is |$10\, 000$| yr, for the mean time between these six events we adopt Δt = 1700 yr. The semimajor axis of S2 is a_S2 ≃ 1000 au, and we take m = 20 M_⊙ and M = 4 × 10⁶ M_⊙, as we have done throughout. (These units imply that G = 4n² ∼ 40.) Then the quantity

and the larger eigenvalue λ ≃ 4. Thus the variations increase by this factor in 1700 yr, or by almost 4 dex in 10⁴ yr. While this estimate is only about half of the logarithmic growth seen in Figs 3 and 6, it only includes half of the events listed in Table 1. The remaining events involve stars other than S2, which is why they were omitted from this estimate, but we show in Section 4.2 how their influence spreads to all stars, S2 included.

The foregoing argument can be criticized on several grounds, and the main one is that it still depends on the numerical integration, as it draws on the distance of the observed two-body encounters. This may be estimated independently as follows. From the initial conditions of the integration we can estimate the central number density n of the S-star cluster from the distance of the sixth nearest neighbour to the black hole (Hénon 1971), which is about 2000 au. Hence n ≃ 2 × 10⁻¹⁰. The median semimajor axis a is about 3000 au, and so we may estimate the typical speed |$v = \sqrt{GM/a}$| as about 200 au yr⁻¹. It follows that a typical S-star i will experience one two-body encounter with another S-star j in 10⁴ yr at a distance less than about r_ij ≃ 60 au. One or two will do so at still smaller distances, consistent with what is seen in the numerical results.

4.2 The role of the indirect acceleration

While Section 4.1 makes a reasonable case for supposing that direct encounters between pairs of S-stars are a major driver of chaos in the system, it is not clear how many stars are affected, as only a fraction of all S-stars appear in Table 1. Of course all stars are subject to two-body direct perturbations, but it might be thought that weak perturbations would drive chaos on a longer time-scale. We now consider the effect of indirect perturbations, and will show that these affect all stars. Furthermore, they keep the growth of chaos in all stars in lockstep.

The indirect perturbation of star 2 in equation (8) is given by the third term on the right, and clearly depends on its distance r₂ from the black hole. Thus we shall consider the pericentre passages of star 2 as the driver of this component of punctuated chaos. In one such event, whose duration is of order q/v₂, where q is the pericentre distance of star 2, the change in energy of star 1 may be estimated as

Much as in Section 4.1, we estimate |$v_2 \sim \sqrt{GM/q}$|⁠, but for the first star we take |$v_1\sim \sqrt{GM/a}$|⁠, where a is the semimajor axis of star 1, since this star can be anywhere on its orbit when the event occurs. Thus

and so the variation in this quantity is

Now suppose that star 2 is one of those that is vigorously affected by direct interactions with other stars, and its variations have Lyapunov exponent λ. Then we can estimate that

where δq₀ is some constant. Next, substituting this estimate into equation (20), and adding the effect of a succession of such events, we readily see that the sum of the corresponding exponentials can be estimated by the effect of the most recent event. Thus the variation in the energy of star 1 at time t may be estimated as

Thus we conclude that chaos in star 1 grows with the same Lyapunov exponent as star 2, and this would be true even if star 1 were not subject to any direct perturbations.

In Fig. 9, we plot the time evolution of the separation in position space for each of the 27 S-stars, as well as for the central black hole. We observe that, on average, all particles diverge exponentially and with roughly the same rate. By fitting linear slopes to each curve in the figure, we measure an average individual Lyapunov time-scale of 462 ± 74 yr. Fig. 9 shows that all stars, and also the central black hole, diverge at the same rate on average. For the black hole the amplitude is smaller than for the bulk of the S-stars by roughly the ratio of their masses, i.e. 2 × 10⁵.

5 CONCLUSIONS

In order to explain short time-scale chaos in self-gravitating democratic dynamical systems, we outline the theory of punctuated chaos. The essence of punctuated chaos is that instantaneous close encounters, or events, drive the exponential growth due to ephemeral perturbations, each of which is followed by a gradual drift in phase space. Our description of punctuated chaos delivers a qualitative and quantitative description of chaos in self-gravitating systems.

We test this description on converged direct N-body simulations of self-gravitating systems under Newton’s equations of motion. We confirm the working of punctuated chaos on the chaotic orbits of the S-stars in orbit around the supermassive black hole in the Galactic Centre. The rather similar values of the mean orbital period of the S-stars and their Lyapunov time-scale then is no coincidence.

It turns out that also comet Halley has a Lyapunov time-scale quite comparable to its orbital period (Boekholt et al. 2016). In this case the chaos is driven by interaction with Venus and Jupiter. We argue that the Lyapunov time-scale under punctuated chaos is proportional to the mean interacting orbital period of the system, which for our selected systems is of the order of a few 100 yr.

According to our arguments and findings (Section 4.2), major bodies are affected on a similar time-scale, but the amplitude δ scales roughly with the ratio of the perturber to the perturbed mass, as one sees for the black hole in Fig. 9.

For the 27 S-stars in the Galactic Centre, we derive a Lyapunov time-scale of 462 ± 74 yr. This is somewhat comparable to the mean orbital period of the S-stars, or 〈P_orb〉 = 269 ± 383 yr, consistent with the prediction from punctuated chaos.

We also qualitatively compare the measured phase-space distance evolution with the theory of punctuated chaos (see Section 3.1, especially Fig. 6).

In the theory (Section 3), the number of interactions (events) and the strength of these interactions are free parameters. In the comparison presented in Fig. 6, the moment and strength of each interaction are measured through an N-body simulation. But we also show (Section 4.1) that these parameters can be estimated in order of magnitude without resort to an N-body simulation, but from estimates of the density and kinematics of the inner S-stars. This procedure allows us to directly compare the theory with the simulations, as well as to identify the events that drive the exponential growth.

If the previous analysis holds, chaos manifests itself from repeated small perturbations that each induce a linear response in the separation between neighbouring solutions of the system. The composition of these linear responses resembles an exponential behaviour and drives chaos in the system. The perturbation, initiated by a close encounter between two (or more) stars, is subsequently communicated to the rest of the system by the perturbed phase-space characteristics of the massive central body. In the case of the S-stars this is the supermassive black hole, and it is the Sun for comet Halley and the rest of the Solar system. In this way the entire system is subject to the same smallest Lyapunov time-scale.

ACKNOWLEDGEMENTS

It is a pleasure to thank Geneva Observatory and in particular George Meylan, Silvia Extröm, and Steven Rieder for their hospitality. This work started with the Masters’ thesis of Jaro Molemkamp at Leiden University under the cosupervision of Veronica Saz Ulibarrena.

SFPZ thanks Deutche Bahn for the availability of green electricity in combination with long delays on trips to from Amsterdam to Geneva, Copenhagen, and München, which enabled part of the writing, further code development, and verification calculations.

Computing resources. The majority of the computational resources were provided by the Academic Leiden Interdisciplinary Cluster Environment (ALICE), and using LGM-II (NWO grant no. 621.016.701). Software we used includes amuse (Portegies Zwart et al. 2018), python (Van Rossum & Drake 2009), pelt (Killick et al. 2012), numpy (Oliphant 2006), scipy (Jones, Oliphant & Peterson 2001), mpfr (Fousse et al. 2007), brutus (Boekholt & Portegies Zwart 2015), matplotlib (Hunter 2007), and linux (see https://github.com/torvalds/linux/releases/tag/v4.1-rc8).

Energy consumption of this calculation. The calculations using brutus are elaborate and took about 7 × 10⁶ CPU seconds. Data processing and analysis require about one tenth, totalling about 3 months of dual CPU usage. Using the tool http://green-algorithms.org/ and Portegies Zwart (2020), we calculated our energy consumption to be about 160 kWh.

DATA AVAILABILITY

Raw and reduced data and processing scripts are available at figshare DOI: 10.6084/m9.figshare.13637672

References

APPENDIX A: CALCULATING THE EVOLUTION OF THE PHASE-SPACE DISTANCE

We calculate the evolution of the phase-space distance, as presented in Fig. 6 (red curve), using the algorithms presented in the form of a python script using the units module from the Astrophysical Multipurpose Software Environment (amuse; Portegies Zwart & McMillan 2018) in the listing below. The events (time T and mutual distance r_ij between the two closest S-stars i and j near the black hole) used to draw the curves are presented in Table 1. These events were automatically detected by identifying the closest encounters between at least two stars near the black hole. Here we included only those encounters with a mutual distance r_ij ≲ 60 au. We introduce the constant k = 0.5 au that relates the effect of the perturbation to the orbital deviation. Looping over time, we calculate the new phase-space distance between the two solutions, given the perturbation introduced in the last close encounter (see equation 3). We empirically determine ω₁ ∝ (k/r_ij)² from the N-body simulations by measuring r_ij and fitting k based on the evolution of the growth of the phase-space distance δ (see equation 6). The stated dependence on r_ij here is justified by the d-dependence in equation (13).

from amuse.lab import units

def phase_space_distance(T, rij, k = 4|units.au):

 delta = [0] | units.m

 t = [0] | units.yr

 dd = [1] | units.m

 while t[-1] < T[-1]:

 d = dd[-1]*(1 + ((t[-1]-T[len(dd)-1]).value_in(units.yr))*(k/rij[len(dd)-1])**2)

 if t[-1] > = T[len(dd)]:

 dd.append(d)

 t.append(t[-1] + dt)

 delta.append(d)

 return t, delta

The input arrays, T and rij are determined from close encounters between the S-stars and the black hole, as explained in Section 3.1 and present in Table 1. Note that these encounters are not all the same as those identified using the pelt algorithm in Fig. 4 because some close encounters do not show up as prominently in the global evolution of the orbital separation (see equation 7), but they do appear in the individual evolution of the semimajor axis of the S-stars. For S6, S21, and S66 these data are presented in Fig. 5. In particular the earlier encounters at 13 and 408 yr, and the late encounter at 9046 yr do not appear from the analysis of the global change in semimajor axis through the pelt algorithm, but they are visible as jumps in the evolution of the individual S-stars.