Confinement for Repulsive-Attractive Kernels

We investigate the confinement properties of solutions of the aggregation equation with repulsive-attractive potentials. We show that solutions remain compactly supported in a large fixed ball depending on the initial data and the potential. The arguments apply to the functional setting of probability measures with mildly singular repulsive-attractive potentials and to the functional setting of smooth solutions with a potential being the sum of the Newtonian repulsion at the origin and a smooth suitably growing at infinity attractive potential.

We say that the nonlocal equation (1) satisfies a confinement property in certain functional setting if every solution ρ(t, ·) to (1) in that setting with compactly supported initial data ρ 0 is compactly supported for all times and its support lies in a fixed ball whose radius only depends on ρ 0 and W .
In most of the mentioned applications, particles/agents repel to each other in a short length scale while there is an overall attraction in larger length scales. Therefore, we can typically concentrate on repulsive-attractive potentials W as in [26,41,22,40,25,36,33,20,21,5,6]. These potentials lead to a rich ensemble of compactly supported steady states whose stability has recently been analyzed [26,41,40,5,6]. While the existence of these compactly supported stable stationary states is a good indication of confinement properties for these repulsive-attractive potentials, it is not equivalent to confinement. Let us finally mention, that repulsive-attractive potentials have also been used in second order models for swarming [18,16,15,12] where exponential decaying at infinity potentials are more suitable from the modelling viewpoint.
Confinement properties were addressed in [10] taking advantage of the wellposedness theory for weak measure solutions of (1) developed in [11]. Using the continuity with respect to initial data in the functional setting of probability measures P(R N ), the authors reduced the confinement of the solutions to (1) to a similar confinement property for solutions of the associated particle ode system:ẋ where Z(i) = {j ∈ {1, . . . , n} : j = i, x j (t) = x i (t)}, x i (0) ∈ R N for all i = 1, . . . , n, 0 ≤ m i ≤ 1, and i m i = M ≥ 0. The authors obtain a confinement property in probability measures assuming that the potential is radial and attractive outside a ball, apart from other technical assumptions related to the well-posedness theory for probability measures. We will improve over the main result in [10] in terms of the assumed attractive strength at infinity. More precisely, we will allow for slower growing at infinity potentials, see Section 2 for the precise hypotheses.
We will later obtain a confinement property for solutions to (1) in a smooth functional setting of compactly supported initial data in W 2,∞ (R N ). The tradeoff is to allow more singular repulsive at the origin potentials. In fact, we will concentrate on the particular case of Newtonian repulsion plus smooth attractive potential with certain growth at infinity. In this functional setting, we can deal with smooth solutions obtained by a slight variation of the arguments in [8,4,9,5]. Section 3 shows that the same ideas used for particles and solutions in the functional setting P(R N ) apply in the continuum model (1) for smooth solutions with this particular potential. The strong Newtonian repulsion at the origin of the potential allows us to derive a priori L ∞ bounds that otherwise are not known. Let us finally mention that confinement for the repulsive Newtonian plus an attractive harmonic potential was obtained in [9].
In the last section of the paper, Section 4, we perform some numerical computations for (2) with different repulsive-attractive potentials and we study the confinement of the stationary states for these cases. The stationary states are found using numerical techniques similar to the ones used in [41,40]. We do a careful study in each case by computing the radius of the support in order to verify numerically the confinement for the solutions. Our numerical studies indicate that the assumptions we impose on the potentials are not sharp and could possibly be improved.

