Secular diffusion in discrete selfgravitating tepid discs II. Accounting for swing amplification via the matrix method^{⋆,}^{⋆⋆}
^{1}
Institut d’Astrophysique de Paris and UPMC, CNRS (UMR 7095),
98bis Boulevard Arago,
75014
Paris,
France
email:
fouvry@iap.fr
^{2}
Institute of Astronomy & KICC, University of
Cambridge, Madingley
Road, Cambridge,
CB3 0HA,
UK
^{3}
Rudolf Peierls Centre for Theoretical Physics, University of
Oxford, Keble Road,
Oxford
OX1 3RH,
UK
^{4}
Laboratoire de Physique Théorique (IRSAMC), CNRS and UPS,
Univ. de Toulouse,
31062
Toulouse,
France
Received: 24 July 2015
Accepted: 25 August 2015
The secular evolution of an infinitely thin tepid isolated galactic disc made of a finite number of particles is investigated using the inhomogeneous BalescuLenard equation expressed in terms of angleaction variables. The matrix method is implemented numerically in order to model the induced gravitational polarisation. Special care is taken to account for the amplification of potential fluctuations of mutually resonant orbits and the unwinding of the induced swing amplified transients. Quantitative comparisons with Nbody simulations yield consistent scalings with the number of particles and with the selfgravity of the disc: the fewer the particles and the colder the disc, the faster the secular evolution. Secular evolution is driven by resonances, but does not depend on the initial phases of the disc. For a Mestel disc with Q ~ 1.5, the polarisation cloud around each star boosts its secular effect by a factor of a thousand or more, accordingly promoting the dynamical relevance of selfinduced collisional secular evolution. The position and shape of the induced resonant ridge are found to be in very good agreement with the prediction of the BalescuLenard equation, which scales with the square of the susceptibility of the disc. In astrophysics, the inhomogeneous BalescuLenard equation may describe the secular diffusion of giant molecular clouds in galactic discs, the secular migration and segregation of planetesimals in protoplanetary discs, or even the longterm evolution of population of stars within the Galactic centre. It could be used as a valuable check of the accuracy of Nbody integrators on secular timescales.
Key words: galaxies: evolution / galaxies: kinematics and dynamics / galaxies: spiral / diffusion / gravitation
Appendices are available in electronic form at http://www.aanda.org
A copy of the linear matrix response code is only available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/qcat?J/A+A/584/A129
© ESO, 2015
1. Introduction
Because most of the stars are born in a stellar disc, galactic astronomy has striven to understand the dynamics and evolution of discs over cosmic time. Stars in cool discs are born on the almost closed orbits of their parent gas clouds, which means that they are confined to a very small volume of phase space. There is therefore a massive opportunity for increasing the entropy of the disc’s phasespace distribution by spreading stars orbits to a wider range of eccentricities and inclinations while conserving the system’s energy. Much empirical evidence of such heating has now been discovered: the random velocities of coeval cohorts of stars increases with the cohort’s age (Wielen 1977; Aumer & Binney 2009); the velocity distribution at the Sun contains streams of stars (Dehnen 1998); and the density of stars in phase space exhibits elongated features (resonance ridges) consistent with a resonant condition (Sellwood 2010; McMillan 2013).
Such longterm evolution of the system’s distribution function (DF) implies that the disc’s gravitational potential cannot have remained static and axisymmetric. Axisymmetric discs support spiral waves that can move through the disc; leading waves unwind into trailing waves and are amplified at corotation. Eventually each wave Landau damps by transferring its energy to stars that resonate with it: it heats the disc by increasing the eccentricity of resonating stars. In a disc made of a finite number of stars, shot noise is a natural excitation mechanism for launching such waves. When the disc is cool, fluctuations of the potential induced by these discrete encounters may be strongly amplified (Kalnajs 1972), so a weak stimulus can produce significant heat. Capturing this secular heating by shot noise requires an apparatus able to follow the dynamics of spiral waves from launch as leading waves through amplification and to thermalisation at resonances to determine their secular dynamical impact. In this paper we demonstrate how the BalescuLenard formalism provides an appropriate framework for taking such interactions into account.
The kinetic theory of stellar systems was initiated by Jeans (1929) and Chandrasekhar (1942) in the context of hot spherical stellar systems, such as elliptical galaxies and globular clusters for which the gravitational susceptibility can safely be neglected, as opposed to cold selfgravitating galactic discs. More generally, in astrophysics, the secular diffusion of giant molecular clouds in galactic discs, the secular migration and segregation of planetesimals in protoplanetary or debris discs, or even the longterm evolution of population of stars within the Galactic centre are all processes for which it is relevant to quantify the dynamical effect of gravitationally amplified potential fluctuations induced by the finite number of elements involved. More than 55 years ago, Balescu (1960) and Lenard (1960) developed a rigorous kinetic theory taking collective effects into account and obtained the corresponding kinetic equation for plasmas, the BalescuLenard equation. More recently Heyvaerts (2010) and Chavanis (2012) have transposed the corresponding nonlinear kinetic equation to the angleaction variables that are the appropriate variables to describe spatially inhomogeneous multiperiodic systems. The corresponding inhomogeneous BalescuLenard equation accounts for selfdriven orbital secular diffusion of a selfgravitating system induced by the intrinsic shot noise due to its discreteness. The formal transposition from positionvelocity to angleaction implies that the secular interactions need not be local in space: they only need to correspond to gravitationally amplified longrange correlations and resonances, which are indeed the driving mechanism for the secular evolution of isolated astrophysical discs via angular momentum redistribution (LyndenBell & Kalnajs 1972).
The BalescuLenard equation is valid at the order 1 /N in an expansion of the dynamics in terms of this small parameter, where N ≫ 1 is the number of stars. Therefore, it takes finiteN effects into account and describes the evolution of the system on a timescale of about Nt_{D}, where t_{D} is the dynamical time. For selfgravitating systems, the collective effects are responsible for an antishielding which tends to increase the effective mass of the stars, hence reducing the relaxation time. When the system is cold, each particle is dressed by the very strong gravitational polarisation it induces. Therefore secular effects can occur on much shorter timescales than one would naively expect, making the system behave as if it were made up of only N_{eff} ~ N/10^{few} particles. The purpose of this paper is to quantify this effect for stable but strongly susceptible galactic discs.
The BalescuLenard formalism has seldom been applied in its prime context, but only in various limits where it reduces to simpler kinetic equations (Landau 1936; Vlasov 1938; Chandrasekhar 1942; Rosenbluth et al. 1957). Weinberg (1993) presents an interesting first implementation, though in a somewhat oversimplified cartesian geometry. This formalism is quite unique in accounting for the nonlinear evolution of discs and galaxies on secular timescales. While potentially probing similar processes, Nbody simulations should be scrutinised carefully in this regime, as shadowing (Quinlan & Tremaine 1992) may, over many orbital times, impact resonant interactions. Nevertheless, Nbody simulations have been shown to more or less reproduce growth rates of discs on dynamical timescales (see e.g. Sellwood & Evans 2001, and references therein, together with Appendix C); qualifying them quantitatively on secular timescales is now within reach of the BalescuLenard formalism.
A companion paper, (Fouvry et al. 2015a, hereafter Paper I) presented a simple and tractable quadrature for the BalescuLenard drift and diffusion coefficients under the assumption that the transient response of the disc can be described by tightly wound spirals via the WKB and epicyclic approximations. These simple expressions have provided insight into the physical processes at work during the secular diffusion of selfgravitating discrete discs. When applied to the secular evolution of an isolated stationary selfgravitating Mestel disc, it identified the importance of the corotation resonance in the inner regions of the disc leading to a regime with both radial migration and heating, in qualitative agreement with numerical simulations.
Yet, the tightly wound approximation is quantitatively questionable when transient spirals unwind. Indeed Paper I found a discrepancy between the predicted secular evolution timescale and the measured one, which might be driven by the incompleteness of the WKB basis, since such a basis can only correctly represent tightly wound spirals. It considered only resonances between orbits that are close to one another in position space, and did not allow remote orbits to resonate, nor wave packets to propagate between such nonlocal resonances. Yet, the seminal works from Goldreich & LyndenBell (1965), Julian & Toomre (1966), Toomre (1981) show that any leading spiral wave undergoes significant amplification during its unwinding to a trailing wave. Because it involves the unwinding of spirals this mechanism is not captured by the WKB formalism of Paper I.
In this paper, we avoid such approximations by using the matrix method (Kalnajs 1976) to estimate the gravitational amplification of the secular response and the role of nonlocal resonances. From this we compute numerically the drift and diffusion coefficients that appear in the BalescuLenard equation. We then compare those predictions to crafted sets of numerical experiments, allowing us to estimate ensemble averaged secular responses of a sizeable number of simulations as a function of the total number of particles N. Such ensemble averaging allows us to make robust predictions for the Nscaling of the secular response and its dependence on halo to disc mass fraction, hence probing the secular importance of gravitational polarisation.
The paper is organised as follows. Section 2 briefly presents the content of the inhomogeneous BalescuLenard equation. Section 3 presents our implementation of the matrix method to compute the diffusion equation for an isolated selfgravitating tapered Mestel disc. Section 4 computes numerically the exact drift and diffusion coefficient in action space for such a truncated Mestel disc, and compares the divergence of the corresponding flux density to the initial measured rate of change of the distribution function. Section 5 presents our Nbody simulations and compares scaling of the flux with the number of particles and the fraction of mass in the disc. Finally, Sect. 6 wraps up. Appendix A presents the relevant biorthogonal basis function. Appendix C validates the response matrix method and the Nbody integrator while matching growth rates and pattern speeds of unstable Mestel discs. Appendix D investigates the roles of selfgravity and basis completeness. Appendix E describes the sampling strategy for the initial distribution. Appendix G presents briefly the available online codes.
2. The inhomogeneous BalescuLenard equation
We intend to describe the longterm evolution of a system made of N particles. We assume that the gravitational background ψ_{0} of the system is stationary and integrable, and associated with the Hamiltonian H_{0}. As a consequence, one can always remap the physical spacecoordinates (x,v) to the angleaction coordinates (θ,J) (Goldstein 1950; Born 1960; Binney & Tremaine 2008). We define the intrinsic frequencies of motions along the action torus as (1)Within these new coordinates, one has that along the unperturbed trajectories the angles θ are 2πperiodic, evolving with the frequencies Ω, whereas the actions J are conserved. We assume that the system is always in a virialised state, so that its distribution function (DF) can be written as a quasistationary DF of the form F = F(J,t), satisfying the normalisation constraint , where M_{tot} is the total mass of the system. On secular timescales, this isolated DF evolves under the effect of stellar “encounters” (finiteN effects). Such a collisional longterm evolution is descrided by the inhomogeneous BalescuLenard equation (Heyvaerts 2010; Chavanis 2012) which reads (2)where are the dressed susceptibility coefficients, d is the dimension of the physical space, μ = M_{tot}/N is the mass of the individual particles, and where we used the shortened notation Ω_{i} = Ω(J_{i},t). Since the r.h.s. is the divergence of a flux, this diffusion equation conserves the number of stars. One should also note the resonance condition encapsulated in the Dirac delta δ_{D}(m_{1}·Ω_{1} − m_{2}·Ω_{2}), with the integration over the dummy variable J_{2} scanning for points where the resonance condition is satisfied. Note importantly right away that Eq. (2) scales like (since μ ∝ 1 /N), so that increasing N or increasing the heat content of the disc have the same effect. For a more detailed discussion on the content of the BalescuLenard equation, see Paper I.
In order to solve the nonlocal Poisson equation, we follow Kalnajs’ matrix method (Kalnajs 1976), so that we introduce a complete biorthonormal basis of potentials and densities ψ^{(p)}(x) and ρ^{(p)}(x) such that (3)The dressed susceptibility coefficients appearing in Eq. (2) are then given by (4)where I is the identity matrix and is the response matrix defined as (5)In the previous expression, we introduced as the Fourier transform in angles of the basis elements ψ^{(p)}(x) using the convention that the Fourier transform of a function X(θ,J) is given by (6)In order to ease the understanding of the BalescuLenard Eq. (2), one may rewrite it as an anisotropic FokkerPlanck equation by introducing the relevant drift and diffusion coefficients. Indeed, Eq. (2) can be put into the form (7)where A_{m1}(J_{1}) and D_{m1}(J_{1}) are respectively the drift and diffusion coefficients associated with a given resonance m. We note that they both depend secularly on the DF F, but we omit this dependence to simplify our notations. The drift coefficients are given by(8)and the diffusion coefficients are given by (9)One can also introduce the total flux of diffusion ℱ_{tot} as (10)so that the BalescuLenard equation from Eqs. (2) and (7) takes the explicitly conservative form (11)
3. The Matrix diffusion equation
When computing the BalescuLenard diffusion and drift coefficients, three main difficulties have to be addressed. First, one must build the mapping (x,v) → (θ,J), because the drift and diffusion coefficients are associated with a diffusion in action space. The second difficulty follows from the nonlocality of Poisson’s equation and the estimation of the response matrix from Eq. (5). Indeed, as noted in Eq. (3), the matrix relies on potential basis elements ψ^{(p)} which must be integrated over the whole action space with functions possessing a pole 1 / (ω − m·Ω). This cumbersome and difficult evaluation has to be performed numerically, along with the matrix inversion needed to estimate the susceptibility coefficients from Eq. (4). Finally, the third difficulty arises from the resonance condition δ_{D}(m_{1}·Ω_{1} − m_{2}·Ω_{2}), which involves determining how orbits may resonate one with another. In Paper I we relied on the epicyclic approximation to build the angleaction mapping and on a WKB basis to treat gravity locally in order to solve these issues while obtaining tractable though approximate expressions.
For a 2D axisymmetric potential, one can define explicitly the actions of the system. Following LyndenBell & Kalnajs (1972), Tremaine & Weinberg (1984), the two natural actions of the system are given by a quadrature and an identity (12)where r_{p} and r_{a} are respectively the pericentre and the apocentre of the trajectory, while E and L are the energy and angular momentum of the star. The first action J_{1} encodes the amount of radial energy of the star, so that J_{r} = 0 corresponds to circular orbits. The second action J_{2} is the angular momentum L of the star. One can then define the two intrinsic frequencies of motion Ω_{1} = κ associated with the radial oscillations and Ω_{2} = Ω_{φ} associated with the azimuthal oscillations. Indeed, one has (13)while the azimuthal frequency Ω_{2} can then be determined via the relation (14)At this stage, one should note that various coordinates can be used to represent the 2D action space. Indeed, once the background potential ψ_{0} is known, one has the bijections (r_{p},r_{a}) ↔ (E,L) ↔ (J_{r},J_{φ}). As a consequence, any orbit can equivalently be represented by the set of the pericentre and apocentre (r_{p},r_{a}) or by its actions (J_{1},J_{2}). However, determining the actions associated with one set (r_{p},r_{a}) only requires the computation of a 1D integral as in Eq. (12), whereas determing the pericentre and apocentre associated with a set of actions (J_{1},J_{2}) requires the inversion of the same nontrivial relation. Because the peri/apocentres are the two roots of the equation 2(E − ψ_{0}(r)) − L^{2}/r^{2} = 0, one also immediately obtains that for a given value of r_{p} and r_{a}, the energy E and the angular momentum L of the orbit are given by (15)where we used the shortening notations ψ_{p/a} = ψ_{0}(r_{p/a}). Therefore, in the upcoming calculations, we use (r_{p},r_{a}) as the representative variables of the 2D action space.
3.1. The basis elements
The expressions (4) and (5) of the susceptibility coefficients and the response matrix require the introduction of 2D potentialdensity basis elements. The 2D potential basis elements ψ^{(p)} that we consider depend on two indices spanning the two degrees of freedom so that one has (16)where is a real radial function and (R,φ) are the usual polar coordinates. The associated surface densities elements are of the form (17)where is a real radial function. The basis elements therefore depend on two indices ℓ ≥ 0 and n ≥ 0. In all the numerical calculations, we used the radial functions from Kalnajs (1976), which are recalled in Appendix A.
The next step is then to determine the Fourier transform with respect to the angles θ of the basis elements. Indeed, to compute the response matrix from Eq. (5), one has to compute , where the resonance vector is given by m = (m_{1},m_{2}). Following Eq. (6), it is given by (18)From LyndenBell & Kalnajs (1972), the angles θ_{1} and θ_{2} associated with the actions from Eq. (12) are given by (19)where is a contour starting from the pericentre r_{p} and going up to the current position r = r(θ_{1}) along the radial oscillation. Following the notations from Tremaine & Weinberg (1984), one can straightforwardly show that Eq. (18) takes the form (20)where is given by (21)In Eq. (21), the boundaries of the integral are given by the pericentre r_{p} and apocentre r_{a} associated with the action J. Such an expression underlines the reason why (r_{p},r_{a}) appear naturally as good coordinates to describe the 2D action space. One can note that Eq. (21) involves an integral over r thanks to the change of variables θ_{1} → r, which satisfies (22)In Eq. (21), θ_{1} [ r ] and (θ_{2} − φ) [ r ] only depend on r via the mappings from Eq. (19). Provided that is a real function, the coefficients are always real. Because these coefficients involve two intricate integrals, they are numerically expensive to compute. However by parity, they obey , which offers a significant reduction of the number of coefficients to compute.
3.2. Computation of the response matrix
We now have all the elements required to compute the response matrix from Eq. (5). In its definition, one should note the presence of an integral over the dummy variable J, which, as discussed previously, is performed in the 2D (r_{p},r_{a})space. The first step is to go from J = (J_{1},J_{2}) to (E,L). The Jacobian of this transformation is given by (23)so that one immediately has dJ_{1}dJ_{2} = dEdL/Ω_{1}. Given the expression (20) of the 2D Fourier transformed basis elements, the response matrix may be written as (24)where the sum on m_{2} has been executed using the Kronecker delta from Eq. (20). Moreover, we dropped the conjugate over since the latter is always real. We may now perform the change of variables (E,L) → (r_{p},r_{a}), so as to rewrite Eq. (24) as (25)where the functions and are respectively given by (26)and (27)It is important to note here that the response matrix is diagonal with respect to the ℓ^{p} and ℓ^{q} indices so that each ℓ may be treated independently. The definition of the function g from Eq. (26) involves the Jacobian ∂(E,L) /∂(r_{p},r_{a}) of the transformation (E,L) → (r_{p},r_{a}) which can be immediately computed from the expressions (15) of E = E(r_{p},r_{a}) and L = L(r_{p},r_{a}). Moreover, in some situations, the DF F = F(J), may also rather be defined as F = F(E,L). It is straightforward to show that one has (28)
3.3. Subregion integration
The next step in the calculation is to perform the remaining integration over (r_{p},r_{a}) from Eq. (25). Because of the presence of the resonant pole , such a numerical integration has to be performed carefully. We cut the integration domain (r_{p},r_{a}) into various subregions indexed by i. The ithregion is centred around the position and corresponds to the square domain such that and , where Δr corresponds to the size of the subregions. The smaller Δr, the better the approximated estimations of the response matrix. Within the ithregion, one can write firstorder Taylor expansions of the functions g and h from Eqs. (26) and (27) around the centre of the region such that (29)where for convenience we shortened the index dependences from Eqs. (26) and (27). The coefficients , and (similarly for h) are given by (30)In the numerical implementation, the coefficients involving partial derivatives are estimated by finite differences, so that one has for instance (31)which allows us to minimise the number of evaluations of g required. The approximated integration on each subregion can then be performed and takes the form (32)where ℵ is an analytical function which only depends on the coefficients obtained in the limited developments from Eq. (29). In order to have a welldefined integral, we added an imaginary part η> 0 to the temporal frequency ω, so that ω = ω_{0} + iη. When looking for unstable modes in a disc, this imaginary part η corresponds to the growth rate of the mode. It is also crucial to note here that one always has , , , and similarly , , . The effective computation of the function ℵ is presented in Appendix B. Thanks to Eq. (32), the expression (25) becomes (33)In the previous expression, in order to effectively compute numerically the sum on m_{1}, we introduce a bound , so that the sum is only reduced to . Because of the requirement to truncate the action space in various subregions as in Eq. (33), the computation of the response matrix still remains a daunting task, to ensure appropriate numerical convergence. In Appendix C, we detail the validation of our implementation of the response matrix calculation by recovering known unstable modes of truncated Mestel discs (Zang 1976; Evans & Read 1998b; Sellwood & Evans 2001). Once the response matrix is known, the determination of the dressed susceptibility coefficients from Eq. (4) involves a straightforward summation^{1}.
3.4. Critical resonant line
The resonance condition encapsulated in the Dirac delta δ_{D}(m_{1}·Ω_{1} − m_{2}·Ω_{2}) generates an additional difficulty in the calculation of the BalescuLenard drift and diffusion coefficients from Eqs. (8) and (9). Recall the definition of the composition of a Dirac delta, and a function (Hörmander 2003), which in a ddimensional setup takes the form (34)where g^{1}(0) = {x  g(x) = 0 } is the hypersurface of dimension (generically) (d − 1) defined by the constraint g(x) = 0, and dσ(x) is the surface measure on g^{1}(0). We have also naturally defined the euclidean norm of the gradient of g as . Here we assumed that the resonance condition associated with the function g(J_{2}) = m_{1}·Ω_{1} − m_{2}·Ω_{2} is nondegenerate, so that ∀x ∈ g^{1}(0),  ∇g(x)  > 0, which also ensures that the dimension of the set g^{1}(0) is (d − 1). One should note that this degeneracy condition is not satisfied by the Keplerian or harmonic potentials. Because we are considering an infinitely thin disc, the dimension of the physical space is given by d = 2, so that the set g^{1}(0) takes the form of a curve γ, that we call the critical resonant curve. Generically, it takes the form of an application of the type (35)One can then immediately compute the r.h.s. of Eq. (34) via (36)where we have naturally defined .
As noted in Eq. (15), using the pericentres and apocentres (r_{p},r_{a}), given the Jacobian from Eq. (23) and proceeding in the same way as in Eq. (25) for the response matrix, one may rewrite the drift and diffusion coefficients from Eqs. (8) and (9) as (37)and (38)where the functions and are respectively defined as(39)and (40)For a given value of J_{1}, m_{1} and m_{2}, and introducing ω_{1} = m_{1}·Ω_{1}, we define the critical curve γ_{m2}(ω_{1}) as (41)The expressions (37) and (38) of the drift and diffusion coefficients immediately become (42)where the resonant contribution  ∇(m_{2}·Ω_{2})  is defined as (43)The derivatives of the intrinsic frequencies with respect to r_{p} and r_{a} appearing in Eq. (43) are computed as in Eq. (31) using finite differences. Once the critical lines of resonances have been determined, the computation of the drift and diffusion coefficients from Eq. (42) is straightforward, so that the full secular diffusion flux ℱ_{tot} from Eq. (10) may be determined.
4. BalescuLenard flux divergences
We may now illustrate how the previous computations of the response matrix and the BalescuLenard drift and diffusion coefficients can be used to recover some results obtained in wellcrafted numerical simulations of galactic discs. Indeed, Sellwood (2012; hereafter S12), studied the longterm evolution of an isolated stable truncated Mestel disc (Mestel 1963). After letting the disc evolve for hundreds of dynamical times, S12 observed a secular diffusion of the disc’s DF in action space, through the spontaneous generation of transient spiral waves. The most striking result of this evolution is given in Fig. 7 of S12, which exhibits the late time formation of a resonant ridge in action space along a specific resonant direction. Such diffusion features observed in the late evolution of an isolated, stable, and discrete system are thought to be signatures of a secular evolution induced by finiteN effects, as described by the BalescuLenard formalism. Because the system is made of a finite number N of pointwise particles, it undergoes (longrange) resonant encounters leading to an irreversible secular evolution. In order to investigate such a collisional evolution, Paper I applied the WKB limit of the BalescuLenard formalism to S12’s simulation. While most of the secular diffusion was qualitatively recovered, there remained a significant timescale discrepancy, since the typical timescale diffusion predicted by this approach was typically 10^{3} times too slow compared to the observations made in S12. The use of a nonlocal basis such as Eq. (16) and the numerical computation of the response matrix from Eq. (24) allows us to incorporate in the present paper these previously ignored contributions from the WKB approach. In the upcoming sections, we therefore present briefly the disc considered by S12 and our determination of the secular diffusion flux predicted by the BalescuLenard formalism.
4.1. Initial setup
We consider the same disc as presented in Sellwood (2012). It is an infinitely thin Mestel disc for which the circular speed v_{φ} is a constant V_{0} independent of the radius. The stationary background potential ψ_{M} and its associated surface density Σ_{M} are given by (44)where R_{max} is a scale parameter of the disc. Following Toomre (1977), Binney & Tremaine (2008), a selfconsistent DF for this system is given by (45)where the exponent q is given by (46)with σ_{r} being the constant radial velocities spread within the disc. In Eq. (45), C_{M} is a normalisation constant given by (47)In order to deal with the central singularity of the Mestel disc along with its infinite extent, we introduce two tapering functions (48)where the indices ν_{t} and μ_{t} control the sharpness of the two tapers, while the radii R_{i} and R_{0} are two scale parameters. These tapers T_{inner} and T_{outer} respectively represent the bulge and the outer truncation of the disc. In addition to these taperings, we also suppose that only a fraction ξ of the stellar disc is selfgravitating, with 0 ≤ ξ ≤ 1, while the rest of the gravitational potential is provided by the static halo. As a consequence, the active distribution function F_{star} is given by (49)We place ourselves in the same units system as in S12, so that we have V_{0} = G = R_{i} = 1. The other numerical factors are given by q = 11.4, ν_{t} = 4, μ_{t} = 5, ξ = 0.5, R_{0} = 11.5, and R_{max} = 20. The contours of the tapered DF F_{star} are illustrated in Fig. 1. At this stage, it is important to note that S12 restricted the perturbations forces to the harmonic sector m_{φ} = 2, so that we may use the same restriction on the considered azimuthal number m_{φ}. As a consequence, in the double resonance sum on m_{1} and m_{2} present in the BalescuLenard flux from Eq. (2), we assume that m_{1} and m_{2} belong to the restricted set {m_{r},m_{φ} } ∈ {(−1,2),(0,2),(1,2) }, where (m_{r},m_{φ}) = (−1,2) corresponds to the Inner Lindblad resonance (ILR), (m_{r},m_{φ}) = (0,2) to the Corotation resonance (COR), and (m_{r},m_{φ}) = (1,2) to the Outer Lindblad Resonance (OLR). All the upcoming calculations have also been performed while taking into account the contributions from the resonances associated with m_{r} = ± 2, which were checked to be largely subdominant.
Fig. 1 Contours of the initial distribution function F_{star} from Eq. (49), in action space (J_{φ},J_{r}). The contours are spaced linearly between 95% and 5% of the distribution function maximum. 

