Secular dynamics of long-range interacting particles on a sphere in the axisymmetric limit

# Secular dynamics of long-range interacting particles on a sphere in the axisymmetric limit

Jean-Baptiste Fouvry111Hubble Fellow Institute for Advanced Study, Princeton, NJ, 08540, USA    Ben Bar-Or Institute for Advanced Study, Princeton, NJ, 08540, USA    Pierre-Henri Chavanis Laboratoire de Physique Théorique, Université de Toulouse, CNRS, UPS, France
###### Abstract

We investigate the secular dynamics of long-range interacting particles moving on a sphere, in the limit of an axisymmetric mean field potential. We show that this system can be described by the general kinetic equation, the inhomogeneous Balescu–Lenard equation. We use this approach to compute long-term diffusion coefficients, that are compared with direct simulations. Finally, we show how the scaling of the system’s relaxation rate with the number of particles fundamentally depends on the underlying frequency profile. This clarifies why systems with a monotonic profile undergo a kinetic blocking and cannot relax as a whole under resonant effects. Because of its general form, this framework can describe the dynamics of globally coupled classical Heisenberg spins, long-range couplings in liquid crystals, or the orbital inclination evolution of stars in nearly Keplerian systems.

## I Introduction

Long-range interacting systems generically undergo an evolution in two stages. First, a fast (collisionless) violent relaxation (Lynden-Bell, 1967) during which the system reaches a quasistationary state (a steady state of the collisionless Boltzmann equation) and the system is dynamically frozen under the mean field dynamics. Then, as a consequence of the finite number of particles, the system undergoes a slow (collisional) relaxation that drives it towards thermodynamical equilibrium. This second stage is generically described by the Balescu–Lenard (BL) equation (Balescu, 1960; Lenard, 1960), recently generalized to inhomogeneous systems (Heyvaerts, 2010; Chavanis, 2012). These formalisms can account simultaneously for inhomogeneity (i.e. non-trivial orbital structures), collective effects (i.e., spontaneous amplification of perturbations) and non-local resonant couplings.

In this letter, we focus our attention on one such long-range interacting system, namely the problem of long-range coupled particles evolving on a sphere. Because of its general form, this system is of relevance in various physical setups ranging from spin dynamics to stellar systems (see Section II). Here, we show how in the axisymmetric limit, the generic methods of inhomogeneous kinetic theory can be applied, and accordingly derive the associated kinetic equation. In addition to allowing for quantitative predictions of the system’s diffusion coefficients, we clarify how this theory predicts the dependence of the system’s relaxation rate with the number of particles, and the important role played by the frequency profile in that respect.

The paper is organized as follows. In Section II, we present the considered model. Placing ourselves within the axisymmetric limit, we derive in Section III the appropriate Balescu–Lenard equation describing the long-term evolution of that system. In Section IV, we present applications of this formalism to recover the system’s diffusion coefficients as well the scaling of the relaxation rate with the number of particles. Finally, we conclude in Section V.

## Ii The model

We consider a set of particles evolving on a sphere of unit radius, and denote the spherical coordinates with . To any location on the sphere, we associate a normal vector . The specific Hamiltonian of the system is

 H=μ∑i

where is the individual mass of the particles, is the pairwise interaction, and is an imposed external potential. The pairwise interaction is developed in Legendre Polynomials as

 U(Li⋅Lj) =−∑ℓαℓPℓ(Li⋅Lj) (2) =−∑ℓ,mαℓbℓYmℓ(Li)Ym∗ℓ(Lj);bℓ=4π2ℓ+1,

where we used the addition theorem, and introduced the spherical harmonics , where is the associated Legendre functions (Arfken et al., 2013), and . The spherical harmonics are normalized as , with the unit volume . The canonical coordinates associated with this two-dimensional phase space are , and the equations of motion for particle are and . We recast these equations as

 dLidt=∑ℓ,mMmℓ(t)Xmℓ(Li)+Xext(Li), (3)

where are the vector spherical harmonics,

 Mmℓ(t)=μαℓbℓ∑jYm∗ℓ(Lj(t)) (4)

are the system’s instantaneous magnetisations and captures the contribution from the external potential.

