The status of isochrony in the formation and evolution of self-gravitating systems

Abstract

In the potential theory, isochrony was introduced by Michel Hénon in 1959 to characterize astrophysical observations of some globular clusters. Today, Michel Hénon’s isochrone potential is mainly used for his integrable property in numerical simulations, but is generally not really known. In a recent paper (Simon-Petit, Perez & Duval 2018), we have presented new fundamental and theoretical results about isochrony that have particular importance in self-gravitating dynamics and that are detailed in this paper. In particular, new characterization of the isochrone state has been proposed, which is investigated in order to analyse the product of the fast relaxation of a self-gravitating system. The general paradigm consists in considering that this product is a lowered isothermal sphere (King Model). By a detailed numerical study, we show that this paradigm fails when the isochrone model succeeds in reproducing the quasi-equilibrium state obtained just after fast relaxation.

gravitation, methods: numerical, globular clusters: general, Galaxy: kinematics and dynamics, galaxies: general

1 INTRODUCTION

The formation and evolution of self-gravitating systems – globular clusters or any kind of galaxies – are a very complex physical process. It takes place on two distinct temporal scales:

• A rapid phase on a few dynamical times of the system during which the systems collapse under its own gravitational field (Binney & Tremaine (2008)). This fast time-scale will be labelled as T_fa ∼  T_dyn.

• A much longer phase, typically of the order of |$T_\mathrm{lo}\sim \frac{N}{\ln N}T_\mathrm{fa}$|⁠, where N is the number of particles in the system, during which the system slowly relaxes due to internal [two-body relaxation; cf. Heggie & Hut (2003), small potential fluctuations arising from their finite number of particles, e.g. Heyvaerts et al. (2017)] or external [tidal shocks, galactic potential, etc., e.g. Gnedin & Ostriker (1997)] processes.

In addition to this time ordering, there also exists a spatial one that distinguishes between:

• Singular systems, in which fast formation process occurs in a single phase in a non-fluctuating external gravitational potential and in which long-time evolution is not affected by merging. Typically, generic globular clusters and low surface brightness (LSB) galaxies are singular systems.

• Composite systems, in which formation process involves successive merging phases in a fluctuating gravitational field [e.g. hierarchical scenario for high surface brightness (HSB) galaxies, e.g. Kauffmann et al. 1999] or numerical simulations with clumps [e.g. van Albada (1982) for the initial study and Trenti, Bertin & van Albada (2005); Iguchi et al. (2005) for more recent works].

Using this nomenclature, a composite system can be viewed as the final result of a complex gravitational interaction between singular systems. The evolution of a composite system over T_lo is a difficult problem and must be addressed in a cosmological framework.

In this work, we are essentially interested in singular systems. Their fast formation process results in a spherical density profile with specific properties: a large homogeneous core containing half of the system mass and surrounded by an r⁻⁴ steep halo (Roy & Perez (2004)). This initial characteristic profile is altered by the long dynamical evolution: the core shrinks and the halo softens. The situation is rather different in composite systems. As a matter of fact, on the first hand, modern simulations with clumps produce highly concentrated core-halo structures, see Trenti et al. (2005). On the other hand, the cosmological hierarchical scenario confers a cusp to HSBs’ profiles.

The evolution of singular systems over T_lo has been studied intensively, especially for globular cluster systems. The relatively small sizes of such stellar clusters (about 1/1000 of the size of the host galaxy) make them stand as points in the galactic potential (Binney & Tremaine 2008). Additionally, their typical dynamical time is about one tenth of that of their hosts. Therefore they can be considered as singular systems.

After decades, the lowered isothermal King (1966) model and its generalizations [e.g. Binney & Tremaine (2008), p. 307–311] appear to be quite universal to describe any state of the evolution of globular clusters. Reasons for such a success in the description of the dynamical properties of globular clusters are physically understandable. After the theory of violent relaxation by Lynden-Bell (1967), one understood that this full process leads to an isothermal sphere with an infinite mass [see e.g. Binney & Tremaine (2008), p.305–307]. To get rid of this singularity, some physical improvements have been proposed. These refinements essentially consist in the introduction of some physical cut-off. When this cut-off is crudely introduced by hands in the isothermal model, it produces the original King (1966) model. It is originally justified by spatial limitation corresponding to the tidal cut-off imposed by the surrounding galaxy in which the globular cluster evolves. When the cut-off is included as a parameter of the statistical problem, it produces an isothermal sphere in a box (ISB), which is a trademark problem of the gravitational statistical physics. This fundamental problem has a centennial history and it finally appears to be closely related to the King model (see Chavanis, Lemou & Mehats 2015, and references therein). Studying the stability properties of this ISB in terms of the conditions imposed on the box, theoretical astrophysicists were gradually able to understand the dynamical history of globular clusters through the more and more refined analysis of the so-called gravothermal catastrophe (Lynden-Bell & Wood 1968; Katz 1978; Padmanabhan 1990; Chavanis 2002). The conclusion of this history is now well established concerning the mass density [see for example a very nice synthesis in section 7.5 of Binney & Tremaine (2008)]: after the initial fast relaxation process an isolated N-body system settles down in a spherical core-halo and quasi-equilibrium state. Considering the sole density profile, this state appears well described by a lowered isothermal sphere (King Model) and evolves adiabatically to a more concentrated lowered isothermal state under the influence of slow relaxation. When the stability of this lowered isothermal state – governed by its density contrast – is no more possible, the core of the system shrinks. For example, in our galaxy, about 20 per cent of the globular clusters possess the steep cusp in their surface-brightness profile (e.g. Djorgovski & King 1986) that is predicted by models of post-collapse evolution. The complementary set is not yet collapsed and possess a core-halo structure with a slope between −4 and −2 in the halo, fixed by the level of relaxation of the system. When the system is assumed to have a non-constant mass function, this mechanism is reinforced by mass segregation: the most massive objects become concentrated in the high-density central core.

Since the ‘King-fit’ scenario is globally accepted for the mass density (at least for globular clusters in our Galaxy¹), the interest of the community progressively moved to the evolution of the globular cluster’s mass function (see Forbes et al. 2018, for the more recent review). This is a much more complicated and tricky problem, which can be attacked only by observations or simulations. This late interest lets the theoretical approach of this problem at rest for a decade.