Confinement for probability measures
In this section, we will work with the theory developed in [11] in the framework of optimal mass transportation theory applied to (1). We remind the reader that the equation (1) can be classically understood as the gradient flow of the interaction potential energy [13,3,14]. Optimal transport techniques allow to construct a well-posedness theory in the space of probability measures with bounded second moments P 2 (R N ) at least for smooth potentials [3]. The regularity assumptions on the potential were relaxed in [11] allowing for potentials attractive at the origin with a at most Lipschitz singularity there, i.e., allowing for local behaviors like W (x) |x| a , with 1 ≤ a < 2.
More precisely, we assume that the potential W (x), see [11], satisfies (NL2) There exists a constant C > 0 such that to derive the well-posedness theory of gradient flow solutions to (1) with initial data in P 2 (R N ).
Under this set of assumptions (NL0)-(NL2), we can derive from [11, Theorems 2.12 and 2.13] that the mean-field limit associated to the model in (1) holds. On one hand, this means that approximating the initial data by atomic measures, we can approximate generic solutions of (1) by particular solutions corresponding to initial data composed by finite number of atoms (particle solutions). On the other hand, this also implies that the solution of (1) given in [11,Theorems 2.12 and 2.13] coincides with the atomic measure constructed by evolving the locations of the atoms through the ODE system (2). In other words, if one is interested in showing a confinement property for (1), it suffices to prove the confinement property for the particle system solving (2) since the solutions of the particle system (2) approximate accurately in finite time intervals the solutions of the partial differential equation (1). All these details are fully explained in [10, Section 3] allowing us to reduce directly to particle solutions.
To show confinement, we need additional assumptions on W (x) as in [10]. Throughout this paper, we assume that W is radially symmetric, and attractive outside some ball, i.e., (NL-RAD) W is radial, i.e. W (x) = w(|x|), and there exists R a ≥ 0 such that w (r) ≥ 0 for r ≥ R a .
It is pointed out in [11, Remark 1.1] that (NL1) guarantees that the repulsive force cannot be too strong, more precisely, is bounded above. Here, we take the convention C W = 0 in case R a = 0.
In order to prove confinement results, we need some other condition to ensure that the attractive strength does not decay too fast at infinity. In addition to (NL0)-(NL3) and (NL-RAD), we assume that W satisfies the following confinement condition: which is less restrictive than the assumption in [10], namely lim r→∞ w (r) √ r = +∞.
Therefore, our goal in this section is to show that if the particles interact under a potential satisfying assumptions (NL0)-(NL3), (NL-RAD), and (NL-CONF), have total mass M = 1, center of mass at 0, and are initially confined in B(0,R 0 ), then they will be confined in some ball B(0,R) for all times, wherē R is independent of the number of particles n but only depending on the kernel W and the initial support of the cloud of particlesR 0 . Note that the zero center of mass assumption is possible due to the translational invariance of (1) and (2). Note also that the solution to the particle system (2) in the sense of [10, Remark 2.1] might lead to a finite number of collision times, in which the solution may lose its regularity. Hence when we study the evolution of some quantities in time, we only take the time derivative in the time intervals in which the solution is regular.
The strategy to get confinement for particles is as follows: we need to control quantities that quantifies how much the distribution spreads in time. In [10, Proposition 4.2] the argument was based in following the particle furthest away from the origin and use some energetic arguments to control the mass of the particles nearby pushing the furthest particle. Here we follow a different idea. We consider other moments of the particle system to control the spread of the distribution of particles in conjunction with the evolution of the furthest particle from the center of mass. More precisely, we couple the evolution of the third absolute moment of the particle system with the evolution of the furthest particle.

Evolution of the third moment
Let M 3 (t) denote the third absolute moment of the particle system (2), namely In this section our main goal is to estimate the time derivative of M 3 (t). As we discussed before, there might be a finite number of collision times in which M 3 (t) becomes non-differentiable. Nevertheless, since all the particles have finite velocity, M 3 (t) is Lipschitz continuous in time even during collision. In all the computation below, the time derivative of M 3 is only taken in the time intervals where M 3 (t) is differentiable; and the continuity of M 3 (t) ensures that the fundamental theorem of calculus still holds for M 3 (t).
Since |x| 3 is a convex function on R N , we know that every pair of attracting particles would give a negative contribution to dM 3 /dt, whereas every repulsing pair gives a positive contribution. We can directly evaluate dM 3 /dt as follows: used. Elementary manipulations yield that T ij can be rewritten as which gives the following upper and lower bound for T ij : Next we will find an upper bound for −w (|x i − x j |) in (4). Let us define the nearest particles set N (i) as the set of indexes of particles that are possibly repelling the i-th particle, more precisely, For j ∈ N (i), (NL-RAD) gives that −w (|x i − x j |) ≤ 0. A better bound can be obtained using (NL-CONF): note that for any fixed constant K 1 > 0 to be specified momentarily, there exists some Let us define the set of furthest particles F (i) as the set of indexes of particles whose distance to the i-th particle are larger than R K 1 , namely The definitions of N (i) and F (i) are illustrated in Figure 1. Then the upper bound for −w (|x i − x j |) can be summarized as following: Figure 1: Illustration of the sets N (i) and F (i). For the i-th particle, N (i) is defined as the set of indexes of particles in the red region, while F (i) is the set of indexes of particles in the blue region.