Equation (3) is the exact evolution equation of this problem. The Hamiltonian from Eq. (1) encompasses a wide class of long-range interacting systems: (i) describes globally coupled classical Heisenberg spins (Gupta and Mukamel, 2011; Barré and Gupta, 2014), (ii) is the Maier-Saupe model for liquid crystals (Maier and Saupe, 1958; Roupas et al., 2017), (iii) captures the process of vector resonant relaxation in galactic nuclei (Kocsis and Tremaine, 2015; Szölgyén and Kocsis, 2018; Takács and Kocsis, 2018) (up to additional conserved quantities).

In the coming section, we consider the most general setup, but place ourselves within the axisymmetric limit, i.e. where the mean field Hamiltonian is invariant w.r.t. . We show how the general kinetic theory of long-range interacting systems (see Appendix A) can straightforwardly be applied to this regime, and accordingly derive the associated evolution equation.222We make the correspondence with Appendix A by noting that plays the role of , the role of the angle , and the role of the action .

## Iii The Balescu–Lenard equation

Let us assume that the system is characterized by a mean distribution function (DF), , normalized so that , with the total mass of the system. Following Eq. (1), the mean specific Hamiltonian of a particle in that system reads

 H0(L) =∫dL′U(L⋅L′)F(L′)+Uext(L) =∑ℓhℓPℓ(u)+Uext(u), (5)

where in the second line, we assumed that the system’s DF and the external potential are axisymmetric, i.e., and , and introduced the coefficients . The associated orbital frequency naturally follows from Eq. (5). For axisymmetric configurations, we have . Therefore, the Poisson bracket satisfies , i.e., any axisymmetric DF is a steady state for the mean field dynamics. In addition, the mean Hamiltonian is integrable, as the action is conserved along the mean motion, while the associated angle , evolves linearly in time with the frequency .

Investigating the long-term evolution of such a quasi-stationary steady amounts to investigating the slow distortion of the system’s mean DF, (assumed to remain linearly stable and axisymmetric throughout its evolution). Following the general kinetic theory of long-range interacting integrable systems (briefly reproduced in Appendix A), deriving the kinetic equation for is immediate. One only needs to proceed by analogies as we detail below.

The interaction potential can be written under the separable form , where the potential basis elements are

 ψ(p)(L)=CℓpYmpℓp(L),Cℓ=√αℓbℓ. (6)

Fourier transform w.r.t. the angle reads

 ψ(p)k(u)=∫dϕ2πe−ikϕCℓpYmpℓp(u,ϕ)=δmpkcmpℓp(u), (7)

with the coefficient . Injected in Eq. (26), the system’s response matrix becomes

 ˆMpq(ω)=2πδmqmp∫dump∂F/∂uω−mpΩ(u)cmp∗ℓp(u)cmqℓq(u). (8)

A system is then said to be linearly unstable if there exists a complex frequency (with ), for which admits an eigenvalue equal to . In that case, the system supports an unstable mode of pattern speed , and growth rate  (see Section 5.3 in Binney and Tremaine, 2008). In present context, Eq. (8) generalizes the stability criteria put forward in Gupta and Mukamel (2011); Barré and Gupta (2014) (see Appendix B).

Following Eq. (25), the system’s dressed susceptibility coefficients read

 (9)

where

 [Ik]ℓℓ′=δℓ′ℓ;[ˆMk(ω)]ℓℓ′=ˆM[ℓ,k],[ℓ′,k](ω). (10)

Assuming that the system is linearly stable, and that the frequency profile is non-degenerate (i.e.  only in isolated points), the long-term evolution of this axisymmetric system is characterized by the inhomogeneous BL equation (see Eq. (20)), that reads here

 ∂F∂t=2π2μ∂∂u ×(∂∂u−∂∂u′)F(u)F(u′)], (11)

where we introduced the total dressed susceptibility coefficients as

 ∣∣ψdtot(u,u′,ω)∣∣2=2∑k≥1k∣∣ψdkk(u,u′,kω)∣∣2. (12)

