Mean-Field Limits: From Particle Descriptions to Macroscopic Equations

We rigorously derive pressureless Euler-type equations with nonlocal dissipative terms in velocity and aggregation equations with nonlocal velocity fields from Newton-type particle descriptions of swarming models with alignment interactions. Crucially, we make use of a discrete version of a modulated kinetic energy together with the bounded Lipschitz distance for measures in order to control terms in its time derivative due to the nonlocal interactions.


Introduction
In this work, we analyse the evolution of an indistinguishable N -point particle system given bẏ subject to the initial data Here x i = x i (t) ∈ R d and v i = v i (t) ∈ R d denote the position and velocity of iparticle at time t, respectively. The coefficient γ 0 represents the strength of linear damping in velocity, ε N > 0 the strength of inertia, V : R d → R + and W : R d → R stand for the confinement and interaction potentials, respectively. ψ : R d → R + is a communication weight function. Throughout this paper, we assume that W and ψ satisfy W (x) = W (−x) and ψ(x) = ψ(−x) for x ∈ R d . They include basic particle models for collective behaviors, see [12,20,25,34,36,46,47,63] and the references therein. Our main goal is to derive the macroscopic collective models rigorously governing the evolution of the particle system (1.1) as the number of particles goes to infinity. On one hand, we will derive hydrodynamic Euler-alignment models given by (1.3) in the mean-field limit: when initial particles are close to a monokinetic distribution ρ 0 (x)δ u 0 (x) (v) in certain sense and ε N = O (1) as N → ∞. On the other hand, we will show that the particle system can be described by aggregation equations of the form (1.5) in the combined mean-field/small inertia limit when initial particles are close to a monokinetic distribution ρ 0 (x)δ u 0 (x) (v), γ > 0 and ε N → 0 as N → ∞. For simplicity of notations when dealing with the mean-field limit, we will take ε N = 1 in the sequel.

Mean-field limits: from particles to continuum
As the number of particles N tends to infinity, microscopic descriptions given by the particle system (1.1) become more and more computationally unbearable. Reducing the complexity of the system is of paramount importance in any practical application. The classical multiscale strategy in kinetic modelling is to introduce the number density function f = f (x, v, t) in phase space (x, v) ∈ R d × R d at time t ∈ R + and study the time evolution of that density function. Then at the formal level, we can derive the following Vlasov-type equation from the particle system (1.1) as N → ∞: (1.6) where ρ f = ρ f (x, t) is the local particle density and F a ( f ) = F a ( f )(x, v, t) represents a nonlocal velocity alignment force given by respectively. Let us briefly recall the reader the basic formalism leading to the kinetic equation (1.6) as the limiting system of (1.1). We first define the empirical measure μ N associated to a solution to the particle system (1.1), that is, ,v i (t)) .
As long as there exists a solution to (1.1), the empirical measure μ N satisfies (1.6) in the sense of distributions. To be more specific, for any ϕ ∈ C 1 0 (R d × R d ), we get d dt (1.7) Notice that the particle velocity can also be rewritten in terms of the empirical This implies that the right-hand side of (1.7) can also be written in terms of the empirical measure μ N as d dt This concludes that μ N is a solution to (1.6) in the sense of distributions as long as particle paths are well defined. In fact, if the interaction potential W and the communication weight function ψ in the classical Cucker-Smale alignment model are regular enough, for instance, bounded Lipschitz regularity, then the global-in-time existence of measure-valued solutions can be obtained by establishing a weak-weak stability estimate for the empirical measure, see [46,Section 5] for more details. The mean-field limit has attracted lots of attention in the last years in different settings depending on the regularity of the involved potentials V, W and communication function ψ. Different approaches to the derivation of the Vlasov-like kinetic equations with alignments/interaction terms or the aggregation equations have been taken leading to a very lively interaction between different communities of researchers in analysis and probability. We refer to [3,4,10,20,30,31,35,44,47,50,[54][55][56]64,67] for the classical references and non-Lipschitz regularity velocity fields in kinetic cases, to [48,49] for very related incompressible fluid problems, and to [7,9,16,17,37,43,45,51,52,61,63,65,66] for results with more emphasis on the singular interaction kernels both at the kinetic and the aggregation-diffusion equation cases.