Open with DEXTER 
4.2. Initial drift and diffusion
As detailed in Eq. (33), the computation of the response matrix requires to consider a grid in the (r_{p},r_{a})space. We considered a grid such that , and Δr = 0.05. The sum on m_{1} appearing in Eq. (33) was reduced to . The basis considered was Kalnajs 2D basis (Kalnajs 1976) with the parameters k_{Ka} = 7 and R_{Ka} = 5. One should note that despite having a disc which extends up to R_{max} = 20, one can still consider a basis truncated at such a small R_{Ka}, so as to be able to efficiently capture the diffusion properties of the system in its inner regions, from where the secular diffusion is known to start. The radial basis elements were restricted to 0 ≤ n ≤ 8. When evaluating the response matrix, as in Eq. (32), one has to add a small imaginary part η to the frequency so as to regularise the resonant denominator. Throughout the calculations presented below, we considered η = 10^{4} and checked that this choice had no impact on our results.
Fig. 2 Illustration of 4 different critical resonant lines in the (r_{p},r_{a})space. As defined in Eq. (41), a critical line is characterised by the resonant vectors m_{1}, m_{2} and a location in action space. Each of the 4 plotted critical lines are associated with the same location , represented by the black dot. The critical lines correspond to various choices of resonant vectors m_{1} and m_{2} among the three inner, outer and corotation Lindblad resonances, respectively ILR, OLR and COR. One should also note that for m_{1} = m_{2}, the critical lines go through the point . 