In Eq. (11), we emphasize the absence of a sum on resonance vectors, owing to the Kronecker symbol in Eq. (9). Collective effects can be switched off by imposing (i.e. replacing the dressed susceptibility coefficients, , by their bare analogs, , see Eq. (22)), which leads to the inhomogeneous Landau equation (Chavanis, 2013a). Finally, we recall that Eq. (11) can be rewritten as a Fokker–Planck equation

 ∂F∂t=−∂∂u[D1(u)F(u)]+12∂2∂u2[D2(u)F(u)], (13)

with the first- and second-order diffusion coefficients

 D2(u) =(2π)2μ∫du′|ψdtot|2δD(Ω(u)−Ω(u′))F(u′), (14) D1(u) =12∂D2∂u+2π2μ∫du′|ψdtot|2δD(Ω(u)−Ω(u′))∂F∂u′.

In practice, for a given value of , one can carry out the integral in Eq. (11) by finding the resonant actions satisfying , which allows for the replacement

Because of its prefactor , the BL equation describes the system’s long-term self-consistent evolution computed at first-order in the effects, accounting for the amplification by collective effects. Here, it is important to note that (i) the orbital space is one dimensional (cf. the one dimensional integral in Eq. (11)), (ii) the symmetry of the interaction imposes resonances (cf. the absence of sums over in Eq. (11)). As a consequence, if the system’s mean frequency profile, , is monotonic, the resonance condition only allows for local resonances, i.e. , leading to zero flux and . In that case, the system cannot relax under effects (Barré and Gupta, 2014; Chavanis, 2013b; Rocha Filho et al., 2014; Lourenço and Rocha Filho, 2015). It undergoes a so-called kinetic blocking (Chavanis and Lemou, 2007), and can only relax under weaker finite- effects associated with higher-order correlations. Still, even if the flux vanishes, the diffusion coefficients and remain non-zero. Conversely, for a non-monotonic frequency profile, non-local resonances, , are allowed, the flux is non-zero, and the system can relax at the order . We illustrate these various effects in the coming section. Finally, we emphasize that the Boltzmann DF is always a stationary solution of the BL equation. Yet, the fact that for kinetically blocked systems any axisymmetric DF is a stationary solution of the BL equation does not imply that these states remain stationary when higher order correlation effects are accounted for.

## Iv Application

We now illustrate the previous formalism, and compare it with direct -body simulations (whose details are presented in Appendix C).

Following Gupta and Mukamel (2011), we first consider a system driven by interactions of the form

 U(x)=−α1P1(x);Uext(u)=Dextu2, (15)

with , , and . In that case, Eq. (5) gives that is a first degree polynomial in , i.e. the frequency profile is monotonic. As for the system’s DF, we consider a waterbag DF

 F(u)=CΘ(sin(a)−|u|), (16)

with the Heaviside function, and a normalisation constant. We pick , for which the system is linearly stable (Gupta and Mukamel, 2011). The gradient of this DF involves Dirac deltas, which makes the computation of the response matrix immediate, as one can get rid of the integral from Eq. (8), and we refer to Eq. (29) for the associated explicit expression. Yet, because of these infinite gradients, the system also supports neutral modes (i.e. modes with zero growth rates (Chavanis et al., 2005)), which lead to localized divergences in the system’s diffusion coefficients, as detailed in Eq. (32). In Fig. 1, we illustrate the BL prediction for such diverging diffusion coefficients as well as measurements from direct -body simulations (using the procedure described in Appendix C).

Keeping the same interactions as in Eq. (15), one can avoid the presence of neutral modes by considering a smooth DF, for example

 F(u)=Ce−(u/σ)4, (17)

with and a normalisation constant. In Appendix D, we present our implementation of the matrix method, and check that such a system is linearly stable (see Fig. 6). In Fig. 2, we illustrate the BL diffusion coefficients and the associated -body measurements for that system.

Glancing at Eq. (11), we argued that a system with a monotonic frequency profile undergoes a kinetic blocking and cannot relax under effects. We illustrate this in Fig. 3 for the waterbag DF from Eq. (16).