Local balanced laws, the mono-kinetic ansatz, and the small inertia limit
The classical procedure in kinetic theory of deriving equations for the first 3 moments of the distribution function f leads to the standard problem of how to close the moment system since the equation for the second moment will depend on higher order moments. Suitable closure assumptions are not known so far even in cases where noise/diffusion is added to the system. However, at the formal level, we can take into account the mono-kinetic ansatz for f , as done in [18,21], leading to where ρ and u are the macroscopic density and the mean velocity of particles, that is, the first two moments of f in velocity variable It is standard to check that the strain tensor and heat flux become zero and the moment system closes becoming the pressureless Euler equations with nonlocal interaction forces (1.3): and on the support of ρ. The last equation coming from the closed equation on the evolution of the second moment is redundant but it gives a nice information about the total energy of the system. Although the monokinetic assumption is not fully rigorously justified and it does not have a direct physical motivation, it is observed by particle simulations that the derived hydrodynamic system shares some qualitative behavior with the particle system, see [12,18,[20][21][22]33]. Note that (1.3) conserves only the total mass in time in this generality. However, the total free energy is dissipated due to the linear damping and the velocity alignment force as pointed out in [19] for weak solutions of this system. The hydrodynamic system (1.9) has a rich variety of phenomena compared to the plain pressureless Euler system. This fact is due to the competition between attraction/repulsion and alignment leading to sharp thresholds for the global existence of strong solutions versus finite time blow-up and decay to equilibrium, see [13][14][15]26,63,68]. We emphasize that the additional alignment, linear damping and attraction/repulsion terms can promote the existence of global solutions depending on the intial data. We will show that these hydrodynamical solutions can be obtained directly from particle descriptions as long as they exist, so their physical relevance is dictated by the time of existence of these solutions. It is worth noticing as in [18] that the mono-kinetic ansatz for f is a measurevalued solution of the kinetic equation (1.6). More precisely, one can show that ρ(x, t)δ u(x,t) (v) is a solution to the kinetic equation (1.6) in the sense of distributions as long as (ρ, u)(x, t) is a strong solution to the hydrodynamic equations Using the continuity equation in (1.3), I 1 can be easily rewritten as By multiplying the velocity equation in (1.3) by ρ and using (∇ v ϕ)(x, u(x, t)) as a test function to the resulting equation yields Then similarly as before, we can rewrite the second and third terms on the right hand side of the equality by using the mono-kinetic ansatz (1.8). This implies Combining all of the above estimates yields d dt This shows that ρ(x, t)δ u(x,t) (v) satisfies the kinetic equation (1.6) in the sense of distributions.
Finally, we will be also dealing with the small inertia limit for both the kinetic equation (1.6) and the hydrodynamic system (1.3) combined with the mean field limit. In the small inertia asymptotic limit, we want to describe the behavior of the scaled kinetic equation (1.10) and the scaled hydrodynamic system (1.11) in the limit of small inertia ε → 0. At the formal level, the equations (1.11) will be replaced by (1.4)-(1.5) as ε → 0. The limiting nonlinearly coupled aggregation equations (1.4)-(1.5) have been recently studied in [39,40]. Several authors have studied particular choices of interactions V, W and comunication functions ψ for some of the connecting asymptotic limits from the kinetic description (1.10) with/without noise to the hydrodynamic system (1.11) in [8,11,42,57], from the hydrodynamic system (1.11) to the aggregation equation (1.4)- (1.5) in [23,59,60], and for the direct limit from the kinetic equation to the aggregation equation (1.4)-(1.5) in [8,53].

