Second order two-species systems with nonlocal interactions: existence and large damping limits

We study the mathematical theory of second order systems with two species, arising in the dynamics of interacting particles subject to linear damping, to nonlocal forces and to external ones, and resulting into a nonlocal version of the compressible Euler system with linear damping. Our results are limited to the $1$ space dimensional case but allow for initial data taken in a Wasserstein space of probability measures. We first consider the case of smooth nonlocal interaction potentials, not subject to any symmetry condition, and prove existence and uniqueness. The concept of solutions relies on a stickiness condition in case of collisions, in the spirit of previous works in the literature. The result uses concepts from classical Hilbert space theory of gradient flows (cf. Brezis [7]) and a trick used in [4]. We then consider a large-time and large-damping scaled version of our system and prove convergence to solutions to the corresponding first order system. Finally, we consider the case of Newtonian potentials -- subject to symmetry of the cross-interaction potentials -- and external convex potentials. After showing existence in the sticky particles framework in the spirit of [4], we prove convergence for large times towards Dirac delta solutions for the two densities. All the results share a common technical framework in that solutions are considered in a Lagrangian framework, which allows to estimate the behavior of solutions via $L^2$ estimates of the pseudo-inverse variables corresponding to the two densities. In particular, due to this technique, the large-damping result holds under a rather weak condition on the initial data, which does not require well-prepared initial velocities. We complement the results with numerical simulations.