Open with DEXTER 
Since the total potential ψ_{M} is known via Eq. (44), the mapping to the angleaction coordinates is completely determined. The two intrinsic frequencies of the system can then be computed on the (r_{p},r_{a})grid via Eqs. (13) and (14). Once these frequencies are known, the critical resonant lines introduced in Eq. (41) can be determined and are illustrated in Fig. 2. It is along these lines that one has to perform the integration present in the definitions of the drift and diffusion coefficients from Eq. (42).
Thanks to this expression, one can then compute the secular diffusion flux ℱ_{tot} defined in Eq. (11). Because the mass of the particles is given by μ = M_{tot}/N, it is natural to consider the quantity Nℱ_{tot} which is independent of N. The vector field −Nℱ_{tot} = −(Nℱ_{Jφ},Nℱ_{Jr}), which represents the direction of diffusion of individual particles, is illustrated in Fig. 3. One can already note in Fig. 3 that the diffusion vector field is along a narrow resonant direction. Along this ridge, one typically has ℱ_{Jφ} ≃ −2ℱ_{Jr}, so the diffusion appears as aligned with the direction of the ILR resonance given by m_{ILR} = (2,−1).
Fig. 3 Map of −Nℱ_{tot}, where the flux has been computed with m_{1},m_{2} ∈ {m_{ILR},m_{COR},m_{OLR}}. Following Eq. (11), −Nℱ_{tot} corresponds to the direction of diffusion of individual particles in actionspace. 

Open with DEXTER 
Fig. 4 Left panel: map of Ndiv(ℱ_{tot}), where the total flux has been computed with m_{1},m_{2} ∈ {m_{ILR},m_{COR},m_{OLR}}. Red contours, for which Ndiv(ℱ_{tot}) < 0 are associated with regions from which the orbits will be depleted, whereas blue contours, for which Ndiv(ℱ_{tot}) > 0 correspond to regions where the value of the DF will be increased during the secular diffusion. The contours are spaced linearly between the minimum and the maximum of Ndiv(ℱ_{tot}). The maximum value for the positive blue contours corresponds to Ndiv(ℱ_{tot}) ≃ 350, while the mininum value for the negative red contours is associated with Ndiv(ℱ_{tot}) ≃ −250. Right panel: from Sellwood (2012) – Fig. 7, contours of the change in the DF between the time t_{S12} = 1400 and t_{S12} = 0, for a run with 50M particles. Similarly to the left panel, red contours correspond to negative differences, i.e. regions emptied from their orbits, while blue contours correspond to positive differences, i.e. regions where the DF has increased during the diffusion. Both of these contours are aligned with the ILR direction of m_{ILR} = (2,−1) in the (J_{φ},J_{r})plane, corresponding to the cyan line. 

Open with DEXTER 
Once the diffusion flux Nℱ_{tot} has been obtained, one can compute the divergence of this flux, to determine the regions for which the DF is expected to change during the secular diffusion. Figure 4 illustrates the contours of Ndiv(ℱ_{tot}). In Fig. 4, we obtained that the BalescuLenard formalism indubitably predicts the formation of a narrow resonant ridge aligned with the ILRdirection, as was observed in S12’s simulation. One also recovers that the stars which will populate the resonant ridge originate from the base of the ridge and diffuse along the direction associated with the ILR resonance. It is most likely that the slight shift in the position of the ridge is due to the fact that the BalescuLenard prediction was carried at t = 0^{+}, while S12’s measurements are at t = 1400, so that we do not expect a perfect match. Other sources of discrepancies might be the use of a softening length in numerical simulations, which modifies the twobody interaction potential, or the difference between an ensemble average (as predicted by the secular formalism) and one specific realisation – our own simulations suggest that there is some variation in the position of the ridge between one run and another. Because we explicitly determined the value of Ndiv(ℱ_{tot}), we may now study the typical timescale of collisional relaxation predicted by this BalescuLenard estimation as detailed in Sect. 4.3. One may also investigate the respective roles of the selfgravitating amplification and the limitation to the tightly wound basis elements as presented in Appendices D.1 and D.2.
4.3. Timescale of diffusion
The most significant disagreement found in Paper I, while applying the WKB approximation of the BalescuLenard equation to S12’s simulation was a discrepancy between the time required to observe the resonant ridge in S12’s simulation and the collisional timescale for which the finiteN effects come into play. As already noted in Paper I, since the BalescuLenard Eq. (2) only depends on N through the mass of the individual particles μ = M_{tot}/N, we may rewrite it in the form (50)where C_{BL} [ F ] = Ndiv(ℱ_{tot}) is the Nindependent BalescuLenard collisional operator, i.e. the r.h.s. of Eq. (2) multiplied by N = M_{tot}/μ. As expected, the larger the number of particles, the slower the secular evolution. This also illustrates the fact that the BalescuLenard equation comes from a kinetic Taylor expansion in the small parameter ε = 1 /N ≪ 1. We introduce the rescaled time τ = t/N, so that Eq. (50)reads (51)letting us express the BalescuLenard equation without any explicit appearance of N. In Paper I, we estimated the time Δτ_{S12} required to observe the ridge as Δτ_{S12} ≃ 3 × 10^{5}. When performing the same measurement thanks to the contours of the diffusion flux div(ℱ_{tot}) computed within the WKB approximation, we obtained Δτ_{WKB} ≃ 3 × 10^{2}, so that Paper I obtained the ratio Δτ_{S12}/ Δτ_{WKB} ≃ 10^{3}. This discrepancy was due to the limitation to tightly wound spirals. Because the estimation of the secular diffusion flux ℱ_{tot} presented in Fig. 4 was made using the matrix method (Kalnajs 1976) with a full basis, it captures the additional swing amplification. Indeed, given the map of Ndivℱ_{tot} obtained in Fig. 4, one may estimate the typical time Δτ_{BL} required for such a flux to lead to the diffusion features observed in S12. The contours presented in Fig. 7 of S12 are separated by an increment equal to , where is the maximum of the normalised DF (via Eq. (49)). In order to observe the resonant ridge, the value of the DF should typically change by an amount of the order . From Fig. 4, one can note that the maximum of the divergence of the diffusion flux is given by  Ndiv(ℱ_{tot})  _{max} ≃ 350. Thanks to Eq. (11), one can immediately write the relation ΔF_{0} ≃ Δτ_{BL}  Ndiv(ℱ_{tot})  _{max}, where Δτ_{BL} is the time during which the BalescuLenard equation has to be evolved in order to develop a ridge. With the previous numerical values, one obtains Δτ_{BL} ≃ 3 × 10^{5}. Comparing the numerically measured time Δτ_{S12} and the time Δτ_{BL} predicted from the BalescuLenard equation, one obtains (52)As expected, the projection of the response over an unbiased basis leads to over a hundredfold increase of the susceptibility of the disc and therefore to a very significant acceleration of secular diffusion. Thanks to this mechanism, we now find a very good agreement between the diffusion timescales observed in numerical simulations and the predictions from the BalescuLenard formalism. This quantitative match is rewarding, both from the point of view of the accuracy of the integrator (symplecticity, timestep size, softening...), and from the relevance of the successive approximations underpinning the BalescuLenard formalism (timescale decoupling, truncation of the BBGKY hierarchy, neglect of the close encounter term ...).
In Appendix D, we show that when considering either Fig. D.1, for which the selfgravity of the system has been turned off, or Fig. D.3 for which the loosely wound basis elements were not taken into account, one does not recover a narrow resonant ridge appearing on timescales compatible with S12’s simulations. Therefore, the main source of secular collisional diffusion oberved in S12 and recovered in Fig. 4 has to be the strong selfgravitating amplification of loosely wound perturbations, i.e. a sequence of uncorrelated swing amplified spirals sourced by finiteN effects is indeed the main driver of secular diffusion. The WKB formalism from Paper I identified correctly the family of orbits involved, but fell short in predicting how narrow the resonant ridge is and how strongly amplified the response is.
5. Comparison to Nbody simulations
In Paper I, we relied on simulations presented in Sellwood (2012) to compare the divergence of the diffusion flux to the WKB prediction. In order to probe the expected scalings with the number of particles or with the active fraction of the disc, we now resort to our own Nbody simulations.
5.1. Nbody integration
The initial sampling of particles is critical when investigating the origin of secular evolution, as one must ensure that the disc is initially in a state of equilibrium. The sampling strategy we implemented is described in some detail in Appendix E.
Once sampled, we evolve the initial conditions using a straightforward particlemesh Nbody code with a singletimestep leapfrog integrator (e.g. Binney & Tremaine 2008, Sect. 3.4.1). We follow S12 and split the potential in which the particles move into two parts: (i) an axisymmetric contribution ψ_{M} from the unperturbed Mestel disc, as in Eq. (44), and (ii) a nonaxisymmetric contribution ψ_{1}(R,φ) that develops as perturbations grow in the disc. This splitting avoids difficulties in the treatment of the rigid component of the potential that is not included in the DF, due to the tapering functions and active fraction introduced in Eq. (49). We calculate ψ_{1} using cloudincell interpolation (e.g. Binney & Tremaine 2008, Sect. 2.9.3) of the particles’ masses onto an N_{mesh} × N_{mesh} mesh of square cells spaced Δx apart, then filtering the resulting density field ρ(x,y) to isolate the disc response (see below), before applying the usual Fourierspace “doubling up” procedure to obtain the potential ψ_{1} at the cell vertices. The contribution of ψ_{1} to each particle’s acceleration is then obtained using the same cloudincell interpolation scheme.
When computing the density mesh, we added a filtering scheme to include only the m_{φ} = 2 disc response, similarly to what was considered in S12. We isolate this m_{φ} = 2 mode by calculating (53)immediately after the cloudincell assignment of mass to the mesh at each timestep, then imposing the new mesh mass distribution (54)with (r_{k},φ_{k}) chosen according to (x_{k},y_{k}) = (r_{k}cos(φ_{k}), r_{k}sin(φ_{k})). To obtain ρ_{2}(r) we use bruteforce computation of Eq. (53) on a series of N_{ring} radial rings with spacing Δr ≪ Δx, using the trapezium rule with N_{φ} = 720 points in φ for the angular integrals. These models are designed to reproduce as closely as possible the essential details of S12’s simulations. There are a couple of deliberate technical differences: S12 uses a polar mesh to obtain ψ_{1}, whereas we use a cartesian mesh with a m_{φ} = 2 prefiltering of the density field; S12 has a block timestep scheme instead of our simpler singletimestep one.
For the results presented here we used a timestep Δt = 10^{3}R_{i}/V_{0} on a mesh that extends to ± R_{max} = 20 R_{i} with N_{mesh} = 120 cells, so that Δx = R_{i}/ 3. The filtering of the potential perturbations to the harmonic sector m_{φ} = 2 was performed with N_{ring} = 1000 radial rings, so that Δr = 2 R_{i}/ 100, and N_{φ} = 720 points in the azimuthal direction. Finally, the computation of the potential from the density via Fourier transform was performed with a softening length ε = R_{i}/ 6, comparable to the value used in Sellwood (2012), which considered a Plummer softening with ε = R_{i}/ 8. The results are not significantly changed when we halve the timestep or the mesh size. In Appendix C, we detail the validation of our Nbody code by recovering known unstable modes of truncated Mestel discs (Zang 1976; Evans & Read 1998b; Sellwood & Evans 2001).
5.2. Scaling with N
In order to rid our measurements of individual fluctuations, we run multiple simulations for the same number of particles and perform an ensemble average of different evolution realisations for the same number of particles. This procedure allows us to estimate only the mean evolution, which is effectively what is described by the BalescuLenard formalism.
In order to study the scaling with N of these numerical simulations, one has to extract from the simulations a quantity on which to test this scaling and compare it with the predictions from the BalescuLenard formalism. The statistical nature of the initial sampling presents an additional difficulty. Indeed, because one only samples N stars as described in Appendix E, the initial effective DF fluctuates around the smooth background DF from Eq. (49) as a Poisson shot noise. These statistical fluctuations originate from the initial sampling and are not as such specific to the physical process captured by the BalescuLenard formalism, so that one should carefully disentangle these two contributions. Hence we introduce the function defined as (55)where the operator ⟨·⟩ corresponds to the ensemble average, approximated here with the arithmetic average over the p = 32 different realisations of simulations for the same number of particles, indexed by i: . In Eq. (55) the function h_{i}(t,N) is a lag function which reads (56)where we defined as F_{i}(t,J,N) the normalised DF of the ith realisation for a number N of particles. Such a quantity is designed to quantify the “distance” between the initial mean DF ⟨F(t = 0)⟩ and the evolved DF F_{i}. We are interested in the early time behaviour of the lag function h from Eq. (55), so that we may perform its Taylor expansion (57)where it is important to note that the coefficients , and depend only on N. Let us now estimate each of these coefficients in turn. Thanks to Eq. (56), one can compute which reads (58)where we used the shortened notations ⟨F_{0}⟩ = ⟨F(t = 0,J,N)⟩ and F = F(t = 0,J,N). We note that this coefficient only depends on the properties of the initial sampling, and not on its dynamics. Because discrete sampling obeys Poisson statistics, one can write (59)where α_{0} is a constant independent of N. One may then compute , which takes the form (60)where we used the shortened notation F′ = [ ∂F/∂t ] (t = 0). One should note that the terms appearing in Eq. (60) have two different physical contents. Indeed, the term [F − ⟨F_{0}⟩] involves the initial sampling, whereas F′ is driven by the dynamics of the system. If we assume that the sampling and the system’s dynamics are uncorrelated, one writes (61)As a consequence, one immediately obtains from Eq. (60) that . One can finally compute the coefficient which reads (62)where we used the shortened notation F′′ = [ ∂^{2}F/∂t^{2} ] (t = 0). Using the same argument as in Eq. (61), we may get rid of the second term in the r.h.s. of Eq. (62). If we now also assume that the variance of [ F′ ] ^{2} is small compared to its expectation, one can write ⟨[ F′ ] ^{2}⟩ = [⟨F′⟩]^{2}, so that Eq. (62) becomes (63)Now the dependence of the term ⟨F′⟩ with N follows from Eq. (50), so that one can write (64)where α_{2} is an amplitude independent of N. This scaling is a prediction from the BalescuLenard formalism. If the secular evolution observed in S12’s simulation was a Vlasovonly evolution, i.e. a collisionless evolution, one would expect a scaling of such that .
Fig. 5 Illustration of the behaviour of the function defined in Eq. (55), for an active fraction ξ = 0.5, when averaged on 32 different realisations for particles numbers N ∈ {8, 12, 16, 24, 32, 48, 64 } × 10^{5}. To compute , we binned the actionspace domain (J_{φ},J_{r}) = [ 0;2.5 ] × [ 0;0.2 ] in 100 × 50 regions. The values of have also been uniformly renormalised so as to clarify this representation. The dots corresponds to the snapshots of the simulations for which was computed, whereas the lines correspond to secondorder fits. As expected, the smaller the number of particles, the noisier the simulation and the larger . 