Purpose, mathematical tools and main novelties
Summarizing the main facts of the mean-field limit and the monokinetic ansatz in Sections 1.1 and 1.2, both the empirical measure μ N (t) associated to the particle system (1.1) and the monokinetic solutions ρ(x, t)δ u(x,t) (v), with (ρ, u)(x, t) satisfying the hydrodynamic equations (1.3) in the strong sense, are distributional solutions of the same kinetic equation (1.6). In order to analyse the convergence of the empirical measure μ N to ρ(x, t)δ u(x,t) (v), the goal is to establish a weakstrong stability estimate where the strong role is played by the distributional solution ρ(x, t)δ u(x,t) (v) associated to the strong solution of the hydrodynamic system (1.3). Our main goal is then to quantify the following convergence in the sense of distributions for both the mean-field and the combined meanfield/small inertia limit for well prepared initial data. Our main mathematical tools are the use of a modulated kinetic energy combined with the bounded Lipschitz distance in order to control terms between the discrete particle system and the hydrodynamic quantities. Let us first introduce the modulated kinetic energy as where f is a solution of kinetic equation (1.6) and u is the velocity field as part of the solution of the pressureless Euler equations (1.3). Modulated kinetic energies were used in conjunction with relative potential energy terms for quasineutral limits of Vlasov like equations [5,6,62] for instance. We would like to emphasize that the quantity (1.12) gives a sharper estimate compared to the classical modulated macroscopic energy. Indeed, the macro energy of the system (1.3) is given by Thus its modulated energy, also often refereed to as relative energy, can be defined as On the other hand, by Hölder inequality, we easily find that This yields In fact, we can easily show that This shows that the convergence of the modulated kinetic energy (1.12) implies the convergence of the modulated macro energy (1.13). We notice that if f is a monoki- , then the second term on the right hand side of (1.14) becomes zero, and the two modulated energies (1.12) and (1.13) coincide. For notational simplicity, we denote by the set of trajectories associated to the particle system (1.1). Then let us define the first important quantity that will allow us to quantify the distance between particles (1.1) and hydrodynamics (1.3), it is just the discrete version of the modulated kinetic energy (1.12) defined as The second quantity that will allow us our quantification goal combined with the discrete modulated energy (1.15) is a classical distance between probability measures, the bounded Lipschitz distance, used already by the pioneers in kinetic theory [4,64,67] in the early works for the mean-field limit. Notice that the pressureless Euler system (1.3) includes the nonlocal position and velocity interaction and alignment forces. Furthermore, its relative energy/entropy has no strict convexity in terms of density variable due to the lack of pressure term. In order to overcome these difficulties, ideas of combining the modulated macro energy and the first or second order Wasserstein distance have been recently proposed in [8,11,32] quantifying the hydrodynamic limit from kinetic equation to the pressureless Euler type system. More recently, in [24], a general theory providing some relation between a modulated macro energy-type function and p-Wasserstein distance is also developed. In particular, in [24,Proposition 3.1], it is discussed that the p-Wasserstein distance with p ∈ [1,2] can be controlled by the modulated macro energy functional.
In the present work, we will employ the bounded Lipschitz distance to provide stability estimates between the empirical particle density ρ N defined as with μ N t be the empirical measure associated to the particle system (1.1), and the hydrodynamic particle density ρ solution to (1.3). More precisely, let M(R d ) be the space of signed Radon measures on R d , which can be considered as nonnegative bounded linear functionals on C 0 (R d ). Let μ, ν ∈ M(R d ) be two Radon measures. Then the bounded Lipschitz distance, which is denoted by d B L : M(R d ) × M(R d ) → R + , between μ and ν is defined by where the admissible set of test functions are given by We also denote by Li p(R d ) the set of Lipschitz functions on R d . In Proposition 2.2 below, we provide a relation between the bounded Lispchitz distance and the discrete version of the modulated kinetic energy (1.15). This key observation allows us to overcome the difficulties mentioned above.