Following Lourenço and Rocha Filho (2015), the dependence of the relaxation rate with is estimated through the quantity , with the average over all the particles of a given realisation. For a given , the time series of is averaged over realisations, as illustrated in the top panel of Fig. 3. Finally, for a given threshold value , we determine the crossing time such that . Should the BL equation (11) have a non-vanishing flux, one expects the scaling . The dependence of for the waterbag DF is illustrated in the bottom panel of Fig. 3. In the range , we measure the scaling , which is expected to converge to for larger values of  (Lourenço and Rocha Filho, 2015). This system indeed suffers from a kinetic blocking because of the impossibility of non-local resonant couplings for a monotonic frequency profile.

In order to recover the scaling predicted by the BL equation (while assuming that as in Eq. (15)), one has to consider a model in which higher harmonics ( or higher) contribute to the pairwise interaction. To illustrate this point, we finally consider a system driven by interactions of the form

 U(x)=−α1P1(x)−α3P3(x);Uext(u)=Dextu2, (18)

with , , and . In that case, Eq. (5) gives that is a non-monotonic second degree polynomial in . We choose the system’s axisymmetric DF to be

 F(u)=Ce−(u−u0)2/(2σ2), (19)

with and , and illustrate it in Fig. 4.

Following Appendix D, we checked that such a system is linearly stable. In Fig. 5, we estimate the scaling of the system’s relaxation with the number of particles.

In the range , we measure a scaling of the form , that is in sensible agreement with the prediction from the BL equation. Because this system can support non-local orbital resonances, it relaxes much more efficiently than kinetically-blocked systems.

## V Conclusion

The inhomogeneous BL equation is being increasingly used to constrain complex dynamical regimes, such as the Hamiltonian Mean Field (HMF) model (Benetti and Marcos, 2017), razor-thin stellar disks (Fouvry et al., 2015), or stellar systems with or without central mass (Bar-Or and Fouvry, 2018; Hamilton et al., 2018).

In the present letter, we illustrated how the same method may be applied to characterize the dynamics of long-range coupled particles on a sphere in the axisymmetric limit. Once one has recognized that this system’s evolution equations are formally identical to the ones of a long-range interacting integrable system, the derivation of the kinetic theory becomes straightforward. In the present case, the reduced number of dimensions of phase space imposes additional geometrical constraints to the system’s dynamics, e.g. allowing only for resonance. We detailed how in the presence of a monotonic frequency profile, the system is submitted to a kinetic blocking and cannot relax under effects, a behavior already encountered in the context of axisymmetrically distributed point vortices (Chavanis and Lemou, 2007).333 Because of the absence of a (quadratic) kinetic energy term in the Hamiltonian from Eq. (1), the present model shares some similarities with point vortices, e.g. the existence of negative temperature statistical equilibria, or a similar BL equation for axisymmetric distributions of point vortices. Depending on the harmonic indices present in the interaction potential, the shape of for the present model can be independent of time, while in the vortex case it depends on time as it is obtained self-consistently from the system’s density profile (Chavanis and Lemou, 2007). In particular, point vortices systems can still undergo a kinetic blocking even if the frequency profile is initially non-monotonic, provided that it becomes monotonic during the evolution. This blocking gets lifted in the presence of a non-monotonic frequency profile, for which non-local resonant couplings are possible.

To emphasize the strength of the BL formalism, we presented quantitative comparisons with direct numerical simulations, recovering both the individual diffusion coefficients, as well as the expected scaling of the relaxation rate with the number of particles.

Despite its recent success, the kinetic theory of long-range interacting systems still asks for more developments, in particular to describe systems with fully degenerate frequency profiles (i.e. , e.g. in the isotropic limit of the present system), or to obtain the kinetic equation for systems undergoing a kinetic blocking (Rocha Filho et al., 2014; Lourenço and Rocha Filho, 2015).

###### Acknowledgements.
JBF acknowledges support from Program number HST-HF2-51374 which was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5–26555. BB is supported by membership from Martin A. and Helen Chooljian at the Institute for Advanced Study.

## Appendix A The inhomogeneous Balescu–Lenard equation

In this Appendix, we repeat the main results regarding the inhomogeneous BL equation, first derived in Heyvaerts (2010); Chavanis (2012). These results are used in the main text to concisely derive the kinetic equation for the problem at hand.