Open with DEXTER 
Fig. 6 Top panel: illustration of the behaviour of the function , where N has been rescaled by a factor 10^{5} so as to simplify the representation. The dots correspond to computed values thanks to Fig. 5, while the line corresponds to a linear fit, which takes the form . The coefficients have been uniformly renormalised so as to clarify this representation. Bottom panel: similar representation for the behaviour of the function , whose linear fit takes the form . Similarly, the coefficients have been uniformly renormalised so as to clarify this representation. 

Open with DEXTER 
One may now compare these predictions to the scalings obtained from Nbody runs. We considered number of particles given by N ∈ {8, 12, 16, 24, 32, 48, 64 } × 10^{5}, and for each of these values of N, we ran 32 different simulations with different initial conditions while using the Nbody code described in Sect. 5.1. For each value of N, one may study the function , as illustrated in Fig. 5. Once the behaviour of the function is known, one can fit these to parabolas as in Eq. (57), so as to determine the behaviour of the functions and . The dependence with N of these coefficients is illustrated in Fig. 6. From the top panel of Fig. 6, we recover the scaling of derived in Eq. (59) due to Poisson shot noise. The bottom panel of Fig. 6 displays the scaling . Given the finite number of simulations considered and the uncertainties in the fits, this is in good agreement with the result presented in Eq. (64). This scaling of with N therefore confirms the relevance of the BalescuLenard formalism to describe the secular evolution of S12’s stable Mestel disc. Specifically, as explained below Eq. (64), if the features observed in S12’s simulation had only been the result of a collisionless mechanism, one would not have observed such a scaling of with N. This scaling confirms that the secular evolution of S12’s stable Mestel disc is the result of a collisional evolution seeded by the discrete nature of the system, and the effect of amplified distant resonant encounters. Another probe of the collisional scaling, which allows us to get rid of Poisson shot noise as present in Eq. (59), is described in Appendix F.
5.3. Scaling with ξ
Since the novelty of the BalescuLenard formalism is to capture the effect of gravitational polarisation, we now further compare qualitatively the prediction from Sect. 4 with the results obtained from numerical simulations, by studying the impact of the active fraction ξ of the disc on the observed properties of the secular diffusion. Indeed, as detailed in Sect. 4.1, the disc considered in S12 had an active fraction of ξ = 0.5, so only one half of the potential was due to the active component. If one increases the active fraction of the disc, one increases the strength of the selfgravitating amplification, and therefore accelerates the secular evolution of the disc, while still remaining in a regime of collisional evolution. Therefore the scaling of with N given by Eq. (64) remains the same, but the prefactor α_{2}(ξ) increases because the secular evolution is amplified via a more efficient polarisation. The dependence of α_{2} with ξ can be both measured from Nbody simulations but also predicted using the BalescuLenard formalism via the calculations presented in Sect. 4.2. Let us consider the same sets of simulations as in Sect. 5.2, so that the number of particles are given by N ∈ {8, 12, 16, 24, 32, 48, 64} × 10^{5}, and for each of these values of N, 32 different simulations with ξ = 0.6 are performed, in order to carry out ensemble averages.
Fig. 7 Top panel: illustration of the behaviour of the function for an active fraction ξ = 0.6, using the same conventions as in Fig. 5. As expected, for an increased active fraction, the secular evolution of the system is fastened leading to a faster increase of the distance . Middle panel: behaviour of the function , for an active fraction ξ = 0.6, using the same conventions as in Fig. 6. Its linear fit takes the form . One recovers the expected scaling of the Poisson shot noise sampling derived in Eq. (59). Bottom panel: behaviour of the function , for an active fraction ξ = 0.6, using the same conventions as in Fig. 6. Its linear fit takes the form . One recovers the expected collisional scaling with N obtained in Eq. (64). 

Open with DEXTER 
The equivalent of Figs. 5 and 6 for ξ = 0.6 is illustrated in Fig. 7. Even for ξ = 0.6, one finds that the function follows a parabola given by Eq. (57). One also recovers the predicted scalings with N of the functions and , respectively representing the initial Poisson shot noise of the sampling and the collisional scaling of the BalescuLenard secular evolution. As expected, when the active fraction of the disc is increased the secular evolution is accelerated. Thanks to these fits, one can study the dependence of the ratio α_{2}(ξ = 0.6) /α_{2}(ξ = 0.5), as defined in Eq. (64), both from numerical simulations as described in Figs. 6 and 7 and from the BalescuLenard equation using the matrix method described in Sect. 4.
From the fits of from Figs. 6 and 7, one can write for ξ = 0.5 and for ξ = 0.6, where we shifted the intercept of the fits to correspond to the centre of the considered region log (N) ∈ [ log (8) ;log (64) ]. One therefore obtains the ratio (65)One may now compare this Nbody measurement to the same measurement performed via the BalescuLenard formalism. Following Eq. (63), one obtains that this ratio is given by (66)where stands for the secular diffusion flux at t = 0^{+} with an active fraction ξ. The value of α_{2}(ξ = 0.5) can be determined via Fig. 4, while the secular diffusion flux determined for ξ = 0.6 is illustrated in Fig. 8. Thanks to the contours presented in Figs. 4 and 8, one can perform the same estimation as in Eq. (65) starting from the BalescuLenard predictions. In order to focus on the contributions associated with the resonant ridge, the integrals on J in Eq. (66) were performed for J_{φ} ∈ [ 0.5; 1.2 ] and J_{r} ∈ [ 0.06; 0.15 ]. We measured (67)Despite the difficulty of this measurement which required to consider a much more sensitive disc with ξ = 0.6, the ratios of α_{2} measured either via direct Nbody simulations as in Eq. (65) or via application of the BalescuLenard formalism in Eq. (67) are within the same order of magnitude. As a consequence, one indeed checks that the BalescuLenard equation is able to correctly capture the relative effect of the disc’s susceptibility on the characteristics of the collisional secular diffusion. The strong consequence of modifying the active fraction, observed both in Eqs. (65) and (67), illustrates the relevance of the selfgravitating amplification in determining the typical timescale of secular diffusion of the system.
Fig. 8 Map of Ndiv(ℱ_{tot}) with an active fraction ξ = 0.6 and using the same conventions as in Fig. 4. The contours are spaced linearly between the minimum and the maximum of Ndiv(ℱ_{tot}). The maximum value for the positive blue contours corresponds to Ndiv(ℱ_{tot}) ≃ 4200, while the minimum value for the negative red contours corresponds to Ndiv(ℱ_{tot}) ≃ −3200. As expected, one recovers that when the active fraction of the disc is increased, the susceptibility of the disc is increased, so that the norm of Ndiv(ℱ_{tot}) gets larger and the secular diffusion is fastened. 