Main results and Plan of the paper
We will first assume that the particle system (1.1), the pressureless Euler-type equations (1.3), and the aggregation equations (1.4)-(1.5) have existence of smooth enough solutions up to a fixed time T > 0. We postpone further discussion at the end of this subsection, although we make precise now the assumptions needed on these solutions for our main results.
Our first main result shows the rigorous passage from Newton's equation (1.1) to pressureless Euler equations (1.3) via the mean-field limit as N → ∞.
be a solution to the particle system (1.1), and let (ρ, u) be the unique classical solution of the pressureless Euler system with nonlocal interaction forces as N → ∞. In fact, we have the following quantitative bound estimate: The main novelty of this first result resides in how to control the alignment terms via the modulated energy combined with the bounded Lipschitz distance.

Remark 1.1. (Singular repulsive interaction)
The previous result also applies to singular repulsive interaction potentials. In particular, it holds for the Coulomb interaction potential on R d given by and for Riesz potentials in a sense to be specified in Section 2.3. Here α d denotes the volume of the unit ball in R d . In order to deal with the singularity on the interaction potential, the diagonal term should be eliminated in the modulated energy functional. This has been recently solved in the recent breakthrough result in [66] by introducing a different relative potential energy avoiding the diagonal terms. The details for singular interaction potentials cases are postponed to Section 2.3, see Theorem 2.1. Section 2 is devoted to the proof of Theorem 1.1 and the generalization to singular repulsive potentials using [66] in its last subsection.
Our second main result is devoted to the asymptotic analysis for the particle system (1.1) under the small inertia regime: ε N → 0 as N → ∞. By Theorem 1.1, we expect that for sufficiently large N 1, the system (1.1) in the mean-field/small inertia limit can be well approximated by At the formal level, since ε N → 0 as N → ∞, it follows from the momentum equations in the above system that the hydrodynamic system (1.3) should be replaced by (1.4)-(1.5) as N → ∞. In order to apply our strategy above, we rewrite the equations (1.4)-(1.5) as We can now state our second main result related to a weak-strong stability estimate in the combined mean-field/small inertia limit.
be a solution to the particle system (1.1), and let (ρ,ū) be the unique classical solution of the aggregation-type equation respectively, and the strength of damping γ > 0 is large enough. If the initial data for (1.1) and (1.4) are chosen such that as N → ∞ (and thus ε N → 0). In fact, we have the following quantitative bound estimate:

Remark 1.2. Theorem 1.2 implies that if the initial data satisfies
where C > 0 is independent of both ε N and N . This further yields that the convergences (1.17) and (1.18) and ∂ tū L ∞ can be bounded from above by some constant, which depends only on ∇ x W W 1,∞ , ψ W 1,∞ , ρ L ∞ (0,T ;L 1 ) , and γ . We refer to [24] for details. For general confinement potentials, we can also deal with general strong solutions for compactly supported initial data since their support remains compact for all times. We refer to [1,15] for particular instances of these results.

Mean-Field Limit: From Newton to Pressureless Euler
In this section, we provide the details of the proof for Theorem 1.1. As mentioned before, one of our main mathematical tools is the discrete version of the modulated kinetic energy E N (Z N (t)|U (t)) defined in (1.15).