The three parameters of the King model essentially control the relative size of its core and the slope of its halo. Monitoring the evolution of globular clusters using the King model seems to be a good idea, but several questions remain posed. What are the physical characteristics of the initial steep core-halo r⁻⁴ structure that is produced by the fast relaxation? In what sense is this initial state of the slow relaxation a lowered isothermal sphere?

Quite recently, a fine analysis by Zocchi, Bertin & Varri (2012) points out that although King models usually offer a good representation of the observed photometric profiles, they often lead to less satisfactory fits to the kinematic profiles, independently of the relaxation condition of the systems.

In his pioneering paper, Michel Hénon (1959) proposed that globular clusters can be isochrone. This proposition was based on the fact that, in the late of the fifties, the observed globular clusters’ mass density distributions looked like the one of the isochrone model. The refinement of the observations has actually revealed a wider diversity! First of all, the mass density of globular clusters changes during their dynamical evolution. Although they are characterized by a spherical core-halo structure almost all along their life, the size of their constant density core and the slope of their power law surrounding halo actually change during their evolution. When they are dynamically young, observed or computed globular clusters are mainly characterized by large cores surrounded by steep r⁻⁴ haloes.

The old idea from Michel Hénon to associate globular clusters to the isochrone model has then been progressively forgotten because of the King model’s ability to reproduce their mass density distribution. It is probably why nobody has verified whether the result of the fast relaxation is truly isochrone, i.e. in a kinematic sense by checking orbital properties of its components.

Excepted an interesting work by Voglis (1994), who remark a quasi-isochrone property in the result of very particular simulations of spherical cold collapse in a cosmological context, no work was directly addressed to test the isochrone status of some stage of the evolution of some kind of self-gravitating system.

Recently Simon-Petit et al. (2018) reopened and widely reconsider the isochrony problem in the steps of Michel Hénon. In this fully detailed paper, the essence of isochrony is presented from mathematical concepts of group theory to the astrophysics of the generalized Kepler’s third law, passing through the new isochrone relativity theory. The aim of this paper is twofold. First, we want to summarize this huge paper in order to highlight its main results for practical uses in the context of galactic dynamics. Secondly, we propose to implement some of these new concepts to address the old idea from Michel Hénon. When translated in our nomenclature, this idea becomes: Could the result of the fast relaxation of singular systems be isochrone?

The paper is arranged as follows. In Section 2 we recall the main results of the general theory of isochrony recently obtained by Simon-Petit et al. (2018). We detail in particular the generalized version of Kepler’s third law, which we will use as an isochrony test in this paper. Section 3 is devoted to the numerical experiment, which we propose in order to check the isochrone character of the product of fast relaxation. We then propose some analysis of the numerical results and conclude.

2 ISOCHRONY IN 3D CENTRAL POTENTIALS

Let us recall the definition of isochrony for spherically symmetric systems introduced in Hénon (1959) and briefly summarize the recent paper Simon-Petit et al. (2018). The latter is a complete and detailed paper about isochrony that contains a new approach to this problem and all proofs and details of some properties we use in this paper. In particular, (i) it characterizes all continuous isochrone potentials, (ii) explains why the essence of isochrony is Keplerian, (iii) proposes the so-called isochrone special relativity theory, and (iv) provides useful applications for our purpose like a generalization of the third Kepler law.

Any test particle with a position-bounded and non-collapsing trajectory in a spherical self-gravitating system has its radial distance |$r(t) = \left|\mathbf {r(t)}\right|$| that oscillates between its value at pericentres and apocentres: r_p = min _tr(t) and r_a = max _tr(t). Its radial distance is governed by the ordinary differential equation of motion

\begin{eqnarray*} \frac{1}{2} \left(\frac{\text{d}r}{\text{d}t} \right)^2 + \frac{\Lambda ^2}{2r^2} - \left(\xi - \psi \right) = 0 \,, \end{eqnarray*}

(1)

where ψ is the radial gravitational potential associated with the system, E = mξ the energy of the test particle of mass m, and L = mΛ the norm of its angular momentum. The oscillation is characterized by its period

\begin{eqnarray*} \tau (\xi ,\Lambda)={2\int _{r_{\mathrm{ p}}}^{r_{\mathrm{ a}}}}\frac{\mathrm{ d}r}{\ \sqrt{2\left[ \xi -\psi \left(r\right) \right] -\dfrac{\Lambda ^{2}}{r^{2}}}} \end{eqnarray*}

(2)

called radial period [see for example Binney & Tremaine (2008), p. 146]. A spherical system is isochrone if all its non-escaping and non-colliding particles admit radial periods that only depend on their energy, i.e. τ = τ(ξ). This concept was first introduced by Michel Hénon in a seminal paper in French [Hénon (1959); for an English translation see Binney (2016)].

In his paper, Michel Hénon exhibited three isochrone central potentials. Completing his analysis, there are actually four isochrone potential families. The evidence of the proof comes from the change of variables

\begin{eqnarray*} x=2 r^2 \ \text{ and } \ Y(x) = x \psi (x) \end{eqnarray*}

(3)

exploited by Michel Hénon. In these variables, Y is related to the potential ψ. One can further show that ψ is isochrone if and only if the graph of Y is a parabola. This parabola property geometrically characterizes an isochrone system. If the origin of the plane (x, y) is defined such that it coincides with the centre of symmetry of the spherical system (r = 0), then the parabola intersects the vertical axis (Oy) at least once. Increasing potentials are associated with positive mass particles and are contained in the convex branches of the parabolas, whereas the concave branches contain decreasing potentials (in parabolas with non-vertical symmetry axes).

A reduction of the isochrone potential search is possible through the action of geometrical transformations. In Michel Hénon’s variables, the radial equation of motion (1) writes

\begin{eqnarray*} \displaystyle \tfrac{1}{16}\left(\frac{\mathrm{ d}x}{\mathrm{ d}t}\right) ^{2}+\Lambda ^{2}=x\xi -Y\left(x\right) \text{.} \end{eqnarray*}

(4)

