Quantified overdamped limit for kinetic Vlasov-Fokker-Planck equations with singular interaction forces

We establish a quantitfied overdamped limit for kinetic Vlasov-Fokker-Planck equations with nonlocal interaction forces. We provide explicit bounds on the error between solutions of that kinetic equation and the limiting equation, which is known under the names of aggregation-diffusion equation or McKean-Vlasov equation. Introducing an intermediate system via a coarse-graining map, we quantitatively estimate the error between the spatial densities of the Vlasov-Fokker-Planck equation and the intermediate system in the Wasserstein distance of order 2. We then derive an evolution-variational-like inequality for Wasserstein gradient flows which allows us to quantify the error between the intermediate system and the corresponding limiting equation. Our strategy only requires weak integrability of the interaction potentials, thus in particular it includes the quantified overdamped limit of the kinetic Vlasov-Poisson-Fokker-Planck system to the aggregation-diffusion equation with either repulsive electrostatic or attractive gravitational interactions.


Introduction
Our starting point is the following kinetic Vlasov-Fokker-Planck equation: is a family of probability measures on R d ×R d with d ≥ 1, and γ > 0 is the strength of the linear damping in velocity and diffusion. Here, ρ γ t = R d µ γ t (· , dv) denotes the x-marginal of µ γ t , representing the positional distribution of points, and F : R d ×P(R d ) → R d is the driving force of the system, which in our case arises from an external and/or interaction potential. More precisely, we assume F to take the form where Φ : R d → R and K : R d → R are given functions, and ∇K(x − y) ρ(dy) .
In this paper, we establish quantitative estimates for the Vlasov-Fokker-Planck equation (1.1) in the overdamped regime, i.e. in the regime when γ ≫ 1. Note that equation (1.1) also includes the classical Vlasov-Poisson-Fokker-Planck system when ∇K = ζx/|x| d , d ≥ 3. Here, the constant ζ can be chosen ζ = ±1 according to applications in either plasma physics or astrophysics.
At the formal level, when γ → ∞, it is expected that the kinetic equation (1.1) will converge to the following well-known aggregation-diffusion equation: Equations of the type (1.2) appear in various contexts, e.g. in biological pattern formation [37], semiconductor equations [38], models for granular flows and self-gravitating particles [4,14], hydrodynamic limits of vortex systems [8,9], and even in the theory of combustion [3]. The best known classical example is the Keller-Segel model [34], with K satisfying ∆K = δ 0 , which describes the collective motion of chemotaxis.
More recently, aggregation-diffusion equations have been studied in the framework of gradient flows with respect to the 2-Wasserstein distance defined below [2,13,33]. Formally, this means that we can recast (1.2) into the form where δ ρ E is the L 2 -derivative of the free energy E , which in this case takes the form Here, Ent(µ|ν) is the relative entropy of µ with respect to ν (cf. (1.11) below). We will indeed exploit the gradient flow structure associated to (1.2) to prove our results.
Formal derivation of the aggregation-diffusion equation. Let us briefly and formally explain why we would expect equation (1.2) to arise from (1.1) as γ → ∞. We begin by writing (1.1) as . Since the right hand side of the above equation can be rewritten as with the standard d-dimensional normal distribution (or Maxwellian) N d (v) = 1 (2π) d/2 e −|v| 2 /2 , and it has the order γ 2 , we infer that On the other hand, if we set m γ t := R d vµ γ t (·, dv), we then find from (1.1) that (1.5) We now use the approximation (1.4) to obtain and this shows that the momentum m γ t roughly satisfies We use again the fact γ ≫ 1 to neglect the lowest-order term in the above equation, and obtain: Inserting this into the continuity equation (1.5) yields which is our limiting equation (1.2).
Motivation of (1.1) from interacting particle systems. The kinetic Vlasov-Fokker-Planck equation (1.1) is closely related to classical Newton dynamics for point particles interacting through the force field F in the presence of noise. More precisely, under suitable assumptions on Φ and K, (1.1) can be derived from the following system of stochastic differential equations [48]: where X γ,i t ∈ R d and V γ,i t ∈ R d , i = 1, . . . , N , denote the position and velocity of i-th particle at time t > 0, respectively, and B 1 t , . . . , B N t are N independent d-dimensional Brownian motions. By considering the mean-field limit N → ∞, the so-called nonlinear McKean process    dX γ t = γ V γ t dt , 1 γ dV γ t = F(X γ t , ρ t ) dt − γV γ t dt + √ 2 dB t , (1.7) replaces the system (1.6), where ρ γ Novelty of the current work. The study of the overdamped limit has been of interest since the seminal work of Kramers [35]. Kramers formally discussed convergence results for the kinetic Fokker-Planck equation (K = 0), now known as Smoluchowski-Kramers limit, by introducing a coarse-graining map-we consider the same coarse-graining map in this work (see (1.8)

below).
In [41], Nelson made Kramer's results rigorous by investigating the corresponding stochastic differential equations (1.6). He showed that a suitable rescaling of the solution to the Langevin equation converges almost surely to the solution of (1.2) with K = 0. Since then, several related results for the case K = 0 have been proven, by either using stochastic and asymptotic techniques [25,28,31,40], or by means of variational techniques [20]. The work [20] is particularly interesting, since it provides both quantitative estimates in the relative entropy and the 2-Wasserstein distance, where the latter is deduced from the former by means of the Talagrand and log-Sobolev inequalities.
In the presence of an interaction potential, a variational technique is proposed in [21] to study the rigorous limit from (1.1) to the equation (1.2) under suitable regularity and integrability assumptions on the interaction potential K, for instance K ∈ C 2 (R d ) ∩ W 1,1 (R d ) having globally bounded first and second derivatives. This strategy is based on duality methods and compactness arguments via a coarse-graining map, which yields the overdamped limit without convergence rates. In [23,27,43], the qualitative analysis of overdamped limit from the Vlasov-Poisson-Fokker-Planck system towards the drift-diffusion equation in the absence of the potential Φ is stuided, i.e. from (1.1) to (1.2) with K given as attractive or repulsive Coulomb potential and Φ ≡ 0. These works are mainly based on the uniform-in-γ entropy estimate combined with the compactness arguments, and thus it is not clear how it can be extended or refined to have quantitative error estimates between solutions to (1.1) and (1.2).
To the authors' best knowledge, a quantified overdamped limit for the Vlasov-Fokker-Planck equation with smooth nonlocal forces has only been established very recently in [45], while the case with singular nonlocal forces remains open. The main purpose of this study is, therefore, to develop a new methodology for the quantified overdamped limit from from (1.1) to (1.2) that works even with irregular or attractive/repulsive singular interaction potentials. We would also like to emphasize that a well-posedness theory for the Vlasov-Fokker-Planck (1.1) with a required energy dissipation in a weight Sobolev space is established-this is necessary due to the lack of a well-posedness theory in the literature for the case when the external potential Φ is unbounded.
The aggregation-diffusion equation (1.2) can also be rigorously derived from compressible Euler-type equations via large friction limits [12,17,18,26]. Despite these fruitful studies on the overdamped limit at the hydrodynamic level, there are only a few works on the rigorous limit from the kinetic equation to the associated continuity-type equation. In the absence of diffusion, the large friction limit of a Vlasov equation with nonlocal forces is rigorously studied in [32], and later in [24], where the assumptions on the potentials used in [32] are relaxed. More recently, a quantified large friction limit of a Vlasov-type equation was established in [10], using a pressureless Euler system as an intermediate system.
Notation. Throughout the paper, we adopt a slight abuse of notation by identifying the probability measure with its Lebesgue density, i.e. µ t (dxdv) = µ t (x, v) L 2d (dxdv) for simplicity of presentation. Here, L m , m ∈ N denotes the Lebesgue measure on R m . By C ∞ c (R m ), we mean the space of smooth functions with compact support; C b (R m ) the space of continuous and bounded functions; and C(I; X) the space of continuous maps from an interval I ⊂ R to a metric space X. By B r (x) ⊂ R d we mean the ball of radius r > 0 around x ∈ R d . For notational simplicity, we set B r := B r (0).
Outline of methodology. In the current work, we provide a rigorous limit from the kinetic equation (1.1) to the aggregation-diffusion equation (1.2) with explicit bounds in the 2-Wasserstein distance defined below. Our methodology hinges upon three main ingredients: (i) An auxiliary continuity equation induced by a coarse-graining map; (ii) A stability estimate in the 2-Wasserstein distance for solutions to the continuity equation with different velocity fields; (iii) A priori γ-independent estimates for solutions to (1.1) in a weighted Sobolev space.
We would like to emphasize that this strategy enables us to take into account the attractive or repulsive Coulomb interaction potential K, i.e. with K satisfying ±∆K = δ 0 . Thus, our result includes the overdamped limit from Vlasov-Poisson-Fokker-Planck system to the classical Keller-Segel equation. In fact, for purely repulsive interactions, we can go beyond the Newtonian singularity, i.e. Riesz potentials K can be considered, see Remark 1.8 below for details.
Let us give a brief outline of our strategy. By introducing the map we find that the coarse-grained quantityρ γ := (π x • Γ γ ) # µ γ satisfies an equation of the form (1.9) Here Notice that the coarse-grained system resembles the aggregation-diffusion equation (1.2)-the essential difference being the drift term-and plays the role of an intermediate system between the kinetic equation (1.1) and the limiting equation (1.2). Therefore, the convergence ofρ γ towards ρ is related to the stability of solutions to Fokker-Planck equations with different drifts.
For the quantitative error estimates of solutions, we use the 2-Wasserstein distance, where for any p ∈ [1, ∞), the p-Wasserstein distance is defined by for any Borel probability measures µ and ν on R m , m ∈ N, where Π(µ, ν) is the set of all probability measures on R m × R m with first and second marginals µ and ν, respectively, i.e. for all ϕ , ψ ∈ C b (R m ): We denote by P 2 (R m ) the set of probability measures in R m with bounded second moment. Then P 2 (R m ) is a complete metric space endowed with the 2-Wassertein distance d 2 (see [2] or [51] for a more detailed discussion). By employing the 2-Wasserstein distance, we first estimate the d 2 -distance between ρ γ and ρ γ . By our definition ofρ γ , we easily obtain the quantitative bound of terms of γ and the kinetic energy of µ γ . Indeed, we show in Lemma 2.6 below that The main difficulty arises from the error estimate betweenρ γ and ρ, the solution to the limiting equation (1.2). For this, we introduce in Section 3.1 a modulated interaction energy and derive in Section 3.2 a stability estimate of the form , for some λ > 0, where e γ t , explicitly given in Lemma 3.9, is related to the error between the drift fields in (1.2) and (1.9).
We notice that the modulated interaction energy D K has a positive sign when the interaction potential K gives rise to positive definite kernels-purely repulsive interactions fall within this class of potentials. Thus for the repulsive interaction potentials, we only need to estimate the error term e γ t L 2 (ρ γ t ) . On the other hand, in the presence of an attractive interaction potential, D K may take negative values, and hence we require |D K (µ, ν)| to be controlled by d 2 2 (µ, ν). This can be achieved for appropriate interaction potentials K as discussed in Section 3.1 (see Theorem 3.3).
In the regular interaction potential case, we do not require higher-order regularity of solutions µ γ to (1.1). However, for taking into account singular interaction potentials, we will need additional regularity of the force field F. This requires some uniform-in-γ estimates of solutions µ γ . For this, we introduce a Sobolev space W k,p H weighted by the exponential of the Hamiltonian H(x, v) := Φ(x) + |v| 2 /2 (cf. Section 1.1). Our careful analysis of solutions µ γ in this space makes the interaction term ∇K ⋆ ρ γ Lipschitz continuous and bounded, uniformly in γ, even if the potential K has a strong singularity at the origin. For details, see Section 5.
Unfortunately, we were unable to apply the above strategy to deduce stability estimates in the 2-Wasserstein distance when ∇K is only assumed to be bounded. The reason lies in our inability to relate the modulated interaction energy D K to d 2 2 . Nevertheless, we are able to obtain stability estimates with respect to the bounded Lipschitz distance d BL by other means. This is addressed in Section 7.
1.1. Main assumption and results. In this part, we list our main results depending on the regularity of the interaction potential K in the force field F.
Concerning the potential function Φ, we assume throughout this paper that 0 ≤ Φ ∈ Lip loc (R d ), and that Here, we denoted by Lip(R d ) the set of Lipschitz functions on R d and Note that the assumptions (A Φ ) above allow us to consider both bounded and unbounded external potentials with at most quadratic growth at infinity, i.e. we can consider the confinement potential Φ(x) = |x| 2 /2.
is called a weak solution to the Cauchy problem for (1.1), with initial datum µ 0 ∈ P 2 (R d × R d ), if µ| t=0 = µ 0 , and (1.1) holds in the following sense: for every s, t ∈ (0, T ), and every ϕ ∈ Similarly, we call a continuous curve ρ ∈ C([0, T ]; P 2 (R d )) a weak solution to the Cauchy problem (1.2), with initial datum ρ 0 ∈ P 2 (R d ), if ρ| t=0 = ρ 0 , and (1.2) holds in the following sense: for every s, t ∈ (0, T ), and every ϕ ∈ C ∞ c (R d ), As for the limiting equation, we consider an additional notion of solution, which allows us to make use of the Wasserstein (or Otto) gradient flow structure of (1.2). Definition 1.2. Let E be the free energy given in (1.3). We call a continuous curve ρ ∈ C([0, T ]; P 2 (R d )) an E -regular solution to the Cauchy problem (1.2), with initial datum ρ 0 ∈ P 2 (R d ) ∩ dom(E ), if it is a weak solution in the sense of Definition 1.1, satisfying additionally ρ ∈ L 1 ((0, T ); W 1,1 (R d )), and For any nonnegative Radon measure ν and probability measure µ in R m , m ∈ N, we set Our first main result concerns the case ∇K ∈ L ∞ (R d ) ∩ Lip(R d ). Under this regularity assumption on K, the global-in-time well-posedness of weak solutions µ γ to (1.1), for each γ > 0, follows immediately from the standard existence theory for the corresponding nonlinear McKean process (1.7) [5,48], and similarly for the weak solution ρ to (1.2) [13]. Theorem 1.3. Let T > 0 and assume that the interaction potential K satisfies ∇K ∈ L ∞ (R d ) ∩ Lip(R d ). Moreover, let the family of initial data (µ γ Then the associated family of unique weak solutions is the unique E -regular solution of (1.2) in the sense of Definition 1.2, with initial datum ρ 0 ∈ P 2 (R d ) ∩ dom(E ), and ρ γ t = π x # µ γ t is the first marginal of µ t , then for some constant C > 0 independent of γ ≥ 1.
The first part of Theorem 1.3 provides uniform-in-time γ-independent moment estimates for the family of solutions (µ γ t ) γ≥1 obtained from standard existence theory mentioned above, which is essential for establishing the second part of the statement-the constant C > 0 in the theorem depends explicitly on the moment bounds (see Section 4).
Next, we consider singular interaction potentials, e.g. ∇K ∈ L q loc (R d ). As mentioned above, we require more regular solutions µ γ of the kinetic equation (1.1). For p ∈ [1, ∞), we denote by L p H := L p H (R d × R d ) the space of measurable functions whose p-th powers are weighted by the exponential of the Hamiltonian H(x, v) = Φ(x) + |v| 2 /2, and equipped with the norm We also denote by W k,p x := W k,p x,0 , i.e.
We now state theorems that hold for singular interaction forces.
Theorem 1.4. Let T > 0 and suppose that the interaction potential K satisfies Let the family of initial data (µ γ There exists a positive time T p ∈ (0, T ], and a family of unique weak solutions We note that Theorem 1.4 includes the Vlasov-Poisson-Fokker-Planck system (i.e. (1.1) with Φ ≡ 0) with either repulsive electrostatic or attractive gravitational interactions: We refer to [15,44,49] for the global-in-time existence of weak solutions for the Vlasov-Poisson-Fokker-Planck system. For interaction potentials that are less singular than Coulomb (see (1.13) below), we refer to [11] for the existence of weak solutions.
On the other hand, our result applies to the case with possibly unbounded external potentials Φ, and provides bounds in the weighted Sobolev space W 1,p x,H . More importantly, these bounds are uniform in γ ≥ 1, which plays an essential role in our quantified error estimate. In a nutshell, the existence and uniqueness result follows from a truncation procedure, for which the truncated equation admits unique global-in-time smooth solutions (cf. [6,19,50]), and by a careful analysis when removing the truncation, while maintaining the γ-independent bounds.
The following result summarizes the quantified error estimate with respect to the 2-Wasserstein distance for various types of singular interaction potentials K. Theorem 1.5. Let the assumptions of Theorem 1.4 be satisfied with p > max{d, q/(q − 1)} and T p > 0 be as given in Theorem 1.4, and (µ γ ) γ≥1 be the family of weak solutions to (1.1) given by that theorem.
Suppose that the interaction potential K additionally satisfies one of the following conditions: , or (iii) (Attractive Newtonian) K is given by the Newtonian potential, i.e. ∆K = δ 0 , where δ 0 denotes the Dirac measure on R d giving unit mass to the origin.
for some constant C > 0 independent of γ ≥ 1.
Several remarks regarding Theorem 1.5 are in order: Remark 1.6. The local-in-time existence of E -regular solutions ρ to (1.2) can be deduced from the uniform-in-γ estimates obtained in Theorem 1.4 via compactness arguments. The uniqueness of such solutions for singular interactions K is, however, a subtle business (see e.g. [16] for the Newtonian potential). Since the well-posedness of the limiting equation is not the main focus of our work, we refer the reader to [13]. That being said, in the absence of uniqueness, the estimate in Theorem 1.5 still holds true, though only for a subsequence of γ ≥ 1.
Remark 1.7. The assumption p > d in Theorem 1.5 is only necessary to obtain uniform L ∞bounds via Sobolev embeddings required for relating the modulated interaction energy D K with the 2-Wasserstein distance d 2 2 , and may be removed if one knows a priori that γ-independent L ∞ -bounds hold true.
Remark 1.8. The assumption on K in Theorem 1.5 allows us to consider the Riesz potentials, for the purely repulsive cases, and for the purely attractive cases (α = d−2 being the Coulomb potential). Indeed, for the repulsive case, we only need to check the condition (1.12). Since we can always find q ∈ (1, ∞] such that (α + 1)q < d when α < d − 1, we have for some R > 0.
It is readily seen that ∇K ∈ W 1,∞ (R d \ B R ). For the attractive case, we need to verify the conditions (1.12) and Theorem 1.5 (i) or (iii). Similarly as before, condition (1.12) imposes for some constant C > 0, where the integrability holds since α + 2 < d, thus validating condition Theorem 1.5 (i). In particular, we can consider interaction potentials K that are made up of combinations of repulsive and attractive potentials, satisfying The final result in the present work is a result concerning the case when ∇K is bounded and does not possess any additional regularity. The stability estimate holds in the bounded Lipschitz norm defined, for any finite measure µ and ν on R d , by for some constant C > 0 independent of γ ≥ 1.
Outline of the paper. The rest of this paper is organized as follows. In Section 2, we introduce the Fokker-Planck-type intermediate system via the coarse-graining map, and provide some quantitative bounds and regularity estimates of solutions to that system related to the kinetic equation (1.1). In particular, we present the uniform-in-γ estimate of second moment of µ γ t in velocity. Section 3 is devoted to the derivation of the stability estimate (1.10). We then consider the regular case, i.e. when ∇K is bounded and Lipschitz, in Section 4, and show the quantitative error estimate between ρ γ t and ρ t in the 2-Wasserstein distance, resulting in Theorem 1.3. For the singular interaction potentials, we will need a higher-order regularity of solutions. Section 5 provides the well-posedness result for the kinetic equation (1.1) and uniform-in-γ estimates of solutions in the weighted Sobolev space W 1,p x,H (cf. Theorem 1.4). In Section 6, we combine the stability estimate (1.10) and the uniform-in-γ estimates of solutions obtained in Section 5 to quantify the overdamped limit for the kinetic equation (1.1) with singular interaction potentials (cf. Theorem 1.5). In Section 7, we provide a different approach to deal with the case when ∇K is only assumed to be bounded, resulting in the proof of Theorem 1.9. Finally, we discuss the solvability for our kinetic equation (1.1), which makes all the a priori estimates for the overdamped limit presented in Section 5 completely rigorous, in Appendix A.
2. An intermediate system 2.1. The coarse-graining map. We introduce a transformation of the system via the following coarse-graining map. For each γ > 0, we set In the following lemma, we show the relative entropy estimate of µ γ with respect to Γ −γ # N 2d , whose proof is postponed to Appendix A.1 (cf. Theorem A.4). The lemma provides, in particular, a uniform in γ bound on the second moment of µ γ in velocity, see Remark 2.3 below.
Lemma 2.2. Let µ γ be a solution of (1.1) with sufficient regularity and integrability. Then for γ ≥ 1, there exist constants α, β, λ > 0, independent of γ ≥ 1, such that then the estimate provided in Lemma 2.2 yields the following uniform-in-γ bound: To justify this claim, we make use of the Donsker-Varadhan dual variational characterization for the relative entropy [22,Proposition 1.4.2], which states that for any bounded measurable function ϕ : R 2d → R, and any ϑ ∈ P(R 2d ), This characterization can be generalized to also hold for any nonnegative measurable functions ϕ : R 2d → [0, ∞] (see [29,Lemma B.1]). Applying the above characterization to ϑ = Γ −γ # N 2d and ϕ(x, v) = (|x + v/γ| 2 + |v| 2 )/4, we deduce a moment estimate for any ν ∈ P(R 2d ) with In a similar fashion, one deduces that Using the identities we deduce the following equalities: Due to the non-increasing nature of relative entropy under push-forward and the joint convexity of the functional (µ, ν) → Ent(µ|ν), we can apply Jensen's inequality to obtain the estimate with M (µ γ 0 , t) as given in Lemma 2.2. Remark 2.5. In a similar fashion as in the proof of Lemma 2.2 (cf. Theorem A.4), one can establish the relative entropy estimate for appropriate constantsβ,λ > 0 independent of γ ≥ 1.
In the following two lemmas, we provide the quantitative bound on the error betweenρ γ and ρ γ in 2-Wasserstein distance and regularity estimates onρ γ and ρ γ . Lemma 2.6 (Relationship betweenρ γ and ρ γ ). Let (ρ γ t ) and (ρ γ t ) be given as above. Then the following estimate holds Proof. We begin by relating Γ γ # µ γ t with µ γ t for each t ≥ 0. By choosing a coupling of the form Now let Π t be an optimal coupling between Γ γ # µ γ t and µ γ t for each t ≥ 0. Then (π x × π x ) # Π t is a coupling betweenρ γ t and ρ γ t . Thus we get which, together with the previous estimate yields the assertion.
Proof. We only prove the statement forρ γ -the argument for ρ γ holds analogously. For We then bound the right-hand side as for each i = 1, . . . , d. As above, we deduce ∇ρ γ ∈ L p (R d ), and consequently,ρ γ ∈ W 1,p (R d ).

Modulated interaction energy and Wasserstein stability estimates
This section contains two essential parts to our strategy: The first part concerns the so-called modulated interaction energy between two measures µ, ν for a given interaction kernel K. This energy appears when considering the differences between the free energies E (µ) and E (ν), which needs to be controlled in terms of the 2-Wasserstein distance d 2 (µ, ν). The second part of this section discusses the temporal derivative of d 2 2 (µ t , ν) along E -regular solutions t → µ t of (1.2) and an arbitrary measure ν. The temporal derivate can be related to differences in free energies, and ultimately to d 2 2 (µ t , ν).

Modulated interaction energy.
For measures µ, ν ∈ P(R d ) we consider the modulated interaction energy and further recall theḢ −1 homogeneous Sobolev norm given by for any signed measure ω on R d . We remark that this norm is only finite for signed measures having zero total mass, i.e. ω(R d ) = 0.
The next proposition relates the homogeneousḢ −1 (R d ) distance and the 2-Wasserstein distance for measures that have essentially bounded Lebesgue densities [36,42].
Due to Remark 3.1 and Proposition 3.2, we can relate D K with d 2 . However, as we will observe in the next sections, when D K take negative values, we will require |D K | to be controlled by d 2 2 . This can be achieved for appropriate interaction kernels as given in the next statement. (1) Smooth interaction: If ∇K is globally Lipschitz, then (2) Weakly singular interaction: In particular, if c ∞ < ∞, with c ∞ given in Proposition 3.2, then (3) Newtonian attractive: When K is the fundamental solution of the Laplacian, i.e. ∆K = δ 0 , D K takes the alternative form (1) We begin by noticing that where π is an arbitrary coupling between µ and ν. Optimizing over all couplings of µ and ν in (3.1) yields Similarly, we find and thus, for all x ∈ R d , This together with (3.2) implies |D K (µ, ν)| ≤ ∇K Lip d 2 1 (µ, ν) ≤ ∇K Lip d 2 2 (µ, ν) due to the monotonicity in p of the p-Wasserstein distances.
(2) An application of Remark 3.1 and Proposition 3.2 yields Thus it remains to show that for some constant C > 0. Notice that where T is the optimal transport map such that T # µ = ν and T θ (y) = (1 − θ)T(y) + θy.
Taking gradient of the above gives By an application of Hölder's inequality, the inner integral of the right-hand side above can be estimated as Consequently, the second term on the right-hand side of the former equation can be further estimated Hence, by using Fubini's theorem, we obtain We finally combine this with (3.3) to conclude the desired inequality. (3) We first notice that D K can be rewritten as Indeed, we get and applying the integration by parts to the right hand side yields (3.4). We next show that and hence, This gives the assertion (3.5). Hence we have , where the last inequality follows from Proposition 3.2.
(ii) As mentioned in the introduction, it is possible to consider linear combinations of interaction potentials satisfying each of the conditions above. For instance, one could consider an interaction kernels satisfying for some R > 0, we can combine Theorem 3.3(1) and 3.3(2) to obtain for an appropriate constant λ ∈ R.
In the case when λ ∈ R is independent of ρ and ν, the EVI may be used to deduce the uniqueness of gradient flow solutions. This follows from [2, Lemma 4.3.4] (cf. [2, Theorem 11.1.4] for the uniqueness argument; see also [16]).
which, by application of the Grönwall's lemma, yields the stability estimate In the following, we establish a similar system of inequalities for the solution to our intermediate system (2.2), which we recall here for convenience: We begin by recalling a known result for the temporal derivative of the 2-Wasserstein distance along the continuity equation.
To deduce the EVI for the evolutions we consider, we will need to compute the subdifferential of E with respect to the 2-Wasserstein distance. While the following proposition is simply a combination of results known in the literature, we include the proof for the sake of completeness.
and let T : R d → R d be the optimal transport map between ρ and ν, i.e. T # ρ = ν. Then Proof. Due to the convexity of θ → Ent(T θ # ρ|L d ) [39], we have that In particular, the inequality holds true also in the limit θ → 0, which yields

Together, this yields the inequality
As for the other terms, we simply perform basic algebraic manipulations to obtain Summing each of the contributions concludes the proof.
By combining Propositions 3.5 and 3.6, we easily obtain the following result, which ultimately allows us to deduce our main stability estimate (1.10).
Due to Theorem 3.7, we find that (ρ γ t ), solving the intermediate problem (2.2), satisfies the evolution-variational-like inequality (3.9) with Consequently, we deduce our main stability estimate from the first equality in (3.7): L 2 (ρt) dt < ∞, (ρ γ t ) is required to satisfy the same regularity assumptions for E -regular solutions (cf. Definition 1.2) for Theorem 3.7 to apply.
Finally, we provide the integrability estimate of e γ t . Lemma 3.9. The following holds.

Stability for bounded and Lipschitz interaction forces
We finally arrive at our first quantitative estimate that applies to globally Lipschitz continuous and bounded interaction forces, i.e., Under this assumption, one easily verifies that (i) for any µ ∈ P(R d ), (ii) and for any µ, ν ∈ P 2 (R d ) and every x, y ∈ R d , where we can easily obtain |(∇Φ)(x) − (∇Φ)(y)| ≤ ∇Φ Lip |x − y|. On the other hand, if we set T as an optimal map giving ν = T # µ, then Combining the above estimates yields (4.1b).
Proof of Theorem 1.3. Due to the assumed regularity on ∇K, one easily establishes the global Lipschitz continuity of x → F(x, ρ t ). The existence of weak solutions µ γ ∈ C([0, T ]; P 2 (R d × R d )) to (1.1) can then be established by, e.g. standard theory of stochastic differential equations applied to the nonlinear McKean equation (1.7). Under the assumption on the initial condition µ γ 0 , Lemma 2.2 together with Remark 2.3 applies, thereby giving the first part of the theorem. As for the second part, we begin by providing an estimate for e γ t L 2 (ρ γ t ) . From the estimate (4.1b), and Lemmas 3.9 and 2.6, we estimate Along with Theorem 3.3, the stability estimate (3.10) yields An application of Grönwall's lemma then gives Making use of the triangle inequality for d 2 and applying Lemma 2.6 twice, we obtain where the constants c, C > 0 are independent of γ ≥ 1, thereby concluding the proof.

Uniform-in-γ estimates in the weighted Sobolev space
As stated in the introduction, in order to consider the singular interaction potentials, we need better regularity of the force field F. This subsequently requires some uniform-in-γ estimates of solutions in the higher-order Sobolev space compared to the regular interaction potential case. The main purpose of this section is, thus, to establish Theorem 1.4.
However, since the proof is rather long and technical, we decided to only prove a priori estimates for sufficiently smooth solutions to the kinetic equation (1.1) in this section, and postpone to complete proof of Theorem 1.4 to Appendix A. In particular, Theorem 1.4 is a restatement of Theorem A.1.
for some constant c K,p > 0, then there exists some T p ∈ (0, T ] such that where M k > 0 is a constant independent of γ > 0.
In particular, we have the estimate To establish Proposition 5.1, we first define Here, N d represents the standard normal distribution on R d . We note that considering η is natural from the physical point of view since the free Hamiltonian for the rescaled kinetic equation (1.1) is given by and η is equal to e −H L 2d up to some constant.
for some g ∈ L p ((0, T ); L p H (R d × R d )). Further, set w := df /dη = f e H and ζ = dg/dη = ge H . Then for any p ∈ [2, ∞), Proof. From the equation for f , we find d dt

=: (I) + (II) + (III) + (IV) .
Note that We use this relation to rewrite (I) as where, in the last equality, we used the fact that In a similar fashion, we use (5.1) to rewrite (II) as which, together with (I) implies (I) + (II) = 0. As for (III), we easily find Finally, to estimate (IV), we apply the Young and Hölder inequalities to obtain Putting the terms (I)-(IV) together and using the identity we then obtain the required estimate.
In the following, we assume that µ γ t ≪ η for all γ > 0 and t ≥ 0. Lemma 5.3. Let µ γ be a solution of (1.1) with sufficient regularity. We assume that Then for any p ∈ [2, ∞), we have for u γ t := dµ γ t /dη = µ γ t e H (the density of µ γ t w.r.t. η), for all t ∈ [0, T ]. In particular, this yields Proof. The proof makes use of Lemma 5.2 with we deduce from Lemma 5.2, for almost every t ∈ (0, T ). We then conclude the proof by applying Grönwall's inequality.
Following similar computations as in Lemma 5.3, we obtain the next statement.
Lemma 5.4. Let µ γ be a solution of (1.1) with sufficient regularity. Further, let p ∈ [2, ∞) such that Then there exists t * p ∈ (0, ∞) such that In particular, for any q ∈ [2, ∞), we have Proof. As in the proof of Lemma 5.3, we apply Lemma 5.2 with . From the assumption on K, we have that ) . Therefore, an application of Lemma 5.2 provides the differential inequality d dt H , which consequently leads to the estimate Hence, for any we then obtain the required estimate. A direct application of Lemmas 5.3 and 2.7 yields the other estimates for q ∈ [2, ∞), q = p.
Remark 5.5. Let q ∈ (1, ∞) and p = max{2, q/(q − 1)}. If R > 0 exists such that then for all x ∈ R d , we obtain The first term on the right-hand side may be further estimated by where c q,2,R > 0 is a constant obtained from the standard (L q , L 2 )-interpolation inequality on the bounded domain B R . Altogether, this gives ∇K ⋆ ρ L ∞ ≤ c K,q 1 + ν L p H for some constant c K,q > 0, i.e. the assumption on K in Lemma 5.4 is satisfied. Therefore, this allows us to consider the case where ∇K is singular at the origin.
Next, we provide an estimate for the spatial derivative of µ γ t . Lemma 5.6. Let µ γ be a solution of (1.1) with sufficient regularity. We assume that and denote ξ γ t := (∇ x µ γ t ) e H for all t ∈ [0, T ]. Then for any p ∈ [2, ∞), we have In particular, this implies where C > 0 depends on c 0 , c 1 , p, and T , but independent of γ > 0.
Proof. For notational simplicity, we drop superscript γ throughout this proof. Set q k := ∂ x k µ for k = 1, . . . , d and consider ξ k := q k e H . Recall that η = e −H L 2d . We notice from (1.1) that q k satisfies with . By the assumptions placed on ∇K and ∂ x k F, we deduce that 3, the second term on the right-hand side of the previous inequality can be bounded from above by This, together with applying Grönwall's lemma yields thereby concluding the proof.
In a similar fashion as in the proof of Lemma 5.4, we obtain the following result.
The following lemma provides a class of examples of interaction potentials K that satisfies the requirements of Lemma 5.7.
Lemma 5.8. Let q ∈ (1, ∞), p = max{2, q/(q − 1)}. Suppose R ≥ 1 exists such that Then there exists a constant c ′ K,p > 0 such that the following estimate holds: Proof. By density argument, it is sufficient to show the estimate for K ∈ C 2 b (R d ). Let χ R be a smooth cut-off function satisfying 1 B R ≤ χ R ≤ 1 B 2R and |∇χ R |(x) ≤ 1/2 for all x ∈ R d . Then, for each k, j = 1, . . . , d, we have where we use the integration by parts to estimate Here q ′ = q/(q−1), and c q,2,R > 0 is a constant obtained from the standard (L q , L 2 )-interpolation inequality on the bounded domain B 2R . Since As for (II), we easily obtain Putting the estimates together and summing up in k, j = 1, . . . , d yields with the constant c ′ K,q = max{1, c q,2,R }c ′ K . Proof of Proposition 5.1. The proof follows from a combination of Lemmas 5.3, 5.4, 5.6, and 5.7.

Stability for singular interaction forces
In this section, we deal with the singular interaction potentials and provide the details of the proofs for Theorem 1.5. We begin by providing a Lipschitz-type estimate of the interaction forces which will be crucially used for the stability estimate. Lemma 6.1. Let q ∈ (1, ∞]. Suppose that there exists R ∈ (0, ∞) such that Then for any µ γ ∈ W 1,q ′ x,H (R d × R d ), q ′ = q/(q − 1), and the following estimate holds: for some constant c 0 > 0, independent of µ γ and γ > 0.
Proof. By density argument, it is sufficient to show the estimate for K ∈ C 2 b (R d ). We begin by writing In the following, we set A r = A θ,γ r (x, v, w) := B r (x+ θ(w − v)/γ) for each x, v, w ∈ R d and r ≥ 1. Now consider a smooth cut-off function χ R satisfying 1 B R ≤ χ R ≤ 1 B 2R for all x ∈ R d . Denoting u γ := µ γ e H and using the shorthandx = x + θ(w − v)/γ, we obtain for every k, j = 1, . . . , d and almost every w ∈ R d , Here, the last term on the right-hand side of the above can be bounded by As for the first and third term, we can simply use Hölder's inequality to bound them from above by We then estimate the second term as where c Φ,q > 0 is the same constant that appears in (A 2 Φ ), and we used Hence, for almost every w ∈ R d , we find that Furthermore, since for any f ∈ L r (η), r > 1, where c N ,r > 0 is a constant depending only on N d and r, and for r = 1, we then obtain for the constant c 0 > 0 independent of u γ and (x, v) ∈ R d × R d . This completes the proof.
We next estimate e γ t L 2 (ρ γ t ) in the stability estimate (3.10). Lemma 6.2. Under the same assumptions in Lemma 6.1, we have where C > 0 is independent of γ.
Proof. This follows from a direct application of Lemmas 3.9 and 6.1.
We can finally provide the details to the proof of Theorem 1.5.
Proof of Theorem 1.5. In the case (i) (respectively (iii)), Theorem 3.3(2) (respectively (3)) provides an estimate for the modulated interaction energy: for a constant c 1 > 0, independent of γ ≥ 1. Clearly, the inequality above also holds true for the case (ii), since in this case, D ≥ 0. On the other hand, Lemma 6.2 gives us a handle on the error term e γ t L 2 (ρ γ t ) . Together with the well-posedness result provided by Theorem 1.4, the general stability inequality (3.10) then gives for some constant c 1 > 0, independent of γ ≥ 1 and t ∈ [0, T p ]. Applying Grönwall's lemma yields for some constant c 2 > 0 is independent of γ ≥ 1. Finally, using a similar argument as in the proof of Theorem 1.3, we conclude that where C > 0 is independent of γ ≥ 1.

Bounded interaction forces
As mentioned in the introduction, we will employ a different approach to obtain stability estimates when ∇K is only assumed to be bounded. In any case, we will need a Lipschitz-type estimate in order to control the error term e γ t . Here, we will make an additional assumption on the external potential Φ.
Then for any , the following estimate holds: for some constant c 0 > 0, independent of µ γ and γ > 0.
Proof. The proof is very similar to that of Lemma 6.1. Again, by a density argument, it is sufficient to show the estimate for K ∈ C 2 b (R d ). Elementary computations give us the following estimate: for some constant c 0 > 0 depending only on d, σ Φ (R d ) and ∇K L ∞ . The rest of the proof follows along the same line as that of Lemma 6.1.
Remark 7.2. It is not difficult to see that the additional assumption on Φ in Lemma 7.1 is satisfied for potentials Φ satisfying Φ(x)/|x| → ∞ as |x| → ∞.
Therefore, we have that for all t ∈ [0, T ] .
An application of [20,Theorem 2.18] then yields the entropy estimate We now provide a more precise estimate of h γ =: (I) + (II) .
As for (I), we use the estimate In the former estimate, we used the well-known Csiszàr-Kullback-Pinsker inequality: for any µ, ν ∈ L 1 (R d ) .
As for the second term, we apply Lemma 7.1 to obtain for some constant c > 0 independent of γ ≥ 1. Consequently, we obtain which altogether, yields, β γ 2 , for appropriate constants λ, β > 0 independent of γ ≥ 1 and t ∈ [0, T ] . Inserting this into the entropy inequality (7.1) and applying Grönwall's inequality then yields the desired result.
Proof Theorem 1.9. We first notice that the bounded Lipschitz distance between µ and ν can be bounded either by the 1-Wasserstein distance between them or the L 1 -norm of µ − ν. This, together with the monotonicity of d p -distance and the Csiszàr-Kullback-Pinsker inequality, yields where C > 0 is independent of γ ≥ 1.
Appendix A. Well-posedness for the kinetic equation In this appendix, we prove the well-posedness of our main equation in the solution space W 1,p x,H (R d × R d ) and simultaneously justify the uniform-in-γ estimates obtained in Section 5. Let us begin by recalling our main equation: where F(·, ρ t ) = −∇Φ − ∇K ⋆ ρ t , with ρ t = π x # µ t being the x-marginal of µ t . As mentioned in the introduction, there is available literature on the existence theory for weak solutions to the equation (A.1) when the interaction potential K is bounded and Lipschitz continuous or it has the form of K(x) = 1/|x| α with α ∈ (0, d − 2]. The existence of classical solutions is also well studied for the equation (A.1) without the potential Φ when ±∆K = δ 0 . On the other hand, we have been unable to find the proper literature that provides the wellposedness for our main equation; the case K(x) = 1/|x| α with α ∈ (0, d − 2] and Φ satisfying the assumptions presented in Section 1.1. For that reason, we provide a well-posedness theory for (A.1), which is summarized in the following theorem.
As for the other terms, we apply the Young's inequality to deduce for appropriate constants c ′ 1 , c ′ 2 , c ′ 3 > 0 that are independent of γ. The moment estimates given in (A.4) by Theorem A.2 allows us to apply the Lebesgue dominated convergence to obtain Notice that the second term in Σ may be integrated by parts to obtain Indeed, for any R > 0, we find since µ r ({|x + v/γ| < R}) ≤ 1 for every r ∈ (0, t). Hence, the previous inequality and for some constants c 4 , c 5 > 0 independent of ε > 0 and γ ≥ 1, provide the upper bound with c ′ 5 = d + c 5 . As for the other terms, we apply Fatou's lemma and the Lebesgue dominated convergence to obtain Putting the inequalities together yields H(x, v) dµ r dr .
In particular, we obtain constants α, β, λ > 0 independent of ε > 0 and γ ≥ 1 such that By Grönwall's lemma, we futher obtain Therefore, integrating (A.6) and using the previous inequality yield thereby concluding the proof.
Remark A.5. Notice that in the case when ∇Φ has linear growth at infinity and ∇K is bounded and Lipschitz, the assumption on F ε in Theorem A.4 is trivially satisfied. In particular, the inequality (A.5) holds also for weak solutions µ to (A.3).
i.e. (µ ε ) ε>0 ⊂ C([0, T p ]; P(R d × R d )) is equi-continuous, which provides the existence of a curve µ ∈ C([0, T p ]; P(R d × R d )) such that [2] µ ε t → µ t narrowly in P(R d × R d ) for all t ∈ [0, T p ] . Due to the lower semicontinuity of L p -norms with respect to narrow convergence, we then deduce that µ t ∈ L p (R d × R d ) for every t ∈ [0, T p ], and consequently µ ∈ L ∞ ([0, T p ]; L p (R d × R d )). It is then possible to extend the previous convergence to the convergence (A.8).
As for the first marginal ρ ε = π x # µ ε of µ ε , we obtain a stronger convergence that is needed to pass to the limit in the bilinear term (K ⋆ ρ ε ) · ∇ϕ, µ ε . For this, we make use of a classical result that characterizes strong convergence in L p -spaces (Kolmogorov-Riesz-Fréchet theorem). For each t ∈ [0, T p ], the family of marginals (ρ ε t ) ε>0 of (µ ε t ) ε>0 is relatively compact in L 1 (R d ). Consequently, there exists a limit curve ρ ∈ C([0, T p ], P(R d )) ∩ L ∞ ([0, T p ]; L p (R d )) such that and obtain where K γ t and ∇ ν K γ t can be explicitly expressed as For any a ∈ R, we use these observations to estimate Together, this yields We then estimate the second term on the right hand side of the above inequality. Note that |F ε (ξ, ρ t )| ≤ c Φ (1 + |ξ|) + c 0 , c 0 = sup t∈[0,Tp] ∇K ε ⋆ ρ ε t L ∞ , and the second moment bound estimate of µ t in x gives (A.10) Furthermore, by Taylor's theorem, we get σ(s) ≃ s and d(s) ≃ s 3 implying σ(s) (2γd(s)) 1/2 ≃ s −1/2 .