Now we claim that
Its validity is one of the main reasons for imposing the requirement (NL-CONF). To prove the claim, recall that we assume the center of mass is at 0 at t = 0 without loss of generality. Due to the conservation of the center of mass, for any time t, we have x j (t) satisfies j m j x j (t) = 0. Let e i ∈ R N denote the unit vector pointing in the direction of x i , then it follows immediately that n j=1 m j x j · e i = 0.
For |x i | > R K 1 , we split the above sum into three parts, and get where the first inequality is due to the fact that x j · e i > 0 for all j ∈ F (i).
Moreover, recall that we find R K 1 , we force it to be bigger than 2R a . This is to guarantee that for all |x i | > R K 1 > 2R a and j ∈ N (i), the angle between the vectors x i and x j is less than π 6 . As a result, we have and by combining it with (9) we obtain the claim (8).
Due to (8), we deduce that for any Note that we can easily bound T 1 by a constant only depending on W (since |x i | ≤ R K 1 and |x j | ≤ R K 1 + R a ), thus we can rewrite the above inequality as inequality as where C 1 , C 2 only depends on W . At this point we will take a pause on the evolution of M 3 ; we will revisit the inequality (10) soon in Section 2.2 to couple it with the evolution of the support.

Coupling with the evolution of the support
For all t ≥ 0, let R(t) denote the distance of the furthest particle from the center of mass (which we assumed to be 0 without loss of generality), namely It is pointed out in [11, Proposition 4.2] that R(t) is Lipschitz in time. Our goal is to prove that lim sup t→∞ R(t) <R, whereR only depends on W .
We begin by reminding a claim proved in [10, Proposition 2.2]: Let e be any unit vector. Then i.e. the green region in Figure 2 contains at least 1/3 of the total mass.
It is argued in the proof of [11,Proposition 4.2] that for all time t > 0, there is a particle index i 0 (t) (here i 0 may depend on t), such that where d + dt stands for the right derivative. This technical point is due again to the lack of regularity of R(t) for all t, see [11,Proposition 4.2]. From now on, the index i 0 refers to any index satisfying the previous properties.
To control d + dt R(t), it suffices to look at the outward velocity of the i 0 -th particle at this time. As argued in [10, Proposition 2.2] and illustrated in Figure  R (11). 3, the red region is possibly pushing it out, but all the green region is pulling it towards the origin. We proceed by estimating the compensation between these two competing effects. Let G(i 0 ) denote the set of indexes of particles in the green region in Figure 3, namely Recall that throughout this section, we assume the total mass is 1 without loss of generality. It then follows from (12) that Using (13), (2), and some simple manipulations, the growth of R(t) is controlled by the following inequality: where θ(−x i 0 , x j − x i 0 ) denotes the angle between the two vectors −x i 0 and x j − x i 0 , which is less than π 3 as shown in Figure 3. Due to (NL-CONF), for any large constant K 2 , which we will fix later, there exists some radius R K 2 , such that w (r) > K 2 /r for all r > R K 2 . Hence for all t satisfying R(t) > R K 2 , (15) becomes where (14) was used to obtain the last inequality.
To ensure the coupling between the growth of M 3 (t) and the growth of R(t) go smoothly, let us go back to (10) and perform some elementary manipulation And if in addition we have R(t) > 2(R K 1 +R a ), then it follows that G(i 0 ) ⊂ F (i) for any i ∈ N (i 0 ), hence we can replace the F (i) in (17) by G(i 0 ) and obtain where we used (14) again to obtain the second inequality. Finally, we set R 1 := max{R K 2 , 2(R K 1 + R a )}, to ensure that both (16) and (18) hold for R(t) > R 1 .
Finally we are ready to couple M 3 (t) with R(t). By putting the estimates on . On the other hand, we will directly prove that M 3 (t 2 ) must be bigger than M 3 (t 1 ), which causes a contradiction.
Let A 1 be a sufficiently large constant which we will determine later. If the particles start in B(0,R 0 ) and eventually touch the boundary of B(0, A 1 R 1 ), then there exist 0 < t 1 < t 2 , such that More precisely, by letting t 2 := min{t ≥ 0 : R(t) = A 1 R 1 } > 0, and t 1 := max{0 ≤ t ≤ t 2 : R(t) = R 1 }, they would satisfy all the requirements.
Since R(t 2 ) > R(t 1 ), we have Using (16), the above inequality implies and by plugging it into the integral version of (18) we obtain On the other hand, if R(t) successfully grows from R 1 to A 1 R 1 , we will show that M 3 indeed has to increase, namely M 3 (t 2 ) > 2M 3 (t 1 ) for A 1 sufficiently large. First, we can bound M 3 (t 1 ) above by the very rough bound Finally, by noticing that we obtain M 3 (t 2 ) (A 1 R 1 ) 2 . Therefore we can choose A 1 sufficiently large such that M 3 (t 2 ) > 2M 3 (t 1 ), which leads to a contradiction with M 3 (t 2 ) < M 3 (t 1 ).
Note that the proof above shows that R(t) can never reachesR := A 1 R 1 , which is a large constant only depend on W , and in particular is independent of the number of particles.