In this equation, ξ and Λ are the two free parameters of the problem. The whole set of isochrones is obtained spanning the allowed (ξ, Λ) space when the graph of Y(x) is a parabola in the (x, y) plane. This span is possible with two simple transformations: (i) adding a constant to Y, which corresponds to adding a constant to the angular momentum or (ii) adding a linear term in x, which corresponds to adding a constant to the energy. Both these two transformations do not affect the type of the potential described by Y. In terms of parabolas, the first transformation named J_{0, λ}|$\colon$| (x, y)↦(x, y + λ) is simply a vertical translation. The second transformation is a transvection J_{ϵ, 0}|$\colon$| (x, y)↦(x, y + ϵx). It swivels the parabolas by inclining its symmetry axis, with a slope between −∞ and +∞, and preserving the intersection of the parabola with the ordinate axis, as well as the abscissa of its vertical tangent, when it exists. The two transformations generate a subgroup

\begin{eqnarray*} \mathbb {A} = \lbrace J_{\epsilon ,\lambda }; \, (\epsilon ,\lambda) \in \mathbb {R}^2 \rbrace \end{eqnarray*}

(5)

of the affine group and reduces the search of isochrone potentials to four classes of parabolas and thus potentials. Each class is generated by one of the four following parabolas under the action of |$\mathbb {A}$|⁠:

• the parabola with a vertical tangent at the origin |$\left(y\propto \pm \sqrt{x}\right)$|⁠; it is associated with a Kepler potential, |$\psi _{\mathrm{ke}}\left(r\right) =-\dfrac{\mu }{r}$|⁠, with μ > 0 related to the central mass;

• the classic parabola (y ∝ x²); it is associated with a harmonic potential, |$\psi _{\mathrm{ha} }\left(r\right) =\frac{1}{2}\omega ^{2}r^{2}$|⁠, with ω real;

• a parabola that intersects (Oy) twice:

  • if it opens to the right, it is associated with a Hénon potential, |$\psi _{\mathrm{he}}\left(r\right) =-\dfrac{\mu }{b+\sqrt{b^{2}+r^{2}}}$|⁠, where b > 0;

  • if it opens to the left, it is associated with an isochrone potential, called bounded or pseudo-Hénon potential, |$\psi _{\mathrm{bo}}(r)=\dfrac{\mu }{b+\sqrt{b^{2}-r^{2}}}$| for r ∈ [0, b], where b > 0.

The physical properties of these potentials are discussed in Simon-Petit et al. (2018).

Therefore, the action of the affine subgroup |$\mathbb {A}$| divides the isochrone potentials, or equivalently the associated parabolas, into four families: Kepler, harmonic, Hénon, and bounded (potentials).

The previous geometrical approach classifies and completes the set of continuous isochrone potentials. However, it splits them into four disjoint families and does not help understanding the nature of isochrony. As a matter of fact, there exists an intrinsic link between all isochrone systems as shown below.

Using Michel Hénon’s variables, a radially periodic orbit r₀(t₀) in a potential ψ₀ can be designed by (ξ₀x₀, Y₀) and (4) with its integral of motion (ξ₀, Λ₀). We consider its image when transforming the affine coordinate system (⁠|$w$|₁ = ξx, |$w$|₂ = y) by |$\mathbf {w_0 } = \left[\xi _0 x_0,y_0\right]^\top \mapsto \mathbf {w_1 } = \left[\xi _1 x_1,y_1\right]^\top$|⁠. We assume Λ₀ = Λ₁ = Λ, i.e. the radial orbits share the same area law, and ξx − Y is conserved by the transformation as initiated by Lynden-Bell (2016). If Y₀ is isochrone, then Y₁ is also isochrone if and only if the transformation is linear: |$\mathbf {w}_{1}=B_{\alpha ,\beta }\left(\mathbf {w}_{0}\right)$| with

\begin{eqnarray*} B_{\alpha ,\beta }=\left[ \begin{array}{cc}\alpha & \beta \\ \alpha -1 & \beta +1 \end{array} \right] ,\ \left(\alpha ,\beta \right) \in \mathbb {R}^{2}. \end{eqnarray*}

(6)

This transformation generalizes the so-called Bohlin or Levi-Civita transformation,² linking all isochrone potentials and orbits together. When B_{α, β} is symmetric, it is named iBolst, for symmetric bohlin boost, and is written as B_γ, where γ = α + β ≠ 0. It is essentially the composition of a homothety with a hyperbolic rotation. The length of vectors is then modified by

\begin{eqnarray*} &&\left(\xi _1 x_1\right)^2 - y_1^2=\gamma \left[\left(\xi _0x_0\right)^2-y_0^2\right] . \end{eqnarray*}

(7)

Under the iBolst group action, the Keplerian potentials generate the isochrone potentials set. More generally, iBolsts provide a characterization of isochrone systems and orbits: an orbit is isochrone if and only if it is the iBolsted image of a Keplerian orbit.

Considering the canonical axis |$Ox=\mathbb {R}\mathbf {i}$| and |$Oy=\mathbb {R}\mathbf {j}$|⁠, a natural reference frame |$(O,\mathbf {u},\mathbf {v})$| is attached to each isochrone parabola. These two natural vectors, respectively, generate the symmetry axis of the parabola for |$\mathbf {u}$| and the tangent to the parabola at the origin O for |$\mathbf {v}$|⁠. For the Keplerian parabola, we have |$\mathbf {u}=\mathbf {i}$| and |$\mathbf {v}=\mathbf {j}$|⁠, and for the harmonic there is an inversion and one has |$\mathbf {u}=\mathbf {j}$| and |$\mathbf {v}=\mathbf {i}$|⁠. This inversion corresponds to the Bohlin–Levi-Civita transformation. Remarkably, it is always possible to map the reference frame of any isochrone parabola to its Keplerian equivalent with an iBolst, such that |$(\mathbf {u},\mathbf {v}) = B_\gamma (\mathbf {i},\mathbf {j})$|⁠, see Simon-Petit et al. (2018) for the general case. Consequently, isochrone orbits in |$(O,\mathbf {i},\mathbf {j})$| are Keplerian in their reference frames |$(O,\mathbf {u},\mathbf {v})$|⁠. Hence, isochrony appears as a relative Keplerian property: any isochrone system is Keplerian when seen in its reference frame, this is isochrone relativity (Simon-Petit et al. 2018).