Modulated kinetic energy estimate
In this part, our main purpose is to give the quantitative bound estimate of the discrete modulated kinetic energy E N (Z N (t)|U (t)).
be a solution to the particle system (1.1), and let (ρ, u) be the unique classical solution of the pressureless Euler system with nonlocal interaction forces (1.3) under the assumptions of Theorem 1.1 up to time T > 0. Suppose that the interaction potential W and the communication weight where C > 0 is independent of N and γ .
Proof. By the notion of our classical solution, we obtain from the momentum Then using this and (1.1), we estimate the discrete modulated kinetic energy functional as Here I 1 can be easily estimated as By definition, we obtain . We next estimate I 3 as On the other hand, the fact ∇ x W ∈ W 1,∞ gives and subsequently this asserts For the estimate of I 4 , we note that Then we rewrite J 2 as This yields Here we can easily estimate I 1 4 as that is, On the other hand, we can estimate Similarly, we also find that Combining all of the above estimates, we have ρ(·, t)).
This completes the proof.
Remark 2.1. We assumed that the communication weight ψ is nonnegative, which takes into account the velocity alignment forces, however a similar bound estimate for the discrete kinetic energy E N to that in Proposition 2.1 can be obtained. Indeed, if ψ can be negative, but bounded, then the third term on the left hand side of (2.1) can be estimated as This yields where C > 0 is independent of N and γ .
In order to close the estimate in Proposition 2.1, we need to estimate the bounded Lipschitz distance between ρ N and ρ. In the proposition below, we provide the relation between the bounded Lipschitz distance and the discrete modulated kinetic energy.

Proposition 2.2.
Let ρ N and ρ be defined as above. Then we have where C > 0 depends only on u L ∞ (0,T ;Li p) and T .
Proof. Consider a forward characteristics η = η(x, t) for the system (1.3) satisfying the following ODEs: subject to the initial data: η(x, 0) = x ∈ R d . The characteristic η is well-defined because of the Lipschitz continuous regularity of u. Note that along the characteristic, the solution ρ can be written in the mild form , s), s) ds , and thus we get t)) .
This together with using the change of variables yields and applying Grönwall's lemma to the above gives where C > 0 depends only on u L ∞ (0,T ;Li p) and T , that is, η is Lipschitz continuous in R d . We also get Here the second term on the right hand side of the above inequality can be estimated as Thus we get and applying Grönwall's lemma to the above deduces where C depends only on u L ∞ (0,T ;Li p) and T . In particular, by taking x = x i (0), we get (2.7) For L 1 , we use the Lipschitz continuity together with (2.6) to obtain For the estimate of L 2 , we notice that Using this identity, the Lipschitz estimate for η in (2.5), and the fact φ ∈ W 1,∞ (R d ), we find where C > 0 is independent of N . We then use Proposition 2.2 to have We finally apply Grönwall's to the above to conclude the desired result.

Convergence estimates
For the convergence estimates, it suffices to prove the following lemma: (ii) Convergence of local energy: (iii) Convergence of empirical measure: Here C > 0 is independent of N .

Proof. (i) For any
(ii) Adding and subtracting, we notice that

Singular interaction potential cases: Coulomb and Riesz potentials
In this part, we discuss the singular interaction potentials. Let d 1 and consider a potential W has the form Note that the case α = d − 2 with d 3 or (2.11) with d = 2 corresponds to the Coulomb potential, and the other cases are called Riesz potentials. With these types of singular potentials, in a recent work [66], the quantitative mean-field limit from the particle system (1.1) to the pressureless Euler-type system when γ = 0, V ≡ 0 and ψ ≡ 0. More precisely, in [66], the following modulated free energy is employed to measure the error between particle and continuum systems: where denotes the diagonal in R d × R d .

12)
where C > 0 is independent of N .

Remark 2.2.
If the interaction potential W is singular at the origin, then the term related to W in (1.1) should be replaced by 1 can not be well defined. This is why the diagonal is excluded in the integration in the modulated potential energy. Remark 2.3. If the right hand side of (2.12) converges to zero as N → ∞, then we also have the same convergence estimates in Theorem 1.1.

Remark 2.4.
Our quantified mean-field limit estimate from (1.1) to (1.3) also apply with a simple combination of Theorems 1.1 and 2.1 for interaction potentials of the form W := W + W with W satisfying ∇W ∈ W 1,∞ (R d ) and W appeared in (2.10) or (2.11).
Proof of Theorem 2.1. For the proof, we only need to reestimate I 3 term in the proof of Proposition 2.1. Although this proof is almost the same with that of [66], we provide the details here for the completeness of our work. Let us denote by On the other hand, we find that Here we used This implies We next use (2.13) to get Thus we obtain This together with the estimates in Proposition 2.1 yields

