Online material
Appendix A: Kalnajs basis
We now detail the properties of the basis introduced in Kalnajs (1976) to describe 2D discs^{5}. The basis depends on two parameters: an index k_{Ka} ∈N and a scale radius r_{Ka} ∈ R^{+}. In all the upcoming formulae of this section, in order to shorten the notations, we write r for the dimensionless quantity r/r_{Ka}. As introduced previously in Eq. (16), the basis elements depend on two indices: the azimuthal number ℓ and the radial index n. One should note that we have ℓ ≥ 0 and n ≥ 0. The radial component of the potential elements are then of the form (A.1)The radial component of the density elements is given by (A.2)In Eqs. (A.1) and (A.2), the coefficients and are defined by (A.3)and (A.4)Finally, in Eqs. (A.1) and (A.2), we have also introduced (A.5)and (A.6)In the two previous expressions, we introduced the rising Pochhammer symbol [ a ] _{i} defined as (A.7)
Appendix B: Calculation of ℵ
We now detail how the analytical function ℵ from Eq. (32) may be computed. In order to ease the implementation of its computation, we rewrite ℵ in an undimensionalised way as follows (B.1)where we assumed that a_{g},a_{h} ≠ 0 and used the change of variables x = x_{p}/ Δr and y = x_{a}/ Δr. We also defined the dimensionless function ℵ_{D} as (B.2)To compute this integral, we may now exhibit a function G(x,y) such that (B.3)One possible choice for G is given by (B.4)One should note in the previous expression the presence of a complex logarithm and a tan^{1}. However, because e,f,η ∈R and η ≠ 0, one can easily show that the arguments of both of these functions never cross the usual branchcut of these functions {Im(z) = 0; Re(z) ≤ 0 }. As a consequence, the expression (B.2) can immediately be computed as (B.5)
Appendix C: Response matrix and Nbody validations
The computation of the response matrix as described in Sect. 3 was validated by recovering the results of the pioneer work of Zang (1976), extended in Evans & Read (1998a,b), and recovered numerically in Sellwood & Evans (2001). These papers predicted the precession rate ω_{0} = m_{φ}Ω_{p} and growth rate η = s of the unstable modes of a truncated Mestel disc similar to the stable one described in Sect. 4.1. To build up an unstable disc similar to the ones considered in these previous works, one has to consider a fully active disc, so that ξ = 1. So as to have Q = 1 (Toomre 1964), the velocity dispersion within the disc is given by q = 6, where the parameter q has been introduced in Eq. (46). Finally, a last parameter one can tune in order to modify the properties of the disc is the truncation index of the inner tapering ν_{t} defined in Eq. (48). While looking only for m_{φ} = 2 modes, we considered three different truncations indices given by ν_{t} = 4, 6, 8. To compute the response matrix, we used the same numerical parameters as described in Sect. 4.2. Looking for unstable modes amounts to looking for complex frequencies ω = ω_{0} + iη such that the response matrix from Eq. (24) possesses an eigenvalue equal to 1. Such a complex frequency is then associated with an unstable mode of pattern speed ω_{0} and growth rate η. To determine the growth rate and pattern speed of the unstable modes, we relied on Nyquist contours similarly to the technique presented in Pichon & Cannon (1997). For a fixed value of η, one can study the continuous complex curve . Because for η → + ∞, one has , the number of windings of this curve around the origin gives a lower bound on the number of unstable modes with a growth rate superior to η. By decreasing the value of η, one can then determine the largest value of η admitting an unstable mode, and therefore the most unstable mode of the disc. The Nyquist contours obtained for the truncation index ν_{t} = 6 are illustrated in Fig. C.1, while the measurements are gathered in Table C.1.
Fig. C.1
Top panel: zoomed Nyquist contours in the complex plane of obtained via the matrix method for a truncated Mestel disc with ν_{t} = 6 and q = 6, looking for m_{φ} = 2 modes. Each contour corresponds to a fixed value of η. For a growth rate of η ≃ 0.20, one can note that the contour crosses the origin, which corresponds to the presence of an unstable mode. Bottom panel: illustration of the function for the same truncated Mestel disc. Each line corresponds to a fixed value of η. Such a representation allows us to determine the pattern speed ω_{0} = m_{φ}Ω_{p} ≃ 0.94 of the unstable mode. 

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

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

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

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

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

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

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

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

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

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

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

The calculation of the 2D response matrix via Eq. (33).
It has been tested for the isochrone and the Mestel discs. A copy of the code is also available at the CDS.
© ESO, 2015