Introduction
Nonlocal aggregation models touch various domains of science and technology such as astroparticle physics, microbiology, population biology, social sciences, artificial intelligence and machine learning.The use of integro-partial differential equations in this context, describing the evolution of a density of individuals ρ(x, t) subject to nonlocal interaction forces, such as has become very popular in the literature.In (1), x is a spatial variable typically ranging in R d , t ≥ 0 is time, W = W (x) is a given interaction potential accounting for attractive or repulsive drift among individuals.The modelling approach in (1) allows to formulate concepts of solutions in a "probability measure landscape", which includes the motion of "pointy particles" as a special case of the "continuum theory", see for example the results in [1,9] framed within the Wasserstein gradient flow theory.
At least on a formal level, a similar dichotomy exists in the case of a second-order approach that takes into account inertial effects, namely ∂ρ ∂t + div(ρv) = 0 the "particle-counterpart" of which is given by the second-order particle system In (3), v = v(x, t) is the Eulerian velocity of the fluid-like ensemble of individuals.System (3) can be considered as a nonlocal version of Euler system for gas-dynamics, in which the classical "pressure term" −∇p(ρ) is replaced by a nonlocal interaction force −ρ∇W * ρ.A variant of (3) includes a linear "friction" -or "damping" -term (with σ > 0 a damping parameter) and an external force −∇V , the full model including nonlinear pressure thus looking like ∂ρ ∂t + div(ρv) = 0 The existence theory for systems of the form ( 5) is a classical challenge of the analysis of PDEs, with strong links with the mathematical theory of systems of conservation laws.Since we are not concerned with pressure terms in this paper, we briefly list a few references on this matter such as [28,3,16] for the multidimensional case and [12,14,20,22] for the one-dimensional case.We refer to [10] for a survey on Euler equations.
In the pressure-less case p = 0 in (5), the density ρ is not forced to be absolutely continue with respect to Lebesgue measure.Hence, "particle" solutions in the spirit of (4) are allowed.When two particles collide, a standard way to continue the solution after collision is the so-called sticky particle condition, which forces particles to stay attached to each other after collision, with a post-collisional velocity that is uniquely determined by the conservation of momentum.The existence and uniqueness of such "sticky particle" solutions has attracted the attention of many researchers for decades.We refer to the recent [26] for a through explanation of the issues related with existence and uniqueness in the multi-dimensional case.A case with W = 0 that is particularly interesting in the applications it the Euler-Posson model, in which W is the solution operator to Poisson equation, see [25].In one space dimension the theory is quite rich of results in the literature, we mention here [5,24,4].In particular, [24,4] first addressed the coupling with general nonlocal forces W = 0 in the context of Wasserstein gradient flows.
At least on a formal level, models of the form ∂ρ ∂t = div(∇p(ρ) + ρ∇(V + W * ρ)) can be obtained by rescaling time in (5) as t = στ and letting σ → +∞.Such a singular limit regime is called "overdamped limit", or "large friction limit", and, in case p = 0, it is well-known in the literature of singular limits for systems of conservation laws as a "diffusive relaxation limit".Relevant examples arise in the case of porous-medium like pressures [21] p(ρ) = ρ γ , γ > 1, see also [23] for more general models, and the more recent result in [8].
In recent years the attention of many researchers in this field turned to systems with many species, motivated for example by opinion formation models [15], chemotaxis models with many species of cells and other aggregation phenomena in biology [11,6], pedestrian movements [2].In many of those phenomena, the second order modelling approach via (3) or similar seems more appropriate in that inertial effects, sometimes referred in these contexts as "persistence" effects, do play a role in the model's dynamics.However, while the mathematical theory of many species systems in the first-order modelling approach has been considered in many papers [], very little attention has been devoted to second order models with many species.In this paper we wish to provide a first contribution in that direction, by restricting for the moment to the one-dimensional case.
In system (7), ρ and η are probability measures modelling two species of agents, or individuals, v and w are the corresponding Eulerian velocities of the two species, σ > 0 is the damping parameter, H ρ , H η , K ρ , K η are smooth (to an extent to be specified later) given space-depending potentials.K ρ and K η are called self-interaction potentials as they describe the interaction between the agents of same species, H ρ and H η are called cross-interaction potentials and model the interaction between the agents of opposing species.The convolutions in (7) are meant with respect to the space variable.All potentials appears in the system with their first derivative.This choice of ours is merely motivated by the fact that all those terms should be considered as gradients of potential energies.System (7) has a natural discrete particle counterpart.Let us consider x 1 , . . ., x N as N particles of the first species with masses m 1 , . . ., m N , and y 1 , . . ., y M as M particles of the second species with masses n 1 , . . ., n M .The dynamics of x i and y j is determined by the following equations with i = 1, . . ., N and j = 1, . . ., M and the following initial data The candidate large-friction of ( 7) is the first order system which was extensively studied in [17], see also [13] for the case with cross-diffusion terms.
Our work contributes to the above line of research in what follows.We stress that our results only deal with the one space dimensional case.
(i) In the case of smooth interaction potentials, we provide a well-posedness result in the 2-Wasserstein space of probability measures.Our main result is contained in Theorem 1. (ii) We then investigate the large damping limit and prove that, under a suitable rescaling, our system converges towards the corresponding first order model.Our result in this framework is stated in Theorem 2. We observe that, at least to our knowledge, the technique used in Theorem 2 was never used in the one-species case.(iii) We then consider the case of Newtonian potentials for the self-interaction part and suitably coercive external potentials and prove a large-time collapse result in Theorem 3. The paper is structured as follows.In section 2 we introduce the main concepts of gradient flows in Wasserstein spaces that we need in our paper, we introduce the large-damping scaling limit, we state our model in a suitable Lagrangian framework, and we state our main results (see subsection 2.6).In section 3 we prove Theorem 1, essentially following the classical strategy of [7], which also used in [4] for the one-species case.In section 4 we prove Theorem 2. In Section 5 we consider the case of Newtonian potentials and prove existence of sticky solutions and the large time collapse to Dirac deltas stated in Theorem 3. Finally, in Section 6 we provide some numerical simulations.

Preliminaries and main results
In what follows we will set the notations, the assumptions, and introduce definitions that will be used throughout the paper, see Subsections 2.1 and 2.2.Subsections 2.3 and 2.5 are devoted to the precise description of system (7) in terms of particles and Lagrangian coordinates respectively.In Subsection 2.4 we provide a formal argument for the large damping limit of system ( 7) towards (10).Finally Subsection 2.6 collects the main results of the paper.