Confinement for kernel with Newtonian repulsion
In this section, we consider the interaction kernel W (x) given by with N ≥ 2. Here N (x) is the Newtonian kernel, namely where c N denotes the volume of a unit ball in R N . Throughout this section we assume that W a (x) satisfies the following assumptions: Our goal in this section is to show that under the above assumptions, if a solution has total mass M = 1, center of mass at 0, and are initially confined in B(0,R 0 ), then it will be confined in some fixed ball centered at 0 for all times, where the radius of the ball only depends on W a ,R 0 , the dimension N , and the L ∞ norm of the initial data. Remark 1. For N = 1, the Newtonian kernel becomes |x|. Note that in this case the confinement result does not hold under the assumptions above, since the repulsive velocity field between two particles will be a constant regardless of the distance between them, while the attraction may vanish as the distance goes to infinity. We can compensate this difficulty by imposing stronger assumption on the attractiveness of W a at infinity. More precisely, by replacing (W-CONF) by lim r→∞ (w (r) − 1)r = +∞, the confinement result will hold with a similar proof as in Section 2 carried over at the continuum level.
Remark 2. (W-CONF) is more restrictive than (NL-CONF), especially for large N . In the proof below, one can see that dM 3 /dt does not cause a problem at all, indeed it satisfies the same inequality as in the non-singular kernel case in Section 2. The problem lies in dR/dt: due to the singular repulsive kernel, we got a worse control of dR/dt than before, see (25).
We point out that slight variations of the arguments in [5,Section 5] and [4,9] give a well-posedness theory for smooth solutions constructed by characteristics. More precisely, for any compactly supported initial data ρ 0 ∈ W 2,∞ (R N ), there exists a unique classical solution ρ ∈ C 1 ([0, T ] × R N ) ∩ W 1,∞ loc (R + , W 1,∞ (R N )) to (1) with W satisfying (W1)-(W2). Moreover, the associated velocity field v(t, x) = −∇W * ρ is Lipschitz continuous in both space and time, hence the characteristics are well defined: and the solution ρ is given by Since the initial data is compactly supported, it remains compactly supported for all time (although the support may grow in time), and its support is obtained through the C 1 -characteristic maps X t .
First we remind a lemma showing L ∞ -bounds of the solution. This is a classical argument that can be seen for instance in [21,9] and [5, Section 5] but we give a short proof for completeness. Lemma 1. Let W be given by (19), with W a satisfying (W1)-(W2). Let ρ be a classical solution to (1) with compactly supported initial data ρ 0 ∈ W 2,∞ (R N ). Then ρ(t, ·) ∞ ≤ M 0 for all t ≥ 0, where M 0 only depends on W a and ρ 0 .
Proof. Due to the assumption (W1), we can find r 0 > 0 sufficiently small, such that Then it follows from (W2) that M W := sup Let us denote byM (t) = max x∈R N ρ(t, x). If the desired result does not hold, then there exists some t 1 > 0, such thatM (t 1 ) > M 0 andM (t) is increasing at t = t 1 . This enables us to find some x 1 ∈ R N , such thatM (t 1 ) = ρ(t 1 , x 1 ), and ∂ ∂t ρ(t 1 , x 1 ) ≥ 0. On the other hand, since ρ is a classical solution, we have Now let us split the integral ρ * ∆W a in B(0, r 0 ) and outside to get Plugging it into (20), we have which contradicts with the assumption that ∂ ∂t ρ(t 1 , x 1 ) ≥ 0.
Next we present a technical lemma which will be used in the proof of confinement.
for all x 0 ∈ R N and all R > 0.
Proof. First note that it suffices to prove the following inequality holds for all v ∈ L 1 (R N ) ∩ L ∞ (R N ): by letting v(y) = χ B(0,R) (y)u(y +x 0 ), where χ Ω is the indicator function on a set Ω ⊂ R N . We point out that one could use Hölder inequality and interpolation inequality on weak L p spaces to obtain a slightly weaker inequality than above, but we will use an easier and more elementary approach instead.
Let w be an indicator function taking value v ∞ on some disk centered at 0 and taking value 0 outside, where the size of the disk is chosen such that w and v have the same L 1 norm. More precisely, w is given by here c N is the volume of the unit ball in R N . Since v 1 = w 1 , it is straightforward to verify that B(0,r) v(y)dy ≤ B(0,r) w(y)dy for all r ≥ 0.
Now we start with the left hand side of (21), and Fubini's theorem yields that where (22) was used.