Open with DEXTER 
5.4. Latetime evolution
The predictions of the BalescuLenard secular diffusion flux presented in Sect. 4 were only applied for the initial time of evolution, i.e. for the estimation of ℱ_{tot}(t = 0^{+}). The Nbody simulations presented in Sect. 5 allowed us to verify the appropriate scaling of the response of the system with the number of particles for the initial time of evolution as illustrated in Fig. 6. Using the BalescuLenard formalism to probe the late secular evolution of the system would require to evolve iteratively Eq. (2) over secular times. Such iterations are clearly beyond the scope of this first paper, however the use of Nbody simulations allows us to start probing now such late times of evolution.
As discussed in Sect. 2, the BalescuLenard equation describes the longterm evolution of a discrete selfgravitating inhomogeneous system. Such a collisional evolution is only relevant for dynamically stable systems, i.e. systems assumed to be stable in the Vlasov sense, and remains valid as the system moves through a series of equlibria, while the drift and diffusion coefficients change with F. Because it has been obtained via a Taylor expansion of the dynamics at the order 1 /N in the number of particles, it remains valid only on secular timescales of the order Nt_{D}, with t_{D} the dynamical time.
On such secular timescales, a BalescuLenard evolution can lead to two different outcomes. On the one hand, if the system remains stable during its entire evolution, the BalescuLenard equation will tend towards a 1 /Nstationary state^{2}. Once such a stationary state of evolution has been reached, the dynamics is then governed by the next order kinetic effects in 1 /N^{2}, which are not captured by the BalescuLenard equation. On the other hand, the BalescuLenard collisional evolution may also lead to a destabilisation of the system. Indeed, the longterm effects of the collisional diffusion, because they lead to an irreversible diffusion of the DF, may change its current state w.r.t. the collisionless (Vlasov) dynamics. After a slow and stable evolution sourced by collisional 1 /N effects, the system may then become dynamically unstable with respect to collisionless dynamics, which becomes the main driver of its latertime evolution, as was suggested by Sellwood (2012). S12 observes an outofequilibrium transition between the 1 /N BalescuLenard collisional evolution and the collisionless Vlasov evolution.
One can illustrate such a transition using the Nbody simulations presented previously. In order to capture the change of regime of evolution within the disc (collisional vs. collisionless), for a given value of the number of particles, we define the quantity Σ_{2}(t,N) as (68)where as in Eq. (55), the operator ⟨·⟩ corresponds to the ensemble average, approximated here with the arithmetic average over the p = 32 different realisations of simulations for the same number of particles N. The radii considered are restricted to the range R ∈ [ R_{inf}; R_{sup} ] = [ 1.2; 5 ], where the active surface density of the disc is little affected by the inner and outer tapers. Finally, to obtain the second equality in Eq. (68), as in Eq. (C.2), we replaced the active surface density of the disc by a discrete sum over all the particles of the system, where the sum on n is restricted to all the particle whose radius lies between R_{inf} and R_{sup}, while their azimuthal phase was written as φ_{n}. Such a quantity allows us to probe easily the presence of strong nonaxisymmetric features within the disc. During the initial BalescuLenard collisional evolution of the system, one expects low values of Σ_{2}. Indeed, during this evolution, one relies on the phase averaging approximation, which assumes that F = F(J,t), so that the DF of the system does not depend on the angles θ. During this collisional phase, Σ_{2} still remains nonzero because the system develops transient spiral waves, which sustain the secular evolution. On the longterm, this collisional evolution, through an irreversible diffusion of the DF, leads to a destabilisation of the system. Eventually, the dynamical drivers of evolution are no longer discrete resonant collisional effects but exponentially growing dynamical instabilities. In this regime of collisionless unstable evolution, one expects much larger values of Σ_{2}, because of the appearance of strong nonaxisymmetric bars within the disc. This bifurcation between these two regimes of diffusion is illustrated in Fig. 9, through the behaviour of the function ^{3}. One can similarly observe this transition directly by looking at the active surface density Σ_{t}(R,φ,t) for these two different regimes. This is illustrated in Fig. 10, where one recovers that in the late time collisionless regime of evolution, the galaxy becomes strongly nonaxisymmetric. S12 found that just after the disc becomes unstable, the pattern of the spiral response is consistent with the ILR frequency corresponding to the ridge^{4}. Hence the phase transition observed in Fig. 10 is driven by all the free energy available in a cold disc, which via spiral transients secularly heats the disc, but only along a very tight resonant direction. This in turn leads the disc towards an orbital instability, transverse to the resonance (via the direct azimuthal analogue to the two stream instability in plasma physics, LyndenBell 1979; Pichon 1994). Qualitatively, one expects that the more massive and the narrower the ridge, the larger the number of orbits trapped in ILR resonance with little relative azimuthal dispersion, and the earlier the instability (Penrose 1960; Pichon & LyndenBell 1993).
Fig. 9 Behaviour of the function as defined in Eq. (68), for various values of the number of particles. The prefactor has been added so as to mask Poisson shot noise, allowing the initial values of to be independent of N. It illustrates the bifurcation between the initial BalescuLenard collisional evolution, for which low values of Σ_{2} are expected and the collisionless Vlasov evolution for which the system is no more axisymmetric leading to larger values of Σ_{2}. As expected, the larger the number of particles, the later the transition. 

Open with DEXTER 
Fig. 10 Illustration of the active surface density Σ_{t} for a Nbody run with N = 8 × 10^{5}, restriced to the range R ≤ 6. Top panel: active surface density Σ_{t} at an early time t = 60, for which the galaxy remains globally axisymmetric. In this regime, the dynamics of the system is collisional and governed by the BalescuLenard Eq. (2). Bottom panel: active surface density Σ_{t} at a much later time t = 2400. The galaxy is then strongly nonaxisymmetric. In this regime, the dynamics of the system is collisionless and governed by the Vlasov equation. 

Open with DEXTER 
In closing, it is quite striking that an isolated galactic disc, fully stable in the mean field sense, will, given time, drive itself through twopoint resonant correlations towards instability, demonstrating the extent to which such cold systems are truly secularly metastable.
6. Conclusion
Most astrophysical discs formed through dissipative processes and have typically evolved over many dynamical times. Even in isolation, the longrange force of gravity allows their components to interact effectively through resonances, which given time may drive them secularly towards more likely equilibria. Such processes are captured by recent extensions of kinetic theories rewritten in angleaction variables (Heyvaerts 2010; Chavanis 2012). Solving these equations provide astronomers with a unique opportunity to quantify the induced secular angular momentum redistribution within these discs (LyndenBell & Kalnajs 1972) on cosmic timescales. While challenging, the numerical computation of the corresponding diffusion and drift coefficients is, as demonstrated, within reach of a relatively straightforward extension of the socalled matrix method (Kalnajs 1976), which computes the orbital response of selfgravitating discs using quadratures and linear algebra.
Paper I presented asymptotic expressions in the tightly wound limit and provided a qualitative insight into the physical processes at work during the secular diffusion of a selfgravitating discrete disc. Conversely, in this paper, we computed numerically the drift and diffusion coefficients of the inhomogeneous BalescuLenard equation for such infinitely thin stellar discs. The selfgravity of the disc was taken into account via the matrix method, validated on unstable Mestel discs. We computed the divergence of the flux density in action space, div(ℱ_{tot}). Swing amplification was shown to provide a significant boost for the diffusion timescale, which now matches the numerically measured one. These computations are the first exact calculation of the BalescuLenard diffusion and drift coefficients in the context of inhomogeneous multiperiodic systems. They capture the essence of selfinduced evolution (nature), which should compete with environmentally induced evolution (nurture). We then compared these predictions to idealised numerical simulations of stable stationary and truncated Mestel discs sampled by pointwise particles, which were evolved for hundreds of dynamical times. Using ensemble averages of our Nbody runs, we also identified a clear signature of the BalescuLenard process in the scaling of the diffusion features with N and ξ, the fraction of the mass within the disc. As originally identified by Goldreich & LyndenBell (1965), Julian & Toomre (1966) in the context of their linear response, the susceptibility of cold selfgravitating discs plays a critical role for their secular evolution as it is squared in the BalescuLenard equation, which boosts considerably the effect of discreteness. Indeed, both the numerical experiments and our computation of the fluxes show that N_{eff} ~ N/10^{4}, which is consistent with the predicted rescaling in (it was shown forty years ago that for a Mestel disc , depending on the exact temperature of the disc, Toomre 1981).
Jointly with Paper I we now have a qualitative and quantitative understanding of the initial secular orbital diffusion process induced by the discreteness of galactic discs. Our qualitative understanding allows us to identify the role played by the square susceptibility in boosting the diffusion. Our quantitative agreement in both amplitude, position, width and scaling of the induced orbital signatures strongly suggests that secular evolution is indeed driven by resonances as captured by the BalescuLenard formalism, and that it does not depend on the initial phases of the disc (since the matching BalescuLenard fluxes are phase averaged). It demonstrates that this equation initially reproduces the observed evolution of selfgravitating discs driven by resonant twopoint correlations beyond the mean field approximation.
The next step will be to evolve iteratively Eq. (2) over a Hubble time, and compare with the result of Nbody simulations. One should also model it jointly with an externally induced orbital diffusion (Fouvry et al. 2015b) arising from for example a (possibly anisotropic) cosmic environment (Codis et al. 2012, 2015) so as to assess which process dominates. At the technical level, the BalescuLenard formalism should be used to (in)validate Nbody integrators accuracy on secular timescales. There are indeed very few analytical predictions on which to calibrate Nbody experiments in this regime. Such an exploration would also allow us to get a better grasp of the impact of the numerical parameters used in the Nbody integration (such as timestep, mesh size or softening length) on the longterm dynamics of the system.
Beyond the application described in this paper, the BalescuLenard formalism may in the future also be numerically implemented to describe for instance the secular diffusion of giant molecular clouds in galactic discs (which in turn could play a role in migrationdriven metallicity gradients and disc thickening), the secular migration of planetesimals in partially selfgravitating protoplanetary debris discs, or even the longterm evolution of population of stars and gas blobs near the Galactic centre. In 3D, assuming spherical symmetry, its implementation could be useful to describe spherical systems dominated by radial orbits, or the secular evolution of tidal debris in our possibly flattened galactic halo using Stäckel potentials.
A similar dynamical phase transition has been observed (Campa et al. 2008) in a toy model of systems with longrange interactions called the Hamiltonian Mean Field (HMF) model. During the slow collisional evolution, because of finiteN effects, the distribution function of the system changes with time. In certain cases, the system may become dynamically (Vlasov) unstable and undergo a rapid phase transition from a homogeneous phase to an inhomogeneous phase. This phase transition can be monitored by the magnetisation (see Fig. 1 of Campa et al. 2008) which is an order parameter playing a role similar to Σ_{2}(t,N).
One could also check that the disc’s distribution function corresponds at that stage to an unstable configuration, using the matrix method described in Appendix C.
See also Earn & Sellwood (1995) for a similar rewriting of the basis normalisations.
Acknowledgments
We thank the referee, James Binney, for numerous valuable comments. J.B.F. thanks the Institute of Astronomy, Cambridge, for hospitality while this investigation was initiated. J.B.F. and C.P. also thank the theoretical physics subdepartment, Oxford, for hospitality and the CNRSOxford exchange program for funding. J.B.F., C.P. and P.H.C. also thank the CNRS Inphyniti program for funding. C.P. thanks Clare and Churchill college, Cambridge, the French embassy and the community of http://mathematica.stackexchange.com for their help. J.M. thanks the Institut d’Astrophysique de Paris for their hospitality and the Ville de Paris for a “Research in Paris” award. We thank Donald LyndenBell, John Papaloizou, Walter Dehnen, Rebekka Bieri, Laura Monk, Gordon Ogilvie, Dmitry Pogosyan and Simon Prunet for stimulating discussions, and Eric Pharabod for his help with Fig. 4. Special thanks to Stephane Rouberol for customising the Horizon cluster for our purposes. This work is partially supported by the Spin(e) grants ANR13BS050005 of the French Agence Nationale de la Recherche (http://cosmicorigin.org) and by the LABEX Institut Lagrange de Paris (under reference ANR10LABX63) which is funded by ANR11IDEX000402.
References
 Aumer, M., & Binney, J. J. 2009, MNRAS, 397, 1286 [NASA ADS] [CrossRef] [Google Scholar]
 Balescu, R. 1960, Phys. Fluids, 3, 52 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Binney, J., & Tremaine, S. 2008, Galactic Dynamics, 2nd edn., Princeton Series in Astrophysics (Princeton University Press) [Google Scholar]
 Born, M. 1960, The Mechanics of the Atom (F. Ungar Pub. Co.) [Google Scholar]
 Campa, A., Chavanis, P.H., Giansanti, A., & Morelli, G. 2008, Phys. Rev. E, 78, 040102 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1942, Principles of Stellar Dynamics (University of Chicago Press) [Google Scholar]
 Chavanis, P.H. 2012, Physica A, 391, 3680 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Chavanis, P.H. 2013, A&A, 556, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Codis, S., Pichon, C., Devriendt, J., et al. 2012, MNRAS, 427, 3320 [NASA ADS] [CrossRef] [Google Scholar]
 Codis, S., Pichon, C., & Pogosyan, D. 2015, MNRAS, 452, 3369 [NASA ADS] [CrossRef] [Google Scholar]
 Dehnen, W. 1998, AJ, 115, 2384 [NASA ADS] [CrossRef] [Google Scholar]
 Earn, D. J. D., & Sellwood, J. A. 1995, ApJ, 451, 533 [NASA ADS] [CrossRef] [Google Scholar]
 Evans, N. W., & Read, J. C. A. 1998a, MNRAS, 300, 83 [NASA ADS] [CrossRef] [Google Scholar]
 Evans, N. W., & Read, J. C. A. 1998b, MNRAS, 300, 106 [NASA ADS] [CrossRef] [Google Scholar]
 Fouvry, J.B., Pichon, C., & Chavanis, P.H. 2015a, A&A, 581, A139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fouvry, J.B., Pichon, C., & Prunet, S. 2015b, MNRAS, 449, 1967 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & LyndenBell, D. 1965, MNRAS, 130, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Goldstein, H. 1950, Classical mechanics (AddisonWesley) [Google Scholar]
 Heyvaerts, J. 2010, MNRAS, 407, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Hörmander, L. 2003, The analysis of linear partial differential operators. I, Classics in Mathematics (SpringerVerlag) [Google Scholar]
 Jeans, J. 1929, Astronomy and Cosmogony (Cambridge Univ. Press) [Google Scholar]
 Julian, W. H., & Toomre, A. 1966, ApJ, 146, 810 [NASA ADS] [CrossRef] [Google Scholar]
 Kalnajs, A. J. 1972, in IAU Colloq. 10: Gravitational NBody Problem, ed. M. Lecar, Astrophys. Space Sci. Lib., 31, 13 [Google Scholar]
 Kalnajs, A. J. 1976, ApJ, 205, 745 [NASA ADS] [CrossRef] [Google Scholar]
 Landau, L. 1936, Phys. Z. Sowj. Union, 10, 154 [Google Scholar]
 Lenard, A. 1960, Annals of Physics, 10, 390 [NASA ADS] [CrossRef] [Google Scholar]
 LyndenBell, D. 1979, MNRAS, 187, 101 [NASA ADS] [CrossRef] [Google Scholar]
 LyndenBell, D., & Kalnajs, A. J. 1972, MNRAS, 157, 1 [NASA ADS] [CrossRef] [Google Scholar]
 McMillan, P. J. 2013, MNRAS, 430, 3276 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Mestel, L. 1963, MNRAS, 126, 553 [NASA ADS] [CrossRef] [Google Scholar]
 Penrose, O. 1960, Phys. Fluids, 3, 258 [NASA ADS] [CrossRef] [Google Scholar]
 Pichon, C. 1994, Ph.D. Thesis, University of Cambridge [Google Scholar]
 Pichon, C., & Cannon, R. C. 1997, MNRAS, 291, 616 [NASA ADS] [CrossRef] [Google Scholar]
 Pichon, C., & LyndenBell, D. 1993, in Statistical Description of Transport in Plasma, Astro and Nuclear Physics, eds. J. Misquich, G. Pelletier, & P. Schuck, 261 [Google Scholar]
 Quinlan, G. D., & Tremaine, S. 1992, MNRAS, 259, 505 [NASA ADS] [CrossRef] [Google Scholar]
 Rosenbluth, M., MacDonald, W., & Judd, D. 1957, Phys. Rev., 107, 1 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Sellwood, J. A. 1983, J. Computat. Phys., 50, 337 [NASA ADS] [CrossRef] [Google Scholar]
 Sellwood, J. A. 2010, MNRAS, 409, 145 [NASA ADS] [CrossRef] [Google Scholar]
 Sellwood, J. A. 2012, ApJ, 751, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Sellwood, J. A., & Evans, N. W. 2001, ApJ, 546, 176 [NASA ADS] [CrossRef] [Google Scholar]
 Toomre, A. 1964, ApJ, 139, 1217 [NASA ADS] [CrossRef] [Google Scholar]
 Toomre, A. 1977, ARA&A, 15, 437 [NASA ADS] [CrossRef] [Google Scholar]
 Toomre, A. 1981, in Structure and Evolution of Normal Galaxies, eds. S. M. Fall, & D. LyndenBell, 111 [Google Scholar]
 Tremaine, S., & Weinberg, M. D. 1984, MNRAS, 209, 729 [NASA ADS] [CrossRef] [Google Scholar]
 Vlasov, A. 1938, Zh. Eksp. i Teor Fiz., 8, 291 [Google Scholar]
 Weinberg, M. D. 1993, ApJ, 410, 543 [NASA ADS] [CrossRef] [Google Scholar]
 Wielen, R. 1977, A&A, 60, 263 [NASA ADS] [Google Scholar]
 Zang, T. A. 1976, Ph.D. Thesis, Massachusetts Institute of Technology [Google Scholar]