The parallel with special relativity provides a direct comprehension of isochrone relativity. The isochrone law of dynamics writes in the same way in all reference frames defined by the iBolsts. The ‘isochrone interval’ ξx − Y is conserved and the energy and potential terms are linearly combined. As in special relativity where the nature of time, light, and space-like vectors is conserved by Lorentz transformations, particular-like cones are preserved by iBolsts and characterize aperiodic, radial, and radially periodic orbits in the isochrone relativity. The correspondence between radially periodic orbits through the iBolst relies on an appropriate choice of affine coordinates system. For any isochrone radially periodic orbit, one can then easily construct its Keplerian or harmonic associated periodic orbit, even though the time does not run identically for the corresponding particles.

Eventually, one of the most noticeable properties of Keplerian systems is the Kepler third law, which proportionally relates the squared radial period τ² of any particle to the cube of its semimajor axis a³. The isochrone relativity enables one to understand that this relation holds for any isochrone potentials after adapting the lengths.

As in the Kepler potential, characteristic semimajor axes a are defined in any isochrone potential by:

• |$a=\frac{1}{2}\left(r_{\mathrm{ a}}+r_{\mathrm{ p}}\right)$| in a Kepler potential;

• |$a=\left(\frac{1}{2}\right) ^{2/3}R$| in a homogeneous bowl of radius R;

• |$a=\frac{1}{2}\left(\sqrt{b^{2}+r_{\mathrm{ a}}^{2}}+\sqrt{b^{2}+r_{\mathrm{ p}}^{2}} \right)$| in a Hénon potential;

• |$a=\frac{1}{2}\left(\sqrt{b^{2}-r_{\mathrm{ a}}^{2}}+\sqrt{b^{2}-r_{\mathrm{ p}}^{2}} \right)$| in a bounded potential.

The Kepler third law is then generalized for any isochrone system:

\begin{eqnarray*} \tau ^2 = \frac{4\pi ^2}{\mu } a^3, \end{eqnarray*}

(8)

where μ is the potential parameter that appears in the definition of ψ_ke, ψ_bo, or ψ_he, and μ = GM = ω²R³ for the homogeneous bowl. This isochrone orbital property holds for each orbit; it is then a key microscopic criterion for the isochrone analysis of systems of a given macroscopic density profile.

3 NUMERICAL EXPERIMENTS

The objective of this section is twofold. First, we want to check the old idea from Michel Hénon, which is in our nomenclature: do the fast collapse of a singular system is an isochrone system. Secondly, we want to check if the relevant King model commonly adopted to describe the result of this fast collapse passes or not the isochrone test.

3.1 Description of the simulations

We have performed two accurate N −body simulations using the latest version of the Gadget-2 code Springel (2005). Only the treecode part of treePM algorithm is employed. These simulations contain N = 3.10⁴ equal mass m particles in a typical radius of the order of unity in the units of the simulations. These units are such that the gravitational constant G = 1 and the total mass M = Nm = 1. As shown by Roy & Perez (2004), this number of particles is sufficient for our purpose of capturing the physics of the problem with reasonable runtimes. These physical self-gravitating systems we have computed are a Hénon sphere and a King model.

The Hénon sphere (e.g. Hénon 1964) is particularly suited for our purpose. On the first hand, it is clearly the simplest idealization of what we have called a singular system. On the second hand, it is known to preserve well its spherical nature during the course of the dynamics even being simulated with a N-body technique (see for example Roy & Perez 2004). In this model, the initial phase-space distribution function is isotropic, spatially homogeneous, Gaussian distributed in velocity space and given by

\begin{eqnarray*} f_{\mathrm{ H}}\left(r,v,j\right) =\left\lbrace \begin{array}{lc}\rho _{0}\left(2\pi \sigma _{v}^{2}\right) ^{-3/2} \exp \left(-\frac{1}{2}\frac{v^{2}+j^{2}/r^{2}}{\sigma _{v}^{2}}\right) , & \text{ for }r\lt R;\nonumber \\ 0 & \text{ if }r\gt R. \end{array} \right. \nonumber \\ \end{eqnarray*}

(9)

