Multiple large-time behavior of nonlocal interaction equations with quadratic diffusion

In this paper we consider a one-dimensional nonlocal interaction equation with quadratic porous-medium type diffusion in which the interaction kernels are attractive, nonnegative, and integrable on the real line. Earlier results in the literature have shown existence of nontrivial steady states if the $L^1$ norm of the kernel $G$ is larger than the diffusion constant $\varepsilon$. In this paper we aim at showing that this equation exhibits a"multiple"behavior, in that solutions can either converge to the nontrivial steady states or decay to zero for large times. We prove the former situation holds in case the initial conditions are concentrated enough and"close"to the steady state in the $\infty$-Wasserstein distance. Moreover, we prove that solutions decay to zero for large times in the diffusion-dominated regime $\varepsilon\geq \|G\|_{L^1}$. Finally, we show two partial results suggesting that the large-time decay also holds in the complementary regime $\varepsilon<\|G\|_{L^1}$ for initial data with large enough second moment. We use numerical simulations both to validate our local asymptotic stability result and to support our conjecture on the large time decay.


Introduction
Several phenomena in biology are governed by the combination of long-range attractive effects with short-range repulsive ones. Examples range from chemotaxis of cells (see e.g. [8]) to swarming phenomena in animal biology [35,9,34,39]. These situations can be modelled at either particle (microscopic) or continuum (macroscopic) level. In the former case one assumes that in a given finite set of particles (or individuals) each pair interact through a "drift" velocity field depending on their distance, consisting of attractive and repulsive forces. What characterizes our approach is that the attractive force maintains a long range scale while passing to the continuum regime, whereas the range of interaction in the repulsive terms degenerates as the number of individuals becomes large. Such approach leads to our prototype model in which the unknown ρ = ρ(x, t) is sought in the set of time dependent curves from t ∈ [0, +∞) onto the space of nonnegative L 1 densities with fixed mass. In (1) the quadratic diffusion term models (local) repulsive effects via a diffusion constant ε > 0, and the nonlocal attractive term is governed by an interaction potential G : R d → [0, +∞) satisfying G(x) = g(|x|) for some C 1 function g : [0, +∞) → [0, +∞) with g > 0. The latter assumption in particular implies that the nonlocal drift term in (1) is attractive, in that it yields a decrease in time of all moments |x| p ρ(x, t) dx of any order p ∈ [1, +∞). As a consequence of that, ρ is expected to concentrate to a Dirac's delta centered at the initial center of mass as t → +∞ in case ε = 0. On the other hand, the diffusion part exerts an opposed effect on ρ, in that it yields a growth of all moments, and it would lead the solution to zero as t → +∞ in case G ≡ 0. Models of the form (1) can be recovered not only in biology as mentioned above, but also material sciences (e.g. granular media models [40] and crystal defects modeling [36]) and in social sciences, see for example pedestrian movements [27] and opinion formation modelling [38]. The latter applied discipline in particular justifies the derivation of the PDE (1) from a discrete set of equations, which can be stochastic or deterministic depending on the applied context. The mathematical theory for (1) heavily relies on the associated energy functional defined on the space of probability measures P(R d ) (the total mass is preserved in time in (1)) with finite L 2 -norm, with the obvious extension to +∞ when ρ ∈ P(R d ) \ L 2 (R d ). The energy functional E can be used to interpret (1) as a Wasserstein gradient flow in the sense of [1,26].
The functional E is λ-displacement convex in the sense of [33,1] provided G is 2λ convex as a function on R d for some λ ≤ 0. The existence (and uniqueness) provided by the so-called JKO scheme [26] is global-in-time in L 2 provided G is smooth enough. Unlike the 2d Keller-Segel model for chemotaxis (see e.g. [25,8]), no concentrations to Dirac's deltas occur in finite time if G has (for example) two continuous and bounded derivatives. The case ε = 0 has been studied by many authors [30,10,11,22,4,5,6].
As solutions cannot blow up in finite time when G is smooth, the dichotomy between the repulsive diffusive effect and the attractive effect induced by the nonlocal term is resolved in terms of the large time behavior: a diffusive behavior would result into the large time decay of solutions as t → +∞, for example w.r.t. the L 2 norm, whereas the aggregative behaviour would yield the formation of a pattern, i.e. a (nontrivial) steady state. Such a problem has been studied mainly in two papers in the past years. In [2] the problem has been addressed at the level of the existence of nontrivial minimizers for the energy functional E[ρ]. Almost at the same time, [12] addressed the existence and uniqueness of nontrivial steady states for (1). The combination of the two results provide the following general picture in the special case in which G is nonnegative and integrable: under the constraint ρ dx = M > 0. Moreover, in case G is supported on R, the equation (1) has a unique nontrivial stationary solution up to multiplications by a constant and translations, which coincides with the global minimizer for E. 2 Partial extensions of the previous results to the case of a porous medium term divρ∇ρ m−1 were proven in [13,19,28]. Let us mention at this stage that an extensive (and very deep) literature has been produced for the existence of global minimizers of the energy E in which the interaction potential is not-integrable and featuring confining property at infinity. For the case of Newtonian potentials we mention [32,31,29,15]. For more general kernels we mention [7,14].
Once the existence or non-existence of nontrivial steady states is clear, a natural question arises on whether or not those steady states are attractors for the semi-group (1). Even in the case of potentials with confining properties, the situation is not completely clear except in the case of "smooth" power laws G(x) = |x| γ , with γ ≥ 1, see [17]. We refer to the introduction of [14] for a very clear and detailed explanation. To our knowledge, the only result which deals with this issue is [3], in which the large time decay of solutions is proven for initial data with large second moment and large mass, with potentials essentially behaving as the Bessel potential in R d and in high enough dimension. The problem of detecting the large time behavior seems more difficult in the case of (1), as the quadratic homogeneity in the energy (2) does not allow to play with the initial mass in order to penalize one of the two terms in the right-hand side of (1). Moreover, having to do with an integrable kernel under pretty general assumptions does not allow to use the homogeneity of the kernel and rescale the equation, as often done in the case of confining interaction potentials.
In our paper we first complement existing results on the existence of steady states for (1) and then address the problem of the large time behavior of (1). Our results are restricted to the onedimensional case. More precisely, • In the case of the attractive Morse potential G(x) = 1 2 e −|x| , we find an explicit formula for the unique nontrivial stationary state to (1) with ε ∈ (0, 1), with fixed mass and center of mass. The result is contained in Section 2. Let us recall that the Morse potential is not included in the results proven in [12,28] because of the discontinuity of its gradient at x = 0.
• We address the large-time behaviour in the case of smooth potentials, i.e. with G ∈ C 2 (R), and with unit mass for simplicity. We first prove the unique steady state provided by [12] in the case ε ∈ (0, 1) is locally asymptotically stable in the 2-Wasserstein distance. The result is stated in Theorem 3.1. • We prove in Theorem 4.1 that all solutions with finite energy E decay to zero locally in L 2 and almost everywhere in x ∈ R as t → +∞ in case ε ≥ G L 1 . • In section 4 we provide some (incomplete) arguments suggesting that the large time decay occurs also in the case ε ∈ (0, G L 1 ) for suitable initial conditions. • In section 5 we support our conjecture that (1) has indeed a multiple behaviour for large times (i.e. that there are more than one attractor for (1) as t → +∞ with fixed ε ∈ (0, G L 1 )) with some numerical simulations.
2. An Explicit Formula for G(x) = 1 2 e −|x| In this section the problem we consider is a special case of (1) in one space dimension where the kernel G is given by the one-dimensional Morse potential and ε ∈ (0, 1) is a fixed constant. Stationary solutions of equation (3) have been studied deeply in [12]. In particular, for ε < 1, the uniqueness of solutions with unit mass and zero center of mass was proven for the stationary version of (3) The unique (up to translation and multiplication) solution detected in [12] is symmetric, monotonically decreasing on x > 0, and its support is a bounded interval in R, namely supp[ρ] = [−L, L] (see [12,Theorem 4.13]) for some L = L(ε) > 0 depending monotonically on ε ∈ (0, 1), with limits L(0) = 0 and L(1) = +∞. One of the most important tools related with the study of the stationary version (5) is the energy functional Our aim in this section is to find an explicit formula for the unique solution of (5) where the aggregation kernel is given by the Morse potential (4). We start by recalling that, for compactly supported ρ, the equation (5) is equivalent to requiring for some constant C ∈ R. Moreover, this constant depends on ρ and satisfies C[ρ] = 2E[ρ] (see [12,Lemma 4.4]). Differentiating twice the aggregation kernel (4) in sense of distributions we get where δ 0 is the Dirac delta at zero. Now differentiate twice the equation (7) w.r.t. x ∈ supp[ρ], and substitute G = G − δ 0 to get Thanks to (7) we can replace the convolution G * ρ by ερ − C, so that (8) becomes We end up with the following boundary value problem where the boundary conditions come from the fact that (due to the results in [12]) we are looking for steady states with connected compact support, i.e. with supp[ρ] = [−L, L], that are C 1 with ρ(−x) = ρ(x), and the latter implies that x = 0 must be a stationary point for ρ. Steady states also have to satisfy that C = C[ρ] = 2E[ρ], ρ ≥ 0, and ρ L 1 = 1. Now solving the ODE in (9) using the variation of constants formula we get the following expression without imposing the boundary conditions: Applying the boundary conditions to (10), we obtain the actual solution Now we have to determine uniquely the values of the two constants L and C so that (11) gives a unique solution. The function ρ in (11) has to satisfy the two extra conditions The first condition (12) gives Computing the integral 2E[ρ], we obtain Taking into account that C = 0, from (13) we get Substituting (14) in (15) we obtain the following equation for L tan(εL) +ε tan(εL) −ε 2 + (tan(εL) +ε)e −2L = 0 Before solving (16), we note that the assumption ρ(x) > 0 on (−L, L) provides a specific interval for L. Indeed, recalling that C < 0, ρ(0) = 0 which is equivalent to cos(εL) < 0, which gives On the other hand, ρ should be positive everywhere on (−L, L), and the only integer that ensures this to happen is m = 0. We therefore obtain the condition π 2ε < L < π ε .

