Relaxation of spherical stellar systems
Jun Yan Lau James Binney



=2pt

    ∂ı iR_ b#1⟨#1⟩
        
        =cmmib10
        =cmmib7
        ==#1
        	#1#1#1#11112131415161718192021222324252627282930313233ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz.⊤#1#2.#1#1#2.#2. v⊤
        ⊤. Nḍv_ c Jk_ BMNRASAnnRA&AApJApJSApJLNatAJA&AA&ASPASPJCAPPASJPhys. Rev. Edf NBODY6ΛCDMLOSIℒS Gyr Myr kpcmixrescutconst km s^-1 km mas yr^-1 mas GeV cm^-3  M_⊙v_∥ pcdeR_cR_LR_d#1#2#1#2HItmr_ sz_ stol_Jrms^ T^ Tθ^ TJ^ TΦ^ T^ T^ Tr^ Tϑ^ Tp^ T tm sgn#1Fig. <ref> sos^2^-1_ pω_ pRΩ
         LRelaxation of spherical stellar systems]
        Relaxation of spherical stellar systemsJun Yan Lau & James Binney]
          Jun Yan Lau^1, James Binney^1E-mail: binney@thphys.ox.ac.uk
        ^1Rudolf Peierls Centre for Theoretical Physics, Clarendon Laboratory,
          Parks Road, Oxford, OX1 3PU, UK
        [
            [
            December 30, 2023
        =====================
        
        
        10 000 simulations of 1000-particle realisations of the same cluster are
        computed by direct force summation. Over three crossing times the original
        Poisson noise is amplified more than tenfold by self-gravity.
        The cluster's fundamental dipole mode is strongly excited by Poisson noise,
        and this mode makes a major contribution to driving diffusion of stars in
        energy. The diffusive flow through action space is computed for the
        simulations and compared with the predictions of both local-scattering
        theory and the Balescu-Lenard (BL) equation. The predictions of local-scattering
        theory are
        qualitatively wrong because the latter neglects self-gravity. These results
        imply that local-scattering theory is of little value. Future work on cluster
        evolution should employ either N-body simulation or the BL
        equation. However, significant code development will be required to make use
        of the BL equation practicable. 
        
          Galaxy:
          kinematics and dynamics – galaxies: kinematics and dynamics – methods:
          numerical
        
        
        § INTRODUCTION
        
        
        The central idea of stellar dynamics is that in a first approximation stars
        move in the `mean-field' gravitational potential that one obtains by
        smearing the mass of each star over a region that is somewhat larger than the
        local inter-star distance. Then in a second step one considers how stars
        drift between orbits in the mean-field potential on account of differences
        between the actual potential of the cluster and the mean-field potential.
        This process of drift between orbits is termed `dynamical relaxation' because
        it drives secular increase in a cluster's entropy. Relaxation causes stars
        with larger masses to be more strongly concentrated towards the cluster
        centre than less massive stars (`equipartition') and the cluster's central
        density to increase while its halo becomes more extended (`core collapse').
        It also causes the mass of the cluster to diminish through stars being
        occasionally accelerated to velocities larger than the local escape speed
        (`evaporation'). 
        
        These general principles have been understood since the seminal work of
        <cit.>, <cit.> and <cit.>. For decades it was not
        feasible to follow the relaxation of a globular cluster by brute-force
        computation. Hence relaxation was investigated, both analytically and
        numerically, by adopting an ansatz for the difference between the
        gravitational field of the cluster and its mean-field model. The ansatz was
        that this consisted of the field of a point mass around each star, the idea
        being that close to each star the mean field of the cluster is negligible
        compared to the Kepler field, while far from each star the reverse condition
        holds <cit.>. Hence relaxation can be modelled by summing large
        numbers of deflections as pairs of stars scatter off each other at relatively
        close quarters. Each such scattering could be computed via a hyperbolic
        Kepler orbit of the reduced particle.
        
        When this `local' relaxation theory is worked out in detail, the rate of
        relaxation emerges as an integral over the impact parameters b of these
        hyperbolic orbits that diverges as ln b at large b<cit.>.
        This weak divergence was mastered by imposing a largest impact parameter to
        be considered, b_ max. <cit.> took b_ max to be the
        local inter-particle distance. <cit.> argues that it should rather be
        the larger of the cluster's core radius and distance to the cluster centre.
        Since the divergence is only logarithmic, the relaxation rate computed using
        either of these choices of b_ max does not differ hugely, even though
        the chosen values of b_ max are typically very different.
        
        Quite recently <cit.> and <cit.> developed a radically
        different approach to relaxation. A rather formal analysis in the framework
        of angle-action variables and using potential-density pairs to solve
        Poisson's equation yields the `Balescu-Lenard equation', hereafter the BL
        equation. This equation states
        that a star diffuses through phase space as a consequence of interacting
        resonantly with other stars in the cluster.  <cit.> does not start
        from the assumption that relaxation proceeds by interaction in pairs, but
        arrives at this conclusion as a consequence of the structure of the
        Boltzmann equation
        
            f t+[f,Φ]=0.
        
        Here f(,,t) is the one-particle distribution function (DF),
        Φ(,t) is the cluster's full (fluctuating) potential and [,] denotes
        a Poisson bracket. When f and
        Φ are decomposed into mean-field and fluctuating parts f_0, f_1,
        Φ_0 and Φ_1, it follows that the evolution of f_0 (i.e.,
        relaxation) is given by
        
            f_0 t=[Φ_1,f_1].
        
        That is, relaxation is driven by correlations between the fluctuating parts
        of the potential and the DF. Two processes drive correlations: (i) by
        Poisson's equation, the potential will be deeper where there are more
        particles, and (ii) there will be more particles where the potential pushes
        particles together. Mechanism (ii) operates even on a population of test
        particles, whereas mechanism (i) operates only when particles bear mass.
        Roughly, mechanism (ii) is responsible for evaporation and mechanism (i) is
        responsible for equipartition and dynamical friction.
        
        Thus the sophisticated BL equation confirms Chandrasekhar's
        ansatz that relaxation occurs through interactions in pairs. However, it does
        not confirm that these interactions are local. In fact, it replaces the local
        criterion with a requirement that the stars can resonate in the sense that
        there exist vectors , ' with integer components such that
        ·+'·'=0, where  and ' are
        (two- or three-dimensional) vectors formed by the stars' fundamental orbital
        frequencies.
        
        Unlike Chandrasekhar and his followers, Heyvaerts and Chavanis included the
        self-gravity of the cluster in their computation of the diffusion rate. Thus
        stars interact with each other not through the vacuum but through the
        polarisable medium formed by the rest of the system. <cit.>
        first demonstrated that polarisation effects could be important by computing
        the wake that a mass M that is on a circular orbit in stellar disc raises
        in the surrounding star field. The temperature of the disc is best quantified
        by the Q parameter <cit.>, which is the ratio of the velocity
        dispersion in the disc to the critical value below which the disc is
        (Jeans) unstable to axisymmetric disturbances. <cit.> showed
        that for a realistic value Q≃1.4 the mass in the wake is several times
        larger than M. Hence the effective mass of a disc star is several times
        larger than its real mass because it polarises the disc around it, and we
        must expect this increase in mass to facilitate exchanges of energy and
        momentum with other disc stars. 
        
        The BL equation tells us that  each star shakes the system at its natural
        frequencies ·. If the system has a natural mode of
        vibration near this frequency, the star excites this mode, and the mode's
        energy  may be absorbed by a star that is quite distant but also happens
        to resonate with the mode.
        
        <cit.> showed that inclusion of self gravity accelerates the
        relaxation of a razor-thin disc by a factor nearer 1000 than 10 on account of
        swing amplification, which <cit.> first worked out for a stellar
        disc following its discovery in a gas disc by <cit.>.  Stars,
        afforced by their wakes, launch running spiral waves, which are
        swing-amplified at their corotation annulus and are finally resonantly
        absorbed at a Lindblad resonance. Thus because disc stars communicate via
        running waves that have an amplifier on the line, their interaction is by no
        means local. The interaction is also ∼1000 times stronger than it would
        be in the absence of self gravity. This strengthening, together with the
        insight that the coupling is resonant, resolved a decades-old puzzle as to
        why N-body simulations of discs that only have stable normal modes develop
        strong spiral structure and then bars <cit.>.
        
        Thus the BL equation has led to a major advance in our
        understanding of stellar discs.   (hereafter H18) asked
        what its implications were for globular clusters, which unlike discs, have
        been thought to be adequately described by  local-scattering theory. After
        all, there is no known analogue of the swing amplifier in a globular cluster.
        
        The BL equation describes the diffusion of stars through action
        space, and provides a prescription for the computation of the relevant
        diffusion coefficients.  Local theory provides an alternative prescription
        for computing the diffusion coefficients, so key questions are (i) do the two
        frameworks yield materially different coefficients? and, if so, (ii) which
        framework's coefficients are more accurate?
        
        A razor-thin disc and a spherical cluster both have an effectively
        two-dimensional action space (H18).  To compute the BL
        coefficients at a location  in this space, one must (i) identify lines
        in action space on which the resonance condition
        ·+'·'=0 is satisfied, and then (ii) integrate a
        certain function along each such line (i.e., vary ' so the condition
        remains satisfied). In the case of a razor-thin disc, it is physically
        plausible that a small number of lines dominate the diffusion coefficients,
        but in a spherical cluster this is not true. Given that it was not feasible
        to evaluate the necessary integrals along all resonant lines, H18
        evaluated them along lines defined by integers |n_i|≤2. Since this scheme
        led them to omit very many resonant lines, including the lines most likely to
        be effectively included by local theory, it would be natural for their values
        of the diffusion coefficients to be smaller than the coefficients computed in
        the local approximation. But their values were  not smaller. By
        repeating their calculation with self-gravity neglected, they were able to
        show that their diffusion coefficients are strongly enhanced by self-gravity,
        which is neglected in the local approximation. This finding cast doubt of the
        applicability of the local approximation to globular clusters, but in light
        of the incompleteness of their treatment H18 had no basis to
        claim that their diffusion coefficients were more accurate than those yielded
        by the local approximation.
        
        In a spherical system, the second component of the vector  in a resonant
        condition is always the angular-momentum quantum number ℓ. Hence, H18's
        sums over n_i≤2 were sums over certain monopole, dipole and quadrupole
        distortions of the cluster. That is, they computed the contributions to
        diffusion coefficients that arise from stars exciting/damping simple global
        distortions of the cluster. In these circumstances it is no surprise that
        neglect of self-gravity much diminishes the diffusion coefficients. For
        example, the simplest  dipole (ℓ=1) mode involves displacing the core
        of the cluster with respect to its envelope. If this is done with
        self-gravity neglected, so in a fixed potential, the core being anchored by
        the potential, can barely move.  By contrast, when self-gravity is included
        the core is free to move because it will take its potential with it. Similar,
        though weaker, arguments apply to motion of the envelope.
        
        H18 showed that their (incomplete) BL coefficients and the
        coefficients from the local approximation predict qualitatively different
        diffusive flux vectors in action space. So the two theories will predict
        different evolutionary tracks.  Given that the work of H18 strongly suggests
        that the cluster's low-order global modes, which are excluded from the local
        approximation, contribute significantly to the diffusion coefficients, one
        can have no confidence in the essential validity of  the local approximation.
        Yet  decades of work by many groups has been underpinned by the local
        approximation. Can it really be seriously in error?
        
        It is now feasible to integrate directly the equations of motion of clusters
        with realistic numbers of stars. By examining orbital diffusion in such
        simulations, one can rigorously test any theory of diffusion.  Surprisingly
        few studies have done this, however. An early attempt was that of
        <cit.>, who measured the second-order diffusion coefficient (̣Δ
        E)^2/ṭ in the energies E of stars in <cit.> models with
        different particle numbers, and compared these measurements with values
        obtained in the local approximation.  He concluded that `overall' the values
        agreed to an `impressive' extent although the measured values were larger
        than predicted by factors 1.5–2 in a more centrally concentrated model and
        smaller than expected by similar factors in a less centrally concentrated
        model. <cit.> used a Fokker-Planck (FP) code, which is based on the
        local approximation to compute the evolution of rotating King models with
        10 240 stars and compared its predictions for the evolution of the central
        density, velocity dispersion, rotation rate and velocity anisotropy with what
        they found by direct integration of N bodies. The FP predictions were
        qualitatively correct but gave rise to quantitative discrepancies.
        
        Given that the local approximation has in the Coulomb logarithm an
        effectively free parameter, it calls for more rigorous validation than
        checking consistency with a few numbers.  A more demanding  comparison of the
        predictions of local theory with direct N-body integration is by comparing
        diffusive fluxes at each point in action space. It is these fluxes that drive
        the evolution in density and velocity profiles studied by <cit.>, and
        diffusion coefficients in actions  are more meaningful than the
        diffusion coefficients in energy studied by <cit.> because changes
        in E, unlike changes in the adiabatic invariants  are not solely
        caused by the fluctuating component of the potential, but have a significant
        contribution from the secular evolution of the mean-field potential Φ_0.
        
        In this paper we use direct N-body integration to determine the diffusive
        flux () in action space for a cluster of 1000 stars that has the
        isochrone DF <cit.>, and we compare this () with the
        fluxes computed by H18 for the same system from (i) the local approximation,
        and (ii) the BL equation. Section <ref> describes the simulations,
        Section <ref> discusses the key issue of choice of cluster centre
        and presents evidence for powerful excitation of a dipole distortion.
        Section <ref> characterises the fluctuating component of a cluster's
        gravitational field and links this to diffusion of stars in energy.
        Section <ref> presents the drift of stars through action space and
        discusses rates of entropy generation. Section <ref> sums up.
        Appendices summarise units for N-body models and the computation of diffusion
        coefficients for energy in the local-scattering approximation.
        
        
        
        § THE SIMULATIONS
        
        
        The isochrone is a spherical system of mass M that has the gravitational potential
        
            Φ(r)=-GM b+√(b^2+r^2),
        
        where b is the system's scale radius.
        We created 10 000 realisations of the isochrone, each with  1000
        equal-mass stars. Each star's  radial coordinate r was chosen by picking a
        random number ξ_1 uniformly on (0,1) and then inverting
        the isochrone's cumulative mass distribution
        
            ξ_1(r)=1 M∫_0^rṛ r^2ρ_ I(r),
        
         where ρ_ I(r) is given by equation (2.49) of <cit.>. 
        Once r had been chosen, a second  random number ξ_2 was picked and the
        particle's speed v was determined by inverting
        
            ξ_2(v)=∫_0^vṿ v^2f_ I[12v^2+Φ(r)]∫_0^∞ṿ v^2f_ I[12v^2+Φ(r)],
        
        where f_ I(E) is the DF of the isochrone <cit.> and
        the function of v on the right of this equation was inferred by
        interpolation on a grid of numerically computed values.  The directions of
        the position and velocity vectors  and  were chosen randomly. This
        sampling procedure was validated by binning the stars in energy and comparing
        the resulting histogram with the analytically determined differential energy
        distribution <cit.>. 
        
         <cit.> was then used to integrate the equations of motion of
        each realisation.  automatically boosts to the zero-momentum frame of
        the initial conditions, and rescales positions and velocities such that
        distances and times are measured in units of a length r_N and time t_N
        that are derived from the potential and kinetic energies of the initial
        conditions (Appendix <ref>). The ratio r_N/b varies slightly from
        realisation to realisation, but in Appendix <ref> we show that
        b=0.2375r_N. The isochrone has a natural timescale
        
            t_ I≡(b^3/GM)^1/2.
        
        In Appendix <ref> we show that the N-body time unit t_N=8.642t_
        I. In the following  distances and times are given in units of b and
        t_ I,
        ignoring the small scatter in r_N and t_N between realisations.
        
        The equations of motion of each
        realisation were integrated for a time 10t_N, with snapshots of the
        phase-space coordinates saved after each elapse of δ t=0.5t_N. Thus
        the data comprise 10 000 sets of 21 1000-particle snapshots.
        
        
        
        
        § CHOICE OF CENTRE
        
        
        We wish to study relaxation via changes in quantities that in the mean-field
        model are constants of motion. To do this we have to choose a centre. Three
        possible choices spring to mind:
        
        
          1. The mean-field centre, i.e., the point r=0 around which each
        realisation was assembled. By symmetry it must coincide with the expectation
        value of any credible estimator of the centre. The mean-field centre is the
        only centre that is defined in the context of the BL equation.
        
        
          2. The barycentre. Momentum conservation ensures that a realisation's
        barycentre defines an inertial frame of reference. Moreover, in our chosen
        zero-momentum frame the barycentre is absolutely stationary. A severe
        disadvantage of the barycentre is that it is sensitive to the small number of
        stars that are at large radii, and the distribution of its offsets from the
        mean-field centre is long-tailed: a Gaussian fit to any component of this
        offset yields a dispersion s_ b∼2.3b yet  ∼34 percent of
        realisations had barycentres that lay more than 10b from the
        mean-field centre. In view of this finding, we did not use the barycentre.
        
        
          3. The potential centre. Our procedure for initialising the
        simulations makes it exceedingly unlikely that a realisation initially
        includes a hard binary, and the timescale for hard binaries to form via
        three-body interactions  is much longer than the
        duration of our integrations <cit.>. So binaries are not a concern. This being so,
        the star with the most negative potential energy is likely to lie close to
        a cluster's physical centre. Consequently, serious consideration should be
        given to computing actions with respect to this centre.  
        
        
        
        The problem with the use of the potential centre is its erratic
        motion. In as much as it is displaced from the barycentre, it must be in
        orbit around that fixed point. This orbit can be considered to be a
        consequence of the force that the somewhat
        lop-sided halo exerts on the core. On top of this smooth motion, the potential
        centre occasionally jumps discontinuously when the distinction of being the most bound
        star switches between stars. These jumps can be made more frequent but
        smaller in magnitude by adopting as the potential centre not the location of
        the most bound particle but the barycentre of the n most bound particles.
        We found n=30 to be a sensible choice.
        
        If we use the mean-field centre to compute actions, those actions will have
        non-zero rates of change because the actual potential  Φ differs from
        the mean-field potential Φ_ mf. For this centre the rate of change of any action
        is 
        
            J̣ṭ=[J,ΔΦ]
        
        where [,] is a Poisson bracket and
        
            ΔΦ≡Φ-Φ_ mf.
        
         If we define actions relative to  the potential centre,
        the actual potential Φ will be much more nearly a function of distance
        fronm the potential centre, so if
        equation (<ref>) applied with ΔΦ equal to the 
        difference between Φ and the closest function of distance from the
        potential centre, |J̇| would be
        small. However, since the potential centre does not define an inertial
        coordinate system, terms describing pseudo-forces would need to be added to the right
        side of equation (<ref>), and typical values of |J̇| would
        probably be comparable to the values encountered when using the mean-field
        centre. 
        
        In the following actions will be defined with respect to the
        mean-field centre.
        
        
        
         §.§ Dynamics of the potential centre
        
        
        The practical difficulties with defining integrals of motion with respect to
        the potential centre do not alter the fact that it has physical significance.
        In fact, understanding its motion is crucial for understanding how clusters
        relax. 
        
        Given that the
        potential centre of each cluster is independently and isotropically scattered
        around the mean-field centre, by the central limit theorem the distribution
        of the distances s between these centres must be Maxwellian.
        Fig. <ref> shows  the variance s^2 of these distances as a
        function of integration time t. The variance starts  small and grows as
         advances the particles. We confirmed that it does  not grow if the
        mean-field potential is used to advance particles. Thus the growth of
        s^2 evident in Fig. <ref> is an effect of self-gravity.
        
        The initial variance is small because it
        is entirely due to Poisson noise in sampling the core of the cluster. As the
        equations of motion are integrated, it grows at a rate that for t t_N
        settles to a constant. Thus the variance tends to linear dependence on
        integration time: s^2∝ t. This is the signature of a
        random walk rather than steady motion, which would imply s∝ t. This
        establishes that the potential centre experiences significant acceleration.
        
        Given that each cluster's barycentre is absolutely fixed, the potential
        centre cannot execute an unconstrained random walk, for that would eventually
        take it arbitrarily far from the barycentre. Hence the potential centre must
        executes a tethered random walk, like a randomly kicked ball that is attached
        by a piece of elastic to a fixed peg. The standard deviation of the distances
        s_ b∼2.3b of the barycentres from the mean-field centre sets the scale of
        the region over which we expect the potential centre to wander.  Since this
        scale is bigger than the standard deviations s^2^1/2 plotted in
        Fig. <ref>, the elastic is slack throughout our integrations
        because they start with the potential centre anomalously close to the
        mean-field centre. The elastic being slack, the potential centre is for the
        moment executing an essentially random walk and s^2∝ t in
        consequence.
        
        
        
        Fig. <ref> indicates that the picture just deduced is an
        oversimplification. It shows estimates of the autocorrelation
        
            ξ(τ)=1 N∑_i=1^N1 T-τ∫_0^T-τṭ  (t)·(t+τ)
        
        of the positions of the potential centres. The bottom curve uses all
        snapshots, while the next curve up discards the first snapshot but uses all
        the other ones, and the second curve up discards the first two snapshots, and
        so on. The more snapshots are discarded, the shorter is the range of delays
        τ that can be probed. Aside from this progressive shortening of the
        curves, they have remarkably similar shapes: a steep
        initial decline is followed by a small bump before a plateau is reached. Note
        that this structure in the bottom curve emerges from snapshots that have no
        overlap with the snapshots that give rise to the structure in the topmost
        curves. Thus the structure is highly reproducible.
        
        The curves' vertical offsets are a straightforward consequence of the steady
        rise in s^2 shown in Fig. <ref>: ξ(0) is simply the
        mean of s^2 over the snapshots used to compute ξ.
        
        The horizontal plateaus of the curves are simple consequences of the random
        walk we inferred above. In this picture (t+τ) is
        the sum of (t) and a series of random increments _i, so
        
            ξ(τ)=(t)·[(t)+∑_i_i].
        
        So long as the steps _i are distributed isotropically, the
        expectation value of (t)·_i=0, so ξ(τ)=|(t)|^2 is
        independent of τ.
        
        
        
        What then is the cause of the initial steep drop in ξ(τ)? Since this
        feature occurs in curves that discard the initial conditions and the
        immediately following snapshots, this steep decline cannot be connected to
        our starting the simulations from anomalously regular initial conditions.
        The argument just given to explain the plateaus, suggests that ξ declines
        initially because the first steps _i satisfy
        (t)·_i<0. This condition will be satisfied if steps have a
        tangential bias. Fig. <ref> illustrates this idea. If the motion
        of the potential centre over short timescales resembles an orbit,
        (t)·_i will be negative for of order half an orbital period.
        Later it will be positive for a similar time. When this type of variation is
        averaged over 10 000 realisations, in each of which the potential centre is
        on a slightly different orbit, the positive and negative contributions to
        ξ(τ) tend to cancel once  τ is comparable to an orbital period,
        and a plateau in ξ ensues.
        
        This discussion leads to the tentative conclusion that the delay
        τ≃10t_ I at which the plateaus start in Fig. <ref> is the
        time over which the motion of the potential centre resembles an orbit, and
        after which it is dominated by noise. The noise is initially highly non-stationary
        because the simulations start from anomalously regular initial conditions
        that reflect Poisson noise in the absence of self gravity. As the equations
        of motion are integrated, self-gravity amplifies this noise by a significant
        factor.  
        
        
        
        
        § GLOBAL POTENTIAL FLUCTUATIONS
        
        
        We now quantify the stochastic perturbing potential ΔΦ, which is
        the difference between the true and mean-field potentials.
         returns the value of the true potential at the locations of each
        particle, so it is trivial to compute the value of ΔΦ at each
        particle's location. However, to make contact with the work of H18 we want to
        decompose ΔΦ() in spherical harmonics around the mean-field
        centre. Unfortunately, there is no straightforward way to compute this
        decomposition from the  values on account of the biased
        distribution of particle locations. Consequently, we obtain the desired
        expansion by using spherical harmonics to solve Poisson's equation afresh.
        That is we compute the true potential from the standard formula <cit.>
            Φ(r,θ,ϕ)=∑_ℓ mΦ_ℓ m(r)Y_l^m(θ,ϕ),
        
        where
        
            Φ_ℓ m(r) =  -4π G/2l+1[   1/r^l+1∫_0^rạ ρ_ℓ m(a) a^l+2        + r^l ∫_r^∞ạ  a^1-lρ_ℓ m(a) ].
        
        Here (r,θ,ϕ) are polar coordinates around the mean-field centre and
        
            ρ_ℓ m(a)≡∫^̣2Ω Y_l^m*(θ,ϕ)ρ(a,θ,ϕ).
        
         Our knowledge of the density  ρ() is limited to the Monte Carlo
        sample   provided by . Given this fact, a natural procesdure is
        to estimate ρ() by binning particles in spherical shells of finite
        thickness and then discretising the integrals in equation (<ref>)
        into sums over bins. It is not hard to show that the statistical noise in the
        resulting estimate of Φ_ℓ m is essentially invariant as we reduce the
        thickness of shells: the increase in the noise level of a each ρ_ℓ m is
        precisely compensated by the increase in the number of shells that contribute
        effectively to Φ_ℓ m(r) at a given radius. Hence, it is safe to reduce
        the width of shells so each shell contains at most one particle. In this
        limit the sum over shells in equation (<ref>) becomes a sum over
        particles, with ρ_ℓ m(a) replaced by the value of Y_l^m* at the
        particle's angular location.
        
        
        
        Each black curve in Fig. <ref> shows the average over all
        realisations of the spectral power
        
            c_1(r)=∑_m=-1^1|Φ_1m(r)|^2
        
        of the dipole component of the potential at a given time. The lowest curve is
        for t=0, i.e., the initial conditions, and as the time of integration
        increases, the corresponding curve moves upwards.  Fig. <ref>
        shows c_1(b) as a function of
        time. The curve strongly resembles the curve of s^2 shown in
        Fig. <ref>. In particular, at later times c_1∝ t.  The
        growth of the dipole in Φ by a factor in excess of 10 puts a lower limit
        on the extent to which the initial conditions are unnaturally symmetric
        because they neglect the long-range correlations that grow rapidly when stars
        move in the self-consistent potential.
        
        Each  red curve in Fig. <ref> plots the ensemble average of
        the magnitude
         
            c_2(r)=∑_m=-2^2|Φ_2m(r)|^2
        
         of the quadrupole
        component of Φ(). The curves are similar to the dipole curves, but
        depressed by a factor of order 10 and displaying weaker growth of quadrupole
        distortions as advances the particles.
        
        
        
        From the perspective of perturbation theory, displacement of the physical
        centre from the mean-field centre causes the perturbing potential
        ΔΦ (eqn <ref>) to be dominated by a dipole term, so the
        results of this section are natural consequences of the random walk we
        inferred in the last section.
        
        The dominance of the dipole component of Φ explains why
        H18 found from the Balescu-Lenard equation that the rate of relaxation
        was dominated by the dipole of the perturbing potential. The much reduced
        values of s^2 when the mean-field potential is used to advance
        particles corresponds to strong reduction in the relaxation rate when
        H18 switched off self gravity. On account of the random nature of
        the motion of the physical centre, the perturbing dipole potential has the
        stochastic nature assumed in the derivation of the Balecu-Lenard equation.
        
        H18 found that the largest contribution to the relaxation rate came from the
        ℓ=1 mode in which the cluster's core and halo move in antiphase. In such
        a mode, the dipoles in the density distribution at small and large r should
        have opposite directions.  The complex numbers ρ_1m define a vector
        
            ≡(√2ρ_11,√2 mρ_11,ρ_10).
        
        Each curve in Fig. <ref> plots the average over all
        realisations of
        
            cos(θ)≡_r^-·_r^+
            |_r^-||_r^+|,
        
        where _r^- is the vector (<ref>) defined by the dipole formed by
        all the particles interior to radius r, while _r^+ is the vector
        associated with the remaining particles. The blue curve shows that in the
        initial conditions the interior and exterior dipoles have random relative
        orientations, regardless of the radius taken to define `interior'. As advances the particles, anti-alignment of the dipoles emerges rapidly. We
        take the radius at which a curve passes its minimum as the measure of the
        boundary between the zones of opposite alignment: if r is decreased from
        this value, |cosθ| decreases because shells with interior
        alignments are classified as exterior.  Conversely, if r is increased
        |cosθ| decreases because shells with exterior alignment are
        classified as interior. Fig. <ref> shows that in the first
        snapshot, the boundary lies at r≃4b while in the last snapshot the
        boundary has reached r≃11b. In the initial snapshot
        |cosθ| is very small for r10b, indicating that shells
        outside 10b have randomly aligned dipoles. As  advances the cluster,
        the point at which |cosθ| becomes small moves rapidly outwards
        as larger and larger shells align their dipole in opposition to the cluster's
        core.
        
        
        
        
         §.§ Impact of global fluctuations
        
        
        
        How much relaxation can be attributed to the large-scale fluctuations in
        Φ that we have just quantified? 
        
        To answer this question we 
        start from the equation for the rate of change of a particle's energy
        E=12 v^2+Φ:
        
            Ẹṭ=Φ t,
         
        which yields the change Δ E in E during an interval τ as
        
            Δ E(τ)=∫_ orbit^τṭ Φ t.
        
        Squaring both sides and taking the expectation value over orbits we find
        
            (Δ E)^2_τ   =∫_0^τṭ∫_0^τṭ' Φ tΦ t'
        
        With the expansion of equation (<ref>) this becomes
        
            (Δ E)^2_τ   =∑_ℓ mℓ'm'∫_0^τṭ∫_0^τṭ'    ×⟨Φ_ℓ m^*
            tΦ_ℓ'm' t'
            Y_l^m*(θ,ϕ)Y_ℓ'^m'(θ',ϕ')⟩.
        
        Since we are taking the expectation over spherical clusters, we are
        effectively integrating over angles. If the primed angles were independent of
        the unprimed angles, the expectation value would vanish for l>0. But when
        t≃ t', the primed and unprimed angles are similar, and the
        expectation value approximates
        
            Y_l^m*(θ,ϕ)Y_ℓ'^m'(θ,ϕ)=(4π)^-1δ_lℓ'δ_mm'.
        
        Let 2τ_ℓ be the time during which the primed and unprimed angles are
        sufficiently nearly equal for (4π)^-1δ_ll'δ_mm' to be a
        useful approximation to the expectation of the product of spherical harmonics
        in equation (<ref>). Then we have
        
            (Δ E)^2_τ≃14π∫_0^τṭ∫_t-τ_ℓ^t+τ_ℓṭ' ∑_ℓ m⟨Φ_ℓ m^* tΦ_ℓ m t'⟩.
        
        If Φ were a stationary random process, the integrand of the integral
        over t would be independent of t and equation (<ref>) would
        predict that (Δ E)^2_τ∝τ as we expect from a
        diffusive process. In our simulations, Φ is not a stationary process
        because the simulations start from anomalously small ΔΦ, but
        4π(Δ E)^2_τ/τ̣ is still given by the t integrand of
        equation (<ref>).
        
        
        
        To estimate the typical size of Φ_ℓ m/ t we first note the
        expression for the fluctuating part of the potential:
        
            δΦ(,t)=∑_ℓ mΦ_ℓ m(||,t)Y_ℓ^m().
        
        This yields the mean-square amplitude of the fluctuations as
        
            δΦ^2(r)≃14π∑_ℓ m|Φ_ℓ m(r)|^2.
        
        Hence in  equation (<ref>) we may use the estimate
        
            14π∑_ℓ m⟨Φ_ℓ m^* tΦ_ℓ m t'⟩≃(δΦ)^2τ^2,
        
         where τ is an appropriate average of the τ_ℓ. With
         this estimate equation (<ref>) becomes
        
            (Δ E)^2_τ≃∫_0^τṭ∫_t-τ^t+τṭ' (δΦ)^2τ^2≃∫_0^τṭ (δΦ)^2τ,
        
        where the second equality requires that the expectation now involves
        averaging over particle radii as well as angles.
        
        Now we estimate (δΦ)^2 by ranking stars by radius in each cluster
        snapshot, and for the particle ranked i in the αth realisation we
        compute the values for ℓ=0,1,2 of
        
            Φ^(ℓ)_α i≡∑_mΦ_ℓ m(|_α i|) Y_l^m(_α i),
        
        i.e., the contributions of the monopole, dipole and quadrupole to the
        potential energy of the ith ranked particle in a particular realisation. 
        Then we compute  the sums over realisations
        
            S^(ℓ)_i≡1∑_α=1^(Φ_α i^(ℓ)-Φ_ tot^(ℓ)(_α i))^2,
        
        where Φ_ tot^(ℓ) is obtained from the potential computed from
        the locations of all 10 million stars at the given time. When we average in
        this way over all realisations we obtain a potential that has negligible
        dipole and quadrupole contributions, but a monopole contribution that is very
        similar to, but significantly different from the analytic isochrone
        potential. S_i^(ℓ) is a measure of the amplitude of the fluctuations
        in Φ contributed by the ℓ multipole at the typical radius of the
        ith ranked particle. Hence the quantity (δΦ)^2 encountered
        above is just S_i^(ℓ) summed over ℓ and averaged over i:
        
            (δΦ)^2=∑_ℓS^(ℓ) where S^(ℓ)≡1 N∑_i S_i^(ℓ).
        
        The full, dashed and dotted curves in
        Fig. <ref> show for ℓ=0,1,2, respectively, the  values of
        S^(ℓ) versus time. Initially, the
        monopole is largest, but it is soon overtaken by the dipole, which continues
        to rise throughout the simulated period, whereas the monopole stagnates. The
        quadrupole continues to rise but is much smaller than either the monopole or
        dipole. 
        
        The dashed curve in Fig. <ref> shows the
        resulting estimate from equations (<ref>) and (<ref>) 
        of the rate of growth of (Δ
        E)^2 when one sets τ=t_N=8.6t_ I
        (Appendix <ref>). The  full curve in Fig. <ref> shows the
        actual rate of growth of (Δ E)^2. Although the shapes of the two
        curves are very different, at the latest time they agree in magnitude better
        than one has any right to expect given the approximations involved in the
        computation of the dashed curve.
        
        
        
        In  Appendix <ref> we quote the formulae for the diffusion
        coefficients in energy that are yielded by local scattering
        theory. When the second-order diffusion coefficient is integrated over the
        isochrone, we find
        
            D_2≡∫Ẹ  g(E)f(E)Ẹ^2ṭ=1.20lnΛ×10^-5σ^5 b,
        
         where σ^2≡ GM/b.  The theory predicts that when we average
        Δ E)^2/ṭ over all particles, we should find
        that
        
            (Δ E)^2ṭ=D_2.
        
        The dotted line in Fig. <ref> shows this
        prediction when we adopt lnΛ=3. 
        
        Interpretation of Fig. <ref> is not straightforward. With our chosen
        value of lnΛ, local scattering theory yields the correct order of
        magnitude for (Δ E)^2/ṭ but it predicts that the latter
        should be constant rather than falling. The hypothesis that relaxation is
        driven by large-scale fluctuations also yields the right order of magnitude
        for (Δ E)^2/ṭ and it correctly predicts that the latter's
        value should be a strong function of time rather than constant.
        Fig. <ref>, however, suggests that it should be a rising rather than a
        falling function of time.
        
        Careful consideration reveals that this prediction is likely be an artifact
        of how we have estimated the right side of equation (<ref>): we have
        taken the time derivatives of Φ_ℓ m to be their values divided by
        the characteristic time t_N. This procedure makes sense once the system has
        settled to an approximate steady state, but initially it is misleading
        because the fluctuations start from anomalously small amplitudes, so we
        estimate that they initially have very small time derivatives.
        Fig. <ref> shows clearly that the dipole's time derivative is
        initially large, presumably because over-dense regions are falling together,
        even though its initial magnitude is small. A better estimate of the initial
        time derivatives of the Φ_ℓ m would be the values they will acquire
        a time t_N later, divided by t_N.  If we liken the large-scale
        fluctuations to harmonic oscillators, they start with zero amplitude and
        maximum velocity, and consequently have large time derivatives.
        
        
        
        
        § DRIFT VELOCITIES IN ACTION SPACE
        
        
        We now use the  data to compute the velocities with which stars drift
        through action space. This is a more sophisticated and direct probe of
        relaxation than (Δ E)^2, which is just a global average of the
        extent to which stars shuffle parallel to a particular direction in action
        space and takes no account of the fact that most shuffles cancel one another:
        the distribution function f() is unchanged when two stars at similar
        action-space locations exchange energy but (Δ E)^2 is
        incremented.
         
        H18 show that in the case of a spherical system, action space, which
        is generically three-dimensional with Cartesian coordinates
        (J_r,J_z,J_ϕ), can be reduced to the two-dimensional space spanned by
        ≡(J_r,L), where L≡ J_z+|J_ϕ| is the magnitude
        of the angular momentum vector. The mass density of stars in this space is
        (2π)^3 times
        
            f_2≡ 2Lf,
        
        where f is the full DF, normalised such that M=(2π)^3∫^̣3 f.
        H18 present plots of the stellar fluxes in
        two-dimensional action space computed both from the BL equation and from
        local scattering theory. Their fluxes were normalised to the case of N=10^5
        stars in a cluster. In our comparisons below we renormalise their fluxes to
        clusters of 1000 stars.
        
        
        
        
        
        The flux  in the  plane is associated with a drift velocity
        
            ()≡() f_2().
        
        From each snapshot it is straightforward to compute each
        particle's value of  from its energy and angular momentum
        because for the isochrone we have an expression for the Hamiltonian
        
            H() = -(GM)^2/2[J_r + 12(L + √(L^2 +4GMb))]^2,
        
        so we can solve for J_r(E,L).
        We divide two-dimensional action space into cells of varying size and we
        compute the mean ⟨Δ⟩ of the
        changes in  experienced by stars in each cell during the last
        half of the simulations, i.e., a time Δ t=5t_N. The drift velocity at
        the centre of the cell is taken to be
        
            =ΔΔ t.
        
        The centre panel of Fig. <ref> shows the resulting velocities,
        while the top and bottom panels show the velocities that correspond to the
        flux vectors computed by H18 from the BL equation and local scattering
        theory, respectively. The length ℓ of an arrow encodes the speed to the flow
        via the non-linear scaling
        
            ℓ(GMb)^1/2=1.25×10^-3  asinh(V/V_0).
        
        where V_0=10^-5 GM/b. The length of an arrow is then a linear function
        of V for V≪ V_0 and a logarithmic function for V≫ V_0.
        
        The three panels differ considerably from one another. The BL vectors in the
        top panel have a very strong tendency to point towards the origin, while
        outside the bottom left corner, the local-scattering vectors in the bottom
        panel point resolutely upward. The  vectors are something of a
        compromise between these two extremes, consistent with the BL and
        local-scattering vectors capturing opposite halves of what actually happens
        in a cluster. A feature of the  vectors that is captured by the BL
        vectors is the tendency to converge on a point at low J_r and non-zero L.
        However, the point of convergence lies at substantially larger L in the BL
        panel than in the  panel. A striking feature of the  panel that
        is captured by neither the other panels is a strong tendency for vectors to
        align with lines of constant H, which run from upper left to lower right.
        The  and the BL vectors have similar lengths at all points, while the
        vectors from the local approximation are shorter, especially at low J_r.
        Thus the local approximation materially under-estimates the speed of
        relaxation.
        
        Given that H18 included only a handful of resonant pairings, it is inevitable
        that the BL vectors differ significantly from the  vectors. The
        conflict between the local-scattering and  vectors clearly confirms
        that important physics is missed by local-scattering theory.  That physics is
        the excitation of system-scale oscillations that is made possible by
        self-gravitation. Thus Fig. <ref> further confirms the essential
        message of H18 that the traditional theory of cluster evolution is highly
        unsatisfactory.
        
        
        
         §.§ Entropy generation
        
        
        In the central panel Fig. <ref>, the  velocity vectors appear
        to converge around the point =(0.4,0.2)√(GMb). This
        convergence implies that in this region f_2 is increasing and consequently
        the entropy density is decreasing. In this subsection we show that
        notwithstanding the local rise of f_2, the cluster's entropy is increasing
        as it must be.
        
        With f normalised to unit integral through action space, the  entropy of a
        cluster in equilibrium is
        
            S = - ∫^̣3   f ln(f h^3),
        
        where h is a constant with the dimension of Planck's constant  to ensure that
        the argument of the logarithm is dimensionless. We need to apply this formula
        to the case of an spherical system, when f(J_r,L). As in H18 we define as a
        third coordinate for action space the cosine β of an orbit's  inclination
        angle, and then we have
        that ^̣3=J̣_rLḶβ̣. Substituting this expression into equation
        (<ref>) and integrating over allowed values of β we obtain
        
            S = - 2∫^̣2  Lf ln(fh^3).
        
        Differentiating with respect to time, we obtain
        
            Ṣ/ṭ   = - 2∫^̣2  L f
            t[ln(fh^3)+1]   = - ∫^̣2  f_2
            t[ln(fh^3)+1],
        
        where the second equality uses equation (<ref>). Differentiation of the
        normalisation condition on f_2 leads to the conclusion that we can drop the +1
        from the square bracket. Then we use the BL equation to eliminate f_2/ t in favour of the flux  in  space to obtain
        
            Ṣṭ=∫^̣2 ∇_· ln(fh^3).
        
        The divergence theorem and the fact that in the case of an isotropic system
        the DF has the form f(H) allow us to obtain from this
        
            Ṣṭ   =-∫^̣2 ·∇_ln(fh^3)   =-∫^̣2 ·ḷṇfḤ.
        
        Given that ḷṇf/Ḥ<0, this equation states that entropy increases if
        on average the action-space flux  crosses contours of constant H
        outwards.  
        
        The upper panel of Fig. <ref> combines velocity vectors in action
        space with contours of constant H (in red). We see that in much of the
        space stars flow approximately along lines of constant H. In the lower left
        of the figure there is strong upward movement of stars and thus entropy
        creation, while in the upper right region stars are moving to lower energy,
        so entropy is being destroyed. The lower panel of Fig. <ref> shows the
        rate of entropy generation in a grey-scale plot, with the boundary between
        entropy creation and destruction marked by a blue contour. Entropy is
        destroyed in most of the depicted part of action space, but the typical
        creation rates far exceed the destruction rates, so overall net entropy is
        created, as is physically mandatory.
        
        
        
        § CONCLUSIONS
        
        
        A suite of 10 000 simulations of different realisations of 1000-particle
        isochrone clusters demonstrates that the neglect of self-gravity by the
        traditional theory of cluster evolution leads to a seriously incomplete
        picture of how clusters relax. Over the short durations (∼3 crossing
        times) of our simulations, self-gravity amplified the original Poisson noise
        by an order-of magnitude. Figs. <ref> and <ref>
        show that the amplification was still some distance from attaining
        statistical equilibrium by the end of our simulations, so our estimates of
        the importance of self-gravity must be on the low side. An important
        implication of this finding is that when a simulations is started by randomly
        sampling a DF, it starts from an anomalously regular configuration that would
        not occur in nature. Several crossing times are required for the noise level
        in the system to settle to a stationary level.
        
        Our experiments were stimulated by the finding of H18 from the BL equation
        that self-gravity substantially increases the relaxation rate of clusters.
        Their conclusion was in large measure based on the demonstration of strong
        excitation by Poisson noise of a cluster's fundamental dipole mode, a result
        presaged by the far-sighted work of <cit.>.
        Our direct-summation N-body simulations confirm that the dipole mode, in
        which the centre moves in opposition to the halo, is powerfully excited. 
        
        We have investigated the extent to which system-scale fluctuations in the
        gravitational potential can account for diffusion in energy, which is
        traditionally ascribed to the cumulative effect of stars scattering off other
        stars. This analysis is tricky because the fluctuations are far from a
        stationary random process during our brief simulations. However, 
        Fig. <ref> provides strong evidence that these fluctuations are at
        least as effective as star-star scattering in causing the energies of stars
        to evolve.
        
        Clusters evolve through stars changing their locations in action space. H18
        showed that the action space of a spherical system can be reduced to the
        two-dimensional space spanned by radial action J_r and the magnitude L of
        the angular-momentum vector. Within cells in this space we have determined the
        mean drift velocity of stars in the N-body simulations and compared these
        with predictions from
        the BL equation and local-scattering theory. Fig. <ref> shows that
        the complex action-space flow recovered from the N-body simulations differs
        from those predicted from the BL equation on the one hand and
        local-scattering theory on the other in ways that suggest that each of the
        competing predictions captures roughly a half of the whole picture. The
        theoretical predictions differ from one another in that the prediction obtained from the BL
        equation by H18 is known to be incomplete, and is in principle capable of
        improvement, while that derived from local-scattering theory is definitive and
        incapable of modification.
        
        Our final conclusion is, therefore, that local-scattering theory is of very
        limited value. With a judiciously chosen value of its free parameter,
        lnΛ, it can be forced to yield a reasonable value for the relaxation
        rate. It also correctly predicts that stars in the cluster's halo gain energy
        at the expense of stars in the core. But it is seriously misleading as
        regards the basic physics of relaxation and predicts the flow of stars
        through action space wrongly. 
        
        It seems then that to continue work on cluster
        evolution based on local-scattering theory would be a mistake. Future work
        should either rely on N-body simulations, or employ the BL equation.  For the
        latter to become a realistic option, one would have to develop a much more
        efficient code for the computation of diffusion coefficients than that used
        by H18. the extent to which this is possible is unclear, but it is surely
        worth a try. 
        
        We employed the direct-summation code  because we were
        anxious to exclude any suggestion that our simulations did not handle
        star-star scattering accurately. If, as now seems likely, relaxation is
        largely driven by system-scale fluctuations, costly direct-summation could be
        replaced by a less costly `collisionless' method of force determination, such as tree
        summation or Fourier transforms. This is another promising direction for
        future work.
        
        
        
        
        
        mn2e
        
        § UNITS
        
        
        The isochrone has natural length- and time-scales, b and (GM/b^3)^-1/2.
        Initial conditions for an N-body simulation of an equilibrium system also
        define natural length- and time-scales, r_N and t_N. Here we determine
        the relations between these two systems when the initial conditions are those
        of an isochrone.
        
        From the initial conditions one can compute the cluster's potential energy
        
            W≡ -12∑_i
            jGm^2|_i-_j|,
        
         where m is the mass of each of the N stars in the cluster.
        For any given value of  r_N, we have
        
            W=-
            GM^22r_N1|-'|,
        
        where ≡/r_N is a dimensionless vector. rescales the coordinates by choosing r_N such
        that the expectation value in equation (<ref>) is unity. Then the
        cluster's potential energy is 
        
            W=-GM^2 2r_N.
        
        
        From the initial conditions we can compute the mean-square speed
        
            σ^2≡1 N∑_i=1^Nv_i^2=v^2.
        
        By the virial theorem, σ^2
        is related to W by
        
            2K=Mσ^2=-W=GM^22r_N.
        
        With t_N any time, and t≡ t/t_N and
        
            ≡%̣ṣ%̣ṣt= t_N r_N,
        
        dimensionless times and velocities, the mean-square of the dimensionless
        velocities is
        
            σ^2=GM 2r_N^3 t_N^2.
         scales times with t_N chosen such that
        σ^2=1/2, so
        
            t_N=(r_N^3 GM)^1/2.
        
        
        The potential energy of the isochrone is
        
            W=-(π2-43)GM^22b.
        
        Comparing this with equation (<ref>) we find
        
            b=0.2375 r_N.
        
        From equation (<ref>) we have
        
            t_N=8.642(b^3 GM)^1/2.
        
        The period of orbits confined to the core is P=2π√2 t_ I≃ t_N.
        
        
        
        § DIFFUSION COEFFICIENTS IN ENERGY
        
        
        Here we list the formulae from <cit.> and <cit.> from which we
        have computed the coefficients that govern diffusion in energy in the
        approximation of local scattering.
        
        In a system of total mass M comprising N stars the density of states is
        
            g(E)=(4π)^2∫_0^r_ max(E)ṛ r^2v(E,r)
        
        where v(E,r)=√(2(E-Φ)) and r_ max(E) is the radius at which
        v vanishes. A related quantity, the density of states weighted by v^2, is
        given by
        
            h(E)=(4π)^2∫_0^r_ max(E)ṛ r^2v^3.
        
        With the DF normalised so that ∫Ẹ g(E)f(E)=1, the diffusion coefficients are 
        
            Ẹṭ   =-(4π GM)^2 NlnΛ[
            1 g(E)∫_-∞^EẸ' g(E')f(E')       -∫_E^∞Ẹ' f(E')
            ]Ẹ^2ṭ   =2(4π GM)^2 NlnΛ[
            1 g(E)∫_-∞^EẸ' h(E')f(E')       +h(E) g(E)∫_E^∞Ẹ' f(E')
            ].
        
        Fig. <ref> shows the diffusion coefficients. The first-order
        coefficient changes sign from positive to negative as the binding energy
        decreases. A useful check on the numerics is provided by the requirement of
        energy conservation that ∫Ẹ g(E)Ẹ/ṭ=0.
        
         
        
         §.§ Derivative of (Δ E)^2_t
        
            (Δ E)^2_t+δ t-(Δ E)^2_t   =
            (∑_i^Nδ_i+∑_j^nδ_j)^2-(∑_i^Nδ_i)^2   =2(∑_i^Nδ_i)(∑_j^nδ_i)+(∑_j^nδ_j)^2
        
        Dividing by δ t and taking the limit δ t→0 and assuming the change in the last interval
        is uncorrelated with the changes up to that instant,
        
            (Δ E)^2_tṭ=2E_t-E_0Ẹṭ+Ẹ^2ṭ.
        
        By E conservation E_t-E_0=0, so 
        
            (Δ E)^2_tṭ=Ẹ^2ṭ.
        
        
        
        
        
        Galaxies and star clusters evolve over cosmic time through random exchanges
        of energy between stars. Fluctuations in the system's gravitational field
        around its mean-field value are the agents of these exchanges. For
        generations it was thought that the fluctuations could be modelled by the
        moving Keplerian fields of individual stars, so stars interacted by
        scattering each other. A new formalism reveals that this picture is
        fundamentally wrong. Stars communicate with each other by broadcasting
        narrow-band signals that are carried by collection oscillations of the whole
        system, and absorbed by stars that happen to resonate at the broadcast
        frequency. The relevant fluctuations comprise large-scale collective modes of
        the system. I'll report on the application of these ideas to galactic discs
        and to star clusters.
        
        
        
        Our initial conditions are
        unnaturally spherically symmetric: in these conditions symmetry around the
        mean-field centre is broken only by Poisson noise. After the the N-body
        equations of motion have been integrated for a few crossing times,
        self-gravity has amplified this noise by a significant factor, and σ_s
        has grown accordingly.
        
        
         This mean If the kinematics of the potential centre were a statistically stationary
        process, the curves would all fall on top of the bottom curve. 
        
        
        In detail, for a grid of radii we compute the autocorrelation of the
        dipole potential
        
            C_1(r,t,τ)≡1 T∫_t^t+Tṭ _1(r,t)·_1(r,t+τ)
        
         where the dot product is between the vectors defined by each dipole. Then we
        average C_1(r,t,τ) over r by averaging the values it takes at the radii of
        particles. Then we compare plots of (Δ E)^2 and
        (τ/τ_0) C_1(t,0).
        
        
        The relaxation
        time is the time required for Δ E to become comparable to E, i.e.
        
            t_ relax=E^2τ_0(Ṗτ_0)^2
        
        
        The full curve in Fig. <ref> shows the mean-square change in
        particle energy as a function of time. Equation (<ref>) predicts
        that this should be a straight line, and over the second half of the
        simulated period the full curve is  reasonably straight. 
        
        The dashed curve in
        Fig. <ref> shows
        
            (ΔΦ_1)^2_τ≡1
            N∑_i=1^N|_1^(i)(τ)-_1^(i)(0)|^2,
        
        where _1^(i)(τ) is the vector defined by the dipole contribution to
        the potential at the radius of the i^ th particle. Given that
        _1 has vanishing expectation value, (ΔΦ)^2 should be a
        useful measure of the amplitude of the fluctuating dipole. 
        
        In equation (<ref>), τ_0 is a measure of the time over which
        Φ/ t evolves coherently, so Ṗτ_0≃ΔΦ is
        essentially the characteristic amplitude of the fluctuations in Φ
        averaged over the locations of particles. All multipoles will contribute to
        (Ṗτ_0)^2, but the contribution from the dipole will be
        ≃(ΔΦ_1)^2. Hence a dimensionless measure of the contribution
        of the dipole to relaxation is
        
            D≡ττ_0(ΔΦ)^2(Δ E)^2_τ.
        
        
        
        Now consider the mean-square of the change in the perturbing potential at the
        location of a star between two times, t,t':
        
            (ΔΦ)^2_t,t'≡[Φ(r,t)-Φ(r',t')]^2=   ⟨|∑_ℓ mΦ_ℓ m(r,t)Y_l^m(θ,ϕ)
            -∑_ℓ mΦ_ℓ m(r',t')Y_ℓ'^m'(θ',ϕ')|^2⟩.
        
        When we multiply out the rhs, we get a term
        
            ∑_ℓ mℓ'm'Φ_ℓ m(r,t)Φ_ℓ'm'(r,t)Y_l^m(θ,ϕ)Y_ℓ'^m'(θ,ϕ)
        
        added to a term that differs only in having primes on t,r,θ,ϕ. Equation
        (<ref>) enables us to simplify these terms to
        
            14π∑_ℓ m(|Φ_ℓ m(r,t)|^2+|Φ_ℓ m(r,t')|^2).
        
        The product in equation (<ref>) also yields
        
            -2∑_ℓ mℓ'm'Φ_ℓ m(r,t)Φ_ℓ'm'(r',t')Y_l^m(θ,ϕ)Y_ℓ'^m'(θ',ϕ')
        
        As in the calculation of (Δ E)^2, we may argue that the
        expectation of the product of spherical harmonics is small unless
        |t-t'|τ_0. If fact, when |t-t'| is small, the expectation of the
        product of spherical harmonics approximates the rhs of equation
        (<ref>) and this term essentially cancels the sum (<ref>),
        as it must given that the lhs of equation (<ref>) vanishes for
        t=t'. In summary, for |t-t'|>τ_0, we have
        
            (ΔΦ)^2_t,t'≃14π∑_ℓ m(|Φ_ℓ m(t)|^2+|Φ_ℓ m(t')|^2).
        
        This equation implies that (ΔΦ)^2_t,t'
        would be independent of t,t'
        if Φ were a stationary random process. In fact it would just be a
        measure of the amplitude of the fluctuating component of the potential.
        Indeed, we can use this estimate in equation (<ref>) to gauge the
        rate at which (Δ E)^2_τ grows when the fluctuations become
        stationary:  with Φ/
        t∼Φ/τ_0, equation (<ref>) yields
        
            (Δ E)^2_τ∼ττ_0(ΔΦ)^2_τ,0.
        
        
        
         for the early evolution of [E(t)-E(0)]^2. At later times the
        quantity on the left of this relation should saturate on account of the
        non-vanishing value of the diffusion coefficient Ẹ/ṭ. The
        relaxation time t_ relax is of order the time required for the linear
        increase predicted by equation (<ref>) to yield the
        value[From the expression (<ref>) for the potential
        energy of the isochrone, we have that the average kinetic energy of a star is
        ∼0.06σ^2.]0.0036σ^4 that is comparable to the saturation
        value of the left side. Thus
        
            t_ relax≈0.0036σ^4D_2=100bσ.
