Stationary States and Asymptotic Behaviour of Aggregation Models with Nonlinear Local Repulsion
We consider a continuum aggregation model with nonlinear local repulsion given by a degenerate power-law diffusion with general exponent. The steady states and their properties in one dimension are studied both analytically and numerically, suggesting that the quadratic diffusion is a critical case. The focus is on finite-size, monotone and compactly supported equilibria. We also investigate numerically the long time asymptotics of the model by simulations of the evolution equation. Issues such as metastability and local/ global stability are studied in connection to the gradient flow formulation of the model.
The derivation and analysis of mathematical models for collective behaviour of cells, animals, or humans have been receiving increasing attention in recent years. In particular, a variety of continuum models based on evolution equations for population densities has been derived and used to describe biological aggregations such as flocks and swarms [20, 28, 29, 32]. A typical aspect of these models is the competition of social interactions (repulsion and attraction) between group individuals, which is also the focus of current research.
In this paper we add novel results on the study of a canonical model for the competition of attraction and repulsion, namely the following one-dimensional aggregation equation for the population density :
Here, is an attractive interaction potential (to be detailed in Section 2.2), is a diffusion coefficient and is a real exponent. Equation (1) falls into the general class of aggregation equations with degenerate diffusion in arbitrary dimension :
The left-hand-side of (2) represents the active transport of the density associated to a non-local velocity field . The potential is assumed to incorporate only attractive interactions among individuals of the group, while repulsive (anti-crowding) interactions are accounted for by the nonlinear diffusion in the right-hand-side. Alternatively, by transferring the nonlinear diffusion to the left-hand-side, one can regard equation (2) as active transport of that corresponds to a velocity field that has a non-local attractive component, , and a local repulsion or dispersal part, . Nonlinear diffusion terms have been suggested in several instances for dispersal and repulsion (cf. e.g.  and [28, 32] in the above context). From a microscopic point of view, the case () can be easily justified, by taking the repulsive force modelled by a potential, similar to the aggregative one, and then performing a scaling limit as the interaction range approaches zero (cf. ). In Section 2 we present a microscopic derivation of the model with arbitrary based on nearest-neighbour interactions for the repulsion, which provides a unified interpretation of the nonlinear diffusion in terms of local forces. Models of type (1) also appear in earlier works on population dynamics [22, 23, 30], but the potentials considered there have different properties than the ones considered in this article.
Regardless of the interpretation, there is a delicate balance between attractive and dispersal effects which results in very interesting (and biologically relevant) dynamics and long-time behaviour of solutions to (2). Well-posedness of solutions to (2) has been studied intensively, we refer for instance to [8, 10, 11]. In particular, the analysis in  takes advantage of the formulation of these models in terms of gradient flows on spaces of probability measures equipped with the Wasserstein metric (cf. ). Also, a wide literature exists in relation to the Keller-Segel model for chemotaxis (see [5, 9] and references therein). As pointed out in such works, existence theory is more delicate in the presence of singular kernels, where finite time blow-up of solutions is possible [5, 24].
where . The special case (1) corresponds to power-law diffusion
Stationary states of (2) are critical points of the energy (3). A recent work by Bedrossian  investigates the existence of global minimizers of (3) using calculus of variations techniques. In particular, the existence of a radially symmetric and non-increasing global minimizer can be inferred for power-law in (4) with . The case is critical and yields a global minimizer only for small enough diffusion coefficient , where the threshold value for is shown to be by Burger et al. in . Energy considerations have also been employed in  to study the large time behaviour of solutions to (2) in one dimension.
This paper considers the power-law diffusion (4), though the results are expected to be true for any strictly increasing function . Previous works considered various specific values of the exponent . In  the case is studied and stationary solutions with compact support in one spatial dimension are characterized. The study in  makes extensive use of the interplay between energy minimizers of (3) and equilibria of (1). Topaz et al.  investigate the case in one and higher dimensions, where the focus, again, is the characterization of the steady states. The authors find compactly supported steady states with steep edges which they call “clumps”.
The results from  and  are the main motivation for this paper. We investigate the existence of finite-size, compactly supported stationary states (“clumps”) for general power exponent . Such biologically relevant equilibria have also been sought and studied in aggregation models with nonlocal repulsion [7, 25]. We show that such equilibria exist for regardless of the size of the diffusion coefficient . Analytical results and numerical experiments suggest that these steady states are the global minimizers of the energy, identified by Bedrossian in . However, the convergence in the aggregation dominated regimes could be arbitrarily slow, leading to metastable dynamics, as already observed in . The case is more subtle, as dynamics depends on the size of , as well as on the spread of the initial data. More precisely, for values of above a certain threshold (which we can identify in a nonlinear eigenvalue problem), diffusion dominates and spreading to a trivial equilibria occurs. For values of below this value, clumps exist, but have only a limited basin of attraction. As shown in the numerical experiments, the spread of the initial data is important in determining whether the solutions approach these stationary states or disperse to infinity.
The paper is organized as follows: in Section 2 we present briefly several derivations of the model (2) and discuss its basic properties. Section 3 is devoted to rigorous analysis of the existence of compactly supported stationary solutions to (1). These stationary states are computed with an iterative scheme in Section 4. In Section 5, the long time asymptotic behaviours of solutions to the evolution equation (1) are investigated numerically.
2 Model Derivation and Basic Properties
2.1 Model Derivation
In the following we discuss two natural derivations of the general model (2).
We consider a system of particles at positions , , in one dimension, with two kinds of interactions: a long-range attraction, which we assume to be in a standard additive form, and a nearest-neighbour repulsion. For simplicity we assume that . Consider the model
where denotes the local repulsion potential and the global attraction potential. The repulsive interaction only concerns nearest neighbours and is scaled via to obtain a meaningful locality. Note that we can also rewrite (5) as a gradient flow for the energy
In order to derive a continuum model, consider a discretization of the pseudo-inverse of the cumulative density distribution function , with and . Then (5) is equivalent to
The limit naturally yields
Now we immediately obtain a microscopic interpretation of (2), by noticing that we can write (6) in this form, provided , or equivalently, . In particular, the choice can be interpreted as the limit of a microscopic nearest-neighbour repulsive potential
or a repulsive force proportional to . It is obvious that higher exponents in the nonlinear diffusion correspond to stronger local repulsion, and we also observe an interesting demarcation at . For the microscopic repulsion is weak, i.e., it has an integral potential, while it becomes non-integrable for . We shall see different properties of the model for and on several instances in our analysis.
For an alternative derivation, which is quite standard for nonlinear diffusions, we consider directly the macroscopic compressible Euler equations with linear friction (friction coefficient scaled to one) and an additional nonlocal force, i.e.
For a diffusive scaling of macroscopic type, i.e. , , and , the left-hand side in the second equation is of order , while the right-hand side is of order . Hence we obtain the leading order terms
and inserting into the first equation yields (2) with the relation .
2.2 Gradient Flow Structure
Assumptions on the kernel . Given the focus of this paper, we list the properties of G only in one dimension. Throughout the paper the interaction kernel is assumed to satisfy
for all ,
for all , .
Kernel having infinite support (Assumption 1) means that the attractive interactions are global. This is an important assumption used to conclude that the support of a steady state is connected. Assumption 2 concerns regularity properties on which are needed to pass derivatives inside the integral operator , as well as in other estimates. In particular, Sobolev embeddings imply . The assumptions on could be potentially relaxed, but providing such sharp conditions is not a purpose of the present paper. Note that pointy potentials such as are included in the present theory. Assumption 3 is symmetry (isotropic interactions) and Assumption 4 means that kernel is purely attractive and interactions decay at infinity.
Using the recently developed techniques on gradient flows in metric spaces, in particular in the space of probability measures equipped with the Wasserstein metric (cf. ), it is straightforward to analyze the well-posedness of the dynamic model. The evolution equation (2) can be written in the form
2.3 Basic Properties of Stationary Solutions
The central issue in this paper is to understand stationary states of (1) and their implications for the long time behaviours of the dynamics. Using the more general form (2) in one dimension, equilibria satisfy a.e. in the equation
Due to the gradient flow formulation of the equation, such equilibria are closely linked to critical points of the energy. In particular, minimizers of the energy functional (3) on the manifold of probability measures are stationary solutions. We present here a few facts about the equilibria of (2) which can be derived relatively easily by generalizing results from . These facts are not directly used in the paper, but they provide strong motivation for our subsequent studies. We do not present their proofs, except for the last one (see the Appendix). Below, denotes the space of non-negative integrable functions of mass one. The nonlinear diffusion function is assumed to be monotone, for , and .
A stationary solution of (2) that has connected support, satisfies
The two facts above may be shown in fact for general dimension. The remaining ones apply to one dimension only.
A stationary solution of (2) in one dimension has connected support.
A minimizer for the energy (3) under the constraint that the center of mass is zero, is necessarily symmetric and monotonically decreasing on by Riesz rearrangement inequality.
Consider a stationary solution of (2) with compact support. Then there exists a symmetric stationary solution such that
which is positive for and negative for .
Consider the possibility of a stationary solution with having infinite measure. For we can conclude from the finiteness of the energy that and thus also . Hence,
which implies . Thus, for stationary solutions with unbounded support. In the case it is not even clear that whether , hence there could be stationary solutions with positive value of the entropy.
Motivated by the facts above, we focus in the rest of the paper on analytical and numerical investigations of equilibria of (1) that have compact and connected support.
3 Compactly Supported Stationary States
Based on previous works [32, 12] and our own numerical investigations, there is a class of steady states which is particularly relevant and important for the dynamics of (1). We refer here to symmetric, monotone and compactly supported steady states. Note also that solutions to (1) preserve mass, that is, remains constant during the time evolution. Hence, equilibria can be considered of fixed mass, given by the initial configuration.
We are interested in finding a symmetric (even) density that has finite support , vanishes on the boundary of the support, and has unit mass. In mathematical terms we look for a solution of the equation (see (8)):
where vanishes outside and . In particular, the continuity condition implies that .
In formulating the integral equation for the steady state, we included the requirement that solutions have unit mass. This is a natural normalization that is inherited from the dynamics of the model. However, in the existence of steady states investigated in this section and for the alternative numerical calculations of the Bessel potential in Section 4.3, it is more convenient to work with the normalization instead. Hence, depending on the context, we may work with any of the two normalizations, and transfer results from one normalization to the other using the following simple scaling argument.
We present the scaling argument for one dimension, but it is straightforward to see that it holds identically in any dimension. Note that if is a solution of the steady state equation (10) that has the desired properties, then
and consequently, satisfies
We can choose to pass from one normalization to another. Note that the support is fixed through the transformation. If in (11) has unit mass, we take , and in (12) has the normalization . On the other hand, if in (11) has the normalization , then by choosing , in (12) has unit mass. Note however that, in making this transformation, the diffusion coefficient changes as well.
3.1 Existence of Stationary Solutions in One Dimension
Using the symmetry of the solution and the fact the vanishes at , equation (10) can be written as
Case (when the problem becomes linear) was investigated in . The approach there was to take a derivative in (13) and transform the equation into an eigenvalue problem for . The main advantage is that the kernel of the resulting integral operator is positive and enables the use of the Krein-Rutman theorem, while the kernel of defined in (14) is not. We follow a similar approach here, in dealing with a general exponent . However, the problem is nonlinear now, and the methodology from  has to be adapted and significantly extended.
Take a derivative in (13) to get:
Although can be infinite at (true for ), is integrable and the equation (15) is still well-defined.
In the following we shall prove the existence of compactly supported stationary states with arbitrary support size by a combination of the Krein-Rutman theorem with fixed point techniques. The first step is to “freeze” and study the eigenvalue problem (15) in . Krein-Rutman theorem (strong version) can be used to argue that a unique nonpositive eigenfunction exists, and hence we can define a map that takes into . The next step is to consider a second map, which maps into its primitive, and show that the composition of the two maps has a fixed point. Such a fixed point is a solution of (15).
We point out the transition at for the study of solutions of (15) performed here. Since vanishes at the boundary, is either zero () or infinity () at . Consequently, since the right-hand-side of (15) is finite at , is either infinite () or zero (). We will treat the two cases separately.
We first state the strong version of the Krein-Rutman theorem used in the arguments.
Theorem 3.1 (Krein–Rutman Theorem, strong version).
Let be a Banach space, let be a solid cone, i.e. such that for all and has a nonempty interior . Let be a compact linear operator, which is strongly positive with respect to , i.e. . Then the spectral radius is strictly positive and is a simple eigenvalue with an eigenvector . There is no other eigenvalue with a corresponding eigenvector .
Another useful tool in the following is a comparison result on the leading eigenvalue, which can be derived easily from the Krein-Rutman theorem:
Let be a compact linear integral operator from with positive continuous kernel and a positive continuous function on such that . Then the leading eigenvalue (from the Krein-Rutman theorem) is greater or equal .
Assume that the maximal eigenvalue is smaller than . Since is an integral operator with positive kernel, its -adjoint operator is a positive operator as well and hence, the Krein-Rutman theorem also implies the existence of a nonnegative function such that
On the other hand, due to the positivity of we have
which is a contradiction. ∎
In the case we rewrite the problem (15) for the derivative in the form
In order to avoid potential issues with dividing by zero (due to ) we start with a regularized problem of the form
Our final goal is to prove the existence of a solution in the following subset of , the set of continuous functions on :
Note that is a bounded set, since any attains its maximum at and hence .
We first define the map that takes to the leading eigenfunction of the eigenvalue problem
The existence and uniqueness of the leading eigenvalue follows ideas from , where the strong version of the Krein-Rutman theorem was used to study the eigenvalue problem.
A major difficulty is that we cannot simply work in the space and the cone of nonpositive continuous functions, i.e. , since always vanishes at and hence it is not in the interior of . To avoid this issue one can work instead with the cone of nonpositive functions in the subspace of defined by functions with vanishing left boundary value. The control of the derivative and the favourable properties of the operator (which yields a definite sign of the derivative at ) help to apply the strong version of the Krein-Rutman theorem. More precisely we use (see also  for further discussion) the space
with the cone
As in  one can see that is now in the set
which is a subset of the interior of the cone . Then we can apply the Krein-Rutman theorem to obtain:
For given and there exists a unique maximal eigenvalue and a unique nonpositive with and satisfying (20). Moreover,
where is the maximal eigenvalue of the positive linear operator .
The existence and uniqueness of a maximal eigenvalue with eigenfunction as above follows from Theorem 3.1 using the cone . The eigenfunction has been normalized to have integral .
In order to obtain a lower bound on we employ Lemma 3.2 with being the principal eigenvalue of , which exists and is positive due to the Krein-Rutman theorem, and an associated nonnegative eigenfunction. Then
Hence, we conclude
Note that with the lower bound on we can estimate the supremum norm of the nonpositive eigenfunction using its normalized -norm. Indeed,
Due to the setup of the Krein-Rutman theorem we have to assume , which in turn affects the setup of the fixed-point argument, as it now needs to be carried out in rather than . Note that is counter-intuitive at a first glance since we expect stationary solutions of the original problem to have infinite derivative at in the case . However, the derivative will be finite for positive used in the fixed point argument, which is also encoded somehow in the fact that as . In the limiting procedure we will subsequently use a limit of solutions in and convergence in .
We construct a fixed-point operator in two steps on the convex set
Note that . For given define
where is the unique nonpositive eigenfunction for the leading eigenvalue of (20), which satisfies . Then consider the map
Therefore, by the normalization of and the above estimate on its supremum norm, solutions of (18) are fixed points of the map .
We now investigate the properties of the two maps and .
The operator is well-defined and continuous. Moreover, is nonpositive with and satisfies .
From Lemma 3.3 we obtain the well-definedness and the properties of . It remains to verify the continuity of . To do so we symmetrize the eigenvalue problem by the transform , which is obviously continuous on . Hence, it suffices to prove the continuity of the map . Then is the unique leading eigenvalue of the self-adjoint operator
which depends on the integral kernel, and consequently on , in a locally Lipschitz continuous way.
For given and , and , eigenpairs of , , respectively, we further have
Since the right-hand side is orthogonal to each element in the null-space of the operator (which consists only of multiples of ), the Fredholm alternative implies the existence of a unique solution , which depends continuously on the right-hand side (also in the supremum norm, cf. ). Hence we have for some constant depending on ,
For the first term on the right-hand-side we already know the local Lipschitz continuity on , while the local Lipschitz continuity of the second term is a straightforward computation. Finally, we can differentiate the eigenvalue equation with respect to and obtain an estimate in the norm of in an analogous way. ∎
The linear operator is well-defined by (23) and compact. Moreover, each function with and for maps to a function .
Putting these properties together we can show:
For each sufficiently small and , there exists a solution and of (18).
A density satisfies (18) if it is a fixed-point of the map , i.e.,
From Lemma 3.4 and 3.5 we conclude that is continuous, compact, and maps the convex and bounded set into itself. Thus, the assumptions of the Schauder fixed-point theorem are satisfied and we conclude the existence of a solution . The existence of then comes from the fact that is a Krein-Rutman eigenfunction of (20) and is obtained from the spectral radius: . Thus, we obtain for each a solution . ∎
The final task remaining is to send to zero. We proceed in several steps. As mentioned above, the derivative approaches infinity at the boundary as and hence is not suitable for performing the limit. For this reason we first integrate (18) to write the eigenvalue problem in terms of only. We write (18) as
and integrate from to to get
where was defined in (14).
Evaluating at , we can get an uniform upper bound on :
We are now ready to prove the main existence result for the case .
Theorem 3.7 (Existence of compactly supported equilibria for ).
For each , there exists and solving (13).
We start with a subsequence such that converges to some , which exists due to the above arguments. We employ Lemma 3.3 to obtain
which yields a uniform lower bound on as tends to zero. Subtracting the equation for (see (27)) evaluated at , respectively , yields
Using the smoothness of and the lower bound for , we estimate
for some positive constant independent of .
Hence, is uniformly Hölder continuous, in particular equicontinuous. With the Arzela-Ascoli theorem we deduce the existence of a convergent subsequence to satisfying (13) (the limiting equation of (27)).
In the case we use the transformation of variables from to before differentiating, i.e., we analyze the problem
For the existence proof we use a very similar strategy based on a fixed-point formulation and the equation for the derivative of interpreted as an eigenvalue problem; hence we only give a sketch of the proof.
The derivative satisfies
with and defined as in (16).
Define as the map from to the eigenfunction of the leading eigenvalue of (31), which is again normalized as a nonpositive function with and . Note that appears in (31) under the integral sign, hence enough regularity is obtained with being just continuous. Solutions of (31) are fixed points of the map , i.e.,
given by (24).
With an analogous proof as for we obtain
The operator is well-defined and continuous. Moreover, is nonpositive with and satisfies .
The map is the concatenation of and an embedding operator. Using the above properties of (see Lemma 3.5) we conclude that is continuous and compact and maps into itself. Thus we conclude from Schauder’s fixed point theorem the existence of and being the derivative of such that (31) is satisfied. Since is continuously differentiable and the exponent is positive, it is straightforward to justify integration of this equation and conclude the existence of and solving (30). Thus we obtain the following result:
Theorem 3.9 (Existence of compactly supported equilibria for ).
For each , there exists and solving (13).
We finally explain why we choose to work in this case with a transformed variable and not with the original formulation. In principle, a fixed point approach could be constructed as in the case , but it is more difficult to justify the way back from the equation for the derivative to the density . For we end up with as a factor of , and since we cannot justify its well-definedness in cases where . This however seems to be rather a technical issue. More severely is the fact that we find in this case. Thus, there is no argument that is positive close to and the fixed-point proof in the variable would only yield the existence of a solution whose support is a subset of . In the formulation with the new variable we know that is decreasing and is positive, so the support is exactly .
3.2 Further Properties of Compactly Supported Stationary Solutions
We finally discuss some finer properties of compactly supported stationary solutions. A first observation concerns the behaviour at . From (13)-(14), using the Lipschitz-continuity of we can find a constant depending on and only, such that
Equation above yields another change at . The solution is Hölder continuous for (as the Barenblatt solutions of the porous medium equation), while it is even differentiable, with vanishing derivative at , for . With , more and more derivatives at become zero.
Another interesting question is the potential concavity of the stationary state. For this sake we compute a formula for the second derivative by differentiating (15) on the support of :
Evaluate at the origin and use , to get
Since is positive for negative argument and vice versa we obtain, also using the negativity of , that . Thus, is locally concave at , which is not surprising since we expect to have its maximum there.
The range of concavity of a stationary solution depends on as well as on the concavity of the kernel . As already argued in , the right-hand-side in (34) is nonpositive if is concave on , i.e., if the support is smaller than half the range of concavity of . If this immediately implies that is negative on the support of , because is positive. In the case we already know that , which does not allow to be concave on the whole support.
4 Numerical Calculation of the Steady States
According to the results of the previous section, steady states exist for any fixed domain size . Now we turn our attention to the qualitative features of these equilibria, in particular to the dependence of the diffusion coefficient on and the limit . To this purpose, we develop an iterative method to solve the eigenvalue problem (15) with general interaction potentials. The ultimate goal is to characterize the long time behaviours of solutions to the evolution equation (1), and this approach is much more effective than solving (1) directly (see Section 5).
4.1 Iterative Scheme
When , the nonlinear eigenvalue problem (15) can still be solved iteratively, though not in one single step as for . The analytical study of the previous section suggests the following iterative scheme for the steady state. Fix and a nonnegative function (the initial guess) on . For , perform the following two steps:
Find by integration and normalization to unit mass:
In the limit when , converges to (depending on ) and converges to a solution of (10).
To implement the algorithm numerically, we set an uniform grid on . We assume that is piecewise constant on the grid, that is, on , . Given the piecewise constant function on , in Step 2 is piecewise linear. To avoid however a possible singularity at (see Section 3.1 for ), in Step 1 we evaluate the density at middle points , that is, . Hence equation (36) evaluated at then reads (with subindex omitted for convenience)
where are approximations (via trapezoidal rule) of With known from the previous iteration, the leading eigenvalue problem (37) to find and can be solved with standard algorithms. We terminate the algorithm when the change in (calculated in the supremum norm) is smaller than .
Uniqueness of equilibria is not guaranteed by the analytical considerations from Section 3.1. However, all numerical investigations we performed showed convergence of the algorithm described above to a unique solution of (10). Finally, we point out that the numerical algorithm also works for compactly supported kernels such as , which violate assumption 4 (strict monotonicity) in Section 2.2. For such kernels though, the steady states of the evolution equation (1) are not unique, as multiple disconnected and non-interacting bumps can exist.
4.2 Qualitative Properties of the Steady States
We start with the dependence of the diffusion coefficient on the domain size . Figure 1 corresponds to the Gaussian kernel , while qualitatively similar plots were obtained for other kernels such as the Bessel potential or the hat function . In the limit , goes to infinity when , while for , approaches a constant value. Both behaviours will be commented in detail below.
The steady states shown in Figure 2 exhibit qualitatively different behaviours too, depending on whether or . We discuss the two cases separately.
When (see Figure 2(a)), the steady states are decreasing and concave down, and spread and decay to zero when increases to infinity. Moreover the numerical solution has a very peculiar form in this limit: it is close to a constant on the support, with a sharp drop near the boundary at . Hence, for large, one can approximate the density profile by a rescaled characteristic function (to preserve unitary mass). Then ,
and the steady state equation (10) yields . This implies the scaling law
A similar result, but regarding the dependence of on large mass (with fixed) was derived in .
To summarize, when , for any diffusion coefficient , there exists a compactly supported steady state that solves (10). Numerical results indicate that this steady state is unique and simulations of the evolution equation (1) (see Section 5) indicate that these equilibria are global attractors for the dynamics. The larger the diffusion coefficient , the larger the size of the support, with scaling given by (38). Equilibria look like rescaled characteristic functions as (and ) become large.
Solutions in this case (see Figure 2(b)) are also decreasing, but change concavity and reach the boundary of the support with zero derivative (results consistent with analytical considerations of Section 3.1).
Also in contrast with the previous case, the steady states approach a fixed profile as goes to infinity. We can identify this limiting profile as follows. As , approaches a constant value (Figure 1) and the constant in (10) approaches zero (see discussion on equilibria of infinite support in Section 2.3). Combing these observations, we infer from (10) that the limiting profiles are governed by the nonlinear eigenvalue problem
In general, although no closed form expressions of the eigenpair are expected, (39) can be solved exactly for some special kernels like the Gaussian