We generically consider an Hamiltonian system in dimensions, and write the phase space canonical coordinates as , respectively the angle and action coordinates (Binney and Tremaine, 2008). The system is assumed to be in an integrable steady state, and following Jeans’ theorem (Jeans, 1915), it can be described by a DF of the form, , that we normalize as . We denote the mean (integrable) potential as , with the associated orbital frequencies . The system comprises particles of mass , coupled one to another through a long-range interaction potential .

As a result of the finite number of particles, the system’s orbital structure gets slowly distorted. To first order in this dynamics is described by the inhomogeneous BL equation

 ∂F(J)∂t=π(2π)dμ∂∂J⋅[∑k,k′k∫dJ′∣∣ψdkk′(J,J′,k⋅Ω(J))∣∣2 ×δD(k⋅Ω(J)−k′⋅Ω(J′))(k⋅∂∂J−k′⋅∂∂J′)F(J)F(J′)]. (20)

Following Kalnajs matrix method (Kalnajs, 1976), one introduces a biorthogonal basis of potentials and densities , with the convention

 ψ(p)(w)=∫dw′U(w,w′)ρ(p)(w′), ∫dwρ(p)(w)ψ(q)∗(w)=−δqp. (21)

The pairwise interaction, , can then be cast under the separable form

 U(w,w′)=−∑pψ(p)(w)ψ(p)∗(w′). (22)

The separable decomposition from Eq. (22) is a fundamental equation that determines what are the natural basis elements appropriate for a given problem, e.g. as in Eq. (6). The bare susceptibility coefficients, , are the Fourier transform of the pairwise interaction w.r.t. , namely

 U(w,w′) =∑p∑k,k′∈Zdψkk′(J,J′)ei(k⋅θ−k′⋅θ′), (23) ψkk′(J,J′) =∫dθ(2π)ddθ′(2π)dU(θ,J,θ′,J′)e−i(k⋅θ−k′⋅θ′),

and Eq. (20) with is the inhomogeneous Landau equation. One can write

 (24)

with standing for the Fourier transform of the basis elements. To account for collective effects (i.e. the spontaneous amplification of fluctuations), one replaces the bare susceptibility coefficients by their dressed analogs

 (25)

as introduced in the BL Eq. (20). Finally, in Eq. (25), the response matrix of a long-range integrable system is given by

 ˆMpq(ω)=(2π)d∑k∫dJk⋅∂F/∂Jω−k⋅Ω(J)ψ(p)∗k(J)ψ(q)k′(J′), (26)

and its effective numerical computation is briefly illustrated in Appendix D. We note that (i) Eq. (24) can be obtained from Eq. (25) by taking (i.e. switching off collective effects); (ii) the substitution of Eq. (22) into Eq. (21) leads to an identity; (iii) Eq. (22) can be obtained from Eqs. (23) and (24) showing the unicity of this expression. As emphasized after Eq. (8), the response matrix, , characterizes the linear stability of the system.

## Appendix B The case of the Heisenberg spins

