Secular diffusion in discrete selfgravitating tepid discs
I. Analytic solution in the tightly wound limit^{⋆}
^{1} Institut d’Astrophysique de Paris and UPMC, CNRS (UMR 7095), 98bis boulevard Arago, 75014 Paris, France
email: fouvry@iap.fr; pichon@iap.fr
^{2} Institute of Astronomy & KICC, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK
^{3} Laboratoire de Physique Théorique (IRSAMC), CNRS and UPS, Univ. de Toulouse, 31062 Toulouse, France
email: chavanis@irsamc.upstlse.fr
Received: 19 February 2015
Accepted: 2 May 2015
The secular evolution of an infinitely thin tepid isolated galactic disc made of a finite number of particles is described using the inhomogeneous BalescuLenard equation. Assuming that only tightly wound transient spirals are present in the disc, a WKB approximation provides a simple and tractable quadrature for the corresponding drift and diffusion coefficients. It provides insight into the physical processes at work during the secular diffusion of a selfgravitating discrete disc and makes quantitative predictions on the initial variations of the distribution function in action space. When applied to the secular evolution of an isolated stationary selfgravitating Mestel disc, this formalism predicts the initial importance of the corotation resonance in the inner regions of the disc leading to a regime involving radial migration and heating. It predicts in particular the formation of a ridgelike feature in action space, in agreement with simulations, but overestimates the timescale involved in its appearance. Swing amplification is likely needed to resolve this discrepancy. In astrophysics, the inhomogeneous BalescuLenard equation and its WKB limit may also 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.
Key words: galaxies: evolution / galaxies: kinematics and dynamics / galaxies: spiral / diffusion / gravitation
Appendices are available in electronic form at http://www.aanda.org
© ESO, 2015
1. Introduction
Understanding the dynamical evolution of galactic discs over cosmic times is a longstanding endeavour. Selfgravitating discs are cold dynamical systems, for which rotation represents an important reservoir of free energy. Fluctuations in the potential induced by discrete (possibly distant) encounters may be strongly amplified, while resonances tend to confine and localise their dissipation: such small stimuli can lead to the spontaneous evolution of discs to distinct equilibria. Quantifying the relative importance of this intrinsically driven evolution with respect to that driven by the environment is of renewed interest now that their cosmological environment is firmly established in the context of the ΛCDM paradigm. The effect of intrinsic susceptibility on secular timescales can be addressed in the context of kinetic theory, which takes explicitly into account such interactions.
The kinetic theory of stellar systems is an old, yet fundamental, topic in astrophysics^{1}. It was initiated by Jeans (1929) and Chandrasekhar (1942) in the context of 3D stellar systems such as elliptical galaxies and globular clusters. The kinetic theory of Coulombian plasmas was developed in parallel by Landau (1936) and Vlasov (1938). When encounters between the particles (stars or electric charges) are neglected, one gets a purely mean field equation (Jeans 1915; Vlasov 1938) called the collisionless Boltzmann equation, or the Vlasov equation. When encounters are taken into account, one gets a kinetic equation that includes a collision term. In early works, the collision term was obtained by assuming that a particle experiences a succession of independent twobody encounters with the other particles. The corresponding kinetic equation can be derived either from the Boltzmann equation by considering a limit of weak deflections (Landau 1936), or directly from the general form of the FokkerPlanck equation by evaluating the diffusion and drift coefficients in a binary collisions approximation (Chandrasekhar 1949; Rosenbluth et al. 1957). In the case of neutral plasmas, the system is spatially homogeneous, so the distribution function depends only on the velocity, hence the name kinetic theory. By contrast, stellar systems are spatially inhomogeneous, so the distribution function depends on position and velocity. In early works on stellar dynamics, spatial inhomogeneity was taken into account in the advection term (Vlasov), but the collisional term was calculated by making a local approximation as if the system were homogeneous. In 3D, the collisional term displays a logarithmic divergence at large scales (Jeans 1929; Landau 1936; Chandrasekhar 1942). In the case of plasmas, this divergence is due to the neglect of collective effects that are responsible for Debye shielding. Landau (1936) phenomenologically introduced a cutoff at the Debye length to regularize the divergence.
Later on, Balescu (1960) and Lenard (1960) developed a rigorous kinetic theory of plasmas, taking collective effects into account, and obtained a kinetic equation, the socalled BalescuLenard equation, that does not present any divergence at large scales. The Debye shielding is taken naturally into account in their treatment through the dielectric function (which is absent from the Landau equation). In the case of stellar systems, the divergence at large scales is solved by the spatial inhomogeneity of the system and its finite extent. One can phenomenologically introduce a cutoff at the Jeans scale (Weinberg 1993), i.e. at the system’s size, which would correspond to the analogue of the Debye length in plasma physics, but this ad hoc treatment is not fully satisfactory. Furthermore, it cannot be applied to cold (centrifugally supported) stellar discs where spatial inhomogeneity is more crucial than in 3D.
A more fruitful procedure is to write the kinetic equation with angleaction variables that are the appropriate variables to describe spatially inhomogeneous multiperiodic systems. When collective effects are neglected, one obtains the inhomogeneous Landau equation (Chavanis 2007, 2013b). When collective effects are accounted for, one gets the inhomogeneous BalescuLenard equation (Heyvaerts 2010; Chavanis 2012a). The collision term in these equations describes the effect of (possibly distant) encounters between stars. For selfgravitating systems, where the interaction is attractive (instead of being repulsive as in Coulombian plasmas), collective effects are responsible for an antishielding which tends to increase the effective mass of the stars, hence reducing the relaxation time. 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 (encounters between stars) into account and describes the evolution of the system on a timescale of the order Nt_{D}, where t_{D} is the dynamical time. For times t ≪ Nt_{D}, or for N → + ∞, it reduces to the Vlasov equation which ignores encounters between stars. Although the kinetic theory was initially developed for 3D stellar systems, the final form of the inhomogeneous BalescuLenard equation also applies to stellar discs such as those considered in this paper.
The BalescuLenard nonlinear equation accounts for selfdriven orbital secular diffusion of a selfgravitating system induced by the intrinsic shot noise due to its discreteness and the corresponding long range correlations. Even though this kinetic equation was first written down more than fifty years ago, it has hardly ever been applied in its prime context, but only in various limits where it reduces to simpler kinetic equations, as discussed above.
In this paper, we will focus on solving explicitly such an equation describing the selfgravitating response of a tepid thin disc to its own stochastic fluctuating potential induced by its finite number of components. In this cool regime, the selfgravity of the disc can be tracked via a local WKBlike response, which in turn allows us to simplify the a priori 2D formalism to an effective (nondegenerate) 1D formalism. We will compare the prediction of the WKB limit to a numerical experiment presented in the literature, and discuss its diagnosis power and possible limitations.
The paper is organised as follows. Section 2 briefly presents the content of the inhomogeneous BalescuLenard equation. Section 3 focuses on razor thin axisymmetric galactic discs within the WKB approximation. Section 4 investigates the formation of a narrow resonant ridge in an isolated selfgravitating Mestel disc. Finally, Sect. 5 wraps up. Appendix A provides a short sketch of the derivation of the BalescuLenard equation. Appendix B considers the inhomogeneous BalescuLenard equation without collective effects. Appendix D compares it to other similar kinetic equations, and in particular its FokkerPlanck limit.
2. The inhomogeneous BalescuLenard equation
We consider a system made of N particles. We suppose that the gravitational background, associated with the Hamiltonian H_{0}, is stationary and integrable, so that we may always remap the physical phasespace coordinates (x,v) to the angleaction coordinates (θ,J) (Goldstein 1950; Born 1960; Binney & Tremaine 2008). We also introduce the intrisinc frequencies of the system Ω defined as (1)Along the unperturbed trajectories, the angles θ are 2πperiodic evolving with the frequency Ω, whereas the actions J are conserved. To describe the longterm evolution of such a system, one assumes that there are two decoupled timescales: a short dynamical timescale and a secular timescale of collisional evolution. We assume that the system is always in a virialized stable state (i.e. is a stable stationary solution of the Vlasov equation), so that the distribution function can be written as a quasistationary distribution F = F(J,t), which satisfies the normalization constraint , where M_{tot} is the total mass of the system. This is a function of the actions only that slowly evolves in time due to stellar encounters (finiteN effects)^{2}. From Heyvaerts (2010) and Chavanis (2012a; see also Appendix A for a short sketch of the derivation), the secular evolution, induced by collisional finiteN effects, of such a quasistationary distribution function F(J,t) is given by the inhomogeneous BalescuLenard equation which reads (2)where 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}). The r.h.s. of Eq. (2)is the BalescuLenard operator which encompasses the secular diffusion due to collisional effects, see Fig. 1. Because it is the divergence of a flux, this writing ensures that the total number of stars is exactly conserved during the secular diffusion. The Dirac delta δ_{D}(m_{1}·Ω_{1} − m_{2}·Ω_{2}) is the sharp resonance condition. One must remark that this condition allows us to describe nontrivial gravitational interactions. Indeed, it can cause nonlocal resonances by coupling different regions of actionspace J_{1} and J_{2}.
Fig. 1 Top: a) a set of two resonant orbits in the inertial frame; b) in the rotating frame in which they are resonant – here through ILRCOR coupling. Bottom: c) fluctuations in the distribution function in actionspace caused by finiteN effects showing overdensities for the blue and red orbits. The dashed lines correspond to 3 contour levels of the intrinsic frequency ω = m·Ω respectively associated with the resonance vector m_{1} (gray lines) and m_{2} (black lines). The two sets of orbits satisfy the resonant condition m_{1}·Ω_{1} = m_{2}·Ω_{2}, and therefore lead to a secular diffusion of the orbital structure of the disc according to Eq. (2). One should emphasize that the resonant orbits need not be caught in the same resonance (m_{1} ≠ m_{2}), be close in position space, nor in action space. 