Evolution of the third absolute moment
Similar to the particle system case, we also start with estimating the time derivative of the third absolute moment M 3 . Here the third absolute moment M 3 is given by and note that in the continuum setting M 3 is indeed differentiable in time for all t ≥ 0, since ρ(t, x) is a classical solution. The same computation as (4) leads to Due to the assumptions (W-RAD) and (W-CONF) on W a , for any A > 0 (which will be fixed at the end of this subsection), there exists some R A > 1, such that the following bound for −w (r) holds, where C N is some constant only depending on N : Using this bound and (5), dM 3 Similar to Section 2, we again claim that T x r ≤ T x a /2 for |x| > 2R A . We start with controlling T x r . It is easy to check that Note that the singularity of the Newtonian kernel is more difficult to treat than in Section 2. We compensate this difficulty by using the fact that ρ(t, x) is uniformly bounded by M 0 from Lemma 1. Hence for N ≥ 2, we are able to apply Lemma 2 to ρ(t, ·), and obtain here C 3 only depends on N and M 0 as obtained in Lemma 1.
To simplify notation, from now on, we define by m(t, x) the mass of ρ within radius R A of x at time t, namely Recall that in the beginning of this section we assume that ρ 0 integrates to 1, which implies that ρ(t, ·) also integrates to 1 for all t ≥ 0. Then, for any t and x, one of the two following scenarios must be true: either m(t, x) ≥ 1 2 , or 1 − m(t, x) > 1 2 . If m(t, x) ≥ 1 2 at some |x| ≥ 2R A , it follows that m(t, x) 2/N is comparable to m(t, x). Hence, using (23), we get ρ(t, y)(|x| + |y|)dy (since |x| ≥ 2R A ) .
Hence by repeating the same argument on the center of mass as in Section 2 (see (9) and the paragraph after it), we can choose A to be sufficiently large, then we would obtain that T x r ≤ T x a /2.
On the other hand, if the opposite scenario is true at some |x| ≥ 2R A , i.e.
then one can directly bound T x r by C(W, N )|x| by applying Lemma 1 and Lemma 2. Meanwhile it follows directly from the definition of T x a that T x a ≥ A 4 |x|, hence by choosing A sufficiently large we obtain that T r ≤ T a /2. Finally, we choose A to be the maximum value needed in the two scenarios. As a result, T r ≤ T a /2 holds for for all |x| ≥ 2R A , implying that where C 4 , C 5 only depends on W a , N and ρ 0 ∞ . Note that this inequality is parallel to the inequality (10) for the discrete case.