Considering |$E_\mathrm{ k}=\sum _{i=1}^{N}\frac{1}{2}m\mathbf {v}_{i}^{2}$| and |$E_\mathrm{ p}=\sum _{i=1,j\gt i}^{N}-\frac{m^2}{\left|\mathbf {x}_{i}-\mathbf {x}_{j}\right|}$|⁠, respectively, the total kinetic and potential energy of the system, when the initial virial ratio κ = 2E_k/E_p is in the interval ]−1, 0[ the system collapses to an equilibrium state in a few dynamical times (e.g. Roy & Perez 2004). Provided that the initial state is not too cold (κ ∈ ]−1, −0.25[) the radial orbit instability (see Maréchal & Perez 2012, for a review of this fundamental process) does not occur and the equilibrium state is spherical. Its mass density is characterized by a core-halo structure. It is very well established that the core of this structure contains roughly half of the mass of the system and the halo is well approximated by a power law: ρ(r) ∝ r⁻⁴ [see e.g. Roy & Perez (2004) or Joyce, Marcos & Sylos Labini (2009)³]. For our purpose, we have studied a Hénon sphere with κ = −0.5, an initial size R = 2, which gives after collapse a typical size R₅₀ = 1 in the units of the simulation.

The King model is a stable spherical isotropic equilibrium state of the Vlasov–Poisson equation. It means that if there is no relaxation its distribution function does not change in time. This distribution function is given by

\begin{eqnarray*} f_{\mathrm{ K}}\left(E\right) =\left\lbrace \begin{array}{lc}\rho _{0}\left(2\pi \sigma _{\epsilon }^{2}\right) ^{-3/2} \exp \left(-\frac{E_{\ell }-E}{\sigma _{\epsilon }^{2}}\right) -1, & \text{ for }E\lt E_{\ell };\\ 0 & \text{ if }E\gt E_{\ell }. \end{array} \right.\nonumber\\ \end{eqnarray*}

(10)

The mean field potential ψ(r) associated with this distribution function is given by the Poisson equation

\begin{eqnarray*} \Delta \phi (r) \!=\! 4\pi m G \rho _0 \left\lbrace \sqrt{\frac{4\phi }{\pi \sigma _{\epsilon }^2}}\left(1+\frac{2\phi }{3\sigma _{\epsilon }^2}\right)+ \mathrm{e}^{\frac{\phi }{\sigma _{\epsilon }^2}}\mathrm{erf}\left(\sqrt{\phi }/\sigma _{\epsilon }\right) \right\rbrace\!, \nonumber \\ \end{eqnarray*}

(11)

where

\begin{eqnarray*} \phi (r)=E_{\ell }-m\psi (r)\,\,\mathrm{and}\,\,\mathrm{erf}(x)=\frac{2}{\sqrt{\pi }}\int _0^x \mathrm{e}^{-u^2}\mathrm{ d}u\, . \end{eqnarray*}

(12)

The King model has three free parameters, which are:

• the liberation energy E_ℓ, which is the cut-off in the energy space introduced to cure the infinite mass problem of the isothermal model;

• the depth of the potential well ψ(0) at the origin;

• the energy variance |$\sigma _{\epsilon }^2$|⁠.

These parameters are gathered in one by introducing |$W_{0}=\frac{\phi (0)}{\sigma _{\epsilon }^2}\gt 0$|⁠. As said before, King models are core-halo structures [see for example Binney & Tremaine (2008), p.308–309]: the size of the core is a non-trivial function of r_c = σ_ϵ(4πmGρ₀)^−1/2 while W₀ controls the steepness of the power law r^α in the halo of the mass density profile. The slope α varies from α ≃ −2 (the isothermal value) when W₀ is large, typically greater than 12, to α ≤ −5 when W₀ becomes smaller than 3. The large energy cut-off of the King model introduces a sharp cut-off in the mass density profile for large values of the radius. The slope of the halo we are speaking about concerns the region surrounding the core where the mass density is not less than typically ρ₀/10³. This value corresponds to observable and then measurable values of the projected luminosities profiles. For our purpose, we have studied a King model with W₀ = 9. This value corresponds to the best-fitting concentration parameter c = 2 proposed in fig. 1 of Djorgovski & King (1986), thanks to the correspondence fig. 4.9 p. 310 in Binney & Tremaine (2008). The typical size of our King model is R₅₀ ≃ 1.

Figure 1.

Evolution of the mean physical characteristics of the simulation: on the left-hand side, the Hénon sphere with κ = −0.5 at t = 0; on the right-hand side, the King model with W₀ = 9 at t = 0. Notice the different time-scales and the adapted y-values for each panel.

The parameters for Gadget runs were adapted to our purpose:

• The softening length of the gravitational force is set to |$\varepsilon =\left(\frac{4\pi }{3N}\right) ^{1/3}R_{50}$| with R₅₀ ≃ 1 for both simulations. This value corresponds to an estimation of the initial mean interparticle distance, it is a bit larger than the usual value for ε in standard simulations. Nevertheless, it is well adapted for our purpose for which we want to minimize the effect of two-body relaxation in order to study the properties of orbits in a frozen collisionless equilibrium gravitational potential. The obtained value is adequate for the softening of the force, according to the criterion proposed by Athanassoula et al. (2000). This value is also sufficient to solve the collapse problem when the initial virial ratio of the Hénon sphere is not too small and the corresponding collapse is not too fast, e.g. κ = −0.5 (see Roy & Perez 2004).

• In Gadget, each particle has its individual time-step bounded by |$\delta t=\min \left[ \ \delta t_{\max },\sqrt{2\eta \varepsilon /\left|\mathbf {a}\right|}\right]$|⁠, where |$\mathbf {a}$| is the acceleration of the particle and η is a control parameter. We choose η = 0.025 and δt_max = 0.01. This ensures an energy conservation of the order of 1 per cent for each run.

• The tolerance parameter controlling the accuracy of the relative cell-opening criterion [parameter designed by ErrTolForceAcc in the documentation of Gadget, see equation 18 of Springel (2005)] is set at α_F = 0.005.

The duration of each simulation is 300 equilibrium dynamical times for the Hénon sphere and 500 equilibrium dynamical times for the King model.

3.2 Orbit monitoring and simulation results

The evolution of the physical parameters of the simulations is presented; in fig. 1, the left-hand side concerns the Hénon sphere and the right-hand side the King model:

• R₉₀, R₅₀, and R₁₀, respectively, represent the radii containing 90, 50, and 10 per cent of the total mass of the system. They are plotted in the top two panels. The King model is a stable equilibrium state, so these three quantities remain constant during the dynamical evolution. The initial state of the Hénon sphere with κ = −0.5 suffers Jeans instability and collapses to a steady state in a few dynamical times. The typical size of each system stays of the order of unity.

• Computing the eigenvalues λ₁ > λ₂ > λ₃ of the inertia matrix of the system, the quantities a₁ = λ₁/λ₂ and a₂ = λ₃/λ₂ are usually called the axial ratios of the system. Both of these quantities are equal to +1 when the system is a sphere. We can see on the next two panels that our simulations remain spherical all over their dynamical evolution.

• The virial ratio κ is equal to −1 when the system is at equilibrium. We can see on the third line panels that the King sphere remains at equilibrium and that the Hénon sphere quickly joins such a state just after a fast relaxation initial phase.

• For such Hamiltonian systems, the total energy E_tot = E_p + E_k is conserved during the time evolution. We observe that in the King model this conservation is properly respected in a mean sense. For the Hénon sphere, this is the same after the warm collapse of the system.

Counting particles in concentric spherical shells, we have computed the radial mass density for our simulation every 10 dynamical times. As expected, the mass density of the King sphere does not evolve at all over the 500 dynamical times computed. As expected too, the Hénon sphere initially collapses in a few dynamical times and reaches a steady state characterized by a core-halo structure: the size of the constant mass density core is roughly the radius containing half of the total mass of the system; the slope of the power-law like surrounding halo is roughly −4. Once it is formed (t ∼ 10T_d), this core-halo structure does not evolve at all until the end of our computation at t = 300 T_d. In Fig. 2, we have plotted these mass densities on the same graph in order to compare them. A rough analysis does not reveal differences between the results of our warm collapse (H) and simulated King model with W₀ = 9 (K). In addition, we note that we can adjust the unique parameter b of the Hénon potential ψ_he defined in Section 2 to get an isochrone mass density (i) that resembles to the one of (K) or (H). In fact, the plots of Fig. 2 clearly show that the density analysis is not conclusive concerning the nature of the equilibrium we obtain after the fast relaxation: it could be fitted by either a King model or an isochrone. Hence, a precise kinetic analysis is needed to reach firm conclusions.

Figure 2.

Mass density comparison: H represents the Hénon sphere at t = 100T_d, K the computed King model at t = 200T_d, i and k are, respectively, the plots of the isochrone (b = 0.36) and the King (W₀ = 9) normalized theoretical mass densities.

The kinematic analysis we propose consists in a fine study of orbits. In each simulation, a subset of n = 200 randomly chosen particles was monitored: three quantities are archived at each time-step t, namely the position |$\mathbf {r}(t)$|⁠, the velocity |$\mathbf {v}(t)$|⁠, and the gravitational potential |$\psi \left(\mathbf {r},t\right)$| imposed by the N bodies of the system at the position |$\mathbf {r}$| at time t. We point out that we have computed this potential using the gadget-2 treecode algorithm, hence with the same softening parameter ε defined above. Analysing these quantities, we can determine, for each monitored particle and when it exists, the period τ of its radial distance and in all cases its total energy |$e(t)=\frac{1}{2}m\mathbf {v}^2(t)+m\psi \left(\mathbf {r},t\right)$|⁠.

The analysis of the period of particles is a tricky job. A priori, we deal with orbits of negative energy particles in a spherical self-gravitating system, and thus in a central potential. In such conditions, the orbit is planar and the radial distance |$r(t)=\left|\mathbf {r}(t)\right|$| is a periodic function of time, i.e. there exists |$\tau \in \mathbb {R}^+$| s.t. r(t + τ) = r(t). To find this τ from numerical data and then the values of r_a and r_p, the simplest way should consist in the use of the FFT algorithm; but as it is could be seen in Fig. 3, in a practical way this method is not precise as it should be expected for several reasons.

Figure 3.

Orbits of four particles extracted from the H simulation. On y-axes, the dotted lines indicate the typical lengths R₉₀, R₅₀, and R₁₀ of the whole system. The time is indicated in the simulation’s T_d units on the x-axes.

Depending on the particle properties, its orbit could be concentrated in the deep centre of the system (see the lower orbit on the right-hand side of Fig. 3), be spanning only the halo (see the upper orbit on the right-hand side of Fig. 3), be spanning all the system or all its core (see the left-hand side of Fig. 3), or be even more special. When the particle experiences the deep core, the two-body effects could influence its dynamics; although the radial distance is periodic, this function is modulated both in phase and in amplitude. The phase modulation due to the high density values is weaker when the orbit stays in the outskirts of the system where only large-scale oscillations of the potential modulate its amplitude. These effects affect both K and H simulations and introduce various and uncontrollable biases when we compute the period using FFT on the whole data. Instead, in order to get the right value of the period with the smallest significative error, we carry out a hand-made analysis of each orbit. We first check the planar property of the orbit: we determine the mean angular momentum, compute the orthogonal plane to this mean vector and reject orbits with an amplitude of azimuthal oscillation δ (see Fig. 4) around this plane less than 20°, such that δ ≥ 0.2r_a.

Figure 4.

3D representation of an orbit extracted from the K simulation: A view from above is presented in the left-hand panel while the side-on view is presented in the right-hand panel.

When the orbit is planar, we determine the coordinates |$(r_{\mathrm{ a},i},t_{\mathrm{ a},i})_{i=1,\cdots ,\mathcal {N}}$| of |$\mathcal {N}$| successive maxima of the function r(t). The value of the integer |$\mathcal {N}$| depends on each orbit considered. We set |$\mathcal {N}\ge \mathcal {N}_{\min }=5$| in order to compute the period at least over |$\mathcal {N}_{\min }$| oscillations. This minimal value allows us to initiate the computation of the period |$\tau =\frac{t_{\mathrm{ a},\mathcal {N}}-t_{\mathrm{ a},1}}{\mathcal {N}-1}$|⁠, the value of |$\mathcal {N}$| is incremented while |$\sigma _{\tau }/\tau \lt 5{{\ \rm per\ cent}}$|⁠, where σ_τ is the standard deviation of the sequence of instantaneous periods τ_i = t_{a, i + 1} − t_{a, i} for |$i=1,\cdots ,\mathcal {N}-1$|⁠. When this algorithm does not converge, the orbit is rejected as not periodic. When it gives a unique value, this result is compared with the data. When several values of the period are possible during several phases of the orbit, it is rejected as multiperiodic.⁴ When the period of the orbit is confirmed, we determine its apocentre |$r_\mathrm{ a}=\frac{1}{\mathcal {N}}\sum _{i=1}^{\mathcal {N}}r_{\mathrm{ a},i}$|⁠. Computing the sequence of the minima of r(t) in the interval |$\left[t_{\mathrm{ a},1},t_{\mathrm{ a},\mathcal {N}}\right]$|⁠, we obtain the pericentre of the orbit in the same mean sense. Computing the mean value of the energy E = 〈e〉_t on this same time interval, we get the mean energy of the particle. The standard deviations of e, τ, r_a, and r_p are used for the uncertainty analysis: the amplitudes of error bars in the plots are twice the standard deviation.

Using this manufactured but precise algorithm, we are able to extract with a good level of confidence the energy, period, apocentre, and pericentre of n_H = 155 orbits for the H simulation and n_K = 172 orbits for the K simulation among the n = 200 monitored for each one.

If each set of orbits is isochrone, it must fulfil the generalized Kepler third law: τ² × a⁻³ = cst, where a is the isochrone length defined in Section 2. We then achieve an analysis in the space |$\mathscr {H}_1=\left[\ln (a),\ln (\tau)\right]$|⁠. As the system is neither Keplerian, harmonic nor bounded, we guess it should be in a Hénon potential; hence |$a=\frac{1}{2}\left(\sqrt{r_\mathrm{ a}^2+b^2}+\sqrt{r_\mathrm{ p}^2+b^2}\right)$|⁠. In this formula, b is a macroscopic positive parameter common to all orbits, while r_a and r_p are microscopic ones, specific to each orbit. For a given value of b, we can plot the set [ln (a), ln (τ)] containing n_X points for each simulation X = K and X = H. We can then determine the weighted linear fit y = sln (a) + c (see Appendix  A) of these plots and determine the residue of this fit, namely

\begin{eqnarray*} \chi ^2_{b,\mathrm{X}}= \frac{1}{\ell ^2} \displaystyle \sum _{i=1}^{n_{\mathrm{X}}} \left\lbrace \ln (\tau _i)-\left[s\ln (a_i) +c \right]\right\rbrace ^2 \,\,\mathrm{with}\,\, \ell ^2= \displaystyle \sum _{i=1}^{n_{\mathrm{X}}} \ln (\tau _i)^2 . \nonumber \\ \end{eqnarray*}

(13)

The optimal value |$\tilde{b}_\mathrm{X}$| of b_X is obtained by a minimization algorithm applied to this residue computation. Explicitly, we compute the residue for discrete values of b in the interval |$\mathcal {B}=\left[0,b_{\max }\right]$|⁠, where |$b_{\max }\simeq R_{50{{\ \rm per\ cent}}}$| as roughly estimated from the isochrone model. We then choose

\begin{eqnarray*} \tilde{b}_{\mathrm{X}}=\min _{b\in \mathcal {B}}\chi ^2_{b,\mathrm{X}} . \end{eqnarray*}

(14)

The plots of |$\chi ^2_{b,\mathrm{K}}$| and |$\chi ^2_{b,\mathrm{H}}$| as functions of b are presented in the small boxes of Fig. 5.

Figure 5.

Isochrone test 1 for X = K (upper panel) and X = H (lower panel). In each panel, the small box represents the residue |$\chi _{b,\mathrm{X}}^2$| and the slope s_X of the weighted linear fit as functions of the isochrone parameter b. The residue plot is the dashed line with values on the left of the small box while the slope is the solid line with values on the right of the small box. The optimal value |$\tilde{b}_{\mathrm{X}}$| of b is identified as well as the corresponding optimal value of the slope |$\tilde{s}_{\mathrm{X}}$|⁠. The weighted best linear fit of the data {ln (a), ln (τ)} for |$b=\tilde{b}_{\mathrm{X}}$| is then plotted.

When the optimal value |$\tilde{b}_{\mathrm{X}}$| is found, it should correspond to the optimal fit of y = ln (τ) by a linear function of x = ln (a) for the considered simulation. The best slope |$\tilde{s}_{\mathrm{X}}$| of this best linear fit corresponds to the best power-law relation between a and τ; it should be 3/2 when the system is isochrone. We gather the statistics in Table 1. The uncertainty of s_X is evaluated by the dispersion parameter σ_b defined in Appendix  A. The amplitude of the error bar is equal to this dispersion parameter.

Table 1.

Statistical analysis for the sets [ln (a), ln (τ)] in the simulations K and H.

X	|$\tilde{b}$|	|$\chi _{b}^{2}\left(\tilde{b}\right)$|	|$\tilde{s}$|	|$\delta _{\tilde{s}}$|
H	3.70 × 10−1	2.87 × 10−2	1.50	8.23 × 10−2
K	1.49 × 10−1	1.07 × 10−2	1.13	4.06 × 10−2

From our orbits analysis, we can investigate another space, namely |$\mathscr {H}_2=\left[\ln (\tau),\ln (-E)\right]$|⁠. This one is more direct than the previous because it does not need another parameter as b to build the isochrone length a. For an isochrone model, the Kepler third law applies and gives τ² × (−E)³ = cst, which implies |$\ln (-E)=-\frac{2}{3}\ln (\tau)+\mathrm{ cst}$|⁠. As we have computed the mean energy E for each monitored orbit, we can additionally check if K or H should be isochrone in this sense. The results are plotted in Fig. 6.

Figure 6.

Isochrone test 2 for X = K (top panel) and X = H (bottom panel). The values of the slopes of the weighted linear fits are indicated near each line.

3.3 Analysis of the results

The analysis of our results is unambiguous. Both spaces |$\mathscr {H}_1$| and |$\mathscr {H}_2$| reveal the isochrone nature of the result of the H simulation and the non-isochrone nature of the K one. On the first hand, the regression of ln (τ) in ln (a) is perfectly linear for both simulations for the optimal value of b, but the value |$\tilde{s}_{\mathrm{H}}=1,50\pm 8,23. 10^{-2}$| is fully compatible with the isochrone one s_{1, iso} = 3/2 whereas |$\tilde{s}_{\mathrm{K}}=1,13 \pm 4,063 . 10^{-2}$| is definitely not. On the other hand, the analysis in |$\mathscr {H}_2$| reveals that there is no unique power-law relation between τ and E for orbits in a W₀ = 9 King model while this relation clearly exists with the right value |${s}_{2,\mathrm{iso}}^2=-2/3$| for the collapsed Hénon sphere with an initial virial ratio κ = −0.5.

Why is the result of fast relaxation isochrone? The answer to this question was certainly proposed 60 yr ago by Michel Hénon in his seminal paper Hénon (1959). During the mixing of the fast relaxation, there is a natural tendency for particles to move towards equipartition in the energy. This is a pillar of statistical mechanics. In this context, resonances enable energy exchanges between particles. If stars with the same radial period τ have different energies, they will exchange energy till they reach the same E. But the definition of isochrony is precisely that all stars with a given τ share the same E.

Why has this old result been progressively forgotten? The response to this question is probably twofold. First of all, due to the slow relaxation process the initial state obtained after the fast relaxation is progressively changed into a more and more concentrated core-halo system. This gravothermal evolution was studied and understood during the last 50 yr. During this process, the mass density of the systems changes and it loses its isochrone property. When we let this possibility occur, decreasing the value of the softening parameter by, at least, an order of magnitude (ε → ε/50), the simulation of the same Hénon sphere clearly shows this mass density evolution over the same duration (see Fig. 7). Our statistical study is no more possible in this context as the orbital periods evolve with the gravitational potential: when it is similar to a King model with W₀ = 9, it should have lost part of its isochrone property.

Figure 7.

Evolution of the density for a Hénon sphere with an initial virial ratio κ = −0.5 and a small softening parameter. The slope s corresponds to the best linear fit of the halo. The value of W₀ is obtained by fitting the halo with a King model.

Another simple reason of the isochrone oblivion is the efficiency of the King model: with three free parameters that control the slope of the halo, the size of the core, and the concentration of the system, this model is perfect to fit all the evolution of globular clusters evolution regarding the mass density or the luminosity profile. As long as nobody asks for detailed kinematical analysis, there is no reason to discredit the King model.

4 CONCLUSION AND PERSPECTIVES

Let us summarize the main results we have presented in this paper:

• The analysis of a self-gravitating system using only its mass density distribution (or luminosity profile) is ambiguous: several models can produce similar mass density profiles. In particular, the King model and its three free parameters can be used to produce a pretty good fit of the mass density profile of a globular cluster at any stage of its dynamical evolution.

• The system produced just after the standard⁵ fast relaxation process is a core-halo structure compatible with both a King or an isochrone mass density. However, when kinematic data are taken into account, the King model fails where the isochrone succeeds in reproducing the equilibrium state.

• The isochrone model is just an initial condition obtained after the formation process of the system. Under slow relaxation processes, the system loses its isochrone character as it is confirmed by density profiles.

Such a result highlights the isochrone model in a new perspective. More than an aesthetic model, useful because it distinguishes itself by its ability for producing analytic formulas for the actions and angles of its orbits; the isochrone is in fact a fundamental potential resulting from a homogeneous fast relaxation process.

This process is the one supposed to occur in the formation of isolated self-gravitating systems (most globular clusters, LSB galaxies). It is then not surprising that they are characterized by a core-halo structure density profile as long as they are not so affected by slow relaxation processes. When the formation process is hierarchical (e.g. HSB galaxies), the continuous merging process is combined with an inhomogeneous fast relaxation. The initial isochrone core-halo structure is no more observable as the core is unstable.⁶ This produces their cuspy profile. This could probably also explain the presence of supermassive black holes in the heart of such structures while they are not expected as a rule for globular clusters or LSB galaxies.

Although the isochrone state is explicitly revealed after the homogeneous fast relaxation process, we have no more physical explanations than the one proposed by Michel Hénon in the sixties (resonant coupling arguments). A special investigation on this subject is tricky because it requires to analyse orbits during a non-stationary phase. It deserves to be discussed in a future work.

ACKNOWLEDGEMENTS

This work is supported by the ‘IDI 2015’ project funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02. The authors thank the referee of the article for valuable comments and fruitful suggestions.

Footnotes

1

For other hosts the situation is not so clear, see e.g. Mackey & Gilmore (2003) for GC in the SMC.

2

In fact, this transformation was initiated by Newton in his Principia, formalized by Goursat (1889) and Darboux (1889) followed by Bohlin (1911) and finally Levi-Civita (1906).

3

The fig. 14 of this paper is particularly explicit about this result.

4

By several we mean more dispersed than |$5{{\ \rm per\ cent}}$| around the mean value. This possibility occurs when there is a strong two-body interaction in which the characteristics of the orbit (period, apocentre, pericentre) are modified.

5

By this precision we want to restrict this affirmation to systems that are initially not sufficiently cold to be influenced by the radial orbit instability.

6

In an inhomogeneous system, the critical value of the density contrast could be very low. The collapse of the core could then appear during the merging and virializing phases.

APPENDIX A: LINEAR FIT WITH WEIGHT

Given a sequence {x_i; y_i ± σ_i}_{i = 1, ⋅⋅⋅, n}, we have to compute values for c and s, which minimize the quantity

\begin{eqnarray*} d^{2}=\sum _{i=1}^{n}\left[ y_{i}-\left(c+sx_{i}\right) \right] ^{2}. \end{eqnarray*}

(A1)

The uncertainty on y_i, namely σ_i, allows us to define the weight |$w_{i}=\sigma _{i}^{-2}$| of the i-th value of y. The d²-minimization problem gives

\begin{eqnarray*} c=\frac{ 1 }{\Delta }\left(\displaystyle \sum _{i=1}^{n}w_{i}x_{i}^{2}\right) \left(\displaystyle \sum _{i=1}^{n}w_{i}y_{i}\right) - \frac{ 1 }{\Delta }\left(\displaystyle \sum _{i=1}^{n}w_{i}x_{i}\right) \left(\displaystyle \sum _{i=1}^{n}w_{i}x_{i}y_{i}\right) \end{eqnarray*}

(A2)

and

\begin{eqnarray*} s= \frac{1 }{\Delta }\left(\displaystyle \sum _{i=1}^{n}w_{i}\right) \left(\displaystyle \sum _{i=1}^{n}w_{i}x_{i}y_{i}\right) -\frac{1 }{\Delta }\left(\displaystyle \sum _{i=1}^{n}w_{i}x_{i}\right) \left(\displaystyle \sum _{i=1}^{n}w_{i}y_{i}\right) \,, \end{eqnarray*}

(A3)

where

\begin{eqnarray*} \Delta =\left(\sum _{i=1}^{n}w_{i}\right) \left(\sum _{i=1}^{n}w_{i}x_{i}^{2}\right) -\left(\sum _{i=1}^{n}w_{i}x_{i}\right) ^{2}. \end{eqnarray*}

(A4)

The uncertainties on c and s are given by

\begin{eqnarray*} \sigma _{c}=\sqrt{\frac{1}{\Delta }\displaystyle \sum _{i=1}^{n}w_{i}x_{i}^{2}}\text{ and }\sigma _{s}=\sqrt{\frac{1}{\Delta }\displaystyle \sum _{i=1}^{n}w_{i}}. \end{eqnarray*}

(A5)