Open with DEXTER 
Even for local resonances (i.e. J_{1} = J_{2}), it can allow for nontrivial coupling of oscillations, as soon as m_{1} and m_{2} have nonzero components. The coefficients represent the dressed susceptibilities of the system, for which collective effects have been taken into account^{3}. To deal with the resolution of the nonlocal Poisson equation, following Kalnajs’ matrix method (Kalnajs 1976), one has to introduce a complete biorthonormal basis of potentials and densities ψ^{(p)}(x) and ρ^{(p)}(x) such that (3)Thanks to this basis, the susceptibility coefficients are given by (4)where I is the identity matrix and is the response matrix given by (5)In this expression, corresponds to the Fourier transform in angles of the basis elements ψ^{(p)}(x), where we used the convention that the Fourier transform of a function X(θ,J) is given by (6)One should note that Eq. (2)is therefore a nonlinear function of F, both explicitly via the products of F on the r.h.s., describing the effects of the binary collisions (as in the Boltzmann or Landau equations), but also implicitly via Eq. (5), which encompasses the collective effects and is specific to the BalescuLenard equation. We also importantly note that – in contrast to its counterpart in plasma physics – Eq. (2)does not assume that local encounters drive secular evolution; resonant interactions of correlated possibly distant dressed orbits are sourcing long term orbital distortions.
2.1. Content of the diffusion equation
One may also rewrite the BalescuLenard Eq. (2)as an anisotropic FokkerPlanck equation introducing the associated drift and diffusion coefficients. It then reads (7)where A_{m1}(J_{1}) and D_{m1}(J_{1}) are respectively the anisotropic drift and diffusion coefficients associated with a given resonance m_{1}, i.e. with a given Fourier mode m_{1} in angles. They both secularly depend on the distribution function F, but this dependence has not been explicitly written out to shorten the notations. The drift coefficients are given by (8)while the diffusion coefficients are given by (9)The rewriting from Eq. (7)allows us to discuss some important properties of such anisotropic diffusion equation. We introduce the total flux, ℱ_{tot}, associated with this diffusion, which reads (10)As a consequence, the BalescuLenard diffusion equation given by the expressions (2)and (7)takes the shortened form (11)We then define as M(t) the mass contained in a volume of the actionspace at time t, so that we have (12)Thanks to the divergence theorem, the variation of mass in due to secular diffusion corresponds to the flux of particles through the boundary of this volume so that (13)where n is the exterior pointing normal vector. One can note in Eq. (13)that the contribution from a given resonance m takes the form of a preferential diffusion in the direction m. This diffusion is therefore anisotropic because it is maximum for n ∝ m and equal to 0 for n·m = 0. To emphasize the anisotropy of the diffusion, one may use the formalism of the slow and fast actions (LyndenBell 1979; Earn & LyndenBell 1996). For simplicity, we consider the 2D case. For a given resonance m = (m_{1},m_{2}), we consider the change of coordinates (14)where and are respectively the slow and fast actions associated with the resonance m. Here m^{⊥} corresponds to the direction perpendicular to the resonance so that m^{⊥} = (m_{2}, − m_{1}), and . Thanks to the chain rule, for any function X(J), one has (15)Introducing the natural vector basis elements and associated with this change of coordinates, the diffusion flux ℱ_{m} associated with a resonance m and defined in Eq. (10)takes the form (16)Such a rewriting illustrates the fact that as soon as only one resonance m dominates the secular evolution, the diffusion flux will be aligned with this resonance. Hence one will observe a narrow monodimensional diffusion in the preferential direction. During this diffusion, particles will conserve their fast action , which can therefore be seen as an adiabatic invariant, whereas their slow action gets to change. This strong anistropy in the diffusion is an essential property of the BalescuLenard Eq. (7).
3. Thin tepid discs and their WKB limit
When implementing the inhomogeneous BalescuLenard equation, one encounters two main difficulties. The first one is the explicit construction of a mapping (x,v) → (θ,J) since the BalescuLenard drift and diffusion coefficients must be computed using angleaction coordinates. Assuming the disc to be tepid, one may rely on the epicylic approximation to construct such a mapping. The second difficulty arises from the nonlocality of Poisson’s equation which requires us to introduce potential basis elements ψ^{(p)} as in Eq. (3). One can then compute the response matrix from Eq. (5), which must subsequently be inverted. The following step is to compute the drift and diffusion coefficients from Eqs. (8)and (9)which requires us to explicitly deal with the resonance constraint m_{1}·Ω_{1} − m_{2}·Ω_{2} = 0. In the case of a 2D axisymmetric disc, one may implement a WKB approximation (Liouville 1837; Toomre 1964; Kalnajs 1965; Lin & Shu 1966; Palmer et al. 1989) which assumes that the diffusion of the system is made through tightly wound spirals. Such an assumption has two main consequences. First of all, Poisson’s equation becomes local, resulting in a diagonal response matrix. Moreover, it also entails that all the resonances becomes exactly local, allowing an explicit calculation of the resonant constraint δ_{D}(m_{1}·Ω_{1} − m_{2}·Ω_{2}). We now detail these two elements: epicyclic approximation and WKB assumption.
3.1. Epicyclic approximation
For a sufficiently cold disc (i.e. a disc where the radial excursions of the stars are small), one can explicitly build up a mapping (x,v) → (θ,J) thanks to the epicyclic approximation. We introduce the polar coordinates (R,φ) to describe the infinitely thin galactic disc, and introduce their associated momenta (p_{R},p_{φ}). As the disc at equilibrium is axisymmetric, the stationary Hamiltonian of the system reads (17)where ψ_{0} is the stationary background potential within the disc. Because ψ_{0} is axisymmetric, it does not depend on φ, so that p_{φ} is a conserved quantity. We may then define the first action of the system, the angular momentum J_{φ}, as (18)For a given value of J_{φ}, the equation of evolution of R is then given by (19)where ψ_{eff} is an effective potential defined as (20)The heart of the epicyclic approximation is to assume that small radial excursions can be approximated as harmonic librations. For a given value of J_{φ}, we implicitly introduce the guiding radius R_{g} as (21)Here R_{g}(J_{φ}) is the radius for which stars with an angular momentum of J_{φ} are on exactly circular orbits. It is important to note that the mapping between R_{g} and J_{φ} is bijective and unambiguous (up to the sign of J_{φ}). We may then define the two frequencies of evolution: Ω_{φ}(R_{g}) the azimuthal frequency and κ(R_{g}) the epicyclic frequency as follows (22)As the radial oscillations are supposed to be small, one may perform a Taylor expansion at first order of the evolution Eq. (19)in the neighbourhood of the minimum R = R_{g} so that R satisfies the differential equation . Hence one can notice that in this limit the evolution of the radius of a star is the one of a harmonic oscillator centred on R_{g}. Up to an initial phase, one has therefore R(t) = R_{g} + Acos(κt), where A is the amplitude of the radial oscillations. The associated radial action J_{r} is then given by (23)For J_{r} = 0, the orbit is circular. Within the epicyclic approximation, the frequencies of motion along the actiontorii, introduced in Eq. (1)are given by Ω(J) = (Ω_{φ}(J_{φ}),κ(J_{φ})). An important dynamical consequence of this approximation is that these two frequencies are only function of J_{φ} and do not depend on J_{r}, so that the resonance constraint m_{1}·Ω_{1} − m_{2}·Ω_{2} = 0 becomes simpler. Finally, one can explicitly construct the mapping between (R,φ,p_{R},p_{φ}) and (θ_{R},θ_{φ},J_{r},J_{φ}) (LyndenBell & Kalnajs 1972; Palmer 1994; Binney & Tremaine 2008), which takes at first order the form (24)Thanks to this mapping and the definitions of the actions from Eqs. (18)and (23), the epicyclic approximation allows us to build up an explicit mapping between the physical phasespace coordinates and the angleaction ones.
Finally, throughout our calculation, we will assume that the stationary distribution function of the disc is a Schwarzschild distribution function (or locally isothermal DF) given by (25)where Σ(R_{g}) is the surface density of the disc and , which varies within the disc, represents the radial velocity dispersion of the stars at a given radius. Increasing values of correspond to hotter discs that are therefore more stable.
3.2. The WKB basis
As we are considering a 2D case, the potential basis elements ψ^{(p)} introduced in Eq. (3)must be written as ψ^{(p)}(R,φ) in the disc polar coordinates and must be orthonormal to the associated surface density Σ^{(p)}(R,φ). Using a WKB approximation amounts to building up local basis elements thanks to which the response matrix will become diagonal.
3.2.1. Definition of the basis elements
We introduce the basis elements (26)where the window function ℬ_{R0}(R) is defined as (27)The basis elements are indexed by three numbers: k_{φ} is an azimuthal number which parametrizes the angular component of the basis elements, R_{0} is the radius position in the disc around which the Gaussian window ℬ_{R0} is centred, and k_{r} is the radial frequency of the basis element. We also introduced an additional parameter σ of scaleseparation, which will ensure the biorthogonality of the basis elements, as detailed later on. Finally, is an amplitude which will be tuned in order to normalize correctly the basis elements. Thanks to a somewhat unsual normalization of ℬ_{R0}, we will ensure that is independent of σ. Figure 2 illustrates the radial dependence of the basis elements.
Fig. 2 Two WKB basis elements. Each Gaussian is centred around a radius R_{0}. The typical extension of the Gaussian is given by the decoupling scale σ, and they are modulated at the radial frequency k_{r}. 

Open with DEXTER 
Figure 3 illustrates the shape of these basis elements in the polar (R,φ)plane. The next steps will be to ensure that these WKB basis elements have all the properties required to allow for the computation of the dressed susceptibility coefficients introduced in Eq. (4). Therefore, we will successively compute the associated surface density elements Σ^{[ kφ,kr,R0 ]}, ensure the biorthogonality of the basis elements and their correct normalization, and finally compute the Fourier transform in angles of the basis elements.
Fig. 3 Two WKB basis elements in the polar (R,φ)plane. Each basis element is located around a central radius R_{0}, on a region of size σ. The winding of the spirals is governed by the radial frequency k_{r}. The number of azimuthal patterns is given by the index k_{φ}, e.g. k_{φ} = 1 for the interior dark gray element, whereas k_{φ} = 2 for the exterior light gray one. 