Local Stability of Steady States for Smooth Attractive Potentials
In this section we study the long-time behaviour of the solution to our one-dimensional evolution equation (3). We do all analysis on a new formulation of the evolution equation (3) obtained by using the Wasserstein metric in one dimensional space.
We consider the equation where the unknown ρ(., t) is a time-dependent probability density on R, ε is a fixed constant in (0, G L 1 ), and G is the aggregation kernel with the following assumptions (1) G ≥ 0, and supp(G) = R, lim r→+∞ g(r) = 0. In the following we review briefly the Wasserstein metric and we shall see how one can reformulate (21) using a simplified expression of the p−Wassertien distance obtained in one space dimension written in terms of pseudo inverses of the cumulative distributions of some probability measures.
Let P(R d ) be the space of the probability measures on R d . We denote by P p (R d ) the space of probability measures µ ∈ P(R d ) having a finite p−moment R d |x| p dµ(x) < +∞. Then for two probability measures µ 1 and µ 2 in P p (R d ), the p−Wasserstein distance between them is defined by where Π(µ 1 , µ 2 ) is the space of all measures π on the product space R d × R d having µ 1 and µ 2 as marginals, i.e.
In the one dimensional space the p−Wasserstein distance can be rewritten by a different expression than (22) which makes a lot of simplifications in the analysis. In particular, the p−Wasserstein distance between µ 1 and µ 2 can be written in one space dimension as the L p − difference between the pseudo inverses of the cumulative distributions of µ 1 and µ 2 . More precisely, let F i : R → [0, 1], i = 1, 2, be the distribution function defined by Then the pseudo-inverse function u i : Using these notations, the p−Wasserstein distance between µ 1 and µ 2 has the following new expression Following the procedure in the seminal papers [18,37], we now rewrite our evolution equation (21) in terms of the pseudo inverse function. Assume for simplicity the solution ρ(t) to (21) is smooth, positive, and has a connected compact support. Let F (t) be its cumulative distribution function. Then one can easily show that its inverse u(t) : [0, 1] → R satisfies the following equation Then, thanks to the expression (25), we can study the p−Wasserstein distance between the generic solution and the steady state of (21) by direct computations on the time evolution of the L p distance of the difference between the solution and the steady state to the pseudo inverse equation (26). We consider the time evolution of the Wasserstein distance (25) between a generic solution and the stationary one of our equation (21). In particular, we will prove that the 2-Wasserstein distance between the generic solution and the steady state to our equation (21) goes exponentially to zero with respect to time under some assumptions on the initial datum and the steady state. Actually with these assumptions we are considering a special case where we take the initial datum to be close to the steady state, in particular we consider "confined" initial data.
Let ρ ∞ (x) denote the unique stationary solution to our equation (21) with unit mass and zero center of mass, and let u ∞ (z) denote the pseudo inverse of its cumulative distribution. Consider our equation (21) with an initial datum ρ(x, 0) = ρ 0 (x) with unit mass and zero center of mass, and let u 0 (z) be the corresponding pseudo inverse. Then let u be the corresponding solution to (26). The time evolution of the Wasserstein distance u(z, t) − u ∞ (z) 2k L 2k ([0,1]) where k ∈ N, is easily found to satisfy where, to simplify notation, we set Then since u ∞ is a stationary solution, it satisfies Using (29), (27) becomes (H(z, ξ)) − G (K(z, ξ)) dξdz with H(z, ξ) defined by (28) and Now, we analyze the two integrals I 1 and I 2 separately which implies that I 1 ≤ 0 since the function f (x) = 1 x 2 is decreasing on x > 0. Note that since the solution has a compact support, the boundary terms vanish in the integration by parts.
For I 2 , we first consider the case where k = 1. Using the same technique used in paper [30], we conclude the following identity in which we have used that for any two solutions u and v of (26) with zero center of mass. Moreover, using (34) we also have the following Next we prove that if we take the support of the steady state to be inside the region where the kernel G is strictly concave, assuming as well that the support of the initial datum is close to the steady state, then for all t > 0 the solution ρ(t) will be very close to the steady state ρ ∞ in the ∞-Wasserstein distance.
Proposition 3.1. Let ρ(x, t) be the solution to (21) having initial density ρ 0 ∈ L 2 (R) ∩ L 1 + (R) with unit mass and compact support. Let ρ ∞ (x) be the stationary solution to (21) with unit mass and same center of mass of ρ 0 . Let u 0 and u ∞ be the pseudo-inverse variables corresponding to ρ 0 and ρ ∞ respectively. Assume for some c > 0 and δ > 0 is arbitrarily small. Then for all t ≥ 0.
Proof. Assume by contradiction that there exists t * > 0 such that Since the L ∞ ([0, 1])-norm is the k → +∞ limit of L 2k ([0, 1])-norm, then for k 1, Now, for such values of k, let Due to (37),t < +∞. As a consequence, we get By (30) and from I 1 ≤ 0 we get Changing the role of ξ and z in I 2 t=t , we get so we have By (38), we get that for some ξ, z ∈ [0, 1] : Which implies that there exists ξ, z s.t. |u(z,t) − u(ξ,t)| > λ and this is because G is strictly concave on [−λ, λ]. This gives So we obtain the following strict inequality On the other hand, sincet is the infimum of all times in which u(t) − u ∞ L 2k ([0,1]) starts increasing, we have Finally, since (42) holds for k 1 and we know that which is a contradiction with (41).
We are now ready to prove our main result of local asymptotic stability of steady states.
Theorem 3.1. Let ρ(x, t) be the solution to (21) having initial density ρ 0 ∈ L 2 (R) ∩ L 1 + (R) with unit mass and compact support. Let ρ ∞ (x) be the stationary solution to (21). Assume G ∈ C 2 and G < 0 on the interval [−λ, λ]. Suppose that there exists δ > 0 such that Then, there exists a positive constant c > 0 such that for all t ≥ 0.
Proof. In view of (30), (32), and (33), we have Now by condition (i) we have that |K| < λ 2 − 2δ. Moreover, the condition (ii) and the result in Proposition Then the assertion follows by the Gronwall lemma.