Online material
Appendix A: Kalnajs basis
We now detail the properties of the basis introduced in Kalnajs (1976) to describe 2D discs^{5}. The basis depends on two parameters: an index k_{Ka} ∈N and a scale radius r_{Ka} ∈ R^{+}. In all the upcoming formulae of this section, in order to shorten the notations, we write r for the dimensionless quantity r/r_{Ka}. As introduced previously in Eq. (16), the basis elements depend on two indices: the azimuthal number ℓ and the radial index n. One should note that we have ℓ ≥ 0 and n ≥ 0. The radial component of the potential elements are then of the form (A.1)The radial component of the density elements is given by (A.2)In Eqs. (A.1) and (A.2), the coefficients and are defined by (A.3)and (A.4)Finally, in Eqs. (A.1) and (A.2), we have also introduced (A.5)and (A.6)In the two previous expressions, we introduced the rising Pochhammer symbol [ a ] _{i} defined as (A.7)
Appendix B: Calculation of ℵ
We now detail how the analytical function ℵ from Eq. (32) may be computed. In order to ease the implementation of its computation, we rewrite ℵ in an undimensionalised way as follows (B.1)where we assumed that a_{g},a_{h} ≠ 0 and used the change of variables x = x_{p}/ Δr and y = x_{a}/ Δr. We also defined the dimensionless function ℵ_{D} as (B.2)To compute this integral, we may now exhibit a function G(x,y) such that (B.3)One possible choice for G is given by (B.4)One should note in the previous expression the presence of a complex logarithm and a tan^{1}. However, because e,f,η ∈R and η ≠ 0, one can easily show that the arguments of both of these functions never cross the usual branchcut of these functions {Im(z) = 0; Re(z) ≤ 0 }. As a consequence, the expression (B.2) can immediately be computed as (B.5)
Appendix C: Response matrix and Nbody validations
The computation of the response matrix as described in Sect. 3 was validated by recovering the results of the pioneer work of Zang (1976), extended in Evans & Read (1998a,b), and recovered numerically in Sellwood & Evans (2001). These papers predicted the precession rate ω_{0} = m_{φ}Ω_{p} and growth rate η = s of the unstable modes of a truncated Mestel disc similar to the stable one described in Sect. 4.1. To build up an unstable disc similar to the ones considered in these previous works, one has to consider a fully active disc, so that ξ = 1. So as to have Q = 1 (Toomre 1964), the velocity dispersion within the disc is given by q = 6, where the parameter q has been introduced in Eq. (46). Finally, a last parameter one can tune in order to modify the properties of the disc is the truncation index of the inner tapering ν_{t} defined in Eq. (48). While looking only for m_{φ} = 2 modes, we considered three different truncations indices given by ν_{t} = 4, 6, 8. To compute the response matrix, we used the same numerical parameters as described in Sect. 4.2. Looking for unstable modes amounts to looking for complex frequencies ω = ω_{0} + iη such that the response matrix from Eq. (24) possesses an eigenvalue equal to 1. Such a complex frequency is then associated with an unstable mode of pattern speed ω_{0} and growth rate η. To determine the growth rate and pattern speed of the unstable modes, we relied on Nyquist contours similarly to the technique presented in Pichon & Cannon (1997). For a fixed value of η, one can study the continuous complex curve . Because for η → + ∞, one has , the number of windings of this curve around the origin gives a lower bound on the number of unstable modes with a growth rate superior to η. By decreasing the value of η, one can then determine the largest value of η admitting an unstable mode, and therefore the most unstable mode of the disc. The Nyquist contours obtained for the truncation index ν_{t} = 6 are illustrated in Fig. C.1, while the measurements are gathered in Table C.1.
Fig. C.1 Top panel: zoomed Nyquist contours in the complex plane of obtained via the matrix method for a truncated Mestel disc with ν_{t} = 6 and q = 6, looking for m_{φ} = 2 modes. Each contour corresponds to a fixed value of η. For a growth rate of η ≃ 0.20, one can note that the contour crosses the origin, which corresponds to the presence of an unstable mode. Bottom panel: illustration of the function for the same truncated Mestel disc. Each line corresponds to a fixed value of η. Such a representation allows us to determine the pattern speed ω_{0} = m_{φ}Ω_{p} ≃ 0.94 of the unstable mode. 

Open with DEXTER 
Once the characteristics (ω_{0},η) of the unstable modes have been determined, one can study in the physical space the shape of the mode. Indeed, for ω = ω_{0} + iη, one can compute , and numerically diagonalise this matrix. One then considers its eigenvector X_{mode} (of size n_{max}, where n_{max} is the number of basis elements considered) associated with the eigenvalue almost equal to 1. The shape of the mode is then immediately given by (C.1)where Σ^{(p)} are the considered surface density basis elements. The shape of the recovered unstable mode for the truncated ν_{t} = 4 Mestel disc is illustrated in Fig. C.2.
Fig. C.2 Unstable mode for the truncated ν_{t} = 4 Mestel disc recovered via the matrix method presented in Sect. 3. Only positive contour levels are shown and are spaced linearly between 10% and 90% of the maximum norm. The radii associated with the resonances ILR, COR and OLR have been represented, as given by ω_{0} = m·Ω(R_{m}), where the intrinsic frequencies Ω(R) = (Ω_{φ}(R),κ(R)) have to be computed within the epicyclic approximation. For a Mestel disc, they are given by Ω_{φ}(R) = V_{0}/R and . 

Open with DEXTER 
Fig. C.3 Measurement of the growth rate η and pattern speed ω_{0} of the m_{φ} = 2 unstable modes of truncated Mestel discs with a random velocity given by q = 6 for various values of the truncation index ν_{t} = 6, 8. The basis coefficient plotted is associated with the basis element (ℓ,n) = (2,0), using the same basis elements as for the matrix method in Sect. 4.2. 

Open with DEXTER 
The same unstable modes were also used to validate the Nbody code presented in Sect. 5. To run these simulations, we used the same sampling technique as described in Appendix E. In order not to be significantly impacted by the absence of a quiet start sampling (Sellwood 1983), for each value of ν_{t}, the measurements were performed with simulations of 20M particles. As observed in Sellwood & Evans (2001), the appropriate setting of the parameters of the Nbody code are crucial to recover correctly the unstable modes of a disc. We considered a grid made with N_{mesh} = 120 grid cells, while using a softening length equal to ε = R_{i}/ 60. As described in Sect. 5, we similarly restricted the perturbing forces only to the harmonic sector m_{φ} = 2, using N_{ring} = 2400 radial rings, with N_{φ} = 720 azimuthal points. In order to extract the properties of the mode present within the disc, one may proceed as follows. For each simulation snapshot, one can estimate the active surface density within the disc via (C.2)where the sum on i is made on all the particles of the simulation and x_{i}(t) is the position of the ith particle at time t. Such a surface density can be decomposed on the basis elements from Eq. (3) in the form (C.3)where the sum on p is made on all the basis elements Σ^{(p)} considered. The effective basis elements used during our measurements are the same as the ones used in the matrix method from Sect. 4.2. Thanks to the biorthogonality property from Eq. (3), the coefficients b_{p}(t) can be immediately determined as (C.4)As we are looking for unstable modes within the disc, we expect to have b_{p}(t) ∝ exp [ −i(ω_{0} + iη)t ], where ω_{0} = m_{φ}Ω_{p} is the pattern speed of the mode and η = s its growth rate. As a consequence, one immediately obtains that (C.5)if one is sufficently careful with the branchcut of the complex logarithm. Such a linear scaling with t of Re(log (b_{p}(t))) and Im(log (b_{p}(t))) is therefeore the appropriate measurement procedure to use in order to estimate the growth rate and pattern speed of the unstable modes of these truncated Mestel discs. The measurements for the various values of the truncation index ν_{t} are illustrated in Fig. C.3. Once the basis coefficients b_{p}(t) have been determined, one can study the shape of the recovered unstable modes in the physical space. Indeed, similarly to Eq. (C.1), the shape of the modes is given by (C.6)In analogy with Fig. C.2, for which the unstable modes have been obtained via the matrix method, Fig. C.4 illustrates the unstable mode of the same truncated ν_{t} = 4 Mestel disc.
Fig. C.4 Unstable mode for the truncated ν_{t} = 4 Mestel disc recovered via direct Nbody simulations as presented in Sect. 5. Only positive contour levels are shown and are spaced linearly between 20% and 80% of the maximum norm. Similarly to Fig. C.2, the radii associated with the resonances ILR, COR and OLR have been represented. 