Coupling with the evolution of the support
Next we will proceed similarly as in Section 2.2, where most of the arguments are parallel. We will quickly go through the similar parts in the proof, and emphasize the differences caused by the Newtonian repulsive kernel.
Due to (W-CONF), for any large constant K 3 to be determined later, there exists some radius R K 3 > 6R A such that w (r) > K 3 /r 1/N for all r > R K 3 . Hence whenever R(t) > R K 3 , the growth of R(t) is now controlled by where the second term on the right hand side is obtained in the same way as the last term in (16), except that the power 1 is replaced by 1/N due to (W-CONF).
To deal with the singularity in the first term on the right hand side, recall that ρ(t, ·) ∞ is bounded above by M 0 for all time due to Lemma 1, which again enables us to apply Lemma 2 to obtain where C 6 only depends on N , W a and ρ 0 ∞ . Similar to Section 2.2, we can find some time interval [t 1 , t 2 ], such that R(t) increases from R K 3 to A 2 R K 3 within [t 1 , t 2 ], andṘ(t 2 ) > 0. Here A 2 is a sufficiently large number to be determined at the end of this subsection. Then we have implying that We apply Hölder's inequality on (27), and obtain that Note that this extra step is needed here but unnecessary in Section 2, due to the different powers in (NL-CONF) and (W-CONF).
Now we are ready to couple the growth of M 3 with (28). Since R(t) > 6R A for all t ∈ [t 1 , t 2 ] (recall that when defining R K 3 we set it to be greater than 6R A ), we could treat (24) in the same way as we did in (17) and (18), and bound the growth of M 3 as follows: Then we integrate (29) in [t 1 , t 2 ], and it becomes By putting the above inequality together with (28), we can fix K 3 to be sufficiently large such that M 3 (t 1 ) > M 3 (t 2 ).
Finally, we prove that if A 2 is sufficiently large, we would have M 3 (t 2 ) > M 3 (t 1 ), hence causing a contradiction. It follows from (26) andṘ(t 2 ) > 0 that which can be made to be greater than M 3 (t 1 ) if A 2 is chosen to be sufficiently large, thus we obtain a contradiction with M 3 (t 1 ) > M 3 (t 2 ). This means that R(t) can never reach A 2 R K 3 , thus implies the confinement of support for all times.
Remark 3. Let us emphasize that for potentials W given by (19), we are only able to prove confinement in the continuum setting, not in the particle setting. The reason is that in the coupling method we use, we need to bound the repulsion part of d dt R(t) using the mass in some neighborhood of the outermost particle. In the continuum setting this is achieved by first obtaining an L ∞ bound on ρ(t, ·) in Lemma 1, then applying Lemma 2 to arrive to (26). However, in the particle setting, we are unable to obtain a bound on the "local density" of the particles that is independent of the particle number, and we are unaware of any such results for repulsive-attractive kernels to the best of our knowledge. Intuitively we do expect the "density" of particles to be bounded, since the singular repulsion would not allow the particles to be densely concentrated. We find it an interesting open problem to prove some non-local version of Lemma 1 for the particle system (2) with W given by (19).