t)|U (t))+F N (Z N (t)|U (t))
We then apply [66, Proposition 1.1] to have that the last term on the right hand side of the above inequality can be bounded from above by Applying the Grönwall's lemma to the resulting inequality concludes the desired quantitative bound estimate. The convergence result can be directly obtained by using Lemma 2.1. This completes the proof.

Proof of Theorem 1.2
We first start with the case of smooth interaction potentials as in previous section and apply a similar strategy to the proof of Proposition 2.1 to the system (1.16). Then we get whereĪ i , i = 1, 2, 3, 4 are the terms I i , i = 1, 2, 3, 4 in (2.2) with replacing (ρ, u) by (ρ,ū), andĪ 5 is given bȳ whereē = ∂ tū +ū · ∇ xū . This can be simply estimated as where C > 0 depends only ē L ∞ , independent of N and ε N . For the rest, we employ almost the same arguments as before to have where C > 0 is independent of N , ε N , and γ > 0. This yields where C > 0 is independent of N , ε N , and γ > 0. On the other hand, by Proposition 2.2, we can bound the first term on the right hand side of the above inequality from above by where C > 0 is independent of N , ε N , and γ > 0. This together with integrating (3.1) in time implies We finally apply Grönwall's lemma to conclude the desired result in Theorem 1.2.

Singular interaction potential cases
Similarly as before, Theorem 1.2 can be also easily extended to the case with Coulomb or Riesz potentials W defined in (2.10) or (2.11). More specifically, we have the following theorem.  × (0, T )). We further assume thatρ ∈ L ∞ (0, T ; C σ (R d )) for some σ > α − d + 1 in the case s d − 1. Then there exists β < 2 such that where C > 0 is independent of ε N and N .

Local Cauchy Problem for Pressureless Euler Equations with Nonlocal Forces
In order to make the analysis for the mean-field limit from the particle system (1.1) to the pressureless Euler-type equations (1.3) fully rigorous, we need to have the existence of solutions for both systems. As mentioned in Introduction, we postpone the existence theory for the particle system (1.1) in Appendix A, and here we provide local-in-time existence and uniqueness of classical solutions for the system (1.3). For the reader's convenience, let us recall our limiting system: with the initial data Here we set the coefficient of linear damping γ = 1. For the one dimensional problem, the well-posedness and singularity formation for the system (4.1) without the linear damping, the confinement and interaction potentials, called Euler-alignment system, are discussed in [13]. To be more precise, the sharp critical threshold which distinguishes the global-in-time regularity of classical solutions and finite-time breakdown of smoothness is analyzed. The sharp critical threshold estimate is also obtained in [15] for the pressureless damped Euler-Poisson system with quadratic confinement potential in one dimension, that is the system (4.1) with replacing W by N , V = |x| 2 /2, and ψ ≡ 0. For the pressureless Euler-Poisson system, the critical threshold is also discussed in [2,38], see also [69] for the case with pressure. More recently, in [27], the local-in-time existence of classical solutions and finite-time singularity formation are taken into account.
We introduce the exact notion of strong solution to the system (4.1) that we will deal with.
Notice that due to the choice of s in the previous definition, these strong solutions are also classical solutions to (4.1). Our main result of this section is the following local Cauchy problem for the system (4.1). Theorem 4.1. Let s > d/2 + 1 and R > 0. Suppose that the confinement potential V is given by V = |x| 2 /2, the interaction potential ∇ x W ∈ (W 1,1 ∩ W 1,∞ )(R d ), and the communication weight function ψ satisfies where B(0, R) ⊂ R d denotes a ball of radius R centered at the origin. For any N < M, there is a positive constant T * depending only on R, N , and M such that if ρ 0 > 0 on R d and then the Cauchy problem (4.1) has a unique strong solution (ρ, u), in the sense of Definition 4.1, satisfying The assumption on the communication weight function (4.2) implies ψ ∈ W 1, p (R d ) for any p ∈ [1, ∞].