Open with DEXTER 
Measurements of the pattern speed ω_{0} = m_{φ}Ω_{p} and growth rate η = s for unstable m_{φ} = 2 modes of tapered Mestel discs.
As a conclusion, the growth rates and pattern speeds obtained either via the matrix method or direct Nbody simulations are gathered in Table C.1. As observed in Sellwood & Evans (2001), the recovery of the unstable modes characteristics from direct Nbody simulations when performed for truncated Mestel discs is a difficult task, for which convergence to the values predicted through linear theory may be difficult.
Appendix D: Why swing matters?
Let us investigate here briefly the importance of selfgravitation and the completeness of the projection basis in capturing the role of swing amplification.
Fig. D.1 Map of , corresponding to the bare secular diffusion flux, using the same conventions as in Fig. 4. The contours are spaced linearly between the minimum and the maximum of . The maximum value for the positive blue contours corresponds to , while the minimum value for the negative red contours is associated with . This figure is qualitatively similar to the one obtained in Fig. 9 of Paper I. 

Open with DEXTER 
Appendix D.1: Turning off the selfgravitating amplification
In order to investigate the role of the selfgravitating amplification, one may perform the same estimation as presented in Fig. 4, while neglecting collective effects. When neglecting collective effects, i.e. when assuming that , one recovers the inhomogeneous Landau equation (Chavanis 2013) which reads (D.1)Equation (D.1) involves the bare susceptibility coefficients  A_{m1,m2}(J_{1},J_{2})  ^{2}, which can be equivalently defined (see Appendix B of Paper I) by (D.2)where u(x) is the binary potential of interaction given by u(x) = −G/ x  for gravity. This estimation therefore does not require to estimate the response matrix from Eq. (24), but one still has to perform integrations along the resonant lines as in Eq. (42). Let us provide a first numerical implementation of this equation in the context of galactic dynamics. The contours of are illustrated in Fig. D.1. Comparing the maps of the dressed diffusion flux Ndiv(ℱ_{tot}) from Fig. 4 and the bare diffusion flux , allows us to assess the strength of the selfgravitating amplification. As expected, when turning off the selfgravity of the system, one reduces significantly the susceptibility of the system and therefore slows down its secular evolution by a factor of about 1000. One may also remark that while the secular appearance of a resonant ridge in the dressed diffusion from Fig. 4 was obvious, the shape of the contours obtained in the bare Fig. D.1 do not emphasise as clearly the appearance of such a narrow resonant ridge. One can still remark that the structure of the bare contours obtained in Fig. D.1 is similar to what was obtained in Fig. 9 of Paper I, through the WKB limit of the BalescuLenard equation. One can finally note that the amplitudes of the bare divergence contours obtained previously are similar to the WKB values obtained in Paper I. As a consequence, the comparison of Figs. 4 and D.1 emphasises that the strong selfgravitating amplification of loosely wound perturbations is indeed responsible for the appearance of a narrow ridge, while also ensuring that this appearance is sufficiently rapid, as observed in the diffusion timescales comparison from Eq. (52).
Appendix D.2: Turning off looselywound contributions
As emphasised in the Introduction, the WKB limit of the BalescuLenard equation presented in Fouvry et al. (2015a) was not able to capture the mechanism of swing amplification, which involves unwinding perturbations. By considering a complete and global basis as in Eq. (16), we have shown in Fig. 4 how the missing amplification from Fouvry et al. (2015a) could be recovered. Using the numerical method of estimation of the secular diffusion flux as presented in Sect. 3, one can try to recover the results obtained within the WKB formalism by carefully choosing the considered basis elements generically introduced in Eq. (16) and chosen to be Kalnajs basis elements as detailed in Appendix A. We recall that each basis element depends on two indices: an azimuthal index ℓ and a radial one n. Because in S12’s simulation perturbations were restricted to the harmonic sector m_{φ} = 2, one only has to consider basis elements associated with ℓ = 2. Moreover, as illustrated in Fig. D.2, the larger n the radial index, the faster the radial variation of the basis elements and therefore the more tightly wound the basis elements. So as to get rid of the looselywound basis elements which are the ones which can get swing amplified, we perform a truncation of the radial indices considered. Therefore, we define the secular diffusion flux computed in the same way than Ndiv(ℱ_{tot}) as presented in Sect. 4.2, except that the basis elements are such that n_{cut} ≤ n ≤ n_{max}, with n_{cut} = 2 and n_{max} = 8. By keeping only the tightly wound basis elements, one can therefore consider the same contribution as the one considered in the WKB limit presented in Paper I. The contours of are illustrated in Fig. D.3. One can note that the values of the contours obtained in the map of illustrated in Fig. D.3 are in the same order of magnitude as the ones which were presented in Fig. 9 of Paper I in the WKB limit. The presence of positive blue contours of is also in agreement with a secular heating of the disc (i.e. an increase of J_{r}). However, these contours do not display a narrow resonant ridge as was observed in S12’s simulation or in Fig. 4.
Fig. D.2 Illustration of the radial basis elements of the ℓ = 2 Kalnajs basis elements for k_{Ka} = 7, defined in Appendix A, which were used in the estimation of the BalescuLenard diffusion flux in Sect. 4.2. As the radial index n increases, the basis elements get more and more wound. 

Open with DEXTER 
Fig. D.3 Map of , corresponding to the dressed secular diffusion flux, using the same conventions as in Fig. 4. In order to limit ourselves only to tightly wound contributions, the basis elements associated with the radial basis index n ∈ {0,1 } have not been taken into account. The contours are spaced linearly between the minimum and the maximum of . The maximum value for the positive blue contours corresponds to , while the minimum for the negative red contours is associated with . This figure is to be compared to Fig. 9 of Paper I. 

Open with DEXTER 
Appendix E: Sampling of the DF
In order to use the Nbody integrator described in Sect. 5, one has to sample the particles according to the DF given by Eq. (49). We introduce the probability distribution function F_{sp}, normalised to 1 and thanks to which the sampling is performed. This probability DF F_{sp} is directly proportional to the active distribution function F_{star} from Eq. (49), so that we may write (E.1)where C_{sp} is a normalisation constant which is determined in the upcoming calculations. Because the mapping (E,L) → (J_{r},J_{φ}) from Eq. (12) is not a trivial one, we do not perform the sampling of the stars in the action space (J_{r},J_{φ}), but rather in the (E,L)space. Moreover, one should pay attention to the fact that the DF F_{sp} from Eq. (E.1) is a probability distribution function in the (x,v)space, so that d^{2}xd^{2}vF_{sp}(x,v) is proportional to the number of particles in the infinitesimal volume d^{2}xd^{2}v around the position (x,v). As we want to sample the particles in the (E,L)space, we introduce the function h_{sp}(E,L) such that dEdLh_{sp}(E,L) is proportional to the number of particles in the volume dEdL around the location (E,L). One can now determine h_{sp}(E,L) as a function of F(E,L). Indeed, we have (E.2)using the fact that the tangential velocity satisfies v_{t} = L/r. The last step is then to perform the change of variable v_{r} → E. One has , so that . Because the radial velocity can be both positive and negative, Eq. (E.2) takes the form (E.3)where we used the definition (13) of the radial intrinsic frequency Ω_{1}. One can then correctly normalise the probability distribution h_{sp} and determine the value of the constant C_{sp} from Eq. (E.1). One should pay attention to the fact that in addition to the tapering functions T_{inner} and T_{outer} from Eqs. (48), we also assume that no stars have orbits that extend beyond R_{max}. As a consequence, the allowed region in the (E,L)space has to satisfy two constraints. First of all, the angular momentum L_{star} has to satisfy (E.4)Then, for a given value of L_{star}, one can show that the energy of the star E_{star} must satisfy the constraint (E.5)These two constraints allow us to completely characterise the (E,L)space on which the sampling (E,L) has to be performed. One finally has to satisfy the constraint (E.6)Given the parameters presented after Eq. (49), one can numerically determine the value of the constant C_{sp} which reads(E.7)We may now proceed to the sampling of the coordinates of the particles. Up to the sign of its radial velocity, one star is characterised by the set {E_{star},L_{star},R_{star},φ_{star} }. Given that the initial state is axisymmetric, the azimuthal angle of the star can be uniformly sampled between 0 and 2π. The next step is then to successively sample (L_{star},E_{star}) and finally R_{star}, using successive rejection samplings as we now detail.
The heart of the rejection sampling is as follows. Let us assume that we want to generate sampling values from a function f(x), from which it is difficult to sample. However, we assume that we have at our disposal another distribution function g(x) from which the sampling is simple, and such that there exists a bound M> 1 satisfying f(x) <Mg(x). The smaller M, the more efficient the sampling. One then has to proceed as follows: sample both a proposal x from g and u uniformly between [ 0;1 ]. One then applies the selection (E.8)In order to have an efficient sampling, one should try to consider a function g close to f.
We may now directly sample (E,L) thanks to this algorithm. The true sampling function is f_{(E,L)} = h_{sp} from Eq. (E.2). The simple sampling function is g_{(E,L)} ∝ 1, defined on the domain characterised by the constraints from Eqs. (E.4) and (E.5). When performing a rejection sampling with such a uniform g_{(E,L)}, in order to determine the bound M_{(E,L)}, one only has to determine a uniform bound for f_{(E,L)}. With the numerical values introduced after Eq. (49), one can check that f_{(E,L)} is such that (E.9)The final element required to be able to perform the rejection sampling with f_{(E,L)} is to be able to draw uniformly candidate (E,L) in the domains defined by the constraints from Eqs. (E.4) and (E.5), which is equivalent as sampling candidates (E,L) from the uniform probability distribution function g_{(E,L)}. To perform this uniform sampling, since the constraints from Eq. (E.5) are expressed for a given value of L_{star}, it is more natural to first draw L_{star} and then E_{star}. The probability distribution according to which L_{star} has to be drawn is of the form f_{L}(L) ∝ (E_{max}(L) − E_{min}(L)). When correctly normalised, it reads (E.10)To sample L from f_{L}, we use another rejection sampling by introducing the additional simple probability distribution function g_{L} defined as (E.11)It is straightforward to check that f_{L}< (3 / 2) g_{L}, so that we may use the bound M_{L} = 3 / 2 to perform the rejection sampling of L_{star}. The final remark is to note that sampling L from g_{L} is simple since its cumulative distribution function can be inverted so as to read (E.12)where W_{1} is the lower branch of the Lambert function W(x) for x ∈ [ −1 / e ;0 ]. With all these elements, the rejection sampling of L_{star} following f_{L} from Eq. (E.10) can be performed.
Once L_{star} has been drawn, it only remains to sample uniformly E_{star} on the interval E_{star} ∈ [ E_{min}(L_{star}) ;E_{max}(L_{star}) ], as given by Eq. (E.4). Thanks to these uniformly drawn candidates (E_{star},L_{star}) and the uniform bound from Eq. (E.9), one can perform the rejection sampling from the probability distribution f_{(E,L)}.
For E_{star} and L_{star} succesfully sampled, one may then sample the radius R_{star} using a similar rejection sampling. The radius has to be sampled according to the probability distribution f_{R} given by (E.13)so that one has f_{R} ∝ v_{r}. However, one should note that for r → r_{p/a}, one has f_{R}(r) → + ∞, so that the rejection sampling cannot be used without considering a probability DF g which also diverges for r → r_{p/a}. In order to get rid of these divergences, instead of sampling the variable r, we sample the angle u ∈ [ −π/2 ;π/2 ], where we defined the mapping r → u(r) as (E.14)so that one naturally has r(−π/2) = r_{p} and r(π/2) = r_{a}. The probability distribution function from which u has to be sampled is immediately given by (E.15)Using the fact that the maximum of f_{u} is reached for u = π/2, one can then sample u from f_{u} using a rejection sampling with a uniform control probability distribution function g_{u}(u) = 1 /π. Once u is known, it only remains to compute R_{star} = r(u_{star}), so that the sampling of all the required quantities for one star has been performed.
The final step of the sampling of the particles is to determine the physical coordinates of the particles (x,v) associated with the set {E_{star},L_{star},R_{star},φ_{star} }. These physical coordinates are the ones which are given to the Nbody integrator. We draw uniformly the sign of the radial velocity ε_{r} ∈ {−1,1 }. Because we are considering a disc made only of prograde stars, one immediately obtains that the radial and tangential velocities v_{r} and v_{t} are given by (E.16)The final step of the transformation to the (x,v)coordinates is then straightforward, since one naturally has (E.17)One should note that the sampling procedure described previously does not correspond to a quiet start procedure (Sellwood 1983), which would lead to a reduction of the initial shot noise within the disc, as briefly discussed in Appendix C, with regard to the validation of the Nbody code.
Appendix F: Another test of the scaling with N
One difficulty with the measurement of the scaling with N presented in Sect. 5.2 is that one has to disentangle the contributions from the initial sampling Poisson shot noise present through h_{0} from Eq. (59) and the effects due to the collisional BalescuLenard diffusion scaling through h_{2} from Eq. (64). Indeed, Poisson shot noise leads to fluctuations of the system DF about its mean value. In order not to be sensitive to such fluctuations, one could only consider fluctuations sufficiently large, i.e. fluctuations caused by an effective secular diffusion rather than caused by inevitable Poisson fluctuations. As a consequence, by restricting ourselves only to large fluctuations, we can get rid of Poisson’s effects. We therefore define the function as (F.1)where we introduced a threshold . Here χ [·] is a characteristic function equal to 1 if , and 0 otherwise. As a consequence, measures the volume in action space of the regions (depleted from particles, since ) for which the mean DF has changed by more than . For a sufficiently large value of the threshold , such a construction allows us not to be polluted by Poisson sampling shot noise. For the initial times, as in Eq. (50), it is straightforward to study the scaling of with t and N.
Fig. F.1 Top panel: illustration of the behaviour of the function from Eq. (F.1), when averaged on 32 different realisations for particles numbers N ∈ {8, 12, 16, 24, 32, 48, 64 } × 10^{5}, along with the associated linear fits. To compute , we used the same binning of the actionspace (J_{φ},J_{r}) as in Fig. 5. As obtained in Eq. (F.3), one recovers that for a fixed value of N, the function is linear. The horizontal dashed line illustrates the threshold value for which the threshold time t_{thold} is determined. Bottom panel: illustration of the behaviour of the function N → t_{thold}(N). As derived in Eq. (F.5), one recovers a linear dependence of t_{thold} with N. 