2.1.
One dimensional Wasserstein distance.We start introducing some preliminaries and definitions on the metric structure; the reader can refer to the classical references [1,27] for further details.Let P(R n ) be the set of Borel probability measures on R n .Given µ ∈ P(R n ) and a Borel map T : R n → R m , we denote by ν := T #µ ∈ P(R m ) the push-forward of µ through T defined by We denote by P 2 (R d ) the set of probability measures on R d with finite second moment, i.e., ´Rd where Π(µ, ν) denotes the class of transport plans between µ and ν, i.e., the probability measures where π i is the projection operator on the i-th component of the product space.By introducing the class of optimal plans between µ and ν, i.e., minimizers of (11), denoted by Π o (µ, ν), the Wasserstein distance can be rewritten as In the one-dimensional case, there exists a unique optimal plan γ ∈ Π o (µ, ν) for which the infimum in (11) is attained, and it can be characterised by the monotone rearrangements of µ and ν: given µ ∈ P(R), its monotone rearrangement is where Ω := (0, 1) and M µ is the cumulative distribution of the measure µ, i.e., The map X µ is right-continuous and nondecreasing and satisfies, by denoting the one-dimensional Lebesgue measure on Ω by m, for all Borel maps ζ : R → R. In particular, µ ∈ P 2 (R) if and only if We further recall that, introducing the closed convex set of right-continuous non-decreasing functions in the Hilbert space L 2 (Ω), i.e., the map is a distance-preserving bijection between the space of probability measures with finite second moments P 2 (R) and the convex cone K of non-decreasing L 2 -functions.Since we are dealing with a two-species system, we will work on the product space P 2 (R) × P 2 (R).For all µ = (µ 1 , µ 2 ), ν = (ν 1 , ν 2 ) ∈ P 2 (R) × P 2 (R), we define the product Wasserstein distance as follows

Main assumptions.
In this subsection we collect the main assumptions we will need in the rest of the paper.We start by specifying the class of interaction potentials we are going to use.
An admissible potential K is said to be sub-quadratic at infinity if there exists a constant C > 0 such that An admissible potential K has a sub-linear gradient if there exists C > 0 such that We call an admissible potential attractive if In Section 5 we will also take into account the action of external potentials in the dynamics.More precisely, we consider A ∈ C 2 (R) and assume that there exist the positive constants λ and α such that A(x) ≥ λ|x| 2 (H1) for all x ∈ R. Denoting with •, • L 2 (Ω) 2 the inner product on the space L 2 (Ω) 2 , that is In particular, throughout the paper, we will usually consider as Hilbert spaces Let I K : L 2 (Ω) → [0, +∞) be the indicator function of the L 2 -convex cone K introduced in (13), that is For a given X ∈ L 2 (Ω), the sub-differential of I K in X is given by or in its alternative form We conclude this subsection with the following definition, which we borrow from [4].
for a.e.m ∈ Ω and all X ∈ K.
An operator F : K → L 2 (Ω) is uniformly continuous if there exists a modulus of continuity ω such that for all X 1 , X 2 ∈ K.
2.3.Particles system.We dedicate this subsection to the study of sticky solutions in the finite dimensional case.Let x = (x 1 , . . ., x N ) ∈ R N and y = (y 1 , . . ., y M ) ∈ R M be the positions of particles of the first and second species respectively.The "sticky" condition clearly preserves the ordering of the particles, therefore their evolution is confined in the closed convex set Setting v = (v 1 , . . ., v N ) ∈ R N and w = (w 1 , . . ., w M ) ∈ R M as the velocity vectors of particles of the first species and second species respectively, we consider the following system for i = 1, . . ., N and j = 1, . . ., M. In system (15), The vector field models the interactions between particles of the first species and the i-th particle of the first species, while the i-th component of the vector field b(x, y) : (x, y) ∈ K N × K M → (b 1 (x, y), . . ., b N (x, y)) ∈ R N describes the interactions between the i-th particle of the first species and particles of the second species.Similarly one can describe the j-th component of the terms Assuming that all the potentials in ( 15) are smooth enough (for example with C 2 regularity), a unique solution to (15) exists as long as particles occupy distinct positions.When two or more particles collide, we apply the concept of sticky particle solution sketched in the introduction.Following [4,24], the precise formalisation of sticky collisions requires the definition of the following normal cones Note that the normal cone N x K N is equal to the sub-differential ∂I K N (x) of the indicator function of K N at the point x.When two particles of the same species collide, an instantaneous force is released and the respective particles velocities evolve as elements of the normal cones N x K N and N y K M respectively.Given these premises, we can consider the second-order system of differential inclusions System ( 16) is justified as follows.Introducing the vector where A(x, y) is the vector in R N +M with components a(x)+b(x, y) and c(y)+d(x, y) respectively.Now, due to the smoothness of the interaction potentials, the vector field A(x, y) can be extended by continuity to the boundary of the cone K N × K M .Therefore, as W and (v, w) only differ by a scalar factor, a suitable modified version of the differential equation for W that keeps the dynamics in which easily yields the last two differential inclusions in (16).
According to [4], if x : [0, ∞) → K N satisfies the global sticky condition, i.e., particles are not allowed to split after colliding, then the following monotonicity property on the family of normal cones N x(t) K N holds: for all s < t.