Remark 4.3.
The L 2 -norm of u on the ball is introduced due to the confinement potential V . In fact, if we ignore the confinement potential V in the momentum equation in (4.1), then under the following assumption on the initial data we have the unique strong solution (ρ, u) to the system (4.1) satisfying

Remark 4.4.
In case of a singular interaction potential beyond the Coulomb case, we refer to [27] for the well-posedness theory for the Euler-Riesz system. More precisely, in [27], the local-in-time existence and uniqueness of classical solutions to the system (4.1) with W defined in (2.10) instead of the regular W , γ = 0, V ≡ 0, and ψ ≡ 0 are discussed. One may extend the arguments used in [27] to study the well-posedness for the system (4.1) with W .

Linearized system
In this part, we construct approximate solutions (ρ n , u n ) for the system (4.1) and provide some uniform bound estimates of it.
Let us first take the initial data as the zeroth approximation: We next suppose that the nth approximation (ρ n , u n ) with n 1 is given. Then we define the (n + 1)th approximation (ρ n+1 , u n+1 ) as a solution to the following linear system.
with the initial data Then by the standard linear solvability theory [58], for any T > 0 we have that the approximation {(ρ n , u n )} ∞ n=0 ⊂ Y s,R (T ) is well-defined. For notational simplicity, in the rest of this section, we drop x-dependence of the differential operator ∇ x . Proposition 4.1. Suppose that the initial data (ρ 0 , u 0 ) satisfies ρ 0 > 0 on R d and and let (ρ n , u n ) be a sequence of the approximate solutions of (4.3) with the initial data (ρ 0 , u 0 ). Then for any N < M, there exists T * > 0 such that sup n 0 sup 0 t T * ρ n (·, t) H s + u n (·, t) L 2 (B(0,R)) + ∇u n (·, t) L ∞ + ∇ 2 u n (·, t) H s−1 M.