Open with DEXTER 
Indeed, one can write (F.2)\newpage
Introducing , one can rewrite Eq. (F.2) as (F.3)Therefore, for a fixed value of N, one expects to observe a linear time dependence of the function , as illustrated in Fig. F.1. In order to test the scaling of Eq. (F.3) with N, one may proceed as follows. Introducing a threshold value , for each value of N, we define the associated threshold time t_{thold}(N) as(F.4)Thanks to the scalings from Eq. (F.3), one immediately obtains that (F.5)Such a linear scaling of t_{thold}(N) with N is a prediction from the BalescuLenard formalism and is nicely recovered in Fig. F.1.
Appendix G: Distributed code description
For the sake of reproducibility, which has been lacking in the context of the linear response of stellar systems, we distribute the linear matrix response code we wrote for this paper both as a Mathematica package^{6}, and a notebook^{7}. The functions therein allow for:

The determination as a function of (r_{p},r_{a}) of the orbits quantites: E, L, J_{r}, Ω_{1} and Ω_{2}.

The construction of the 2D basis from Kalnajs (1976).

The computation of the Fourier transform w.r.t. the angles, i.e. the computation of from Eq. (21).

The calculation of the 2D response matrix via Eq. (33).
It has been tested for the isochrone and the Mestel discs. A copy of the code is also available at the CDS.
All Tables
Measurements of the pattern speed ω_{0} = m_{φ}Ω_{p} and growth rate η = s for unstable m_{φ} = 2 modes of tapered Mestel discs.
All Figures
Fig. 1 Contours of the initial distribution function F_{star} from Eq. (49), in action space (J_{φ},J_{r}). The contours are spaced linearly between 95% and 5% of the distribution function maximum. 

Open with DEXTER  
In the text 
Fig. 2 Illustration of 4 different critical resonant lines in the (r_{p},r_{a})space. As defined in Eq. (41), a critical line is characterised by the resonant vectors m_{1}, m_{2} and a location in action space. Each of the 4 plotted critical lines are associated with the same location , represented by the black dot. The critical lines correspond to various choices of resonant vectors m_{1} and m_{2} among the three inner, outer and corotation Lindblad resonances, respectively ILR, OLR and COR. One should also note that for m_{1} = m_{2}, the critical lines go through the point . 

Open with DEXTER  
In the text 
Fig. 3 Map of −Nℱ_{tot}, where the flux has been computed with m_{1},m_{2} ∈ {m_{ILR},m_{COR},m_{OLR}}. Following Eq. (11), −Nℱ_{tot} corresponds to the direction of diffusion of individual particles in actionspace. 

Open with DEXTER  
In the text 
Fig. 4 Left panel: map of Ndiv(ℱ_{tot}), where the total flux has been computed with m_{1},m_{2} ∈ {m_{ILR},m_{COR},m_{OLR}}. Red contours, for which Ndiv(ℱ_{tot}) < 0 are associated with regions from which the orbits will be depleted, whereas blue contours, for which Ndiv(ℱ_{tot}) > 0 correspond to regions where the value of the DF will be increased during the secular diffusion. The contours are spaced linearly between the minimum and the maximum of Ndiv(ℱ_{tot}). The maximum value for the positive blue contours corresponds to Ndiv(ℱ_{tot}) ≃ 350, while the mininum value for the negative red contours is associated with Ndiv(ℱ_{tot}) ≃ −250. Right panel: from Sellwood (2012) – Fig. 7, contours of the change in the DF between the time t_{S12} = 1400 and t_{S12} = 0, for a run with 50M particles. Similarly to the left panel, red contours correspond to negative differences, i.e. regions emptied from their orbits, while blue contours correspond to positive differences, i.e. regions where the DF has increased during the diffusion. Both of these contours are aligned with the ILR direction of m_{ILR} = (2,−1) in the (J_{φ},J_{r})plane, corresponding to the cyan line. 

Open with DEXTER  
In the text 
Fig. 5 Illustration of the behaviour of the function defined in Eq. (55), for an active fraction ξ = 0.5, when averaged on 32 different realisations for particles numbers N ∈ {8, 12, 16, 24, 32, 48, 64 } × 10^{5}. To compute , we binned the actionspace domain (J_{φ},J_{r}) = [ 0;2.5 ] × [ 0;0.2 ] in 100 × 50 regions. The values of have also been uniformly renormalised so as to clarify this representation. The dots corresponds to the snapshots of the simulations for which was computed, whereas the lines correspond to secondorder fits. As expected, the smaller the number of particles, the noisier the simulation and the larger . 

Open with DEXTER  
In the text 
Fig. 6 Top panel: illustration of the behaviour of the function , where N has been rescaled by a factor 10^{5} so as to simplify the representation. The dots correspond to computed values thanks to Fig. 5, while the line corresponds to a linear fit, which takes the form . The coefficients have been uniformly renormalised so as to clarify this representation. Bottom panel: similar representation for the behaviour of the function , whose linear fit takes the form . Similarly, the coefficients have been uniformly renormalised so as to clarify this representation. 

Open with DEXTER  
In the text 
Fig. 7 Top panel: illustration of the behaviour of the function for an active fraction ξ = 0.6, using the same conventions as in Fig. 5. As expected, for an increased active fraction, the secular evolution of the system is fastened leading to a faster increase of the distance . Middle panel: behaviour of the function , for an active fraction ξ = 0.6, using the same conventions as in Fig. 6. Its linear fit takes the form . One recovers the expected scaling of the Poisson shot noise sampling derived in Eq. (59). Bottom panel: behaviour of the function , for an active fraction ξ = 0.6, using the same conventions as in Fig. 6. Its linear fit takes the form . One recovers the expected collisional scaling with N obtained in Eq. (64). 

Open with DEXTER  
In the text 
Fig. 8 Map of Ndiv(ℱ_{tot}) with an active fraction ξ = 0.6 and using the same conventions as in Fig. 4. The contours are spaced linearly between the minimum and the maximum of Ndiv(ℱ_{tot}). The maximum value for the positive blue contours corresponds to Ndiv(ℱ_{tot}) ≃ 4200, while the minimum value for the negative red contours corresponds to Ndiv(ℱ_{tot}) ≃ −3200. As expected, one recovers that when the active fraction of the disc is increased, the susceptibility of the disc is increased, so that the norm of Ndiv(ℱ_{tot}) gets larger and the secular diffusion is fastened. 

Open with DEXTER  
In the text 
Fig. 9 Behaviour of the function as defined in Eq. (68), for various values of the number of particles. The prefactor has been added so as to mask Poisson shot noise, allowing the initial values of to be independent of N. It illustrates the bifurcation between the initial BalescuLenard collisional evolution, for which low values of Σ_{2} are expected and the collisionless Vlasov evolution for which the system is no more axisymmetric leading to larger values of Σ_{2}. As expected, the larger the number of particles, the later the transition. 

Open with DEXTER  
In the text 
Fig. 10 Illustration of the active surface density Σ_{t} for a Nbody run with N = 8 × 10^{5}, restriced to the range R ≤ 6. Top panel: active surface density Σ_{t} at an early time t = 60, for which the galaxy remains globally axisymmetric. In this regime, the dynamics of the system is collisional and governed by the BalescuLenard Eq. (2). Bottom panel: active surface density Σ_{t} at a much later time t = 2400. The galaxy is then strongly nonaxisymmetric. In this regime, the dynamics of the system is collisionless and governed by the Vlasov equation. 

Open with DEXTER  
In the text 
Fig. C.1 Top panel: zoomed Nyquist contours in the complex plane of obtained via the matrix method for a truncated Mestel disc with ν_{t} = 6 and q = 6, looking for m_{φ} = 2 modes. Each contour corresponds to a fixed value of η. For a growth rate of η ≃ 0.20, one can note that the contour crosses the origin, which corresponds to the presence of an unstable mode. Bottom panel: illustration of the function for the same truncated Mestel disc. Each line corresponds to a fixed value of η. Such a representation allows us to determine the pattern speed ω_{0} = m_{φ}Ω_{p} ≃ 0.94 of the unstable mode. 

Open with DEXTER  
In the text 
Fig. C.2 Unstable mode for the truncated ν_{t} = 4 Mestel disc recovered via the matrix method presented in Sect. 3. Only positive contour levels are shown and are spaced linearly between 10% and 90% of the maximum norm. The radii associated with the resonances ILR, COR and OLR have been represented, as given by ω_{0} = m·Ω(R_{m}), where the intrinsic frequencies Ω(R) = (Ω_{φ}(R),κ(R)) have to be computed within the epicyclic approximation. For a Mestel disc, they are given by Ω_{φ}(R) = V_{0}/R and . 

Open with DEXTER  
In the text 
Fig. C.3 Measurement of the growth rate η and pattern speed ω_{0} of the m_{φ} = 2 unstable modes of truncated Mestel discs with a random velocity given by q = 6 for various values of the truncation index ν_{t} = 6, 8. The basis coefficient plotted is associated with the basis element (ℓ,n) = (2,0), using the same basis elements as for the matrix method in Sect. 4.2. 

Open with DEXTER  
In the text 
Fig. C.4 Unstable mode for the truncated ν_{t} = 4 Mestel disc recovered via direct Nbody simulations as presented in Sect. 5. Only positive contour levels are shown and are spaced linearly between 20% and 80% of the maximum norm. Similarly to Fig. C.2, the radii associated with the resonances ILR, COR and OLR have been represented. 

Open with DEXTER  
In the text 
Fig. D.1 Map of , corresponding to the bare secular diffusion flux, using the same conventions as in Fig. 4. The contours are spaced linearly between the minimum and the maximum of . The maximum value for the positive blue contours corresponds to , while the minimum value for the negative red contours is associated with . This figure is qualitatively similar to the one obtained in Fig. 9 of Paper I. 

Open with DEXTER  
In the text 
Fig. D.2 Illustration of the radial basis elements of the ℓ = 2 Kalnajs basis elements for k_{Ka} = 7, defined in Appendix A, which were used in the estimation of the BalescuLenard diffusion flux in Sect. 4.2. As the radial index n increases, the basis elements get more and more wound. 

Open with DEXTER  
In the text 
Fig. D.3 Map of , corresponding to the dressed secular diffusion flux, using the same conventions as in Fig. 4. In order to limit ourselves only to tightly wound contributions, the basis elements associated with the radial basis index n ∈ {0,1 } have not been taken into account. The contours are spaced linearly between the minimum and the maximum of . The maximum value for the positive blue contours corresponds to , while the minimum for the negative red contours is associated with . This figure is to be compared to Fig. 9 of Paper I. 

Open with DEXTER  
In the text 
Fig. F.1 Top panel: illustration of the behaviour of the function from Eq. (F.1), when averaged on 32 different realisations for particles numbers N ∈ {8, 12, 16, 24, 32, 48, 64 } × 10^{5}, along with the associated linear fits. To compute , we used the same binning of the actionspace (J_{φ},J_{r}) as in Fig. 5. As obtained in Eq. (F.3), one recovers that for a fixed value of N, the function is linear. The horizontal dashed line illustrates the threshold value for which the threshold time t_{thold} is determined. Bottom panel: illustration of the behaviour of the function N → t_{thold}(N). As derived in Eq. (F.5), one recovers a linear dependence of t_{thold} with N. 

Open with DEXTER  
In the text 