Open with DEXTER 
3.2.2. Associated surface densities
In order to ensure the biorthogonality of the basis, we will first build up the surface densities associated with the potential elements introduced in Eq. (26). We extend the WKB potential in the zdirection using the Ansatz (28)Poisson’s equation in vacuum Δψ^{[ kφ,kr,R0 ]} = 0 leads to (29)We now explicitly introduce our WKB assumptions. We assume that the spirals are tightly wound so that (30)Introducing the typical size of the system R_{sys}, we also additionally suppose that we have (31)Assuming that k_{φ} is of the order of unity, Eq. (29)becomes (32)Hence within the WKB limit, the extended potential from Eq. (28)takes the form (33)where we ensured that for z → ± ∞ the potential tends to 0, therefore introducing a discontinuity for ∂ψ/∂z in z = 0. Thanks to Gauss theorem, the associated surface density satisfies (34)so that we have (35)
3.2.3. Biorthogonality condition and normalization
One must now ensure that the basis elements introduced in Eqs. (26)and (35)form a biorthogonal basis as assumed in Eq. (3). Indeed, it has to satisfy the property (36)The r.h.s. of this expression takes the form (37)The integration on φ is straightforward and is equal to . In order to perform the integration on R, we have to introduce additional assumptions to ensure the biorthogonality of the basis. The peaks of the Gaussians in Eq. (37)can be considered as separated if satisfies the condition (38)Under this assumption^{4}, the term from Eq. (37)can be assumed to be nonzero only for . The r.h.s. of Eq. (36)then takes the form (39)The integration on R takes the form of a radial Fourier transform of a Gaussian of spread σ at the frequency . It is therefore proportional to exp [ − (Δk_{r})^{2}/ (4 /σ)^{2} ]. Hence we will suppose that the frequency spread Δk_{r} satisfies (40)Under this assumption, the term from Eq. (39)is nonzero only for . In order to have a biorthogonal basis, one must therefore consider a spread σ, central radii R_{0}, and radial frequencies k_{r} such that (41)With these constraints, one must necessarily have , and in order to have a nonnegligible term in Eq. (36). It then only remains to explicitly estimate the amplitude of the basis elements. Equation (36)gives (42)Thanks to the WKB assumption from Eq. (41), the integration can be straightforwardly computed and leads to (43)
3.2.4. Fourier transform in angles
In order to estimate the susceptibility coefficients and the response matrix from Eqs. (4)and (5), one has to be able to calculate for the WKB basis elements. Thanks to the explicit mapping from Eq. (24), we have to compute (44)The integration on θ_{φ} is straighforward and equal to . Regarding the dependence on θ_{R} in the complex exponential, we may write (45)where the amplitude H_{kφ}(k_{r}) and the phase shift are given by (46)For typical galaxies, we have 1 / 2 ≤ Ω_{φ}/κ ≤ 1 (Binney & Tremaine 2008). Assuming that the azimuthal number k_{φ} is of the order of unity, one can use the WKB hypothesis introduced in Eq. (30), so that Eqs. (46)can be approximated as (47)Because we have assumed that the disc is tepid, the radial oscillations are small so that A ≪ R_{g}. We may then get rid of the dependences on A in ℬ_{R0}(R_{g} + Acos(θ_{R})) by replacing it with ℬ_{R0}(R_{g}), so that the only remaining dependence on A will be in the complex exponentials. To be able to explicitly perform the remaining integration on θ_{R} in Eq. (44), we introduce the Bessel functions of the first kind which satisfy the relation (48)We then finally obtain the expression of the Fourier transform in angles of the WKB basis elements which reads (49)
3.3. Estimation of the response matrix
Thanks to the WKB basis introduced in Eq. (26), one can now explicitly compute the response matrix from Eq. (5). Indeed, we use the expression (49)of the Fourier transform of the basis elements, and after simplification of the phaseshift terms thanks to the approximation from Eq. (47), one has to evaluate (50)The first step of the calculation is to show that in the WKB limit, the response matrix becomes diagonal. One should note that the previous expression (50)is similar to Eq. (37), where we discussed the biorthogonality of the WKB basis. In Eq. (50), the azimuthal Kronecker symbols ensure that . Moreover, because of our WKB assumptions from Eq. (41)on the step distances of the basis elements, the product of the two radial Gaussians in R_{g} imposes that in order to have a nonzero contribution. In order to shorten the notations, we temporarily introduce the function h(R_{g}) defined as which encompasses all the additional radial dependences appearing in Eq. (50). Thanks to the change of variables J_{φ} → R_{g}, the integral on J_{φ} which has to be evaluated in Eq. (50), when estimated for , is qualitatively of the form (51)This expression corresponds to a radial Fourier transform ℱ at the frequency Δk_{r}. It can be rewritten as a convolution of two radial Fourier transforms so that it becomes (52)where . We now rely on the WKB assumption from Eq. (40). If one has Δk_{r} ≠ 0, because of the Gaussian from Eq. (52), the contribution from ℱ [ h ] will come from the region k′ ~ Δk_{r} ≫ 1 /σ. We assume that the function h is such that its Fourier Transform is limited to the frequency region  k′  ≲ 1 /σ. This is consistent with assuming that the properties of the disc are radially slowly varying, and this implies that nonzero contributions to the response matrix can only be obtained when . Therefore, we have shown that within our WKB formalism, the response matrix from Eq. (5)is diagonal.
In order to shorten the notations, we will denote the matrix eigenvalues as (53)For these diagonal coefficients, the last step is to explicitly compute the integrals over J_{φ} and J_{r} in Eq. (50)to obtain the expression of the response matrix eigenvalues. We now detail this calculation. Thanks to our scaledecoupling approach, we may replace the radial Gaussian from Eq. (51)by a Dirac delta while paying a careful attention to the correct normalization of the Gaussian. Hence we have to evaluate (54)Because of the presence of the azimuthal Kronecker symbol, we may drop the sum on m_{φ}. The intrinsic frequencies from Eqs. (22)allow us to compute (55)Moreover, we assume that the galactic disc is tepid so that  ∂F/∂J_{φ}  ≪  ∂F/∂J_{r} . We may then only keep the term corresponding to a gradient with respect to the radial action J_{r}. Thanks to the expression of the Schwarzschild distribution function from Eq. (25)and the expression of the basis ampitude from Eq. (43), Eq. (54)becomes after some simple algebra (56)We may now use the following integration formula (see formula (6.615) from Gradshteyn & Ryzhik 2007) (57)where α> 0, β> 0, m_{r} ∈Z and ℐ_{mr} are the modified Bessel functions of the first kind. We apply this formula with and . We also introduce the notation (58)so that Eq. (56)becomes (59)We now define the dimensionless shifted frequency s as (60)Because we have ℐ_{− mr} [ χ ] = ℐ_{mr} [ χ ], we may rewrite Eq. (59)using the reduction factor (Kalnajs 1965; Lin & Shu 1966) defined as (61)As a conclusion, we obtain that within our WKB formalism the response matrix becomes diagonal and in the limit of tepid discs reads (62)This eigenvalue recovered using the WKB basis introduced in Eq. (26)is in full agreement with the seminal results from Kalnajs (1965) and Lin & Shu (1966). In order to handle the singularity of the eigenvalue appearing for s = n ∈Z, one adds a small imaginary part to the frequency of evaluation, so that s = n + iη. Indeed, as long as η is small compared to the imaginary part of the least damped mode of the disc, adding this complex part makes a negligible contribution to the expression of Re(λ).
3.4. Estimation of the susceptibility coefficients
One can now estimate the dressed susceptibility coefficients from Eq. (4). In order to shorten the notations, we will write the WKB basis elements introduced in Eq. (26)as (63)We have shown previously in Eq. (62)that within our WKB basis, the response matrix is diagonal. Its eigenvalues will be noted as λ_{p} so that we have . Hence the expression (4)of the susceptibility coefficients takes the form Using the expression of the Fourier transformed basis elements obtained in Eq. (49), we obtain (64)where we used the shortened notations κ_{i} = κ(J_{i}), and R_{i} = R_{g}(J_{i}). We also used the approximation introduced in Eq. (47)for the values at which the Bessel functions have to be evaluated. Thanks to the Kronecker symbols in m_{φ} and , we necessarily have (65)so that the sum on can be dropped.
3.5. Restriction on the loci of resonance
Before proceeding with the evaluation of the susceptibility coefficients obtained in Eq. (64), one may first emphasize a crucial consequence of the WKB basis from Eq. (26), which is the restriction to only exactly local resonances. Indeed, one can note that the expressions (8)and (9)of the drift and diffusion coefficients all involve an integration over the mute variable J_{2}. For a given value of J_{1}, m_{1} and m_{2}, this integration should be seen as a scan of the entire actionspace, searching for resonant regions where the constraint m_{1}·Ω_{1} − m_{2}·Ω_{2} = 0 is satisfied. We first recall the rule for the composition of a Dirac delta and a function which reads (66)where Z_{f} = { y  f(y) = 0 }, and we have supposed that all the poles of f are simple. As noted in Eq. (22), within the epicyclic approximation, the intrinsic frequencies Ω = (Ω_{φ},κ) only depend on R_{g} = R_{g}(J_{φ}) and are independent of J_{r}. Hence, the resonance condition m_{1}·Ω_{1} − m_{2}·Ω_{2} = 0 only depends on and is independent of . Hence if we consider fixed J_{1}, m_{1} and m_{2}, the resonant Dirac delta which has to be studied takes the form (67)where the resonance condition is given by (68)The radii therefore correspond to the resonant radii for which the resonance condition is satisfied. When writing Eq. (67), we assumed that the zeros of the resonance function are simple, which corresponds to the assumption that for any resonant radius , we have . Assuming that , this condition can be rewritten as (69)Resonance poles are therefore simple as long as the rates of change of the two intrinsic frequencies are not in a rational ratio. One can note that the Keplerian case for which κ = Ω_{φ} and the harmonic case for which κ = 2Ω_{φ} are in this sense degenerate. They can lead to resonant poles of higher multiplicity and would therefore require a more involved evaluation of the BalescuLenard collision operator. In what follows we assume that the potential is nondegenerate.
One can now use the properties of the WKB basis to restrict the range of resonant radii, . The expression (64)of the susceptibility coefficients, thanks to the two Gaussians, imposes that the relevant resonant radius must necessarily be close to R_{1}. As noted in Eq. (65), in order to have a nonzero susceptibility, one also has to satisfy the constraint . The resonant condition which has to be satisfied is therefore given by (70)where the distance is such that  ΔR  ≤ (few) σ. Because the scaledecoupling parameter σ is supposed to be small compared to the size of the system, we may approximate the previous resonant condition as (71)On the l.h.s. of Eq. (71), the term within bracket is nonzero because we assumed in Eq. (69)that the resonant poles are simple. Moreover, ΔR is small because of our scaledecoupling approach. The r.h.s. of Eq. (71)is discrete: it is either zero or at least of the order of κ(R_{1}). Because the l.h.s. is necessarily small, we must have (72)This result is a crucial consequence of our WKB tightly wound spiral assumption because it implies that only local resonances are allowed. In particular this implies that the WKB limit does not allow for distant orbits to resonate (through e.g. propagation of swing amplified wave packets, see below). Then the sum from Eq. (67)can be limited to the evaluation in . Hence within this WKB limit, the susceptibility coefficients from Eq. (64)have to be evaluated only for m_{2} = m_{1} and R_{2} = R_{1}, so that we have to deal with the expression (73)
3.6. Asymptotic continuous limit
One can note that in Eq. (73)the susceptibility coefficients are still expressed as a discrete sum on the basis index and . Our next step is to replace these sums by continuous integrals. The discrete basis elements are separated by the step distances ΔR_{0} and Δk_{r}, which must satisfy the WKB hypothesis detailed in Eq. (41). We use the Riemann sum formula , where Δx controls the distance between the basis elements. This transformation is a subtle stage of the calculation because one has to consider step distances ΔR_{0} and Δk_{r}, which have to simultaneously be large to comply with the WKB assumption from Eq. (41)and small to allow the use of the Riemann sum formula. As we are going to transform both the sums on and , the exact value of the susceptibility coefficients will depend on our choice for ΔR_{0} Δk_{r}. One has to consider the case (74)This sampling corresponds to a critical sampling condition (Gabor 1946; Daubechies 1990; see also Fouvry et al. 2015b). Equation (73)then takes the form (75)One can now assume that the radial Gaussian present in Eq. (75)is sufficiently peaked. Because it is correctly normalized, we may in this limit replace it by δ_{D}(R_{1} − R_{0}). The integration on R_{0} can then be immediately performed to give (76)where only λ_{kr} depends on the frequency of evaluation ω. One may note that in Eq. (76), all the dependencies in σ have disappeared, so that the value of the susceptibility coefficients is independent of the precise choice of the WKB basis. The square of the susceptibility coefficients which is required to estimate the drift and diffusion coefficients from Eqs. (8)and (9)is therefore given by (77)
where we introduced a cutoff at 1 /σ_{k} for the integration on k_{r}. This bound is justified by the WKB constraint from Eq. (41), which imposes that the probed radial frequency region is bounded from below. It is also important to note that these susceptibility coefficients should be evaluated at R_{2} = R_{1}, since we proved in Eq. (72)that, consistently with our WKB approximation, exactly local resonances are the only ones which have to be considered.
At this stage, there is an arbitration to make between two possible behaviours depending on the physical properties of the underlying disc. First of all, if the amplification function k_{r} → 1 / (1 − λ_{kr}) is asymptotically a sharp function reaching a maximum value λ_{max} for k_{r} = k_{max}, one can assume that the susceptibility coefficients are dominated by the contribution from the peak in λ_{kr}. In this situation, we can perform an approximation of the small denominators. The second possible behaviour arises if the function k_{r} → 1 / (1 − λ_{kr}) is asymptotically flat, so that there is no characteristic scale of blowup of the amplification eigenvalues. In such a situation, the susceptibility coefficients are mostly dominated by the behaviour at the boundaries of integration from Eq. (77)where k_{r} → 1 /σ_{k}. The detailed response structure of the selfgravitating disc then does not play a significant role.
We place ourselves within the approximation of the small denominators, assuming that the biggest contribution to the susceptibility coefficients comes from waves which yield the largest λ_{kr}. Therefore, one has to suppose that the function k_{r} → 1 / (1 − λ_{kr}) is a sharp function reaching a maximum value λ_{max}(R_{1},ω), for k_{r} = k_{max}(R_{1},ω), with a characteristic spread Δk_{λ}(R_{1},ω). The expression (77)then becomes (78)While still focusing on the contribution to the susceptibility coefficients due to the waves with the largest amplification λ(k_{r}), one can improve the approximation of the small denominators from Eq. (78). Indeed, starting from Eq. (77), one can instead perform the k_{r}integration for k_{r} ∈ [ k_{inf}; k_{sup} ], where these bounds are given by λ(k_{inf / sup}) = λ_{max}/ 2. This approach is numerically more demanding but does not alter the principal conclusions drawn in this paper, while allowing a more precise determination of the secular flux structure. All the calculations presented in Sect. 4 were performed with this improved approximation. Finally, in Appendix B, we detail how this same WKB formalism may be applied to the inhomogeneous BalescuLenard equation without collective effects (Chavanis 2013b).
3.7. Estimation of the drift and diffusion coefficients
The drift and diffusion coefficients are given by Eqs. (8)and (9). Within the WKB approximation, we have shown in Eqs. (65)and (72)that the susceptibility coefficients have to be evaluated only for m_{1} = m_{2}, so that the sum on m_{2} in the expressions of the drift and diffusion coefficients may be dropped. As the resonances are exactly local, using the formula from Eq. (67)adds a prefactor of the form 1 /  ∂(m_{1}·Ω) /∂J_{φ} , so that the drift coefficients from Eq. (8)become (79)Similarly, the diffusion coefficients are given by (80)In both Eqs. (79)and (80), the susceptibility coefficients are given by Eq. (77)or, within the approximation of the small denominators, by Eq. (78).
Such simple expressions of the drift and diffusion coefficients along with the required expression of the susceptibility coefficients constitute a main result of this paper. One should also note that this WKB formalism is self contained and is obtained without any ad hoc assumptions or fittings. Except for the explicit recovery of the amplification eigenvalues from Eq. (62), the calculations presented previously are not limited to the Schwarzschild distribution function from Eq. (25). Indeed, the drift and diffusion coefficients from Eqs. (79)and (80)are valid for any tepid disc, as long as the epicyclic anglesactions mapping from Eq. (24)may be used. In Appendix C, we show how the drift and diffusion coefficients from Eqs. (79)and (80)can be explicitly computed, when one considers a Schwarzschild distribution function as in Eq. (25)and when the susceptibility coefficients are estimated thanks to the approximation of the small denominators from Eq. (78). Finally, in Appendix D, we compare this 2D WKB BalescuLenard equation with other similar kinetic equations.
4. Application
We may now illustrate how this WKB approximation of the inhomogeneous BalescuLenard equation can be implemented to recover some results obtained via wellcrafted numerical simulations. Indeed, Sellwood (2012; hereafter S12) using careful numerical simulations studied the longterm evolution of an isolated stable and stationary Mestel disc (Mestel 1963) sampled by pointwise particles. When evolved for hundreds of dynamical times, such a disc would secularly diffuse in actionspace through the spontaneous generation of transient spiral structures. Figure 7 of S12 shows for instance the late time formation of resonant ridges along very specific resonant directions. Such features are possible signatures of secular evolution, which results in a longterm aperiodic evolution of a selfgravitating system, during which small resonant and cumulative effects can add up in a coherent way. These small effects, which are then amplified through the selfgravity of the system originate from finiteN effects. Indeed, the distribution function of the system is made of a finite number N of pointwise particles. Even with a perfect numerical integrator, the system would necessarily undergo encounters during which orbits feel the discreteness of the joint DF through its two point correlation. One should remark that these interactions need not be local but assume that the potential fluctuations are resonant so as to build up a secular evolution of the system. This effect, which is still present even in the absence of any numerical noise, is the effect captured by the BalescuLenard equation (see Fig. 1).
4.1. The disc model
The disc considered by S12 is an infinitely thin Mestel disc, for which the circular speed v_{φ} is a constant V_{0} independent of the radius. Such a model represents fairly well the observed rotation curve of real galaxies. The stationary background potential ψ_{M} of such a disc and its associated surface density Σ_{M} are given by (81)where R_{i} is a scale parameter. Because of this scale invariance, the relationship between the angular momentum J_{φ} and the guiding radius R_{g} is straightforwardly given by (82)Within the epicyclic approximation, the intrinsic frequencies Ω_{φ} and κ can be computed thanks to Eqs. (22)and read (83)We note that , so that the Mestel disc could be seen as an intermediate case between the Keplerian case for which κ/ Ω_{φ} = 1 and the harmonic case for which κ/ Ω_{φ} = 2. The ratio of the intrinsic frequencies is a important parameter for the system since it will determine the location of the resonances and a constant ratio may introduce dynamical degeneracies. This is the case for the Keplerian and harmonic discs for which κ/ Ω_{φ} is a rational number, as discussed below Eq. (69). By contrast, for the Mestel disc, the nonrational ratio ensures that the potential is nondegenerate. Using the epicyclic approximation, the DF considered by S12 takes, as in Eq. (25), the form of a Schwarzschild DF, where the intrinsic frequencies are given by Eq. (83), the velocity dispersion σ_{r} is constant throughout the entire disc, and the surface density is given by Σ_{t}, i.e. the active surface density of the disc. Indeed, in order to accommodate the central singularity and the infinite extent of the Mestel disc, one introduces tapering functions T_{inner} and T_{outer} to damp out the inner and outer regions, which read (84)where ν_{t} and μ_{t} control the sharpness of the two tapers and R_{0} is an additional scale parameter. The two tapers are physically motivated by the presence of a bulge and an outer truncation for the disc. Moreover, in order to reduce the susceptibility of the disc, we also suppose that only a fraction ξ of the disc is selfgravitating, with 0 ≤ ξ ≤ 1, so that the rest of the gravitational field is provided by the static halo. As a conclusion, the active surface density Σ_{t} is given by (85)where Σ_{M} is the surface density of the Mestel disc from Eq. (81). 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 σ_{r} = 0.284, ν_{t} = 4, μ_{t} = 5, ξ = 0.5 and R_{0} = 11.5. The shape of the active surface density is illustrated in Fig. 4. The initial contours of the Schwarzschild DF from Eq. (25)are shown in Fig. 5. One can also straightforwardly estimate the total active mass of such a disc, given by M_{tot} ≃ 5.4. For such an almost scale invariant disc, the local Toomre parameter, Q (Toomre 1964) (86)which for Q> 1 ensures the stability of the disc with respect to local axisymmetric disturbances, becomes almost independent of the radius, especially in the intermediate regions of the disc. As illustrated in Fig. 6, Q ≃ 1.5 between the tapers and increases strongly in the tapered regions.
Fig. 4 Surface density Σ_{t} of the tapered Mestel disc. The unit system has been chosen so that V_{0} = G = R_{i} = 1. Because of the tapering functions, the selfgravity of the disc is turned off in the inner and outer regions. 

Open with DEXTER 
Fig. 5 Contours of the initial distribution function in actionspace (J_{φ},J_{r}), within the epicyclic approximation. The contours are spaced linearly between 95% and 5% of the distribution function maximum. 

Open with DEXTER 
The expression (10)of the secular diffusion flux requires summing on all the resonances m. S12 restricted pertubations forces to m_{φ} = 2, so that we may impose this same restriction on the considered azimuthal number m_{φ}. Throughout our numerical calculations, we will restrict ourselves to only three different resonances which are: the inner Lindblad resonance (ILR) corresponding to , the outer Lindblad resonance (OLR) given by and the corotation resonance (COR) for which . Moreover, all the calculations in the upcoming sections have also been performed while taking into account the contributions from the resonances with m_{r} = ± 2, which were checked to be subdominant. Being able to perform such a restriction to the relevant resonances appearing in the secular flux ℱ_{tot} from Eq. (10)is an important step of the calculation.
Returning to the fast and slow coordinates from Eq. (14), one can note that the diffusion associated with the COR resonance amounts to diffusion along the J_{φ}axis. Such diffusion brings stars from one quasicircular orbit to another of a different radius and is called radial migration. Conversely, the diffusion associated with the ILR and OLR resonances exhibits a nonzero diffusion component in the J_{r}direction. It therefore increases the velocity dispersion within the disc so as to heat it, while either decreasing (ILR) or increasing (OLR) star’s angular momentum.
4.2. Disc amplification
One may now study the behaviour of the amplification eigenvalues λ_{kr} given by Eq. (62), thanks to which one can perform the improved approximation of the small denominators. For a given resonance m and angular momentum J_{φ}, the amplification function k_{r} → λ_{kr} is presented in Fig. 7.
Fig. 6 Dependence of the local Q Toomre parameter with the angular momentum. It is scale invariant except in the inner/outer regions because of the presence of the tapering functions T_{inner} and T_{outer}. The unit system has been chosen so that V_{0} = G = R_{i} = 1. 

Open with DEXTER 
Fig. 7 Variations of the response matrix eigenvalues λ with the WKBfrequency k_{r}, for m = m_{COR} and two values of J_{φ}. The curve that peaks at large k_{r} is for the smaller value of J_{φ}. 

Open with DEXTER 
As Eq. (62)only depend on s^{2}, the ILR and OLR resonances always have the same response matrix eigenvalues. One can also note that the eigenvalues λ(k_{r}) are maximum for a frequency k_{max}(J_{φ}), where λ(k_{max}) = λ_{max}, in a region whose size is given by the width at half maximum Δk_{λ}. Because of the scaleinvariance property of the Mestel disc, it is straighforward to show that Δk_{λ} ∝ 1 /J_{φ}, k_{max} ∝ 1 /J_{φ} and k_{inf / sup} ∝ 1 /J_{φ}. One can then consider the behaviour of the amplification factor 1 / (1 − λ_{max}), which encodes the strength of the selfgravitating amplification, as shown in Fig. 8.
Fig. 8 Dependence of the amplification factor 1 / (1 − λ_{max}) with the position J_{φ} in the disc. The amplification associated with the COR is always larger than the one associated with the ILR and OLR. 

Open with DEXTER 
We note that the COR resonance is always more amplified than the ILR and OLR resonances, but the maximum amplification (~3 for the COR and ~1.5 for the ILR and OLR) remains sufficiently small, so that the susceptibility coefficients from Eq. (4)are not dominated only by the selfgravitating amplification.
Fig. 9 Map of N div(ℱ_{tot}), where the total flux ℱ_{tot} has been summed over the three resonances (ILR, COR and OLR). Red contours, for which N div(ℱ_{tot}) < 0, correspond to regions from which the orbits will be depleted, whereas blue contours, for which N div(ℱ_{tot}) > 0, correspond to regions where secular diffusion will tend to increase the value of the DF. The net fluxes involve simultaneously radial migration near (J_{φ},J_{r}) ~ (1.8,0), and heating near (J_{φ},J_{r}) ~ (1,0.1). 

Open with DEXTER 
4.3. Computing the diffusion flux
Given the knowledge of the eigenvalues, it is now straightforward to compute the susceptibility coefficients within the improved approximation of the small denominators thanks to Eq. (77), where the integration on k_{r} is performed for k_{r} ∈ [ k_{inf} ; k_{sup} ]. One can then compute the associated drift and diffusion coefficients respectively given by Eqs. (79)and (80). The diffusion flux ℱ_{tot} defined in Eq. (10)immediately follows, where the sum on m is restricted only to the three resonances ILR, COR and OLR. In Appendices E and F, we discuss two specific properties of such a truncated Mestel disc, namely the cancellation between the radial components of the diffusion and drift elements (the Schwarzshild conspiracy, Appendix E) and the natural and intrinsic presence of a temporal frequency bias (Appendix F), both of which enlighten the subtle arbitrations between the different resonances.
Fig. 10 Overlay of the WKB predictions for the divergence of the diffusion flux N div(ℱ_{tot}) and the differences between the initial and the evolved DF in S12’s simulation. The opaque contours correspond to the differences in the actionspace for the DF in S12 between the time t_{S12} = 1000 and t_{S12} = 0 (see the upper panel of S12’s Fig. 10). The red opaque contours correspond to negative differences, so that these regions are emptied from their orbits, whereas blue opaque contours correspond to positive differences, i.e. regions where the DF has increased through diffusion. The transparent contours correspond to the predicted values of N div(ℱ_{tot}) using the same conventions as in Fig. 9. We note the overlap of the predicted transparent red and blue contours with the measured solid ones. 

Open with DEXTER 
Fig. 11 Overlay of the WKB predictions for the divergence of the diffusion flux N div(ℱ_{tot}) on top of the contours of the DF in actionspace measured in the S12 simulation. The black background contours are the levels contours of the DF at time t_{S12} = 1400 (see the lower panel of Fig. 7 of S12). These contours are spaced linearly between 95% and 5% of the DF maximum and exhibit clearly the appearance of a resonant ridge. The coloured transparent contours correspond to the predicted values of N div(ℱ_{tot}) using the same conventions as in Fig. 9. We note that the developed late time ridge is consistent with the predicted depletion (red) and enrichment (blue) of orbits. 

Open with DEXTER 
Finally, we may compute the divergence of this flux, div(ℱ_{tot}) given by Eq. (11)in order to compare quantitatively the WKB predictions with the results from S12’s simulations. Because the mass of the individual particles is given by μ = M_{tot}/N, we will consider the quantity N div(ℱ_{tot}), which is independent of N. Figure 9 represents the contours of N div(ℱ_{tot}). Comparisons of our WKB predictions with the results from S12’s simulation are illustrated in Figs. 10 and 11. In Fig. 9, red contours correspond to regions for which N div(ℱ_{tot}) < 0, so that thanks to Eq. (11)they are associated with actionspace regions where the WKB BalescuLenard equation predicts a decrement of the DF during secular evolution. In contrast, blue contours are associated with regions for which N div(ℱ_{tot}) > 0, so that the DF will increase there. The overall picture involves two competing processes: i) the beginning of a ridge forming towards (J_{φ},J_{r}) ~ (1,0.1); and ii) the formation of an over density near (J_{φ},J_{r}) ~ (1.8,0). Point i) is in fact consistent with both the early time measurement of S12 as shown in Fig. 10 and the late time measurement of S12 as shown in Fig. 11. These qualitative agreements are in fact surprisingly good, given that the WKB theory is approximate and was only estimated for t = 0^{+}. Interestingly, the early time measurement from Fig. 10 also displays a hint of an overdensity on the J_{r} = 0 axis, in agreement with point ii), while the late time measurement suggests that the over density has split, with a hint of a second ridge forming. The time evolution of Eq. (2) is likely to explain why this over density on the axis seems to split, and why the ridge gets amplified.
In Fig. 9, we explicitly computed N div(ℱ_{tot}), so that we may now study the typical timescale of collisional relaxation predicted by the WKB BalescuLenard equation. This is the purpose of the next section.
4.4. Physical timescales
Given estimates of the diffusion flux, one can explicitly compute the collisional timescale, i.e. the timescale for which the finiteN effects become significant. Indeed, the larger the number N of particles, the later these effects will come into play. One should remark that our writing of the BalescuLenard Eq. (2)depends on N only through the mass of the individual particles given by μ = M_{tot}/N. Indeed, one can rewrite this kinetic equation under the form (87)where L is the operator of pure advection, and C_{BL} = N div(ℱ_{tot}) is the Nindependent BalescuLenard collisional operator, i.e. the r.h.s. of Eq. (2)multiplied by N = M_{tot}/μ. Equation (87)underlines the fact that the collisional term is associated with a kinetic Taylor expansion in the parameter ε = 1 /N ≪ 1. Within the angleaction coordinates, the advection operator is immediately given by (88)Because we have assumed that F is always quasistationary, so that F = F(J,t) (adiabatic approximation), one has L [ F ] = 0. We now introduce the rescaled time (89)so that Eq. (87)immediately becomes (90)Equation (90)corresponds to a rewriting of the BalescuLenard equation, where N is not present anymore. This will allow us to quantitatively compare the time during which the S12 simulation was run to the diffusion time predicted by our WKB BalescuLenard formalism. In order to ease this comparison, we place ourselves in the same units system as the one used by S12. Figure 7 of S12 for which the ridge was observed was obtained with the parameters N = 50 × 10^{6} and Δt_{S12} = 1400. Using the rescaled time introduced in Eq. (89), one obtains that S12 observed the resonant ridge after a time Δτ_{S12} = Δt_{S12}/N ≃ 3 × 10^{5}. One can then compare this time with the typical time required to obtain a resonant ridge within our WKB formalism. Given the map of N div(ℱ_{tot}) described in Sect. 4.3, one can estimate the typical time for which this flux could lead to the features observed in S12. The contours presented in the Fig. 7 of S12 are separated by a value of , where corresponds to the maximum of the normalized DF from Eq. (25). As a consequence, to observe the resonant ridge, the DF should typically change by a value of the order of . From Fig. 9, one can note that the maximum value of the divergence of the flux is given by  N div(ℱ_{tot})  ≃ 0.4. Finally, thanks to Eq. (11), one may write the relation ΔF_{0} ≃ Δτ_{WKB}  N div(ℱ_{tot}) , where Δτ_{WKB} is the minimal time during which the WKB BalescuLenard equation has to be considered in order to develop a ridge. Thanks to the previous typical numerical values, one obtains that Δτ_{WKB} ≃ 3 × 10^{2}. When comparing these two typical times, Δτ_{S12} the duration during which S12 simulation was performed and Δτ_{WKB} the duration required to observe secular diffusion in the WKB BalescuLenard formalism, one obtains the order of magnitude (91)Hence the direct application of the WKBlimited BalescuLenard equation does not allow us to predict the observed timescale for the diffusion features in simulations. Indeed, the timescale of collisional diffusion predicted by this WKB formalism seems much larger than the time during which the numerical simulation was performed. This discrepancy is also strengthened by the use of a softening length in numerical simulations, which induces an effective thickening of the disc, so as to slow down the collisional relaxation. A possible explanation for this timescale discrepancy is discussed in the next section.
4.5. Interpretation
In order to interpret S12’s simulation under the light of a collisional secular diffusion equation, such as the BalescuLenard equation and its WKB limit, one should first note the undisputed presence of collisional effects in S12’s simulation. Indeed, Fig. 2 of S12 shows that when the number of particles of the simulation is increased, the strength of the density fluctuations is delayed, which in turn is likely to be related to the amplitude of the secular features. The larger the number of particles, the later the effect of secular diffusion. Such dependence illustrates the fact that discreteness effects do play a role in the secular diffusion observed in S12.
Sellwood & Kahn (1991) have argued that a sequence of causally connected swing amplified transients could occur subject to a (possibly nonlocal) resonant condition between successive spiral waves. The BalescuLenard formalism captures precisely such sequences – in as much as it integrates over dressed correlated potential fluctuations subject to relative resonant conditions, but does not preserve causality nor resolve them on dynamical timescales. The exact initial phases are not relevant in the BalescuLenard formalism: see Appendix A for a sketch of a full derivation which makes this point clear.
The timescale discrepancy observed in Eq. (91)might be driven by the incompleteness of the WKB basis. Indeed, Eq. (26)’s basis – thanks to which the susceptibility coefficients were evaluated in Eq. (77)– does not form a complete set, as it can only represent correctly tightly wound spirals. It also enforces local resonances, and does not allow for remote orbits to resonate, or wave packets to propagate between such nonlocal resonances. The seminal works from Goldreich & LyndenBell (1965), Julian & Toomre (1966), Toomre (1981) showed that any leading spiral wave during its unwinding to a trailing wave undergoes a significant amplification, coined swing amplification. Because it involves open spirals this mechanism is not captured by the WKB formalism. This additional amplification is expected to increase the susceptibility of the disc and therefore accelerate secular diffusion (both drift and diffusion), so that the timescales discrepancy from Eq. (91)will become less restrictive. Following the notations from Toomre (1981), the truncated Mestel disc considered in the S12 simulation corresponds to Q ≃ 1.5 and X = 2, so that Fig. 7 from Toomre (1981) shows that significant swing amplification (of order 10) may be expected. It has also been claimed (Toomre & Kalnajs 1991) that swing amplified shot noise in the shearing sheet approximation would behave like significantly heavier macroparticles. Such an amplification would keep a dependence of the secular response with the total number N of particles, but would reduce significantly the effective number of particles.
The BalescuLenard WKB limit seems to capture qualitatively the main features of the initial diffusion process in action space (as discussed in Sect. 4.3), but falls short in predicting the relevant timescale. The remaining questions are therefore: what is the exact impact of swing amplification? Can it explain the timescale discrepancy?
5. Conclusion
We implemented the inhomogeneous BalescuLenard Eq. (2)for an infinitely thin galactic disc using two main approximations. We first assumed the disc to be tepid. We could then use the epicyclic approximation which allowed for an explicit mapping between the physical coordinates (x,v) and the angleaction coordinates (θ,J) via Eq. (24). Our second approximation relied on the introduction of the tightly wound basis elements from Eq. (26). Because of the corresponding WKB approximation, we obtained in Eq. (62)a diagonal response matrix, so that gravity is effectively treated locally. The associated scaledecoupling hypothesis yields a crucial restriction to only local resonances, as shown in Eq. (72). We then derived in Eq. (77)a simple quadrature for the susceptibility coefficients, given by Eqs. (B.7)and (B.8)for the bare ones. Thanks to this restriction to local resonances, we were also able to write the drift and diffusion coefficients as simple quadratures in Eqs. (79)and (80).
These simple expressions derived within the WKB formalism yield, to our knowledge, a first nontrivial explicit expression for the BalescuLenard diffusion and drift coefficients. They are certainly useful to provide insight into the physical processes at work during the secular diffusion of a selfgravitating discrete disc. Moreover, modulo the restriction to the three physically motivated resonances ILR, COR and OLR, our WKB formalism can be used for quantitative comparisons to numerical experiments such as the one presented in Sellwood (2012). It considered a stable isolated Mestel disc sampled by pointwise particles, whose secular evolution is induced by finiteN effects ideally captured by the BalescuLenard equation.
The straightforward calculation in the WKB limit of the divergence of the full diffusion flux, N div(ℱ_{tot}), illustrated in Figs. 9–11, recovered most of the secular features observed in S12. This qualitative agreement is impressive given the level of approximation involved in the WKB limit. The hints for the formation of a ridge – depletion and enrichment of orbits along a preferred direction – is qualitatively consistent with the findings of S12 and Fouvry & Pichon (2015), Fouvry et al. (2015a), without postulating additional assumptions about the source of fluctuations^{5}.
The comparison of the collisional time predicted in the WKB limit (Eq. (91)) to the diffusion time of S12 simulation, highlights nonetheless a significant quantitative overestimation. We provided a possible explanation which relies on the intrinsic limitations of the WKB formalism, as it cannot account for swing amplification, during which unwinding transient spirals are strongly amplified. This additional amplification, which involves explicitly nonlocal wave absorption and emission, may be the missing contribution required to reconcile quantitatively our predictions and the simulation. One venue will be to compute numerically exactly Eqs. (8)and (9)in action space – without assuming tightly wound spirals or epicyclic orbits – with a complete basis^{6}. This is the topic of an upcoming numerical investigation in Fouvry et al. (2015c).
Should this complementary investigation explain the timescale mismatch, one would be in a stronger position to validate the accuracy of Nbody schemes to correctly capture secular evolution of discrete selfgravitating cold discs over very long timescales. This would clearly be a worthy assessment of such schemes relying on the BalescuLenard theory. Once the above described conundrum is resolved, we will also be able to evolve over secular times the BalescuLenard equation and predict the full cosmic time evolution of such discrete discs. This may also contribute to solving the timescale discrepancy.
In closing, beyond the application presented in Sect. 4, the above developed tightly wound BalescuLenard formalism may for instance describe the secular diffusion of giant molecular clouds in galactic discs (which in turn could play a role in migration driven metallicity gradients and disc thickening), the secular migration of planetesimals in partially selfgravitating protoplanetary discs, or even the longterm evolution of population of stars, gas blobs and debris near the Galactic centre. Such topics will be subject to further investigations.
For a historical account of the development of kinetic theories in astrophysics, plasma physics, and for systems with longrange interactions, see the Introduction in Chavanis (2013a,b).
In this paper, we are not interested in the initial complex mechanism of violent relaxation (LyndenBell 1967), during which the system gets virialized, since we intend to describe the longterm evolution of an already and continuously virialized system.
In the secular timescale limit, the amplification through the propagation of waves between resonances is assumed to be instantaneous, see Appendix A.
In contrast, the formalism of secular forcing presented in Fouvry & Pichon (2015), Fouvry et al. (2015a) postulated a partially ad hoc shape of the perturbation power spectrum, see Eq. (F.2).
An alternative middle ground would be to account for nonlocal interferences of the WKB wave packets, given by Eq. (26), which in turn would allow for nonlocal resonances to come into play.
This also explains why the factorization assumption A_{m1,m2}(J_{1},J_{2}) = A_{m1}(J_{1}) A_{m2}(J_{2}) used in Chavanis (2007) does not hold.
Acknowledgments
J.B.F. thanks the Institute of Astronomy, Cambridge, for hospitality while this investigation was completed. J.B.F. and C.P. also thank the theoretical physics subdepartment, Oxford, for hospitality and the CNRSOxford exchange program for funding. We thank Donald LyndenBell, James Binney, Simon Prunet, Walter Dehnen, John Magorrian and Mir Abbas Jalali for their feedback. 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
 Balescu, R. 1960, Phys. Fluids, 3, 52 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Binney, J., & Tremaine, S. 2008, Galactic Dynamics: (second Edition), Princeton Series in Astrophysics (Princeton University Press) [Google Scholar]
 Born, M. 1960, The Mechanics of the Atom (F. Ungar Pub. Co.) [Google Scholar]
 Chandrasekhar, S. 1942, Principles of Stellar Dynamics (University of Chicago Press) [Google Scholar]
 Chandrasekhar, S. 1949, Rev. Mod. Phys., 21, 383 [NASA ADS] [CrossRef] [Google Scholar]
 Chavanis, P.H. 2007, Physica A, 377, 469 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Chavanis, P.H. 2012a, Physica A, 391, 3680 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Chavanis, P.H. 2012b, Physica A, 391, 3657 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Chavanis, P.H. 2012c, Eur. Phys. J. Plus, 127, 19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chavanis, P.H. 2013a, Eur. Phys. J. Plus, 128, 126 [CrossRef] [EDP Sciences] [Google Scholar]
 Chavanis, P.H. 2013b, A&A, 556, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chavanis, P.H., & Lemou, M. 2007, Eur. Phys. J. B, 59, 217 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
 Daubechies, I. 1990, Information Theory, IEEE Transactions on, 36, 961 [Google Scholar]
 Earn, D. J. D., & LyndenBell, D. 1996, MNRAS, 278, 395 [NASA ADS] [CrossRef] [Google Scholar]
 Fouvry, J.B., & Pichon, C. 2015, MNRAS, 449, 1982 [NASA ADS] [CrossRef] [Google Scholar]
 Fouvry, J.B., Binney, J., & Pichon, C. 2015a, ApJ, 806, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Fouvry, J.B., Pichon, C., & Prunet, S. 2015b, MNRAS, 449, 1967 [NASA ADS] [CrossRef] [Google Scholar]
 Fouvry, J. B., Pichon, C., Magorrian, J., & Chavanis, P.H. 2015c, A&A, accepted [arXiv:1507.06887] [Google Scholar]
 Gabor, D. 1946, Electrical Engineers – Part III: Radio and Communication Engineering, 93, 429 [Google Scholar]
 Goldreich, P., & LyndenBell, D. 1965, MNRAS, 130, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Goldstein, H. 1950, Classical mechanics (AddisonWesley) [Google Scholar]
 Gradshteyn, I. S., & Ryzhik, I. M. 2007, Table of integrals, series, and products (Amsterdam: Elsevier/Academic Press) [Google Scholar]
 Heyvaerts, J. 2010, MNRAS, 407, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Jeans, J. 1915, MNRAS, 76, 70 [NASA ADS] [CrossRef] [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. 1965, Ph.D. Thesis (Harvard University) [Google Scholar]
 Kalnajs, A. J. 1976, ApJ, 205, 745 [NASA ADS] [CrossRef] [Google Scholar]
 Klimontovich, I. 1967, The statistical theory of nonequilibrium processes in a plasma (M.I.T. Press) [Google Scholar]
 Landau, L. 1936, Phys. Z. Sowj. Union, 10, 154 [Google Scholar]
 Lenard, A. 1960, Ann. Phys., 10, 390 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, C. C., & Shu, F. H. 1966, Proc. National Academy of Science, 55, 229 [NASA ADS] [CrossRef] [Google Scholar]
 Liouville, J. 1837, J. Math. Pures Appl., 1, 16 [Google Scholar]
 LyndenBell, D. 1967, MNRAS, 136, 101 [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]
 Mestel, L. 1963, MNRAS, 126, 553 [NASA ADS] [CrossRef] [Google Scholar]
 Palmer, P. 1994, Stability of Collisionless Stellar Systems: Mechanisms for the Dynamical Structure of Galaxies, Astrophys. Space Sci. Lib. (Netherlands: Springer) [Google Scholar]
 Palmer, P. L., Papaloizou, J., & Allen, A. J. 1989, MNRAS, 238, 1281 [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]
 Rosenbluth, M., MacDonald, W., & Judd, D. 1957, Phys. Rev., 107, 1 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Sellwood, J. A. 2012, ApJ, 751, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Sellwood, J. A., & Kahn, F. D. 1991, MNRAS, 250, 278 [NASA ADS] [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]
 Toomre, A., & Kalnajs, A. J. 1991, in Dynamics of Disc Galaxies, ed. B. Sundelius, 341 [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]
Online material
Appendix A: Sketch of BalescuLenard derivation
Two derivations of the inhomogeneous BalescuLenard equation have been presented in the literature. The first one (Heyvaerts 2010) is based on the appropriate truncation at the order 1 /N of the BBGKY hierarchy. The second (Chavanis 2012a) relies on the Klimontovich equation, using a quasilinear approximation. We now briefly sketch the derivation presented in Chavanis (2012a). We consider an isolated system of N particles in interaction, of mass μ, in a physical space of dimension d. Their dynamics is entirely described by Hamilton’s equations (A.1)where the Hamiltonian of the system is given by (A.2)Here u(  x_{i} − x_{j} ) is the binary potential of interaction. In the gravitational case, it satisfies u(x) = − G/  x . One can now introduce the discrete distribution function F_{d}(x,v,t) defined as (A.3)along with the corresponding potential (A.4)One can show that F_{d} satisfies the Klimontovich equation (Klimontovich 1967) (A.5)where we have defined the Hamiltonian H_{d} as (A.6)At this stage, it is important to note that the Klimontovich Eq. (A.5)contains exactly the same information as Hamilton’s Eqs. (A.1). We now introduce the smooth distribution function F(x,v,t) = ⟨ F_{d}(x,v,t) ⟩, corresponding to an average of F_{d} over a large number of initial conditions. One can then write F_{d} = F + δF, where δF denotes fluctuations around the smooth distribution. In a similar way, we introduce ψ(x,v,t) = ⟨ ψ_{d}(x,v,t) ⟩, so that ψ_{d} = ψ + δψ. We have therefore decomposed the discrete distribution function F_{d} into a smooth component F that evolves slowly with time, whereas the fluctuating component δF evolves more rapidly. As a consequence, when considering the evolution of the fluctuations, one can assume the smooth distribution to be frozen. Using this timescaledecoupling approach, one can use the angleaction coordinates (θ_{1},J_{1}) associated with the quasistationary smooth potential ψ to describe the fast evolution of the fluctuations. Using these decompositions and this change of coordinates, Eq. (A.5)takes the form of two evolution equations (A.7)and (A.8)where we have introduced the intrinsic frequencies of the system Ω_{1} = Ω(J_{1}), as in Eq. (1), and where ⟨ . ⟩ denotes an angle average. Because of our timescaledecoupling approach, we may neglect the time variation of F(J_{1},t) in the calculation of the collision term (adiabatic approximation). It implies that F may be treated as a constant in Eq. (A.7)because F evolves on a (relaxation) timescale much larger than the (dynamical) time corresponding to the evolution of δF. In order to be valid, this approximation requires that N ≫ 1. Finally, we also assume that the distribution F remains Vlasov stable, so that its evolution is only governed by correlations and not by dynamical instabilities.
The first step of the derivation of the BalescuLenard equation is then to study the short timescale evolution Eq. (A.7). One defines the FourierLaplace transform of the fluctuation δF as (A.9)valid for Im(ω_{1}) sufficiently large. Similarly to Eq. (6), one also defines the spatial Fourier tranform of the initial value as (A.10)Thanks to these transformations, Eq. (A.7)may be rewritten under the form (A.11)We now use the basis elements introduced in Eq. (3), so that we may decompose the potential fluctuations under the form (A.12)We introduce the Laplace transform of a_{p}(t) as (A.13)One can then take the inverse Fourier transform of Eq. (A.11), multiply by and integrate over θ_{1} and J_{1} (using the property that dxdv = dθ_{1}dJ_{1}). One gets (A.14)where the response matrix is given by Eq. (5). Equation (A.14)can be rewritten under the form (A.15)where the susceptibility coefficients have been introduced in Eq. (4). One can now compute the collision term appearing in the r.h.s. of Eq. (A.8). It requires us to evaluate (A.16)Using Eq. (A.11), one immediately obtains that (A.17)In Eq. (A.17), the first term corresponds to the selfcorrelation of the potential whereas the second term corresponds to the correlations between the potential fluctuations and the distribution function at time t = 0. Each of these terms must then be considered one at a time. Assuming that there is no correlation in the initial phases, one can show (see Appendix C from Chavanis 2012a) that (A.18)Hence, given Eq. (A.15), one can rewrite the first term of Eq. (A.17)under the form (A.19)If we consider only the contributions that do not decay in time, one can perform the substitution
Starting from Eq. (A.19), thanks to the previous substitution and using the fact that , one can show that the first contribution from Eq. (A.16)takes the form (A.20)The last step of the calculation is to use the Landau prescription ω → ω + i0^{+} along with Plemelj formula (A.21)where is the principal value. One can then rewrite Eq. (A.20)under the form (A.22)Similarly, one can rewrite the second term of Eq. (A.17)under the form (A.23)Using symmetries of the matrix , starting from Eq. (A.23), one can show that the second contribution from Eq. (A.16)finally takes the form (A.24)Combining the two contributions obtained in Eqs. (A.22)and (A.24), one can rewrite Eq. (A.16)under the form (A.25)Hence, using the slow evolution Eq. (A.8), we recover the BalescuLenard equation introduced in Eq. (2). As a final remark, one may notice that on the short dynamical timescale, the evolution is governed by Eq. (A.7), which involves the fluctuating components δF and δψ of the distribution function and the potential. In contrast, on the long secular timescale, the evolution is governed by Eq. (A.8), which after an angleaverage only involves the mean distribution function F. Indeed, thanks to the ensemble average, all the crosscorrelations between the fluctuations δF and δψ, as in Eqs. (A.18), (A.19)and (A.23)can be expressed in terms of the underlying smooth distribution function F only, so that the fluctuating components are absent from the secular BalescuLenard collision operator from Eq. (A.25).
Appendix B: Turning off collective effects
The reference Chavanis (2013b; see also Appendix A of Chavanis 2012a) considers the inhomogeneous BalescuLenard equation without collective effects. This collisional kinetic equation is the equivalent of the Landau equation for inhomogeneous systems. It can be straightforwardly obtained as an approximation of the BalescuLenard Eq. (2). Indeed, one has to make the substitution from the dressed susceptibility coefficients to the bare ones given by A_{m1,m2}(J_{1},J_{2}), so that the inhomogeneous BalescuLenard equation without collective effects (i.e. the inhomogeneous Landau equation) is given by (B.1)where the coefficients A_{m1,m2} are associated with the Fourier transform in angles of the interaction potential (Pichon 1994; Chavanis 2013b) and read (B.2)where u(x) is the binary potential of interaction given by u(x) = − G/  x  in the gravitational case. A key remark at this stage is that the expression (B.2)does not require the introduction of biorthogonal basis elements as in Eq. (3), whereas in order to estimate the dressed susceptibility coefficients from Eq. (4), one must necessarily rely on Kalnajs’ matrix method (Kalnajs 1976). It is however possible to express A_{m1,m2} using the potential basis. Indeed, for a fixed value of x_{2}, we consider the function x_{1} → u(x_{1} − x_{2}). One can then decompose this function on the basis elements ψ^{(p)}(x_{1}), so that we may write (B.3)where it is important to note that the basis coefficients a_{p}(x_{2}) are functions of x_{2}. Thanks to the biorthogonality property of the basis detailed in Eq. (3), one can obtain the expression of the coefficients a_{p}(x_{2}) which reads (B.4)where we used the fact that u is a real function. Hence we finally obtain (B.5)Taking appropriately a Fourier transform with respect to the angles as in Eq. (B.2), we obtain (B.6)This fairly simple relation allows us to express the bare susceptibility coefficients A_{m1,m2} using the potential basis. One does not need anymore to perform a Fourier transform in angles of the interaction potential because the resolution of Poisson’s equation has been implicitly hidden in the effective construction of the basis elements^{7}. Neglecting the collective effects in the expression (4)of the dressed susceptibility coefficients amounts to taking , so that we obtain (B.7)The negative sign in Eq. (B.7)plays no significant role in the kinetic equations since one has to make the substitution of the square modulus in the BalescuLenard Eq. (2)to obtain the inhomogeneous Landau Eq. (B.1). Using our WKB basis, one can proceed similarly as in the dressed case, by limiting oneself only to local resonances. Equation (76)therefore becomes in the bare case (B.8)This expression of the bare susceptibility coefficients can be straightforwardly obtained from the dressed ones by imposing λ = 0. One should note that because of the absence of the amplification term 1 / (1 − λ_{kr}), the approximation of the small denominators cannot be used. Hence a physically motivated evaluation of the expression (B.8)becomes more subtle to perform, especially because of the possibly important role that the Coulomb logarithm 1 /k_{r} might play. Finally, one obtains that even for exactly local resonances, i.e. m_{1} = m_{2} and R_{1} = R_{2}, the use of our WKB formalism allowed us to obtain nondiverging bare susceptibility coefficients. On the other hand, one could try to estimate the bare susceptibility coefficients starting from Eq. (B.2), i.e. without using any potential basis. Using the polar coordinates (R,φ), one has to compute (B.9)By azimuthal symmetry, it is straightforward to show (Pichon & Cannon 1997) that (B.10)Hence the different m_{φ}modes of the A_{m1,m2} coefficients are independent. This result is identical to what was obtained in Eq. (65)for the dressed case. In order to illustrate the regularizing role of the WKB basis from Eq. (26), we will place ourselves in the context of an extremely tepid disc, and therefore assume . Thanks to the epicyclic mapping from Eq. (24), one can drop all the dependences on θ_{R} appearing in R and φ. Equation (B.9)then immediately implies that (B.11)and the bare susceptibility coefficients are given by (B.12)In this illustrative limit, we have by construction no contributions from the ILR and OLR resonances for which m_{r} ≠ 0. After an immediate change of variables, using , it becomes (B.13)The resonance condition from Eq. (68)takes the form m_{φ}Ω_{φ}(R_{1}) = m_{φ}Ω_{φ}(R_{2}) because we restricted ourselves in Eq. (B.11)to the case m_{r} = 0. Hence assuming that the function R_{g} → Ω_{φ}(R_{g}) is a monotonic function, one has to satisfy the constraint R_{1} = R_{2}, so that we are restricting ourselves only to local resonances as in Eq. (72). For such resonances, Eq. (B.13)becomes (B.14)At this stage, one must note that this integral is divergent. Indeed, for Δ → 0, one has . Hence the expression of the bare susceptibility coefficients derived from the Fourier transform of the interaction potential as in Eq. (B.2), when restricted to exactly local resonances becomes logarithmically divergent. This divergence, which is observed in the case of local resonances, is induced by the interaction of singular orbits. It is important to remark that this divergence was not present in Eq. (B.8)when computing the bare susceptibility coefficients using the WKB basis from Eq. (26). This implies that the WKB basis is not complete. For a complete biorthogonal basis, the expression (B.6)of the A_{m1,m2} coefficients is an exact expression. However, our calculation shows that the scaledecoupled WKB basis we used, possesses a subtle regularizing incompleteness which gets rid of the diverging contributions to the coefficients A_{m1,m2} in the limit of exactly local resonances.
Appendix C: The Schwarzschild DF case
When considering a Schwarzschild distribution function as in Eq. (25), while relying on the approximation of the small denominators from Eq. (78), one can explicitly perform the remaining integration on the radial action in the expressions (79)and (80)of the drift and diffusion coefficients. We now detail this explicit calculation. For such a Schwarzschild distribution function, it is straightforward to check that the gradients of the distribution function with respect to the actions are given by (C.1)Using the expression of the susceptibility coefficients from Eq. (78), after some simple algebra, one can rewrite the drift coefficients from Eq. (79)under the form (C.2)Similarly, the diffusion coefficients from Eq. (80)take the form (C.3)In Eqs. (C.2)and (C.3), in order to shorten the notations, we introduced the function defined as where we used the shortened notation (C.4)for the term appearing as a prefactor in the expressions (79)and (80)of the drift and diffusion coefficients. In addition to the integration formula (57), we may also rely on the additional identity (C.5)where α> 0, β> 0, and m_{r} ∈Z. In analogy with the definition from Eq. (58), we also introduce χ_{max} as (C.6)One can then immediately perform the integration on from Eq. (C.3), so that the diffusion coefficients are given by (C.7)where the function is defined as After some algebra, the drift coefficients from Eq. (C.2)are given by (C.8)where the function is defined as (C.9)These explicit expressions of the diffusion and drift coefficients obtained in Eqs. (C.7)and (C.8)allow us to estimate in a simple way the secular flux in the entire action space J = (J_{φ},J_{r}), once we assume that the distribution function is a Schwarzschild DF given by Eq. (25)and that the susceptibility coefficients can be approximated by Eq. (78).
Appendix D: Relation to other kinetic equations
The kinetic equation governing the collisional evolution of a system of N stars at the order 1/N is the inhomogeneous BalescuLenard Eq. (2). This equation conserves the total number of stars and the energy, and monotonically increases the Boltzmann entropy (HTheorem). We note that the collisional evolution of the system is due to a condition of resonance encapsulated in the term δ_{D}(m_{1}·Ω_{1} − m_{2}·Ω_{2}). In general, this condition can allow for local and nonlocal resonances. In the case of tepid discs considered in the present paper, assuming that only tightly wound spirals are sustained by the disc, we justified in Eq. (72)the fact that the resonances are purely local, so that m_{1} = m_{2} and . Furthermore, because of the epicyclic approximation, the intrinsic frequencies of the system, given by Eq. (22), depend only on J_{φ}, so that Ω = Ω(J_{φ}). Under these conditions, the BalescuLenard equation giving the collisional evolution of F = F(J_{φ},J_{r},t) may be rewritten as (D.1)The integration on is straightforward because of the δ_{D}function (local resonance), and we are left with (D.2)where the susceptibility coefficients are generally given by Eq. (77).
The kinetic Eq. (D.2)is an integrodifferential equation that governs the evolution of the system as a whole. It describes the effects of encounters between any test particle characterized by the angleaction coordinates (J_{φ},J_{r}) and the field particles characterized by the (running) angleaction coordinates . Actually, there is no distinction between test and field particles, so that they are characterized by the same distribution function F(·,t) that evolves in a selfconsistent manner, hence the integrodifferential character of the kinetic equation. This is a characteristic of the BalescuLenard equation describing the evolution of the system as a whole.
Appendix D.1: FokkerPlanck limit
We can also use this formalism to directly obtain the FokkerPlanck equation governing the relaxation of a test star in a bath of field stars, assumed to be in a steady state with a distribution function . Proceeding as in Chavanis (2012a), we just have to replace in Eq. (D.2)the distribution function of the field particles by the static distribution , while the time evolving distribution function of the test particle is rewritten as P(J_{φ},J_{r},t) for clarity. This heuristic procedure is justified in Chavanis (2012a) by an explicit calculation of the diffusion and drift coefficients of the FokkerPlanck equation. It transforms the integrodifferential Eq. (D.2)into a differential equation (D.3)which can be interpreted as a FokkerPlanck equation. If we assume that the field particles are at statistical equilibrium (thermal bath), described by the Boltzmann distribution (D.4)then, using the relation ∂F_{0}/∂J = − βF_{0}Ω (see the definition of Ω in Eq. (1)), we can reduce the FokkerPlanck Eq. (D.3)to the form (D.5)with the diffusion coefficient (D.6)We note that the friction term in the FokkerPlanck Eq. (D.5)is proportional and opposite to the intrinsic frequency Ω, and that the friction coefficient ξ satisfies a (generalized) Einstein relation ξ = Dβ (for each resonance) (Chavanis 2012a).
This FokkerPlanck formalism may have other applications. For example, if we consider a system with two species of particles (e.g. characterized by different masses m_{A} and m_{B}), and if species B has reached an equilibrium state with a distribution , we can use Eq. (D.3)to describe the relaxation of particles of species A due to the encounters with particles of species B (but neglecting the encounters between particles of species A). This approach is further discussed in Appendix F of Chavanis (2013b). On the other hand, for a single species system, we may describe the early dynamics of the system as a whole with a good approximation by replacing in the BalescuLenard Eq. (D.2)by the initial distribution function , leading again to an equation of the form (D.3).
Appendix D.2: Other kinetic equations
We now compare the previous results to other related kinetic equations. For spatially homogeneous systems with longrange interactions, the BalescuLenard equation reads (see Chavanis 2012c): (D.7)where (D.8)is the dielectric function. For onedimensional systems, it reduces to the trivial form (D.9)In d> 1, there are always resonances between particles with different velocities, implying that the BalescuLenard equation relaxes towards the Boltzmann distribution on a timescale of the order Nt_{D}, where t_{D} is the dynamical time. By contrast, in d = 1, the resonances become local (in velocity space) and since the term in brackets is antisymmetric with respect to the interchange of v and v′, the BalescuLenard diffusion current vanishes exactly. This implies that the relaxation time is larger than Nt_{D}, presumably of the order N^{2}t_{D}, corresponding to the next order term in the expansion of the dynamics in powers of 1 /N. The FokkerPlanck equation describing the evolution of a test particle in a bath of field particles is discussed in Chavanis (2012c, 2013a).
The kinetic equation governing the collisional evolution of a system of N point vortices in twodimensional hydrodynamics at the order 1 /N can be written, in the axisymmetric case, as (see Chavanis 2012b): (D.10)where ω(r,t) is the profile of vorticity, Ω(r,t) the profile of angular velocity, and χ(r,r′,Ω(r,t)) is related to the dressed potential of interaction between the point vortices (see Chavanis 2012b for more details). This equation conserves the total number of point vortices and the energy, and monotonically increases the Boltzmann entropy (Htheorem). If the profile of angular velocity is monotonic, the kinetic equation reduces to the form (D.11)For nonmonotonic profile of angular velocity, one can have nonlocal resonances (i.e. distant collisions between point vortices), as studied in Chavanis & Lemou (2007). This produces a diffusion current. If the profile of angular velocity is, or becomes, monotonic, the resonances are purely local and, since the term in brackets is antisymmetric with respect to the interchange of r and r′, the diffusion current also vanishes. This implies that the relaxation time is larger than Nt_{D} as discussed above. The FokkerPlanck equation describing the evolution of a test vortex in a sea of field vortices is discussed in Chavanis (2012b).
If we focus on purely local resonances, we note that the inhomogeneous BalescuLenard Eq. (D.1)is different from Eqs. (D.9)and (D.11)because it is twodimensional in J_{r} and J_{φ}, and the resonances act only on J_{φ}. Therefore, purely local resonances do not yield a zero flux, contrary to Eqs. (D.9)and (D.11). This really is an effect of the twodimensionality of the system. Indeed, for local resonances, the 1D inhomogeneous equation also yields a zero flux: (D.12)
Appendix E: The Schwarzschild conspiracy
The Schwarzschild distribution function introduced in Eq. (25)and considered in S12 simulation has the specificity to be exponential in the J_{r}direction, so that F is Boltzmannian with respect to the J_{r} variable. It is known (e.g. Chavanis 2012a) that the (complete) Boltzmann distribution is the steady state of the BalescuLenard Eq. (2). Therefore, we can expect that the (partial) exponential behaviour of the Schwarzschild distribution will induce simplifications that we now detail. In analogy with Eq. (10), the flux associated with a given resonance m is defined as (E.1)where the nonbold ℱ_{m} is a scalar. Using shortened notations and forgetting numerical prefactors, we may rewrite the drift and diffusion coefficients from Eqs. (79)and (80)under the form (E.2)where we used the same shortened notation for the resonant factor 1 / (m·Ω)′ as in Eq. (C.4). The flux can then be decomposed as with (E.3)and (E.4)We are interested in the value of the flux at the initial time, where the distribution function is given by the Schwarzschild distribution. Because of the exponential dependence in J_{r} of the Schwarzschild distribution function, one has . As a result, the radial component (E.3)of the flux cancels out and the flux is simply given by Eq. (E.4). This coincidence could be called the Schwarzschild conspiracy and has important consequences on the properties of the collisional diffusion. Indeed, for a tepid disc, one has  ∂F/∂J_{r}  ≫  ∂F/∂J_{φ} . Hence one would expect the gradients with respect to J_{r} to be the major contributors to the diffusion. When considering independently the drift and diffusion coefficients, the gradients in ∂F/∂J_{r} dominate the diffusion current. Thus for m_{r} ≠ 0, the diffusiononly flux can be approximated by (E.5)As the ILR and the OLR have a nonzero m_{r} compared to the COR, these resonances should dominate independently the drift and diffusion components. However, when considering the full flux made of the contributions from the drift and diffusion coefficients, because of the Schwarzschild conspiracy, there is a simplification of the dominant terms in ∂F/∂J_{r}, so that one recovers as in Eq. (E.4)that only the smaller gradients ∂F/∂J_{φ} remain present. The Schwarzschild conspiracy between drift and diffusion will therefore tend to slightly reduce the magnitude of the full diffusion flux, so as to moderately slow down the collisional relaxation. More importantly, the Schwarzschild conspiracy will favour the COR resonance (radial migration) over the ILR resonance (J_{r}heating). One should note that the DF which is effectively sampled in S12 is of the form , with (Toomre 1977; Binney & Tremaine 2008). It is only within the epicyclic approximation that this DF takes the form of the Schwarzschild DF from Eq. (25). As a consequence, in S12 simulation, the Schwarzschild conspiracy is not exactly satisfied as observed in Eq. (E.3), but the residual difference driving the secular diffusion is likely to be subdominant, as illustrated in Fig. E.1.
Fig. E.1 Contours in actionspace (J_{φ},J_{r}) of the difference between the sampled anisotropic DF of S12 simulation and its Schwarzschild epicyclic approximation F_{Sch} from Eq. (25). The plotted quantity is , and contours labels are expressed in percentages. 

Open with DEXTER 
Appendix F: Temporal frequency selection
An important feature of the diffusion Eq. (7)is that the diffusion takes place along specific resonance directions associated with the vectors m as discussed in Eq. (16). Hence being able to determine the dominant resonance is crucial in order to estimate the direction of the secular diffusion in action space. The temporal frequency associated with a given resonance m in a location J of actionspace is given by ω = m·Ω. Thanks to the expression (83), one immediately obtains that for a Mestel disc, one has (F.1)In Fouvry & Pichon (2015), Fouvry et al. (2015a), we studied the same S12 simulation using the WKB limit of the secular diffusion equation, which intends to describe the secular forcing of a collisionless selfgravitating system perturbed by an external source. An essential assumption of this approach was to consider the external perturbation as originating from numerical Poisson shot noise, and therefore assume it to be proportional to the local active surface density. The autocorrelation of the external perturbation ψ^{ext} was taken to be equal to (F.2)One can remark that this crude assumption on the noise properties has no ω dependence, so that all resonances are equally favoured by the Poisson shot noise and are perturbed similarly whatever their associated intrinsic frequencies ω = m·Ω. This ad hoc and simple noise approximation is one of the limitations of the formalism presented in Fouvry & Pichon (2015), Fouvry et al. (2015a). In contrast, in the WKB BalescuLenard equation described in this paper, this preferential selection of the resonances based on their intrinsic frequency is naturally present. Indeed, one can note in the expression (E.4)of the flux associated with a resonance m, the presence of the prefactor 1 / (m·Ω)′
which arose in Eq. (67)when handling the resonance condition. For the Mestel disc, whose intrinsic frequencies are given by Eq. (83), this term can be straightforwardly computed and reads (F.3)Comparing the ILR resonance to the OLR and COR, one immediately obtains that (F.4)Hence because the ILR resonance is associated with lower intrinsic temporal frequency ω = m·Ω, the resonant factor 1 / (m·Ω)′ naturally tends to favour the ILR resonance with respect to the OLR and COR, and therefore performs natively a temporal frequency biasing which was absent from the ad hoc assumption of Eq. (F.2)describing the external forcing considered in Fouvry & Pichon (2015), Fouvry et al. (2015a).
One can even be more specific when comparing the ILR and OLR resonances. The only difference between these two resonances is the sign of m_{r}. The expression (62)of the amplification eigenvalues shows that its value only depends on s^{2}, so that λ_{ILR} = λ_{OLR}. The expression of the susceptibility coefficients from Eq. (77)is also independent of the sign of m_{r} so that . Hence, when considering the flux ℱ_{m} given by Eq. (E.4), one notes that between the ILR and OLR resonances, ℱ_{m} only changes through the factor 1 / (m·Ω)′. Thanks to Eq. (F.4), one immediately obtains (F.5)The secular diffusion flux associated with the OLR resonance is therefore always much smaller than the one associated with the ILR resonance because of this effect of temporal frequency biasing.
As a conclusion, the temporal frequency selection effect described in Eq. (F.4)will tend to favour the ILR resonance because it is associated with a smaller intrinsic frequency. However, one can observe in Fig. 8 that the COR resonance is always more amplified than the ILR resonance. Finally, the crucial remark is to note that the susceptibility coefficients from Eq. (77)involve Bessel functions , which are such that if m_{r} = 0, or = 0 otherwise. As a consequence, close to the J_{r} = 0 axis, the COR resonance will always tend to become the dominant resonance. There is therefore a nontrivial arbitration between these opposite effects when considering the respective contributions of the various resonances to the full secular diffusion flux.
All Figures
Fig. 1 Top: a) a set of two resonant orbits in the inertial frame; b) in the rotating frame in which they are resonant – here through ILRCOR coupling. Bottom: c) fluctuations in the distribution function in actionspace caused by finiteN effects showing overdensities for the blue and red orbits. The dashed lines correspond to 3 contour levels of the intrinsic frequency ω = m·Ω respectively associated with the resonance vector m_{1} (gray lines) and m_{2} (black lines). The two sets of orbits satisfy the resonant condition m_{1}·Ω_{1} = m_{2}·Ω_{2}, and therefore lead to a secular diffusion of the orbital structure of the disc according to Eq. (2). One should emphasize that the resonant orbits need not be caught in the same resonance (m_{1} ≠ m_{2}), be close in position space, nor in action space. 

Open with DEXTER  
In the text 
Fig. 2 Two WKB basis elements. Each Gaussian is centred around a radius R_{0}. The typical extension of the Gaussian is given by the decoupling scale σ, and they are modulated at the radial frequency k_{r}. 

Open with DEXTER  
In the text 
Fig. 3 Two WKB basis elements in the polar (R,φ)plane. Each basis element is located around a central radius R_{0}, on a region of size σ. The winding of the spirals is governed by the radial frequency k_{r}. The number of azimuthal patterns is given by the index k_{φ}, e.g. k_{φ} = 1 for the interior dark gray element, whereas k_{φ} = 2 for the exterior light gray one. 

Open with DEXTER  
In the text 
Fig. 4 Surface density Σ_{t} of the tapered Mestel disc. The unit system has been chosen so that V_{0} = G = R_{i} = 1. Because of the tapering functions, the selfgravity of the disc is turned off in the inner and outer regions. 

Open with DEXTER  
In the text 
Fig. 5 Contours of the initial distribution function in actionspace (J_{φ},J_{r}), within the epicyclic approximation. The contours are spaced linearly between 95% and 5% of the distribution function maximum. 

Open with DEXTER  
In the text 
Fig. 6 Dependence of the local Q Toomre parameter with the angular momentum. It is scale invariant except in the inner/outer regions because of the presence of the tapering functions T_{inner} and T_{outer}. The unit system has been chosen so that V_{0} = G = R_{i} = 1. 

Open with DEXTER  
In the text 
Fig. 7 Variations of the response matrix eigenvalues λ with the WKBfrequency k_{r}, for m = m_{COR} and two values of J_{φ}. The curve that peaks at large k_{r} is for the smaller value of J_{φ}. 

Open with DEXTER  
In the text 
Fig. 8 Dependence of the amplification factor 1 / (1 − λ_{max}) with the position J_{φ} in the disc. The amplification associated with the COR is always larger than the one associated with the ILR and OLR. 

Open with DEXTER  
In the text 
Fig. 9 Map of N div(ℱ_{tot}), where the total flux ℱ_{tot} has been summed over the three resonances (ILR, COR and OLR). Red contours, for which N div(ℱ_{tot}) < 0, correspond to regions from which the orbits will be depleted, whereas blue contours, for which N div(ℱ_{tot}) > 0, correspond to regions where secular diffusion will tend to increase the value of the DF. The net fluxes involve simultaneously radial migration near (J_{φ},J_{r}) ~ (1.8,0), and heating near (J_{φ},J_{r}) ~ (1,0.1). 

Open with DEXTER  
In the text 
Fig. 10 Overlay of the WKB predictions for the divergence of the diffusion flux N div(ℱ_{tot}) and the differences between the initial and the evolved DF in S12’s simulation. The opaque contours correspond to the differences in the actionspace for the DF in S12 between the time t_{S12} = 1000 and t_{S12} = 0 (see the upper panel of S12’s Fig. 10). The red opaque contours correspond to negative differences, so that these regions are emptied from their orbits, whereas blue opaque contours correspond to positive differences, i.e. regions where the DF has increased through diffusion. The transparent contours correspond to the predicted values of N div(ℱ_{tot}) using the same conventions as in Fig. 9. We note the overlap of the predicted transparent red and blue contours with the measured solid ones. 

Open with DEXTER  
In the text 
Fig. 11 Overlay of the WKB predictions for the divergence of the diffusion flux N div(ℱ_{tot}) on top of the contours of the DF in actionspace measured in the S12 simulation. The black background contours are the levels contours of the DF at time t_{S12} = 1400 (see the lower panel of Fig. 7 of S12). These contours are spaced linearly between 95% and 5% of the DF maximum and exhibit clearly the appearance of a resonant ridge. The coloured transparent contours correspond to the predicted values of N div(ℱ_{tot}) using the same conventions as in Fig. 9. We note that the developed late time ridge is consistent with the predicted depletion (red) and enrichment (blue) of orbits. 

Open with DEXTER  
In the text 
Fig. E.1 Contours in actionspace (J_{φ},J_{r}) of the difference between the sampled anisotropic DF of S12 simulation and its Schwarzschild epicyclic approximation F_{Sch} from Eq. (25). The plotted quantity is , and contours labels are expressed in percentages. 

Open with DEXTER  
In the text 