Proof.
For the proof, we use the inductive argument. Since we take the initial data for the first iteration step, it is clear to find We now suppose that for some T 0 > 0. In the rest of the proof, upon mollifying if necessary we may assume that the communication weight function ψ is smooth. Since this proof is a rather lengthy, we divide it into four steps: • In Step A, we provide the positivity and H s (R d )-estimate of ρ n+1 : for t T 0 , where C > 0 is independent of n, and E : [0, Step C, we estimate the higher order derivative of u n+1 : for t T 0 , where C > 0 is independent of n, and E satisfies the same property as in Step B.
• In Step D, we finally combine all of the estimates in Steps A, B, & C to conclude our desired result.
Step A.-We first show the positivity of ρ n+1 . Consider the following characteristic flow η n+1 associated to the fluid velocity u n by with the initial data η n+1 (x, 0) = x ∈ R d . Since u n is globally Lipschitz, the characteristic equations (4.4) are well-defined. Then by using the method of characteristics, we obtain and applying Grönwall's lemma yields We next estimate H s -norm of ρ n+1 . We first easily find for 2 k s. Here we use Moser-type inequality to estimate I i , i = 1, · · · , 4 as Combining all of the above estimates entails (4.5) for t T 0 , where C > 0 is independent of n.
Step B.-Due to the positivity of ρ n+1 , it follows from the momentum equation in (4.3) that u n+1 satisfies (4.6) Taking the differential operator ∇ to (4.6) gives where we used ∇V = x and I d denotes the n × n identity matrix. Note that We also estimate the last terms on the right hand side of (4.7) as These estimates together with integrating (4.7) along the characteristic flow η n+1 implies By using Grönwall's lemma, we obtain This together with (4.5) asserts where For the L 2 -estimate of u n+1 on B(0, R), we multiply (4.6) by u n+1 and integrate it over B(0, R) to yield (B(0,R)) .
Here we used Thus we obtain d dt u n+1 where C > 0 depends only on R and ψ L 2 . Integrating this over [0, t] with t T 0 and using the estimates (4.5) and (4.8) imply u n+1 where Step C.-For 2 k s + 1, we find where J 1 and J 2 can be estimated as For the estimate of J 4 , we use the fact that W is the Coulombian potential to deduce We next divide J 5 into two terms: Note that Thus for = k − 1 we get and for 0 k − 2 we obtain This asserts Similarly, by integration by parts, we notice that On the other hand, we find that Moreover, for 0 k − 3 we obtain Thus we have and subsequently we get We finally combine all of the above estimate to have and applying Grönwall's lemma gives where we used the estimates in Steps B & C and E 3 : Step D.-We now combine (4.5), (4.8), (4.9), and (4.10) to have ρ n+1 (·, t) H s + ∇u n+1 (·, t) L ∞ + u n+1 (·, t) L 2 (B(0,R)) + ∇ 2 u n+1 for t T 0 , where C > 0 is independent of n, and E : [0, On the other hand, the right hand side of (4.11) converges to ρ 0 H s + u 0 L 2 (B(0,R)) + ∇u 0 L ∞ + ∇ 2 u 0 H s−1 as t → 0 + and that is strictly less than N . This asserts that there exists T * T 0 such that This completes the proof.
On the other hand, if k = 1, we obtain We now combine all of the above estimates to have and subsequently this yields where C > 0 is independent of n. This together with (4.13) asserts that (ρ n , u n ) is We finally provide the uniqueness of strong solutions. Let (ρ, u) and (ρ,ũ) be the strong solutions obtained above with the same initial data (ρ 0 , u 0 ). Set (t) a difference between two strong solutions: (t) := ρ(·, t) −ρ(·, t) L 2 + u(·, t) −ũ(·, t) H 1 .
Then by using almost the same argument as above, we have This concludes that (t) ≡ 0 on [0, T * ] and completes the proof.
. Similarly, we also find Combining all of the above estimates, we obtain for t ∈ [0, T 0 ), where F N (x, v) denotes the discrete free energy given by If d = 2, then we have either where α ∈ (0, 2). On the other hand, if d 3, we obtain for all t ∈ [0, T 0 ), where α ∈ (d − 2, d). Since the right hand side of the above inequality is uniformly bounded in t, we conclude T 0 = ∞ for the case d 2. An upper bound estimate of the distance between particles is a simple consequence of the uniform-in-time bound estimate of the free energy F N due to the confinement potential whenever is present. If V = 0, one can obtain that particles cannot escape to infinity in finite time as soon as ∇ x V has linear growth as |x| → ∞.
Remark A.1. If the interaction and confinement potentials W and V are regular enough, that is ∇ x W ∈ W 1,∞ (R d ) and ∇ x V ∈ W 1,∞ (R d ), we have global-in-time existence and uniqueness of solutions by the standard Cauchy-Lipschitz theory. Moreover, an uniformin-time bound of the distance between particles can also obtained due to the confinement potential if V → +∞ as |x| → ∞.
Let us finally comment on the one dimensional case. If d = 1 and the interaction potential W is given by (2.11), then we apply Theorem A.1 to get the global unique classical solution and uniform-in-time bound estimate. If W is given by the Coulomb potential, that is Thus the interaction force −W is discontinuous, but bounded. In this sense, it is not so singular compared to the other cases. Since the velocity alignment force is regular, we can use a similar argument as in [50,Proposition 1.2], see also [10,41], to have the following proposition.
Proposition A.1. Let d = 1. For any initial configuration Z N (0), there exists at least one global-in-time solution to the system of (1.1) with (A.1) in the sense that (x i (t), v i (t)) satisfies the integral system: Even though Proposition A.1 does not provide the uniqueness of solutions, it is not necessary for the analysis of mean-field limit or mean-field/small inertia limit from the particle system (1.1) to the pressureless Euler system (1.3) or the aggregation equation (1.4).