{'fraction_non_alphanumeric': 0.035233824471492634, 'fraction_numerical': 0.016655989750160152, 'mean_word_length': 4.405966968566862, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Detecting Irregular Orbits in Gravitational N-body Simulations

We present a qualitative diagnostic based on the continuous wavelet transform, for the detection of irregular behavior in time series of particle simulations. We apply the method to three qualitatively different gravitational 3-body encounters. The intrinsic irregular behavior of these encounters is well reproduced by the presented method, and we show that the method accurately identifies the irregular regime in these encounters. We also provide an instantaneous quantification for the degree of irregularity in these simulations. Furthermore we demonstrate how the method can be used to analyze larger systems by applying it to simulations with 100-particles. It turns out that the number of stars on irregular orbits is systematically larger for clusters in which all stars have the same mass compared to a multimass system. The proposed method provides a quick and sufficiently accurate diagnostic for identifying stars on irregular orbits in large scale N-body simulations.


Introduction
Gravitational N-body systems are intrinsically chaotic, at least for N > 2 [1]. For N = 2 Newtonian systems are regular, but post Newtonian dynamics [2,3,4] can already reveal chaotic behaviour [5,6] by N = 2. Chaoticity in a self gravitating systems with N > 2 can reveal itself on a very short time scales, of the order of orbital periods. On the other hand, in some systems like the planetary system around the Sun, chaoticity only reveals itself on a time scale of billions of orbits [7,8]. In a stellar cluster which does not contain a dominant central mass, orbits can be chaotic on a much smaller time scale. In such an environment seemingly regular orbits can suddenly become highly irregular, to return later to regular again [9,10].
In this study we focus on the characterisation of orbits which show irregular behaviour on a short (dynamical) time scale, and not on a long (relaxation) time scale. Star clusters, with contain between about 100 and O(10 7 ) stars, are governed by microscopic few body interactions [11]. Analysis of such systems is often hampered by this internal microscopic physics and a qualitative indicator of chaotic behaviour could enormously assist in the understanding of large scale gravitational N-body simulations. In particular, since it is thought that stars on irregular orbits in self gravitating N-body systems have an important effect of the bulk properties of such systems [12].
We present a transparent diagnostic to qualify chaotic behaviour of individual trajectories in N-body systems. In addition, the method has some quantitative qualities. The application of this method ranges from 3-body interactions to simulations of entire star clusters and galaxies. For clarity, we define irregular orbits as orbits which are deterministic though sensitive to the initial conditions and which cannot be described as a sum of periodic motions.
The motion of a star in an N -body simulation can change from regular to irregular (and vice versa) due to gravitational interactions with other stars. Irregularity of an orbit then is a local quantity, and this forces us to deviate from the conventional methods based on Lyapunov [13] numbers such as explained in [14,15,16]. In addition, the calculation of Lyapunov exponents is costly, may require reruns and are therefore less suited for a direct diagnostic, whereas we are predominantly interested in a diagnostic that can be evaluated at runtime. The Lyapunov indicator, however, is well suited for discriminating between ordered and weak chaotic motion, like planetary systems, whereas we are in particular interested in very chaotic systems [14]. Others methods for detecting chaos in Hamiltonian systems, such as SALI [17,18], Fourier Transform [19,20,21], Poincaré Section [22], the zero-one method [23] and the geometric indicator [24] are also less suited for analysing gravitational N-body simulations during run time, since they also are not practical in providing an instantaneous quantification of the chaoticity of the system we are interested in here.
In § 2 we discuss the basics of the Continuous Wavelet Transform (CWT), which we present in a package called CWaT. We then proceed by applying this method to several N -body simulations in § 3, to conclude in § 4.

The continuous wavelet transform
The continuous wavelet transform [25] (CWT) gives the time-frequency representation of a time series f (t) by fitting a wavelet Ψ to it at subsequent points in time t. We have used the Morlet-Grossmann wavelet [26], which is described by Here σ is a measure of the spread in time and ω 0 the base frequency of the wavelet. The wavelet is a periodic sinusoidal signal enveloped by a Gaussian, from which we can construct a family of wavelets using a scaling factor α and a time shift β We define the CWT as Here Ψ is the complex conjugate of Ψ and f (t) a time series. Fitting the wavelet to f (t) at discrete moments in time (β ∈ [0, t]) yields T (α, t). We express T (α, t) in terms of frequency ω rather than scale α (since we are more accustomed to the former) using the relation This leads to a time-frequency representation T (ω, t). For an extensive discussion see [27].

Analysing the CWT
The main features of the time series can be described by its instantaneous frequencies ω(t). These can be extracted from T (ω, t) by selecting its local maxima. We only consider all local maxima larger than a threshold τ = 0.95 of the global maximum. In this we take a different approach than [28,29] who explicitly consider the global maximum. Using a lower threshold enables us to filter most of the noise from the CWT without losing the relevant information on which the χ can be determined. We also tried τ = 0.93 and τ = 0.97, which gave very similar results in the resulting value of χ.
We find the extrema by simulated annealing [30,31], following [32]. Connecting the local maxima results in a number of curves in the time-frequency domain, which we call ridges. The time series is subsequently analysed by studying the ridges, enabling a qualitative analysis of the behaviour of the time series. This is considerably more practical than using the CWT directly, an example is presented in fig. 2.
A periodic time series (i.e. with a constant frequency) results in a uninterrupted horizontal ridge indicating the fundamental frequency P and possibly overtones spaced at periods P/n where n is integer. A quasi-periodic timeseries can be described by a sum of periodic time series, each with its own ridge at a specific frequency. Such a time series is then represented by multiple horizontal ridges. Irregular time series on the other hand lead to curved ridges of limited length [33].

A quantitative measure of chaos
To bring the chaos detecting qualities of the method one step further we quantify the number of ridges and their curvature in the irregularity indicator χ. We define χ at any moment in time by: Here n(t) is the number of ridges at time t and s i a measure for the slope of ridge i. Using a least squares fit we calculate the slope of the ridge. Because χ would be dominated by large values of the slope we use an empirical cutoff. The slope is normalised to unity to assure that s i ∈ [0, 1].
The ridges of periodic and quasi-periodic time series are represented by horizontal ridges, they imply χ(t) ≃ 0 (see § 2.2). We note however that the proposed indicator (see Eq. 5), cannot give a yes/no definition about the occurrence of local exponential instability, but it is capable of giving a measure of the degree of instability.

Application
To demonstrate the effectiveness of CWaT we apply it to three gravitational 3-body interactions. Such interactions are irregular [34], with some regular characteristics (see [35,36] for an interesting earlier analysis of this sort).
In each of the three following examples we inject a single object (star) in a system consisting of two objects (stars) which are bound in a binary. A time series is constructed by projecting the orbit of each of the stars on three orthogonal axes (see also [19]). The stellar orbits are computed using the starlab software environment [37,38] which is publicly available from http://www.ids.ias.edu/~starlab. This N-body integration package calculates the orbits of the stars with individual time steps [39] using a fourthorder Hermite predictor-corrector scheme [40]. A regular time series is constructed by Hermite interpolation at evenly spaced time intervals. The time step adopted for CWT is 1/1024 th of the total integration time. Further analysis is carried out with the S-WAVE software package [32], to ultimately calculate the CWT and the ridges of the time series.
In the following two § § we discuss the chaotic behaviour in one non-resonant preserving encounter and one resonant encounter. In the first encounter the incoming star escapes after a close encounter with a binary. In the following example, in § 3.2 (see [41] for more examples, and even some very peculiar orbits in [42]), the incoming star has multiple close interactions with the binary, we call such interactions 'resonant'. The example shown is a democratic resonance which is characterised by all three stars interacting strongly with multiple close passages; each star generally interacts with comparable frequency with any of the other stars.
We use the following convention: star #1 and #2 form the initial binary which is encountered by star #3. The time series give the position of the star with respect to the centre of mass of the 3-body system. The units are given in standard N -body units [43].

Irregular motion in a non-resonant gravitational interaction
In the first example we select an interaction between a binary, consisting of star #1 and #2 to interact with a single incoming star, #3. After interacting with the binary, star #3 escapes without disrupting the binary. The incomer, however, may change the binding energy of the binary. Precise initial conditions are presented in Tab. .1.
In fig. . 1 we present a schematic diagram of the interaction. In this example the binary has a close encounter with the incomer at t ≃ 150. After the encounter the incoming star is deflected along a hyperbolic path.
The analysis using CWaT of one of the nine (three stars and three axes) time series is given in fig. 2 The ridge at P ≃ 6.3 in the third panel indicates the period of the binary system, which after t ≃ 150 has reduced to P ≃ 6.0. This ridge is joined by another on the interval t = [150, 330] at a shorter period of P ≃ 3. The lower period ridge is an overtone caused by an increase of the eccentricity of the binary from e = 0 to e ≃ 0.38.
As we discussed in Section 2.3, the bottom panel of fig. 2, shows that χ ≃ 0 during the entire interaction in accordance with the periodic behaviour of the system. The short fly-by interaction has not let to a qualitative indication of irregular behaviour in the system.
A similar analysis can be performed for any star and for each projection. The result will be quite similar, supporting the argument that the presented method is independent of the selection of a fundamental plane.

Irregular motion in a democratic resonant interaction
In a democratic resonance encounter the incoming star interacts closely with both components of the initial binary, until one of the stars eventually escapes. The initial conditions used for this example are presented in Tab. 2 and a schematic representation of the interaction is shown in fig. 3, and the analyse using CWaT is presented in fig. 4.
This interaction is considerably more complicated that the previous two ex- amples. In fig. 3 and 4 we see how the incomer, star #3, is captured by the binary (consisting of star #1 and #2) at t ≃ 400 and interacts several times strongly with the other two stars before being ejected again at t ≃ 1700. The initial binary has a period P ≃ 6.3, but after t ≃ 400 the CWT becomes very The strong democratic resonance interaction in this § provides an excellent example to compare the CWaT method with Poincaré sections, which we will do in the subsequent §.

Poincaré section of a democratic resonant interaction
Poincaré sections have, like the CWaT diagnostic, the advantage that a qualitative impression of the irregularity in time series can be obtained instantaneously, id est; at run time. The main disadvantage of Poincaré sections is the   6. Evolution of the radii which contains 5% (bottom thick solid curves), 50% (middle solid curve) and 75% (thin solid curve) of the mass for simulation model EqM (all masses are equal). For model S, with a Salpeter initial mass distribution for the individual particles, the same data is presented as the thick dotted curve (5% Lagrangian radius) and the two thin dotted curves representing the 50% and the 75% Lagrangian radii.

Irregular motion in a 100-body system
After having confirmed that the CWaT method enables a reliable detection of irregular orbits in small N-body systems it is time to focus on the larger systems. Here we will apply CWaT to self gravitating N-body systems with 100 particles. In particular we report here the result for several of such simulations. The initial realizations of the particle clusters ware generated as follows.
We make a distinction between two basic models, model EqM and model S. Model EqM was constructed by distributing 100 objects of the same mass in a Plummer [45] Spherical density distribution. The velocity of all point masses are subsequently scaled such that the entire system is in virial equilibrium, and to N-body units [43]. We run the system for 100 N-body time units, which corresponds to about 35 dynamical time units.
In fig. 6 we present the evolution of several fundamental radii for simulation model EqM which contain 5%, 50% and 75% of the entire mass of the system. Although it is not our intention to discuss the dynamical evolution of this system it is worth mentioning that the core experiences a gravothermal runaway [46] at about t ≃ 35, noticeable as the deepest point reached by the lower solid curve (5% Lagrangian radius), whereas the outer parts of the cluster continues to expand, mainly due to binary heating.
A further analysis of simulation EqM using CWaT, is presented in fig. 7. For each star we compute the CWT for the duration of the simulations. In fig. 7 we then show the fraction of stars (f EqM ) which has a χ > 0.1 (top curve), 0.2, 0.4 and 0.6 (bottom curve). The various curves are rather smooth and flat. Some caution is well taken here, though, as near the edges (within the first and last ∼ 20 time units) the CWaT becomes unreliable. The range between t = 20 and t = 80, where the method is accurate, we notice that the fraction of stars with any degree of chaoticity is roughly constant. It is interesting to note that > ∼ 20% of the stars have χ > 0.6 and that 30-40% of the stars have a χ < 0.2. In § 3.2, by the way, we demonstrated that a value of χ > ∼ 0.6 is a strong indication that the underlying orbits change on a very short time scale. The core collapse around t ≃ 35 is not clearly noticeable. Core collapse has apparently no direct influence in the number of irregular orbits in the equal-mass N-body system.
We run a second simulation with the same number of particles, distributed initially using the same 3 dimensional density profile, in which we assign randomly masses to each of the stars from a power-law distribution. We opted for a [47] mass function which has the form of N(m) ∝ m −2.35 . Here N(m) is the number of stars formed per unit mass range. We apply this mass function over a range of two orders of magnitude, i.e, the most massive star is a hundred times more massive than the least massive stars. We call this model S.
We perform several simulations using this mass function but with new initial realizations of the system. The run-to-run variations can be quite large, mainly due to the small number of stars. A consistent trend, however, is evident.
The dotted curves in fig. 6 present the evolution of the fundamental radii containing 5%, 50% and 75% of the mass of simulation model S. The evolution of the structure of model S is quite different than for the simulation with equal masses (model EqM). In the former case core collapse occurs much earlier (at about t = 14), which is consistent with simulations of [48]. However, model S also shows a second core collapse near t ∼ 63. Noticeable is the much more dramatic expansion of the outer parts of the cluster, which is driven by the heating of the cluster core by close bound pairs which form during the core collapse (exempli gratia see [42,49]).
We compute the CWT for each star in model S, for which we show the results of two simulations. For a better comparison we present the ratio q S/EqM ≡ f S /f EqM in fig. 8. Here f S and f EqM are computed for the same value of χ, being either 0.1, 0.2, 0.4 or 0.6. Here the bottom curve represents the ratio in the number of stars for χ > 0.6, the two curves give q S/EqM for χ > 0.1. The ratio q S/EqM < 1 for all values of χ > 0.1 and over the entire time span of the simulations.
The result of two simulations with a mass spectrum (solid and dotted curves, in fig. 8 show the same behaviour. The structure of the curves for χ > 0.4 and χ > 0.6 are also quite interesting. They continue to drop at a constant rate until model S reaches a second core collapse around t = 60-70. The comparison carried out here, though not the main objective of this paper, is quite interesting from an astrophysical perspective. As the multi-mass simulation systematically appears to have fewer stars on irregular orbits. This deficiency can be as small as q = 0.05 during core collapse (or there about). We are currently extending this analysis to larger N-body systems on order to study core collapse in large N-body simulations.

Conclusions
We present a qualitative diagnostic to analyze irregular behaviour in particle based gravitational N-body simulations. The method, dubbed CWaT, is based on continuous wavelet transforms and provides a direct and instantaneous diagnostic for irregular behaviour.
In our analysis we demonstrate that CWaT is suitable for analyzing self gravitating N-body simulations, as it does not require a long time series for the analysis and it provides an instantaneous indicator for the degree of chaos, allowing a quantification of the number of irregular orbits at any moment in time in a large N-body simulation.
We applied CWaT successfully to analyse gravitational 3-body interactions. The results of the CWaT method for 3-body systems compares well with Poincaré sections. Here CWaT provides a qualification of the irregularity of the orbit, whereas for Poincaré sections such an analysis is technically much harder to perform.
Eventually we applied the CWaT method to several N-body systems with 100 particles in which all stars had the same mass, and another set of simulations in which a mass spectrum was adopted. The multi mass system systematically produced fewer stars on irregular orbits, compared to the equal mass system. In addition, we noticed that during a deep core collapse the number of irregular orbits dropped even further, whereas in the relatively low density inter-core collapse stages the number of irregular orbits tends to increases again.

Appendix: Initial conditions
The calculations for testing the method are performed with the kira integrator which is part of the starlab gravitational N-body environment [37] (see also http://www.manybody.org/manybody/starlab.html). The initial conditions used for the 3-body encounters used as examples in sections 3.1 and 3.2 are presented in the tables 1 and 2, respectively. The numbers in the tables are in standard N-body units [43]. Exact initial conditions for each of the three encounters studied is available via http://modesta.science.uva.nl. The CWaT source code is also upon request. Table .1 Initial conditions of the non-resonant preserving 3-body interaction. All numbers are in N -body units The initial orbital parameters are: m 1 = m 2 = m 3 = 0.5, a = 1, e = 0, b ≃ 1.34, v ∞ ≃ 1.24v c where the critical velocity v c corresponds to a kinetic energy equal to the binding energy of the binary. [43].