Hence, for any function
Consequently, integrating the last two equations in ( 16) on a time interval [s, t], one obtains System ( 16), together with ( 17) and ( 18), can be rewritten in a more compact form in the new variables (x, y, p, q) where p and q are defined by yielding the following first-order system of differential inclusions with the additional characterisation of v and w in terms of p and q 2.4.Time scaling and formal large damping limit.One of the purposes of the present work is to study system (7) in the large time / large damping regime, namely we aim to send σ → +∞ in (7) after having suitably rescaled the time variable.We start performing the scaling at the level of particles, namely for system (9).Consider the new time variable τ defined by and introduce the scaled particle trajectories as follows: Notice that we can scale the initial velocities accordingly as := ω j = σw j .
Hence, system (9) becomes A formal limit σ → +∞ leads to the following first-order system of differential equations for particle positions A similar time scaling can be performed at the level of (7).Using the definition of τ in (19) and considering ( ρ, v, η, w) solution to we can introduce the rescaled densities and velocities as Then the quadruple (ρ, v, η, w) solves and formally, as σ → ∞, we have 2.5.Lagrangian description of the continuum model.We now transpose the considerations above in terms of a Lagrangian description for system (7).For any X ∈ K, where K denotes the convex cone introduced in (13), we define the set and the closed subspace A crucial quantity in the following analysis is the projection P HX : L 2 → H X given by for all U ∈ L 2 (Ω).The proof of the following Lemma is an easy consequence of Jensen's inequality, see [4,Lemma 2.2].
Lemma 1 (H X -contraction).Let ψ : R → [0, ∞) be a convex l.s.c.function.Then P HX is dominated by X, namely and we write P HX ≺ X.
Consider a quadruple (ρ, η, v, w) solution to (7) and define the maps where Ψ is the isometry defined in ( 14) that associates to a probability measure its monotone rearrangement.In the new unknowns (X, Y, V, W ), system (7) can be (formally) rephrased as Similarly to Section 2.3, one can show that the previous system can be reformulated in terms of differential inclusions to incorporate particles collisions.Moreover, since we will investigate on the large-damping limit, through the paper we consider the Lagrangian counterpart of the rescaled system (20).Then, according to the previous calculations, we get the system with ε := σ −2 and where we have denoted by We observe that if K ρ , H ρ , K η , H η are C 1 functions that satisfy (A) and (SL) then the two operator F and G defined in (26) and ( 27) are uniformly continuous and bounded according to Definition 3.
In order to consider the case of Newtonian potentials, we introduce the following notion of generalised Lagrangian solutions for system (25) under globally sticky dynamics, see [4].
and, similarly, where F [X(t), Y (t)] and G[X(t), Y (t)] are the operators defined in ( 26) and ( 27).(2) Semigroup property: for all t ≥ t 1 ≥ 0, the right derivatives (3) Projection formula: ] with F and G as in ( 26) and ( 27) and the interaction potentials K ρ , K η , H ρ and H η satisfying (A) and (SL), then any Lagrangian solution is a generalised Lagrangian solution.
In the following we will make use of the auxiliary variables and that allow to rephrase system (25) in the equivalent form 2.6.Main results.We collect in this subsection the main results presented in the paper.The first result concerns the well-posedness of system (7) in the 2-Wasserstein space of probability measures and in the sense of sticky solutions, under smoothness assumptions on the interaction kernels.
Theorem 1.Let T > 0 and suppose that the kernels H ρ , K ρ , H η , K η ∈ C 1 (R) satisfy (A) and (SL).Let ρ, η ∈ P 2 (R) and v ∈ L 2 (dρ) and w ∈ L 2 (dη).Then, there exists a unique quadruple that is a distributional solutions to system (7) such that We then address the σ → ∞ limit of ( 7) towards (10) using the rescaling in (20), making rigorous the formal argument presented in Section 2.4.This task is performed at the level of the Lagrangian system (25) sending the parameter ε = σ −2 → 0, coming back to the Eulerian description through the isometry (14).The following result is proved in Section 4.
Theorem 2. Let T > 0 and suppose that the kernels H ρ , K ρ , H η , K η ∈ C 1 (R) satisfy (A) and (SL).Let (ρ ε , η ε , v ε , w ε ) be solution to system (20) with ε = σ −2 under the initial condition (ρ ε , η ε , v ε , w ε ) and let (ρ, η) be solution to system (21) with initial data (ρ, η).Furthermore, assume that (i) Remark 1 (Initial data are not well-prepared in the velocity variable).In Theorem 2, recalling that dρ) and w ∈ L 2 (dη) are given and independent of ε.Therefore, assumption (ii) is quite general in the context of singular limits.Assumption (i) instead imposes that the initial density should converge to the one of the limiting first order system.
Lastly, under the action of Newtonian self-interaction kernels, K ρ (x) = K η (x) = N (x) := |x|, symmetric and attractive cross-interactions, H ρ (x) = H η (x) = H(x) and suitably coercive external potentials, we focus on a different aspect, that is the convergence to stationary solutions of (7).More precisely, we will consider the following system and its Lagrangian counterpart Stationary solutions in this case are (ρ s , η s ) = (δ 0 , δ 0 ) where δ is the Dirac measure, which corresponds to (X s , Y s ) = (0, 0) in terms of the Lagrangian description.The last result we present in the paper shows that solutions to (37) converge to the stable stationary solution in the 2−Wasserstein distance.