In this Appendix, we consider the case where the system’s interaction is limited to only one harmonic , so that . This is of particular relevance for classical Heisenberg spins ((Gupta and Mukamel, 2011; Barré and Gupta, 2014), and the Maier-Saupe model for liquid crystals ((Maier and Saupe, 1958; Roupas et al., 2017). In that limit, the dressed susceptibility coefficients from Eq. (9) become

 ψdkk′(u,u′,ω)=−δk′kckℓ(u)ck∗ℓ(u′)εkℓ(ω), (27)

with , and the susceptibility coefficient, , follows from Eq. (8) reading

 εkℓ(ω)=1−2π∫duk∂F/∂uω−kΩ(u)|ckℓ(u)|2. (28)

Let us now follow Eq. (15) and restrict ourselves to the pairwise interaction , with . In that context, Eq. (28) reduces to

 ε±11(ω)=1∓π∫du(1−u2)∂F/∂uω∓Ω(u), (29)

with (see Eq. (5) for the definition of ). As such, in Eq. (29), we immediately recover the susceptibility coefficients obtained by a more complex method in Barré and Gupta (2014) (see Eq. (33) therein). For a waterbag DF as in Eq. (16), the expression of the susceptibility coefficients can be further simplified to become

 ε±11(ω)=1+Dext(1−sin2(a))ω2−(2Dextsin(a))2≡ε(ω). (30)

Such a system is linearly stable if there exists no (with ) for which . The constraint , naturally imposes . As for the constraint , and introducing the energy of the system as , one concludes that the system admits no unstable modes if

 κ=Dext−3ϵ12Dextϵ (31)

satisfies . Introducing the critical energy , the system is therefore stable if , and we recover the criterion put forward in Gupta and Mukamel (2011); Barré and Gupta (2014).

For a waterbag DF from Eq. (16), one can finally compute explicitly the diffusion coefficient from Eq. (14). The frequency, , being monotonic, the resonance condition is straighforwardly solved, and Eq. (16) gives

 D2(u)=(2π)2μDext|c1(u)|4|ε(Ω(u))|2F(u), (32)

with . Introducing , with , the inverse of the susceptibility coefficient reads

 1ε(ω)=1−˜ω2(1−κ)−˜ω2, (33)

where we recall that for a linearly stable system. As a consequence, for , the inverse of the susceptibility coefficient becomes infinite, i.e. the system supports a neutral mode (Chavanis et al., 2005). This leads in particular to a divergence of the diffusion coefficient in , as highlighted in Fig. 1. We finally note that in the absence of collective effects (i.e. in the Landau limit), these divergences vanish.

## Appendix C The N-body implementation

In this Appendix, we present our -body implementation of the problem at hand. Glancing back at the equations of motion Eq. (3), one notes that the velocity vector, , is expressed only as a function of the current particle’s location, , and the instantaneous values of the magnetisations, . There are such evolution equations, but since magnetisations are shared by all particles, their computation can be done only once per timestep. As a result, the overall complexity of advancing the particles for one timestep scales like , with the maximum harmonic number appearing in the considered pairwise interaction in Eq. (2).

The heart of the -body implementation is then (i) to compute efficiently the spherical harmonics (and the vector ones) at the location of the particles, (ii) to compute the magnetisations in Eq. (4), and (iii) to compute the velocity fields in Eq. (3). To compute the spherical harmonics, it is convenient to work with real spherical harmonics. These are computed following (Press et al., 2007, see Eq. (6.7.9)) for the renormalized associated Legendre polynomials, and the second-order recurrence relation (similarly for ) for the azimuthal component. For the (real) vector spherical harmonics, we follow the recurrences presented in (Mignard and Klioner, 2012, see Appendix (B.2)), adapted to the renormalized associated Legendre polynomials. Once all velocity vectors are determined, particles are advanced for a timestep ( in all our applications) using a fourth-order Runge-Kutta integrator (Press et al., 2007, see Eq. (17.1.3)).

To measure diffusion coefficients in -body simulations, we proceed similarly to Benetti and Marcos (2017) in the case of the HMF model. We note that the present model shares some similarities with the HMF model (Antoni and Ruffo, 1995), in particular the property of being a decoupled -body problem, i.e. it can be integrated in operations per timestep, rather than . Here, to measure diffusion coefficients, we perform different realisations, with particles. At the initial time, particles are divided among action bins of size . For each realisation and action bin, we determine the time series of the mean square variation of , averaged over all the particles initially within the bin. For every action bin, these time series are then averaged over all realisations, and considered up to the time where . The diffusion coefficients are obtained finally through a linear fit of these ensemble-averaged series, as in Figs. 1 and 2.

## Appendix D The matrix method

The generic expression of the response matrix is given by Eq. (8). It asks to compute an expression of the form

 ∫1−1dug(u)h(u)+iη ≃∑i∫δu2−δu2dxaig+bigxaih+bihx+iη (34)

where are real, and is an imaginary part added to the frequency . To get the r.h.s. of Eq. (34), we followed Fouvry et al. (2015), truncated the integration domain into regions, introducing , so that the center of each region is with , and performed a linear expansion of the numerator and denominator, so that and (similarly for ). To get the second line, we assumed , and introduced the dimensionless function

 ℵD[b,c,η]=∫12−12dx1+bx1+cx+iη=G[12]−G[−12], (35)

where for the primitive , one can choose

 G[x]= bxc+bη+i(c−b)2c2 (36)

We illustrate this method in Fig. 6, by representing the Nyquist contours associated with the system from Eq. (17). This shows that this particular system is linearly stable.

Throughout the applications presented in the main text, we truncated the orbital space in elements, and added to the frequency a small imaginary part to regularize the resonant denominator. We checked that these choices had no impact on our results.

## References

• Lynden-Bell (1967) D. Lynden-Bell, MNRAS 136, 101 (1967).
• Balescu (1960) R. Balescu, Phys. Fluids 3, 52 (1960).
• Lenard (1960) A. Lenard, Ann. Phys. 10, 390 (1960).
• Heyvaerts (2010) J. Heyvaerts, MNRAS 407, 355 (2010).
• Chavanis (2012) P.-H. Chavanis, Physica A 391, 3680 (2012).
• Arfken et al. (2013) G. Arfken, H. Weber, and F. Harris, Mathematical Methods for Physicists (Elsevier Science, 2013).
• Gupta and Mukamel (2011) S. Gupta and D. Mukamel, J. Stat. Mech. 3, 03015 (2011).
• Barré and Gupta (2014) J. Barré and S. Gupta, J. Stat. Mech. 2, 02017 (2014).
• Maier and Saupe (1958) W. Maier and A. Saupe, ZNatA 13, 564 (1958).
• Roupas et al. (2017) Z. Roupas, B. Kocsis, and S. Tremaine, ApJ 842, 90 (2017).
• Kocsis and Tremaine (2015) B. Kocsis and S. Tremaine, MNRAS 448, 3265 (2015).
• Szölgyén and Kocsis (2018) Á. Szölgyén and B. Kocsis, Phys. Rev. Lett. 121, 101101 (2018).
• Takács and Kocsis (2018) Á. Takács and B. Kocsis, ApJ 856, 113 (2018).
• Binney and Tremaine (2008) J. Binney and S. Tremaine, Galactic Dynamics: Second Edition (Princeton University Press, 2008).
• Chavanis (2013a) P.-H. Chavanis, A&A 556, A93 (2013a).
• Chavanis (2013b) P.-H. Chavanis, Eur. Phys. J. Plus 128, 126 (2013b).
• Rocha Filho et al. (2014) T. M. Rocha Filho, A. E. Santana, M. A. Amato, and A. Figueiredo, Phys. Rev. E 90, 032133 (2014).
• Lourenço and Rocha Filho (2015) C. R. Lourenço and T. M. Rocha Filho, Phys. Rev. E 92, 012117 (2015).
• Chavanis and Lemou (2007) P.-H. Chavanis and M. Lemou, Eur. Phys. J. B 59, 217 (2007).
• Chavanis et al. (2005) P.-H. Chavanis, J. Vatteville, and F. Bouchet, Eur. Phys. J. B 46, 61 (2005).
• Benetti and Marcos (2017) F. P. C. Benetti and B. Marcos, Phys. Rev. E 95, 022111 (2017).
• Fouvry et al. (2015) J.-B. Fouvry, C. Pichon, J. Magorrian, and P.-H. Chavanis, A&A 584, A129 (2015).
• Bar-Or and Fouvry (2018) B. Bar-Or and J.-B. Fouvry, ApJ 860, L23 (2018).
• Hamilton et al. (2018) C. Hamilton, J.-B. Fouvry, J. Binney, and C. Pichon, MNRAS 481, 2041 (2018).
• Jeans (1915) J. Jeans, MNRAS 76, 70 (1915).
• Kalnajs (1976) A. J. Kalnajs, ApJ 205, 745 (1976).
• Press et al. (2007) W. Press et al., Numerical Recipes 3rd Edition (Cambridge University Press, 2007).
• Mignard and Klioner (2012) F. Mignard and S. Klioner, A&A 547, A59 (2012).
• Antoni and Ruffo (1995) M. Antoni and S. Ruffo, Phys. Rev. E 52, 2361 (1995).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters