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.
© ESO, 2015