Large time decay
The nonexistence of a nontrivial steady state for ε ≥ G L 1 is reasonably seen as a consequence of a stronger impact of the diffusion term on the dynamics of (21) compared to the case 0 < ε < G L 1 . For this reason, in the case ε ≥ G L 1 we expect solutions to decay to zero for large times in a diffusive fashion, similarly to what happens to the solution to the Cauchy problem of the porous medium equation, see [41]. This task is one the goals of this subsection.
On the other hand, there is another question which naturally arises in the case 0 < ε < G L 1 , that is whether or not solutions may exhibit a diffusive behaviour also in the case 0 < ε < G L 1 . As specified in the introduction, this question is of difficult solution even in classical problems such as the (modified) Keller-Segel system with linear diffusion, see [3]. As we will show later on, numerical simulations suggest that solutions may decay to zero provided the initial condition is "spread" enough, i.e. has large enough second moment. In this section we will produce some incomplete mathematical arguments that support such a conjecture.
We start by proving the decay of the solution to our equation (21) in case ε ≥ G L 1 .
Therefore, up to a subsequence ρ k := ρ(t k ) we know that I[ρ k ] → 0 as k → +∞. This means that there exists some constant C 1 > 0 s.t. I[ρ k ] ≤ C 1 . Expanding the integral I[ρ], using ρ ∈ P(R), and integrating by parts we get From the assumptions on G, (49) implies that there exists a constant C such that uniformly w.r.t. k. Using Nash inequality in one space dimension, namely with f = ρ 3 2 k in (51), and using standard L p interpolation inequalities we obtain uniformly w.r.t. k for some constant C 2 > 0. This means that ρ 3 2 k is uniformly bounded in H 1 . By Sobolev embedding, up to a subsequence ρ This shows that I[ ρ] = 0, which implies that ρ(ερ − G * ρ) 2 x dx = 0 and so (ερ − G * ρ) x = 0 on the support ofρ. This means thatρ is a steady state. Since we are in the case ε > G L 1 where we have no steady state but zero (see [12]), we have ρ = 0. As a trivial consequence of all the above procedure, every family ρ(·, t k ) with t k → +∞ has a subsequence that converges to zero almost w.r.t. x ∈ R. By a.e. uniqueness of the a.e. limit, the whole family {ρ(·, t)} t≥0 converges to zero as t → +∞.
We now focus on the case 0 < ε < G L 1 . We consider two arguments suggesting that the large time decay holds also in this case for suitable initial conditions.
It is well known that the one dimensional porous medium equation Figure 1. Two different asymptotic behaviours of the solution depending on the initial datum.
can be approximated as N → +∞ by the empirical measure of the N -particle system see for instance [23]. As a toy model, we therefore consider the following two particle approximation of (21), Assuming the two particles occupy symmetric positions w.r.t. zero, i.e. X 1 = −X 2 = −X, then (54) becomes the one ODE Assuming that −G is zero at x = 0, it increases on an interval [0, ] and it decreases to zero on ( , +∞), for small ε it is easy to show the existence of two equilibrium points a, b for the ODE in (55), a stable and b unstable. Hence, for the Cauchy problem (55) one can show the following behaviour: =⇒ lim t→+∞ X(t) = a. This means that the behaviour of the discrete density varies for different choices of the initial datum. If we take X 0 > b then X(t) increases, which means that the support of the density ρ(t) will increase to +∞ as t → +∞ i.e. the density ρ(t) decays to zero, see Figure 1.
As a second argument in favour of the large time decay for suitable initial conditions, we get back to the continuum equation (21) and consider initial conditions of the form We compute the time derivative of the second moment as follows: By substituting the explicit initial condition ρ R , after few computations we get we recall that zG (z) ≤ 0 on z ≥ 0. Hence, under the simple assumption that there exists a range of values R ≥ R 0 , for some R 0 > 0, for which the second moment of the corresponding solution has a positive derivative at t = 0. Clearly, this computation does not solve the problem, but it shows nevertheless that there are initial data (with large enough second moment) producing a dominant repulsive behaviour for short times.