Existence and uniqueness for smooth potentials
In this section we prove Theorem 1, namely existence and uniqueness of solution to system (7).To perform this task, we pass through existence of solutions to the Lagrangian system (25), where we apply the theory of Maximal Monotone Operators subject to Lipschitz perturbations in the spirit of [7, Theorem 3.17], see Proposition 1 below.The result in the original variables is then proved using the properties contained in Proposition 2 below.
We start proving the following Lemma.
Lemma 2. Let (X, Y ),( X, Y ) ∈ K×K be given.Consider the interaction kernels H ρ , K ρ , H η , K η under assumptions (A) and (SL) and let F and G be the operators defined in (26) and ( 27) respectively.Then there exist two positive constants C 1 and C 2 depending on the Lipschitz constants of the kernels, such that Proof.We only prove (i) since (ii) follows from a similar argument.By the definition of F in (26) we have Using the fact that |x + y| 2 ≤ 2(|x| 2 + |y| 2 ), the right hand side of (39) can be controlled by Let L(K ′ ρ ) and L(H ′ ρ ) be the Lipschitz constants of K ′ ρ and H ′ ρ respectively, then, using Jensen's inequality, the right hand side of (40) is bounded by Thus, there exists a positive constant ) .Analogously, one can prove the inequality (ii), we omit the details.
We are now ready to state existence result for Lagrangian solution to system (25).
Proof.According to the discussion in Section 2.5, system (25) can be rewritten in the following equivalent form where P and Q are definded in (34) and ( 35) respectively.In order to prove the result we will follow the strategy in [7, Theorem 3.17].Consider the operator 2 .Note that A is convex and bounded from below.Consider the iterative sequence defined as follows: fix U 0 := (X, Y , P , Q) ≡ (X, Y , εV + X, εW + Y ) and, for n ≥ 1 construct U n+1 (t) := (X n+1 (t), Y n+1 (t), P n+1 (t), Q n+1 (t)) recursively as the weak solution to the implicit-explicit system the previous system (42) can be rewritten in the following compact form Since the functional A is convex, its sub-differential is a maximal monotone operator in the sense of [7] and R can be seen as a Lipschitz perturbation of it, see [7,Lemma 3.1].A direct computation shows that 1 2 , then proceeding as in [7, Lemma A.5], we have that Invoking Lemma 2 and the definitions for P and Q in (34) and ( 35) respectively, we can say that there exists a positive constant C depending on T , ε and on the Lipschitz constants of the kernels An easy iterative procedure implies that thus, U n uniformly converges on [0, T ] to some U .Due to the Lemma 2, R is continuous in L 2 in each component.Moreover, since the subdifferential of A is closed, we can pass to the limit in (43) and obtain that U is a weak solution to the system (41).Concerning uniqueness, let U 1 = (X 1 , Y 1 , P 1 , Q 1 ) and U 2 = (X 2 , Y 2 , P 2 , Q 2 ) be two solutions to system (41) with the same initial condition U 1 = U 2 = U .Proceeding in an analogous way as before, we can argue that where the positive constant that proves the uniqueness.
The following Proposition collects some properties of Lagrangian solution.
Proposition 2. Let F, G : K × K → L 2 (Ω) be uniformly continuous operators in (26) and (27) and let (X, Y ) be the Lagrangian solution to (25).Then, the following properties hold: exist for all t ≥ 0. (ii) V and W are the unique elements of minimal norm in the closed convex sets and respectively.In particular, by replacing Ẋ by V and Ẏ by W , (25) and (36) hold for all t ≥ 0. (iii) The functions t → V (t) and t → W (t) are right-continuous for all t ≥ 0. (iv) If T 0 X ⊂ (0, ∞) and T 0 Y ⊂ (0, ∞) denote the subsets of all times at which the maps s → V (s) L 2 (Ω) and s → W (s) L 2 (Ω) respectively are continuous, then (0, ∞) \ T 0 X and (0, ∞) \ T 0 Y are negligible, V and W are continuous, X and Y are differentiable in L 2 (Ω) at every point of T 0 X and T 0 Y respectively.(v) Setting ρ(t, •) := Ψ −1 (X(t, •)) and η(t, •) := Ψ −1 (Y (t, •)) where Ψ is the isometry introduced in (14), there exist a unique map v(t, •) ∈ L 2 (R, ρ) and a unique map w(t, for every t ∈ T 0 X , and for every t ∈ T 0 Y .Proof.The results in (i), (ii), (iii) are consequences of the general theory of [7,Theorem 3.5].Concerning (iv) and (v), we follow [4,Theorem 3.5].We prove only (47), since the proof of (48) is similar.By applying [7, Remark 3.9], one can see that if t is a point of differentiability of X, the derivative with respect to time of X in t is the projection of 0 onto the affine space generated by P (t) − ∂I K (X(t)) − X(t), i.e., the orthogonal projection of P (t) − X(t) onto the orthogonal complement of the space generated by ∂I K (X(t)).By using [4, Lemma 2.5], we obtain (47).Since any element of H X(t) can be written as v • X, where v ∈ L 2 (Ω) is a suitable Borel map, we have that there exists a Borel map We are now in the position of proving the main result of this Section, namely Theorem 1, that concerns existence and uniqueness of the solution to system (7).
Proof of Theorem 1.Let ρ, η ∈ P 2 (R) and v ∈ L 2 (dρ), w ∈ L 2 (dη) be given initial conditions.Define the L 2 (Ω)-functions X = Ψ(ρ) and Y = Ψ(η) and the compositions V = v • X and W = w•Y .Then (X, Y , V , W ) is an admissible initial condition for system (25), thus Proposition 1 ensures existence and uniqueness of a couple (X, Y ) that is the Lagrangian solution to (25).According to Proposition 2 we can define the right-continuous functions V and W such that (44) holds for all t ≥ 0 and introduce ρ(t, •) := Ψ −1 (X(t, •)) and η(t, •) := Ψ −1 (Y (t, •)).Let v(t, •) be the map given by Proposition 2 and ϕ be a test function on (0, T ) × R, then Using (47) and integrating by parts, the r.h.s. of ( 49) is equal to As proved in Proposition 2 we have that Ẋ(t, m) = V (t, m) and from the definition of the operator P (t, m) in (34), one obtains that (50) equals that is the distributional formulation of the momentum equation in (7).Similarly, for the continuity equation we have Concerning the initial conditions, since lim t↓0 X(t) = X in L 2 (Ω) for Proposition 1 and X = Ψ(ρ), we have that ρ → ρ in P 2 (R) as t → 0.Moreover, A similar argument can be used for the pair (η, w).

Large-damping limit
In this section we study the large-damping limit of system (7) for the damping parameter σ → ∞ as stated in Theorem 2. In particular, we aim at making the formal argument introduced in Section 2.4 rigorous, and showing that solutions to system (20) converge to the ones of the first-order system In what follows we will assume that the potentials H ρ , H η K ρ , K η are under assumptions (A) and (SL).
Recalling the definition of F [X, Y ](m) and G[X, Y ](m) in ( 26) and ( 27), we introduce the operator . 25) can be rewritten in the following compact form We are now in the position of proving Theorem 2.

