Efficient numerical calculation of drift and diffusion coefficients in the diffusion approximation of kinetic equations

In this paper we study the diffusion approximation of a swarming model given by a system of interacting Langevin equations with nonlinear friction. The diffusion approximation requires the calculation of the drift and diffusion coefficients that are given as averages of solutions to appropriate Poisson equations. We present a new numerical method for computing these coefficients that is based on the calculation of the eigenvalues and eigenfunctions of a Schr\"odinger operator. These theoretical results are supported by numerical simulations showcasing the efficiency of the method.


Introduction
In this paper we consider models from statistical physics that have the general form (1.1) right-hand side, Q is a linear operator, which can be either an integral or a differential operator with respect to the variable v; it is intended to describe some collision or friction-dissipation mechanisms. Of course, the asymptotics is driven by the properties of the leading operator Q: • We assume that Q is conservative in the sense that Q f ε dv = 0 holds. Consequently, integrating (1.1) with respect to the velocity variable, we obtain the local conservation law ∂ t ρ ε + ∇ r · J ε = 0, where we denote ρ ε (t, r) = f ε (t, r, v) (t, r, v) dv.
• We also suppose that the kernel of Q is spanned by a positive normalized function v → F(v). As we shall detail below on specific examples, Q also satisfies some dissipation properties: roughly speaking, these properties encode the fact that Q forces the distribution function to become proportional to the equilibrium F. In turn, the dissipation permits us to establish estimates which lead to the ansatz f ε (t, r, v) = ρ ε (t, r)F(v) + √ εg ε (t, r, v). In order to have a current of order 1, we should assume vF dv = 0, so that J ε = g ε dv.
Finally, let us suppose that ρ ε , J ε and g ε have well-defined limits as ε → 0, denoted ρ, J and g, respectively. Multiplying (1.1) by √ ε and letting ε go to 0 yields which provides us with a closure for the equation obtained by passing to the limit in the conservation law that the numerical simulation of an equation like (1.2) is far less costly than the one of (1.1); on the one hand, we have eliminated the velocity variable, on the other hand, when ε → 0, (1.1) contains stiff terms that induce prohibitive stability conditions. However, we are left with the difficulty of calculating the effective diffusion and drift coefficients D and K which relies on being able to calculate the inverse of the operator Q, i.e. on solving equations of the form Qφ = h. We address these questions in the specific case of kinetic models for swarming for which the auxiliary equations that determine D and K do not have explicit solutions. In recent years, increasing efforts have been devoted to the development of mathematical models able to describe the self-organization of a large set of living organisms (fish, birds, bacteria, . . .), after the pioneering work of Vicsek et al. (1995). Based on simple rules of information exchange between the individuals about their close environment, the whole population organizes itself in remarkable patterns. Starting from individual-based descriptions, kinetic equations can be derived by means of mean-field regimes, as is usual in statistical physics (Dobrušin, 1979;Spohn, 1991;Sznitman, 1991;Golse, 2003); we refer the reader to Ha & Tadmor (2008), Carrillo et al. (2010Carrillo et al. ( , 2011, Bolley et al. (2011), Motsch & Tadmor (2011 and Vecil et al. (2013) for a thorough study of such models. For the models we are interested in, the interaction operator Q has the form that involves a quite complicated potential function v → W (v). The role of this potential function in swarming models is to fix an asymptotic speed of the individuals and several possibilities of repulsive at the origin and attractive at infinity potentials have been used in the swarming literature (Levine et al., 2000;D'Orsogna et al., 2006). These potentials in swarming models have been one of the motivations of this work. However, we will treat this problem in general since it has potential applications in other areas where the diffusion approximation or overdamped limit with nonexplicit eigensystems is important. The equilibrium function simply reads F(v) = (1/Z) e −W (v) , where Z denotes the normalization constant, but inverting Q is not that simple. In contrast to the standard Fokker-Planck operator corresponding to the Langevin dynamics with linear friction, i.e. when W (v) = v 2 /2, in general v → vF(v) is not an eigenfunction of Q, and we do not have explicit formulas for the effective coefficients. Our approach for computing the coefficients is based on the spectral properties of the operator. In fact, a unitary transformation enables us to reformulate the Poisson equation Qφ = h in the form Hu = μ, where H is a Schrödinger operator associated to a certain potential Φ, that depends on W ; see equation (3.1). In turn, the operator admits a spectral decomposition. Then, our method is based on expanding the functions on the eigenbasis, and then computing the coefficients of φ in this basis. The latter step does not present any difficulty and the computational effort is concentrated on the determination of the eigenelements. In practice, we work on the discrete form of the equations, and we expect that only a few eigenmodes are necessary to capture the effective coefficients. This is the strategy we shall discuss, adopting high-order finite elements discretizations of the Schrödinger equations, which allows us to make use of performing computational tools (Martin, 2010).
We mention now some related problems, to which the numerical method developed in this paper can also be applied, and alternative approaches for the calculation of the coefficients that appear in macroscopic equations. First, it should be noted that a similar formalism, based on the calculation of the eigenfunctions and eigenvalues of the linearized collision operator, has been developed for the calculation of transport coefficients using the linearized Boltzmann equation (Resibois & De Leener, 1977). Second, the method can be adapted to compute the effective, possibly space-dependent, coefficients that come from the homogenization process of advection-diffusion equations with periodic coefficients (e.g. Bensoussan et al., 1978 andPoupaud, 2004), the diffusion approximation for fast/slow systems of stochastic differential equations (Pardoux & Veretennikov, 2001, 2003, 2005, and in connection to functional central limit theorems for additive functionals of Markov processes (Komorowski et al., 2012). In all these problems, the drift and diffusion coefficients that appear in the macroscopic equation can be expressed in terms of the solution of an appropriate Poisson equation. These effective coefficients can be alternatively calculated using either Monte Carlo simulations, e.g. Pavliotis et al. (2006) and Boozer & Kuopetravic (1981), the heterogeneous multiscale method (Weinan et al., 2005), numerical solution of the Poisson equation using spectral methods (e.g. Majda &McLaughlin, 1993 andPavliotis, 2002). Other numerical approaches include the use of linear response theory and of the Green-Kubo theory (Joubaud & Stoltz, 2012) and the expansion of the solution to the Poisson equation in appropriate orthonormal basis functions, e.g. Hermite polynomials. This technique, which is related to the continued fraction expansion (Risken, 1989) has been applied to the calculation of the diffusion coefficient for the Langevin dynamics in a periodic potential (Pavliotis & Vogiannou, 2008) and to simple models for Brownian motors (Latorre et al., 2014).
The rest of the paper is organized as follows. In Section 2, we introduce the kinetic model for swarming we are interested in. By using the Hilbert expansion, we detail the diffusion asymptotics, and we discuss some properties of the effective coefficients. The convergence can be rigorously justified, and a complete proof is given in the appendix. In Section 3, we switch to the formalism of Schrödinger equations and we explain in detail our numerical strategy. We also provide estimates to justify the approximation by truncating the Fourier series. Section 4 is devoted to the numerical illustration of the approach, dealing with a relevant set of potentials W . In particular we bring out the role of the scaling coefficients that appear in the potential, with difficulties related to the so-called tunnelling effect (Helffer & Sjöstrand, 1984;Robert, 1987). Section 5 is reserved for conclusions. The proof of Theorem 2.1 on the mean-field limit approximation can be found in the appendix.

Motivation: swarming models
The analysis of interacting particle systems with random noise finds applications in collective behavior and self-organization of ensembles of large individuals. These models lead typically to systems with nonstandard friction terms. In particular, the following example was proposed in the literature in D' Orsogna et al. (2006), which includes the effect of self-propulsion and a Rayleigh-type friction to impose an asymptotic cruising speed for individuals. This system, with noise and linear Stokes friction, was considered in Carrillo et al. (2009). For N particles, it reads where Γ i (t) are N independent Brownian motions with values in R d and σ > 0 is the noise strength.
Here, α > 0 is the self-propulsion strength generated by the organisms, β is the friction coefficient and γ = α/β is the squared asymptotic cruise speed of the individuals. Note that α − β|v| 2 v = −α∇ v W (v) with W (v) = 1 4γ |v| 4 − 1 2 |v| 2 . Different models for friction can lead to asymptotic fixed speeds for individuals; see for instance Levine et al. (2000). It is interesting, therefore, to study properties of more general potentials in velocity W (v), with the basic behaviour of having a confinement when the speed is large and that zero speed is a source fixed point of the speed dynamics. The class of potentials for the velocity that we will consider in this paper includes, for instance, potentials of the form The mean-field limit of the stochastic particle system (2.1) above under suitable assumptions on the interaction potential U is given by the following kinetic Fokker-Planck equation (see Bolley et al., 2011): and stands for the usual convolution product with respect to the space variable.
Let us first remark that T F = 1/α is the natural relaxation time for particles to travel at asymptotic speed √ α/β. We introduce the time and length scales T and L, which are determined by the time/length scales of observation. They define the speed unit U = L/T that will be compared to V, the typical particle speed, and V th = √ σ/α, the typical value of fluctuations in particle velocity, called the thermal speed. Then we can define dimensionless variables, denoted by primed quantities, as Note that C has the dimension of velocity. With this rescaling, we obtain the dimensionless kinetic equation where primes have been eliminated for notational simplicity and where the operator Q is defined as with W (v) = 1 4γ |v| 4 − 1 2 |v| 2 for the problem corresponding to (2.1). Here, η 1 , η 2 , γ , θ and ε are dimensionless parameters given by The operator Q defined in (2.3) is the Fokker-Planck operator corresponding to the stochastic differential equation where Γ (t) denotes standard d-dimensional Brownian motion. The generator of this diffusion process is Now, assume the following relation between the dimensionless parameters, with a finite asymptotic dimensionless speed,

With these assumptions, equation (2.2) becomes
In this regime, the dominant mechanisms are the noise and the nonlinear friction. This scaling is the so-called diffusion scaling for kinetic equations; see Bouchut et al. (2000) and Degond et al. (2000) and the references therein. We remark that for the particular application of swarming, other distinguished limits can also be considered; see Carrillo et al. (2009) for details. Different interesting features of the model arise, depending on the relative magnitudes of the scaling parameters γ , θ .
In order to obtain a closed macroscopic equation for the density ρ in the limit ε → 0, we use the standard Hilbert expansion method. Inserting the Hilbert expansions into (2.2) and identifying terms with equal power of √ ε, we get is the Maxwellian distribution associated to the Fokker-Planck operator Q, and Z denotes the normalization constant. This is particularly clear by rewriting the Fokker-Planck operator Q as (2.7) • ε −1/2 terms: To invert this equation, we need to solve the following problems: Note that h dv = 0 is clearly a necessary condition for the problem Q(f ) = h to admit a solution.
Assuming that the potential W increases sufficiently fast as |v| → +∞ and given that the right-hand side of (2.8b) is equal to −∂M (v)/∂v i , from the divergence theorem we deduce that the solvability condition is satisfied for (2.8b). The solvability condition is satisfied for (2.8a) under, for example, the assumption that the velocity potential W is spherically symmetric. In Section 4.3, we will see how the case of nonsymmetric potentials where the compatibility condition (2.8a) is not fulfilled can be dealt with.
Existence of a solution for equations of the form Q(f ) = h relies on the possibility to apply the Fredholm alternative. For this it is sufficient to show that Q has compact resolvent in the space L 2 (R d ; M −1 (v)). This follows, for example, by assuming that the potential satisfies See for example, Villani (2009, Theorem A.19). Under this assumption we can apply the Fredholm alternative to obtain • ε 0 terms: However, using again the compatibility condition, we conclude that where the diffusion and drift matrices are given by Therefore, in the ε → 0 limit regime we expect the macroscopic density to be well approximated by the solution to the aggregation-diffusion equation (2.10). We remark that by settingχ = χ M andκ = κM , the Poisson equations (2.8) become where L, defined by (2.5), denotes the generator of the diffusion process t → v(t) defined in (2.4). The drift and diffusion coefficients (2.11) are also given by the formulas This is the form of the Poisson equation that appears in functional central limit theorems and in the diffusion approximation for fast/slow systems of stochastic differential equations (Ethier & Kurtz, 1986;Komorowski et al., 2012). Furthermore, we note that the diffusion matrix is positive definite and thus, we can talk properly about a diffusion matrix. To show this, we use (2.7) to deduce (2.13) Now, given any vector ξ ∈ R d \ {0}, we can compute The strict equality follows from the fact that χ · ξ/M | = const, as we can immediately deduce from (2.8a).
The analysis of the asymptotic regime remains technically close to the derivation of the diffusion regimes for the Vlasov-Poisson-Fokker-Planck equation (Poupaud & Soler, 2000;Goudon, 2005;Masmoudi & Tayeb, 2007;El Ghani & Masmoudi, 2010). We detail in the appendix the proof of the following statement.
Theorem 2.1 Let us consider a sequence of initial data f ε Init 0 that satisfies We suppose that the potentials U and W satisfy the technical requirements listed in the appendix A. Let 0 < T < ∞. Then up to a subsequence, still labelled ε, the solution f ε associated to (2.6) converges , with ρ being the solution to (2.10) and initial data ρ(t = 0, r) given by the weak limit in L 1 (R d ) of f ε Init dv. For the time being, let us discuss the numerical evaluation of the effective coefficient matrices D, K. In the particular case considered here, the right-hand side of both equations in (2.8) (or equivalently (2.12)) is of the form v times a radial function; then we can simplify the diffusion and drift matrices by taking into account the symmetries of the problem. We leave the reader to check the following result.
Lemma 2.2 Given v ∈ R d , let us define for any indices i, j the linear operators that exchange v i and v j and change v i by −v i , respectively. Then, the following relations hold: As a consequence, we deduce that there exist reals D > 0 and κ such that and the macroscopic aggregation-diffusion equation (2.10) becomes (2.14) Therefore, in order to compute the effective macroscopic equation (2.14), we need to find the solutions to one component of (2.8) (or (2.12)). They are given by explicit formulas in very few cases. A particular example is given by the quadratic potential W (v) = |v| 2 /2: then This is due to the fact that −vM (v) is an eigenfunction of the Fokker-Planck operator Q with quadratic potential. However, for more general, nonquadratic potentials such as the one used in the swarming model, it is not possible to obtain explicit formulas for the coefficients of the limiting equation (2.14).
In the next section we will study this problem by eigenfunction expansion of the Fokker-Planck operator and we will discuss how to accurately approximate those coefficients.

Approximation of the diffusion and drift coefficients
We start by recalling a well-known (unitary) equivalence between Fokker-Planck (in L 2 (R d ; M −1 dv)) and Schrödinger operators (in L 2 (R d )). By setting it is easy to check that H reduces to the Schrödinger operator with the potential We remark that Φ is precisely the potential that appears in assumption (2.9). The operator We point out that this also defines the domain of the operator Q. When working with the operator H, the Lebesgue space L 2 (R d ) is a natural framework; in what follows, we shall denote by (· | ·) the standard inner product in L 2 (R d ). Using classical results for Schrödinger operators (see for instance Reed & Simon, 1978, Theorem XIII.67), we have a spectral decomposition of the operator H under suitable confining assumptions.
We remark that the spectral gap Λ of the Schrödinger operator H is the Poincaré constant in the Poincaré inequality associated to the Fokker-Planck operator Q, i.e. the Poincaré inequality for the probability measure M (v) dv = Z −1 e −W (v)/θ dv.
Note that property (2.13) implies that the operator H is positive and that the kernel is spanned by √ M . We wish to solve the equations or equivalently, The diffusion and drift coefficients are defined by the quadratic quantities The Schrödinger operator is self-adjoint in L 2 (R d ) and under assumption (2.9) it has compact resolvent. We denote its eigenvalues and eigenfunctions by {λ n , Ψ n } +∞ n=1 . A lot of information on the properties of the eigenvalues and eigenfunctions of Schrödinger operators is available (Reed & Simon, 1978). Using the spectral decomposition of H we can obtain a formula for the effective drift and diffusion coefficients. These formulas are similar to the Kipnis-Varadhan formula for the diffusion coefficient in the functional central limit theorem for additive functionals of reversible Markov processes (Kipnis & Varadhan, 1986). We can use these formulas to develop a numerical scheme for the approximate calculation of the drift and diffusion coefficients. Indeed, if we develop H(u χ ), u χ , H(u κ ) and u κ in the eigenbasis, we get Substituting in (3.2) and (3.3), we obtain the following formulas for the diffusion and drift matrices: We can obtain approximate formulas for the drift and diffusion coefficients by truncating the data: given h ∈ L 2 (R d ) and denoting We denote by u N the solution of H(u N ) = h N , and compare it with u, the solution of H(u) = h. We use the Sobolev and Cauchy-Schwarz inequalities as follows: Hence, we get The accuracy of the approximation is therefore driven by the accuracy of the approximation of the data h by its truncated Fourier series; we need information on the behaviour of the eigenvalues λ n for large n and on the accuracy of the spectral projection of h. We proceed by analogy to the standard theory of Fourier series, where the behaviour of the Fourier coefficients is related to the regularity of the function. More precisely, the estimate Now, we estimate the difference as where we used that H(Ψ n ) = λ n Ψ n and the fact that H is self-adjoint. Therefore, we deduce that since the eigenvalues are in increasing order, leading to the desired estimate (3.5). A similar argument A direct application of the strategy above to which satisfy H k (h) ∈ L 2 (R d ) for all k > 0, together with the symmetry of the potential in Lemma 2.2, leads to the main result of this section: estimating the error due to the truncation in (3.4).
Theorem 3.2 Given D N and K N the truncated diffusion and drift coefficients defined by then the following error estimate holds: for all k > 0 and all N ∈ N, there exists C k > 0 (depending on k but not on N) such that

Numerical approximation and simulations
The numerical method in practice works as follows: • Step 1. H → H R : we consider the problem set on [−R, +R] d , with R 1 and completed with homogeneous Dirichlet boundary conditions; owing to the functional framework, we expect that the eigenfunctions are localized around the origin and are exponentially decreasing far away from the wells of the potential; see Agmon (1985). In particular, it holds under our assumptions on the behaviour of the potential as |v| → +∞, from Hislop & Sigal (1996, Theorems 3.4, 3.10). Then, we choose R large enough to reduce the truncation error. In our examples, we fix R = 10.
• Step 2. H R → H R,μ : we use finite element methods to discretize the operator H R . In practice, we use the library Mélina (Martin, 2010), a uniform mesh on [−R, R] with 1000 uniform P 10 elements and quadrature of degree 21. We have chosen this method since we need an approximation of the solution of the partial differential equation with high-order accuracy; see Remark 4.2.
Here and in what follows we denote by μ > 0 a measure of the accuracy of the underlying discretization method. It thus contains information on both the refinement of the mesh, and the degree of the piecewise reconstruction.
• Step 3. λ n , Ψ n → λ R,μ n , Ψ R,μ n : having at hand the discrete operator on the truncated domain, denoted H R,μ , we determine its first N eigenelements λ In the finite elements framework, the eigenvectors Ψ R,μ n are piecewise polynomials functions approximating the nth eigenfunctions. We recall that the eigenvectors form an orthonormal family.
• Step 4. (η n , ω n ) → η R,μ n , ω R,μ n : given the data we compute the corresponding N Fourier coefficients by using an appropriate quadrature formula, depending on the approximation framework, for the discrete analogue (η R,μ n , ω R,μ n ) of Our results here are computed using a simple composite rectangular rule.
• Step 5. (D, K) → D R,μ,N , K R,μ,N : we now approximate the diffusion and drift coefficients using (3.6) to conclude (4.1) Once the eigenelements are known, the computational cost of the evaluation of the coefficients η R,μ n , ω R,μ n is linear with respect to the size of the linear problem to be solved (that depends directly on μ). Hence, the main source of the computational cost relies on the determination of the N eigenpairs.
For the potentials that we consider in this paper, Lanczos-like algorithms can be used. As iterative methods, their computational cost cannot be estimated a priori. Nevertheless, we expect that only a few eigenpairs are sufficient to provide us with an accurate result (for the quadratic case, the problem is exactly solved with the first eigenpair associated with a positive eigenvalue), so that the resolution would be far less costly than solving the linear system, a problem that, for small μ's, would also require iterative methods. Here we use standard Lanczos techniques; we refer the reader to Lehoucq & Sorensen (1996), Sorensen (2002) and Saad (2011) for further information on these methods and to Bonnaillie-Noël et al. (2007) for the computation of the first few eigenpairs of complicated Schrödinger operators, based on finite element approximations.
We will show numerical simulations for three different potentials in one dimension given by • case A: the symmetric smooth potential given by W (v) = 1 4γ v 4 − 1 2 v 2 with γ > 0; • case B: the symmetric singular potential given by W (v) = 1 4γ v 4 − 1 3 |v| 3 with γ > 0; • case C: the tilted smooth potential given by W (v) = 1 4γ v 4 − 1 2 v 2 − δv with γ , δ > 0. We can gather all of them in a single potential: Note that in case C (or σ = 0, δ > 0 in (4.2)), the potential W is not symmetric and the compatibility condition for solving the auxiliary equation (2.8a) is not satisfied. We shall see how the theory can be adapted to this case (see Section 4.3).
Remark 4.1 Note that in the one-dimensional case, the drift coefficient can be expressed in a simpler form. This is due to the fact that we can solve explicitly the one-dimensional Poisson equation, up to quadratures (Pavliotis & Stuart, 2008, Section 13.6). Indeed, according to (2.12) and the expression of the operator Q in (2.7), we note that Then, direct integration yields Therefore the drift coefficient defined in (2.11) becomes In cases A and B, we have R vM (v) dv = 0 by symmetry and K is given by This explicit formula will be used to check the accuracy of the method.
In order to reduce the number of free parameters, we rescale the velocity by defining v = √ γṽ into the Fokker-Planck operator in (2.7), to get where we have dropped the tildes in the velocity variable for notational simplicity, with the rescaled potential In this way, θ → 0 and γ → ∞ play the same role. In fact, for the symmetric smooth potential of case A (σ = 0), all terms involving γ disappear in the rescaled potential and we can remove one parameter by setting θ = 1. In our simulations, we consider the Schrödinger operator with the potential in each of the different cases above and denote by λ k (γ ) the kth positive eigenvalue. The first eigenvalue for H is simple, equal to 0. In cases A and B, the potentials are spherically symmetric. Then the second eigenvalue λ 1 (γ ) tends exponentially to 0 as γ tends to ∞ (Helffer & Sjöstrand, 1984;Simon, 1984, Theorem 1.5;Robert, 1987). In fact, a careful reading of these references gives that λ 1 (γ ) e −cγ for some positive constant c and λ 2 (γ ) O(1). This is a manifestation of the tunnelling effect. This behaviour leads to numerical difficulties. Indeed, we have to capture two simple eigenvalues (0, λ 1 (γ )) but with λ 1 (γ ) → 0 exponentially fast as γ → ∞. The first eigenfunction is symmetric and the second one antisymmetric. Numerically, when γ is very large, the gap between 0 and λ 1 (γ ) becomes negligible compared with the order of the accuracy of the method or even compared with machine precision. Then numerically it appears as if the problem has a double eigenvalue. Then the computation breaks down. Similar difficulties appear for the magnetic tunnelling effect; see Bonnaillie-Noël et al. (2014). A good way to determine whether or not the computation is accurate is to look at the eigenfunction: as soon as the symmetry is broken for the first two eigenfunctions, the computation is wrong. Remark 4.2 We have also tested the method by using the standard finite difference discretization. We roughly obtain similar results for small values of γ and with the same number of numerical unknowns as for the finite element algorithm (which means with a very refined grid for the finite discretization method). Discrepancies appear as γ increases: the loss of symmetry of the eigenfunctions is sensitive earlier. This is reminiscent of the well-known fact for similar problems, that increasing the degree of polynomials involved in the approximation (p-extension) is more efficient than refining the mesh (h-extension), see Ainsworth (2004) and Bonnaillie-Noël et al. (2007).

Case A
In Fig. 1 we present the first two eigenfunctions for γ = 1, 10, 50, 100, 120 on [−R/2, R/2]. We can observe the localization of the eigenfunction and the exponential decay far away from the wells of the potential. Looking at the symmetry of the eigenfunction, we see that the computation becomes problematic for γ > 100 (the eigenfunctions for γ = 120 are neither symmetric nor antisymmetric). In Fig. 2, we plot the convergence of the positive eigenvalues: γ → λ j (γ ) and log γ → log 10 λ j (γ ). We clearly observe the tunnelling effect, with the first eigenvalue λ 1 (γ ) converging exponentially fast to 0 while the higher eigenvalues remain of order 1.
In Fig. 3, we present the behaviour of the diffusion and drift coefficients as functions of γ . The first two columns give the result of our algorithm Steps 1-5. In the last column, we compare the numerical drift coefficient with the value K * given by formula (4.3). More precisely, we display the relative error V. BONNAILLIE-NOËL ET AL. |K * (γ ) − K(γ )|/K * (γ ) as a function of γ . It validates the accuracy of the algorithm. As expected from the results on the exponential decay of the second eigenvalue, the existence of a spectral gap and the formulas for D and K, the diffusion and the drift coefficients grow exponentially fast as γ → ∞. Figure 4 presents the convergence of the algorithm with respect to the number of eigenmodes for γ = 1, 10, 50. Using relations (4.1), we represent where we use the shorthand notation D(N) = D R,μ,N for fixed R and μ respectively, κ(N) = κ R,μ,N and D * is a reference value for the diffusion coefficient obtained for N = 50 and K * is as above a numerical approximation of (4.3) with a composite rectangular rule. We observe that only 10 eigenmodes are enough to calculate accurately the diffusion coefficient and 15 eigenmodes for the drift coefficient. We can conclude that our numerical method leads to an efficient and accurate calculation of the drift and diffusion coefficients.

Case B-nonsmooth potential
For the case B, we recall that   We show similar computations to those presented in the previous subsection. Figure 5 shows that the computations are no longer accurate beyond γ 7: the numerical eigenfunctions have lost their symmetry. For the nonsmooth potential, we cannot take a larger value of γ .
In Fig. 6, we present the convergence of the positive eigenvalues. As for the previous potential, we observe the exponential decay of λ 1 (γ ) to 0 as γ increases. For large values of γ we note a change of slope in the log-log graph and the loss of monotonicity for the eigenvalues λ 3 , λ 4 , λ 5 . A further investigation would be necessary to decide whether this behaviour is a numerical artifact or whether it is the actual behaviour of the eigenvalues.  In Fig. 7, we display the behaviour of the diffusion and drift coefficients as functions of γ showing the exponential growth as γ → ∞. The comparison with (4.3) justifies the quality of the approximation. In Fig. 8, we present the convergence of the drift and diffusion coefficients with respect to the number of modes for γ = 1 and γ = 5: as in case A, the coefficients are well captured with a few eigenmodes.

Case C-smooth potential with a linear drift
Now we present our computations for the nonsymmetric potential The compatibility condition, necessary so that we can apply Fredholm's alternative, does not hold: here However, we can adapt the asymptotics in order to handle this situation where the flux of the equilibrium state does not vanish. Starting from (2.6), we set with f ε the solution of (2.6). We check that withρ ε = f ε dv. We perform the Hilbert expansion onf ε as in Section 2: we still have Q f (0) = 0 and thusf (0) (t, r, v) = ρ(t, r)M (v) at leading order, whilef (1) is defined by inverting Therefore, equations (2.8) become The compatibility condition is now satisfied and we getf (1) = χ · ∇ r ρ + ψ · (∇ r U ρ)ρ. Considering the ε 0 terms in the expansion, we finally arrive at where the diffusion and drift matrices are given by the following analog of (2.11):   . Diffusion and drift coefficients as functions of δ for the nonsymmetric potential (in Case C), with γ = 1, 5, 10, 25 (with legend as in Fig. 10).
The analysis of this case can be performed as in the appendix; see Goudon & Mellet (2003) for a similar problem. In Fig. 9, we present the first two eigenfunctions for γ = 1, 10 and δ = 1, 5, 10. We observe that the eigenfunctions are no longer symmetric. In Fig. 10, we plot the first positive eigenvalue λ 1 (γ , δ) as a function of δ for γ = 1, 5, 10, 25. As in the previous cases, Fig. 11 gives the approximation of the diffusion and drift coefficients using our algorithm and the comparison with (4.3). We observe, even for this nonsymmetric case, the good approximation of the drift coefficient. Figure 12 illustrates the convergence of the algorithm for γ = 1 and δ = 1, 5, 10. As expected, we need more eigenmodes than for the smooth symmetric potentials in order to obtain an accurate approximation of the drift and diffusion coefficients.

Conclusions
A new method for calculating the drift and diffusion coefficients in the diffusion approximation for a swarming model was presented in this paper. Our method is based on the calculation of the eigenvalues and eigenfunctions of an appropriate Schrödinger operator. This operator was obtained after a unitary transformation of the Markov generator that appears in the Poisson equation arising in the definition of the coefficients of the limiting problem. The eigenvalue problem for this Schrödinger operator was solved using a high-order finite element approximation. Our numerical method was tested on a few simple potentials and the effects of a tilt and of the lack of smoothness of the potential on the drift and diffusion coefficients were investigated. We also investigated the difficulties related to the tunnelling effect that appears in the 'semiclassical' limit γ 1. We believe that the numerical method developed in this paper can be applied to the calculation of effective coefficients in a wide variety of diffusion approximations, coarse-grained and mean-field models that appear in kinetic theory or in homogenization theory. The crucial observation is that in many 1557 EFFICIENT NUMERICAL CALCULATION OF DRIFT AND DIFFUSION COEFFICIENTS different settings effective coefficients are given in terms of the solution to an appropriate linear Poisson equation (Resibois & De Leener, 1977;Bensoussan et al., 1978;Goudon & Poupaud, 2004;Pavliotis & Stuart, 2008;Komorowski et al., 2012). Thus, we believe that the spectral approach advocated in this paper can be of more general interest. Moreover, this approach is also open to relevant perspectives: (1) It would be interesting to develop a detailed analysis of rates of convergence, and a careful study of the computational cost, depending on the numerical parameters of the approximation (R, μ, etc). Depending on the underlying operators, the method might also benefit from the use of appropriate preconditioners.
(2) It would be interesting to compare the performance of the proposed numerical method with alternative techniques such as Monte Carlo methods (Pavliotis et al., 2006), expansions in orthogonal polynomials (Pavliotis & Vogiannou, 2008), etc.
(3) The method can be extended to more general types of Poisson equations (and the calculation of the corresponding effective coefficients) including hypoelliptic operators of Schrödinger type that appear in, e.g. Hairer & Pavliotis (2004, or auxiliary equations, possibly depending on the space variable, that appear in homogenization theory (Bensoussan et al., 1978;Goudon & Poupaud, 2004). holds. This is crucial to the analysis. As mentioned above, the discussion is strongly inspired by the study of the Vlasov-Poisson-Fokker-Planck system (Masmoudi & Tayeb, 2007;El Ghani & Masmoudi, 2010). In this appendix, we present a proof of Theorem 2.1. Naturally, the proof of this theorem consists of two parts: the derivation of a priori estimates presented in Section A.2 and the passage to the limit, presented in Section A.3.

A.2 A priori estimates
We start by observing that d dt f ε dv dr = 0 holds, which means that the total mass is conserved. Here and below, we always assume sup ε>0 f ε Init L 1 < ∞.
Owing to the regularity of U, we observe that Next, we compute The last term can be rewritten as |∇ v W f ε + 2∇ v f ε | 2 dv dr, bearing in mind (A.2).