Numerical Simulations
In this section we shall present two numerical simulations. In the first simulation we validate our result in the section 3 using the finite volume method introduced in [16]. The second simulation shows the large time decay to zero of the solution to (21) using the particle method already sketched in the previous section.
We begin with the first simulation where we apply a 1D positive preserving finite-volume method for (21), see the paper [16]. Divide the computational domain into finite-volume cells denote the averages of the solution ρ computed at each cell U i . Then integrating the equation (21) over each cell U i , we obtain a semi-discrete finite-volume scheme given by the following system of ODEs for ρ i where the numerical flux F i+ 1 2 is an approximation for our continuous flux −ρ(ερ − G * ρ) x . We obtain the following expression for F i+ 1 and where the minmod limiter is given by Finally, we integrate the semi-discrete scheme (56), which is a system of ODEs, numerically using the third-order strong preserving Runge-Kutta (SSP-RK) ODE solver used in [24]. We apply this method to our equation (21) where we take the steady state inside the interval where G is concave and we are able to do that by experimenting our equation (21) for different values of ε and actually we get that for small ε. Then we take the initial density close to the steady state and we get our result, see Figure 2.
In the second simulation we consider a particle method (see [23]) in which we approximate the PDE (21) by a system of N particles X 1 (t), . . . , X N (t) with equal masses m i = 1 where Then we solve the particle system (61) using the Runge-Kutta MATLAB solver ODE23. Note that the initial mesh sizes are automatically determined by the total number of particles N and the initial density values. We take the initial positions X(0) = X 0 = (X 1 0 , X 2 0 , . . . , X N 0 ) s.t.
We apply this method to our equation (21) where we take the initial density ρ 0 with large second moment and we get the decay to zero for large times, see Figure 3. This behaviour supports our conjecture that a multiple behaviour holds for (21) in the aggregation dominated regime, namely 0 < ε < G L 1 .