Newtonian potentials
This section is devoted to study existence of solutions and asymptotic property of system (7) when self-attractive forces are driven by Newtonian potentials, i.e., K ρ (x) = K η (x) = |x|.In order to proceed, we need to restrict the analysis to the case of equal cross potentials, namely H ρ = H η =: H.We also consider two uniformly convex external potentials A ρ and A η acting on the system.More precisely, we assume A ρ , A η ∈ C 2 (R) under assumptions (H1) and (H2).
These additional terms don't affect the study of existence of solutions, in the generalised sense specified in Definition 5, but are only required in the study of asymptotic behaviour in Theorem 3. The system we are dealing with can be expressed in Lagrangian coordinates as follows We can associate to the system (56) the following functional In particular, we write where As shown in [18,9], it is easy to prove that the self-interaction contributions in F are linear when restricted to K.
Since X ∈ K, X is non-decreasing, then the set {X(m) ≥ X(s)} can be characterised as follows The first result in this Section consists in proving the existence of a map t → (X(t), Y (t)) that is a generalised Lagrangian solution to (25) with respect to the choice Θ = P HX (F 1 )(t, m) and Ξ = P HY (F 2 )(t, m), i.e., the system (56) can be written as follows where and are the force operators and describe the external and interaction forces that act on the system.The following proposition ensures that a generalised Lagrangian solution exists.
defined on the intervals L i := [l i−1 , l i ) and Z j := [z j−1 , z j ), for i = 1, . . ., N − 1 and j = 1, . . ., M − 1.Consider the finite dimensional Hilbert set and its closed convex cone Note that the projected forces are well defined and Lipschitz continuous according to the definitions in ( 59)-( 60) and assumptions (A) and (SL).Now, assume that the initial condition then, we obtain (61)-( 62) by solving We have that H m = H X(t) and H n = H Y (t) in [0, t 1 ), thus the projection onto the set H m yields piecewise constant functions defined on the same intervals as (X, V ), and similarly the projection onto H n .Taking t 1 as the new initial time, we can consider a new initial condition < N and M ′ < M and, proceeding in the same fashion, we can define t 2 > t 1 and consider the evolution in the time interval [t 1 , t 2 ).
Iterating the procedure, we obtain a sequence of collision times 0 =: t 0 < t 1 < • • • < t K := ∞ and the quadruple (X, Y, V, W ) such that When an inelastic collision occurs, we have that In order to prove inclusion (30), it is not restrictive to assume t 1 = 0. We proceed by induction on the collision times.In the first time interval [0, t 1 ), inclusion (30) holds by considering the empty set for the subdifferential ∂I K (X(t)).Now, suppose that (30) is satisfied in [t k−1 , t k ).Hence, by induction assumption, with ξ ∈ ∂I K (X(t k )).By (64), for any t ∈ [t k , t k+1 ).Combining equations ( 67) and (68) we get Invoking again (64), we have hence using (66), we derive Applying [4, Lemma 2.6], we find that , and using the monotonicity property of the sub-differential, one obtains that for all t ∈ [t k , t k+1 ).Therefore inclusion (30) is satisfied.Now, let us prove that (32) holds.
Consider system (36) with P replaced by Thus, we have that for any t ≥ s ≥ 0, where we used the monotonicity of the sub-differential.Integrating on s ∈ [0, t] we obtain for a.e.t ≥ 0. Since the following property holds (cfr.[4]) we derive A similar proof holds for the equations (31) and (33).Finally, since the construction above starts form descrete initial data in the form of the piecewise constant functions as in (61)-(62), and since these functions are dense in L 2 (Ω), we can approximate any given initial data and then combine the procedure into the proof with the stability Theorem 4.4 in [4].
Thanks to the non-negativity of the cross-potential H, assumption (H1) and the fact that which holds since m 2 − m ≤ 0 for m ∈ (0, 1) and ∂ m X + ∂ m Y ≥ 0 for X, Y ∈ K, we obtain that where C 1 is a constant depending on initial data.Thus Computing the temporal derivative of the L 2 -distance between (X, Y ) and (X s , Y s ), we derive In order to control the last term in the chain of equality above we compute Using the definitions of F 1 and F 2 in (59) and (60) and the property for the projection operator we have Using assumption (AT) we can bound the terms involving the cross-interaction potential H as follows thus, using assumption (H2) and ( 73), (77) can be bounded from above by Note that for any A > 0 we have −XV ≤ X 2 A 2 + V 2 4A 2 .Then, applying this inequality to −σXV and −σY W , we obtain the following inequality holding for any A 1 , A 2 > 0: By taking sufficiently small A 1 and A 2 , we have that (78) is bounded from above by for some constants C 1 , C 2 > 0. Putting together estimates (75) and (80), we have that Integrating in time inequality (81), for all T > 0 we obtain where C 2 is a constant which depends on initial data.Proceeding as in (79) and using the bound in (74), we have that Combining estimates (74) and (82) we get Since the operator F defined in (57) is a monotone operator, then and ℓ is unique.Moreover, Lemma 2 ensures that the operator F is continuous, thus Using the coercivity of the external potentials A ρ and A η and (83), we have that ℓ is necessarily zero, hence the statement holds.