Numerics
In this section we numerically check the confinement properties of several potentials together with the long-time behavior of the corresponding particle systems. Let us remark that in all the cases we have simulated, for which confinement holds, the long time behavior of the system seems to converge toward a compactly supported stationary state. In some of the potentials below, this has not been rigorously proved. This is an interesting theoretical question that will be treated elsewhere. Our objective in this section is to check if the conditions under which confinement has been shown in previous sections are sharp or not. With this aim, we remind, as it was said in Section 2, that equation (1) is a gradient flow of the interaction energy with respect to the Wasserstein distance. Thus, stable stationary states of (1) are local minimizers of the interaction energy.
In Section 2 we have shown that the radius R(t) defined by (11) is bounded by a constant R that depends only on the potential W and the initial data. Moreover, R is independent of the number of particles and under certain additional assumptions, see Section 2, we know that the particle systems are indeed good approximations of the solutions to the continuum model (1). For this reason, we have chosen a particle framework to perform our numerical investigation. We also follow the idea of decreasing the energy since stationary states are local minimizers of the energy. Given n particles located at x 1 , . . . , x n ∈ R N with masses m 1 = m 2 = · · · = m n = 1/n, their discrete interaction energy is given by The simulations are done by an explicit Euler scheme leading to a trivial gradient descent method as long as the energy is decreasing at each time step. This method allows to efficiently solve for stationary states of (2). In stiffer situations, as for the Morse potentials below, an explicit Runge-Kutta method is used instead. These methods are essentially the same as the ones used in [41,40] for finding stationary states of different repulsive-attractive potentials. Our stopping criterion is to achieve a numerical steady state. For us, a numerical steady state is a particle distribution for which the discrete l ∞ -norm of the velocity field in (2) is below some predetermined threshold, which we impose to be 0.001/n.
The section is divided into three subsections, each one showing the results for a particular chosen potential. The limit growth for the attractiveness of the potential at infinity under condition (NL-CONF) is log(r). For this reason in Section 4.1 we constructed a piecewise potential with exact logarithmic attraction at infinity. This selection has been done to check the sharpness of condition (NL-CONF). In Subsection 4.2 we go further to take a piecewise potential growing at infinity exactly like log(log(r)), which grows even slower at infinity compared to log(r). This potential does not satisfy the condition (NL-CONF). At the end of this section, in Subsection 4.3, we also analyze the case of the Morse potential. This potential is known to be a repulsive-attractive potential under certain choices of the parameters with negligible attractive strength at infinity, i.e., W (x) → 0 as |x| → ∞. These potentials are more interesting in terms of biological relevance as discussed in [18,12].
As a final remark, we point out that all the used potentials are not singular at the origin and simulations are performed in dimension N = 2.

Logarithmic attraction at infinity
We show confinement when the potential has exact logarithmic attraction at infinity. The chosen potential is It can be checked that it is a repulsive-attractive satisfying w(0) = w(1) = 0, w(r) ∈ C 3 (0, +∞), and the repulsion at the origin is −r. The stationary states are shown in Table 1. We have chosen initial data in such a way that n = 100 n = 200 n = 500 n = 1000 Table 1: Stationary states and radius of their support as a function of the number of particles for the potential w(r) given in equation (30).
the particles feel the logarithmic interaction by randomly placing n particles in a centered square in such a way that |x i − x j | > 1 for some values of i, j ∈ {1, 2, . . . , n}. We have run simulations varying the number of particles n and the initial data. For each simulation the center of mass C n for the particle system was computed and then It is observed that r n (t) < 1 in all the cases for large times and it converges to some asymptotic value. This fact can be explained because when all the particles are out of the range of the log(r) part then the radius and the behavior depends only on the polynomial part of the potential. Simulations indicate that there is confinement for this potential and thus, condition (NL-CONF) is not sharp. Figure 4 shows the evolution of the radius as a function of the particle number n and as a function of time for a particular initial data.

Log-log attraction at infinity
In this case we consider a potential behaving like log(log(r)) at infinity This potential satisfies w(0) = w(e) = 0, it is C 3 (0, +∞), and the repulsion at the origin is −r. The numerical experiments suggest that by increasing the number of particles, the radius of the support of the stationary state increases and stabilizes. Table 2 shows the stationary states as a function of the number of particles n for this potential. These numerical simulations, together with the evolution of the radius of the support both in time and as a function of the number of particles not shown here, indicate that even if the growth at infinity of the potential is less than log(r) there is still confinement. n = 100 n = 200 n = 500 n = 1000 r ∼ 0.7838 r ∼ 0.7954 r ∼ 0.8052 r ∼ 0.8085 Table 2: Stationary states and radius of their support as a function of the number of particles for the potential w(r) given in equation (31).

Morse potential
The usual form of this potential is the following where constants C A and C R are the attraction and repulsion strength respectively and the constants l A and l R are their respective length scales. For our simulations we will take the scaling shown in [18,12]. That is, where V (r) = − exp(−r/l A ) and C = C R /C A and l = l R /l A . It is known for this potential [18,12] that for C > 1 and l < 1 the potential U (r) is shortrange repulsive and long-range attractive with a unique minimum defining a typical distance between particles. Also, in this regime, the condition Cl N = 1 distinguishes between the so-called H-stable and catastrophic regimes.
H-stable case: In our simulations we fix the parameters as C A = l A = 1, C R = 1.9, and l R = 0.8 leading to C = 1.9, l = 0.8, and Cl 2 = 1.216 > 1.  Table 3: Stationary states and radius of their support as a function of the number of particles for the potential U (r) given in equation (31). Table 3 demonstrate that the radius increases by increasing the number of particles, but with a slower rate. As clearly visualized in Figure 5(a), the radius appears to grow like a square root function as the number of particles increases. This observation is further supported by numerical evidence in Figure 5(b), where we plot the square of the radius versus the number of particles, and the linear regression provides a good fit to the data.

Numerical experiments in
It would be interesting to study the H-stable case in more details, although it is outside of the scope of this paper. From these numerical results, we can extract a conjecture that for the continuum system the support of the density would go unbounded over time, i.e., the confinement result should not hold for the H-stable potential.
Catastrophic case: The parameters we choose for the experiments are C A = l A = 1, C R = 1.3 and l R = 0.2 so that Cl 2 = 0.052 < 1.
The results for this case are shown in Table 4. In contrast to the H-stable case, the radius of the support converges to a limiting value. In Figure 6(a) we observe how the radius of the support decreases in time to a limiting value with n = 1000 particles, and in Figure 6(b) we show how the radius of the support of the stationary state increases and converges to a certain value as a function of the number of particles. We conclude that there should be confinement properties for the Morse potential in the catastrophic case.
The final goal would be to find replacements for the condition (NL-CONF) in order to include the cases where w(r) → 0 as r → ∞. One possibility is  Table 4: Stationary states and radius of their support as a function of the number of particles for the potential U (r) given in equation (31).
to invoke scaling limits for integrable potentials. More precisely, we scale the potential in (1) as ρ t = ∇ · (ρ(∇W * ρ)), x ∈ R N , t > 0, (32) in such a way that W (x) = −N W (x/ ) approximates a Dirac Delta at 0 with certain weight as → 0. Now, if the potential is such that then equation (32) is formally approaching ρ t = α∇ · (ρ∇ρ), x ∈ R N , t > 0. In the H-stable case, the Morse potential satisfying α > 0 leads to a limiting nonlinear diffusive equation, which is coherent with the no confinement property. In the catastrophic case, the Morse potential satisfying α < 0 leads to a limiting anti-diffusive nonlinear equation, which might also be coherent with the confinement property of the potential. We conjecture these integrability conditions might have some implications for confinement properties of potentials.