Simulations
This last section is devoted to provide some numerical example on the behaviour of solutions to system (7).Numerical simulations will be performed by using the discrete particle counterpart of ( 7), namely solving numerically (9).We recall that the system of ODEs we are dealing with is the following where x i and y j denotes the particles positions of first and second species respectively, v i and w j their velocities and m i and n j their masses, for i = 1, . . ., N and j = 1, . . ., M .By a normalisation in the masses the total number of particles for each species, N and M respectively, will be modified in each of the examples below in order to highlights possible different changes in the solutions.System (84) will be coupled with an uniform distributed set of particles in the space interval [0, 1] and a randomly distrubution for the velocities.We then let the particles evolve by using an explicit second-order three steps Runge-Kutta method, (cf.[19]) up to the first collision.In order to detect collisions between particles we fix a tolerance parameter toll and we assume that it occurs when the distance between two consecutive particles of the same species, for instance x i and x i+1 , is smaller than toll.Once two consecutive particles collide they are replaced by a single particle with new position and velocity given by and doubled mass, and we let the system evolving again with this new set of particels.In all the simulations below we fix toll = 0.002.
We study numerical solutions to the system (84) both in case of smooth potentials and in case of Newtonian self-potentials.Several examples are presented in the smooth case, where we highlight the possibility of a sticky dynamics, both in attractive and repulsive regime.Furthermore, we will compare solutions to second-order system with solution to first-order one as the increasing values of the damping parameter σ, also comparing the Wasserstein distance between the solution to the second-order system and the solution to the first-order system as σ varies.Wasserstein distance is computed using its one-dimensional equivalence with the L 2 −norm at the level of monotone rearrangements.
The first examples we provide concern the evolution of particles subject to the action of radial smooth potentials.Figure 1 displays the sticky particle dynamics when all the potentials are smooth and attractive.Instead, in Figure 2 the self-potentials are attractive and the crosspotentials are repulsive, while in Figure 3 the self-potentials are repulsive and the cross-potentials are attractive.In particular, we highlight how the behaviour is strongly different by comparing two simulations performed with the same potentials, number of particles and initial position, but different set of initial (random) velocities.
We then show a couple of simulations in which the self-potentials are attractive Newtonian, while the cross-potentials are symmetric, radial and smooth.In particular, in Figure 4, the cross-potentials are attractive, indeed the particles collide, while in Figure 5, they are repulsive and not all the particles collide.According to results in Section 5 also the effect of external potentials is taken into account.
We then focus on the numerical investigation of the large damping regime.Figures 6 and  7 show a comparison between the particle evolution associated to the second order system and the ones associated to the first order system (51), for various choices of potentials.We highlight numerically the relevance of the dumping parameter σ in the evolution: increasing the value of σ solutions of the two different problems becomes indistinguishable.
Finally in Figure 8, considering the same potentials in Figure 4, we display the Wasserstein distance between the solution to the second-order system and the ones to the first-order system for different values of σ.For small values of σ, the Wasserstein distance grows initially, and then decays in time.When σ is bigger, the distance remains controlled for all times.Email address: valeria.iorio1@graduate.univaq.it

Figure 3 .
Figure 3. Two possible outcomes (top and bottom) for the evolution of the system under the action of self-repulsive potentials K ρ (x) = 2e −|x| 2 and K η (x) = e −|x| 3 and attractive cross-potentials H ρ (x) = |x| 2 and H η (x) = −e −3|x| 2 .In both the simulations the numbers of particles are fixed as N = 170 and M = 160, but initial velocities change (randomly).

Figure 8 .
Figure 8. Behaviour of the Wasserstein distance between solutions of the first order system and solutions of the second order system.The self-potentials are Newtonian attractive potentials, while the crosspotentials are given by H ρ (x) = H η (x) = −e −|x| 2 .Increasing the damping parameter Wasserstein distance remain controlled.
2, we recall below the notion of Fréchet subdifferential for a generic operator F on a general Hilbert space.