Biological aggregation driven by social and environmental factors: A nonlocal model and its degenerate Cahn-Hilliard approximation

Biological aggregations such as insect swarms and bird flocks may arise from a combination of social interactions and environmental cues. We focus on nonlocal continuum equations, which are often used to model aggregations, and yet which pose significant analytical and computational challenges. Beginning with a particular nonlocal aggregation model [Topaz et al., Bull. Math. Bio., 2006], we derive the minimal well-posed long-wave approximation, which is a degenerate Cahn-Hilliard equation. Energy minimizers of this reduced, local model retain many salient features of those of the nonlocal model, especially for large populations and away from an aggregation's boundaries. Using the Cahn-Hilliard model as a testbed, we investigate the degree to which an external potential modeling food sources can be used to suppress peak population density, which is essential for controlling locust outbreaks. A random distribution of food sources tends to increase peak density above its intrinsic value, while a periodic pattern of food sources can decrease it.

1. Introduction. This paper has three primary aims. First, beginning with a partial integrodifferential equation model of biological aggregation, we show that a long-wave approximation yields a degenerate Cahn-Hilliard equation akin to models describing phase separation in materials science and evolution of thin fluid films. Second, we study solutions of the degenerate Cahn-Hillard model and demonstrate that they compare well with those of the original model. Using the Cahn-Hilliard model is advantageous because it eliminates many of the analytical and computational challenges of nonlocal equations, particularly in higher dimensions. Our methodology is to study the energy minimizers of this local equation, rather than studying the time evolution of the original nonlocal equation. Finally, we use the Cahn-Hilliard model as a testbed to assess whether imposing an external potential modeling the environment, e.g., food sources, can be used to control peak density in aggregations, which is of interest for locust swarms.
Biological aggregations such as bird flocks, fish schools, and insect swarms are driven by social interactions between group members. Typical social forces include attraction, repulsion, and alignment [33,37] which are preprogrammed by evolution and which are facilitated by sight, sound, smell, touch, or some combination of these senses. The accepted biological paradigm is that repulsion operates over shorter distances than attraction, but that it is much stronger over those distances. In the context of many mathematical models, these are necessary conditions to achieve cohesive groups with finite densities [51,53].
Because sensing operates over finite rather than infinitesimal distances, it is most naturally described via nonlocal operators. The most common mathematical descriptions assume that social interactions take place in an additive pairwise manner, and that the strength of interaction of organisms depends on the distance between them. We will use a kinematic model, where an individual's velocity is determined by social interactions dependent on the positions of other organisms; we will eventually augment this model to allow taxis due to environmental effects.
For a population that is well-described by a continuum density u(x, t), the social velocity v, can be modeled as v(x, t) = f (x − y)u(y, t) dy ≡ f * u. (1) Because the interactions are pairwise, additive and solely a function of position, the velocity can be written as the convolution of a social kernel, f , that encodes the influence of organisms at location y on those at location x. Such convolution-type continuum aggregation models date back to the seminal work of [49], which examineṡ in one spatial dimension. The left hand side of this equation is the material derivative of the population density u, that is, the rate of change of the population density moving at the velocity arising from social interactions. On the right hand side, D is a diffusion constant. The diffusion term was intended to model social repulsion operating over short distances, driving flux down population gradients. This model conserves mass, preserves positivity of the density, and can display an instability from a homogeneous state to patterned states. However, a key feature of biological swarms is that they are compactly supported with sharp edges, and (2) cannot reproduce this property. This is due to the linear diffusion; even compactly supported initial state are instantaneously smoothed and densities become positive throughout the domain. More recently, a population model that has received much attention is the aggregation equation,u + ∇ · (uv) = 0, v = −∇Q * u, where the velocity, v, can incorporate both attractive and repulsive social forces. The social kernel is written as minus the gradient of a social potential, −∇Q. One typically assumes that the social potential is directionally isotropic and symmetric (the social force between two organisms is equal in magnitude and opposite in orientation) which allows us to write the social potential as Q = Q(r), solely a function of the distance between organisms, r. Typically Q(r) is decreasing for small r to model repulsion, increasing over larger r to model attraction, and asymptotically flat for large r as organisms cannot sense each other at large distances.
The aggregation equation has recently been the subject of a rich and rapidly growing literature, including [10, 12, 13, 15-18, 22, 23, 39, 40, 46-48, 63, 68] and has also been adapted to describe specific biological organisms, such as locusts [66]. Crucially, (3) minimizes an energy, which can be used to analyze the existence and stability of equilibria [8,9,30,39,61]. Despite its popularity in the literature, one shortfall of this model is that it does not produce so-called well-spaced groups. In well-spaced biological swarms, individuals tend to pack only up to a maximum density, and no further. In contrast, for (3), if there is a minimizer for a particular mass M , doubling the mass yields a minimizer with double the density. This is unbiological, and results from the attraction and the repulsion scaling equally with density. An alternative aggregation model that does form well-spaced groups iṡ initially proposed in [65]. While (3) models attraction and repulsion nonlocally, (5) models attraction nonlocally via a kernel Q a , but models short-range repulsion via the (local) nonlinear diffusive velocity −u∇u. This model is similar to (2) but incorporates diffusion that is density-dependent and degenerate. That is, the diffusive flux tends to zero as the density tends to zero, which allows for compactly supported solutions. Eq. (5), like (3), allows the formation of compactly supported solutions [3,24,65]. It minimizes the energy, The nonlinear diffusion has two notable effects. First, because the repulsion dominates at higher densities, the density is bounded as the mass M increases, leading to aggregations of approximately constant density at large mass [24,65]. Second, the introduction of the gradient as part of the repulsion in the social velocity smooths the density profile at the support's boundary. Solutions are continuous with a discountinuous gradient, often referred to as a non-zero contact angle. Further work on (5), and on a generalization incorporating repulsion of the form −u p ∇u, has examined existence and uniqueness [4,5,21], regularity [32], global attractors [41], and properties of steady states as a function of the power p of the nonlinear repulsion [24]. Despite the promise of (5) as a model, it poses significant challenges. Beyond the analytical difficulties associated with the nonlocal operator, even numerical simulation is problematic. First, for 2-d simulation on an n × n grid, the convolution requires O(n 4 ) operations for quadrature or O(n 2 log n) for pseudospectral evaluation. Second, the steady states of (5) can have steep edges, like some steady states of (3). These require a high degree of spatial resolution. Third, standard methods produce oscillations emanating from the contact points, and such artifacts must be eliminated. Fourth, and relatedly, a numerical scheme must preserve the positivity of the model. It is notable but perhaps not surprising that fully two dimensional simulations of (5) are rare in the literature [29,65]. At present, most published simulations are either in 1-d or in a radially symmetric 2-d geometry. This is true of (3) as well, excepting special cases where the computation reduces to a boundary integral problem describing the motion of a self-deforming curve [17,64]. Given the attractive features of (5) from a modeling perspective, and yet given the analytical and computational challenges, we investigate whether there exists a simpler, local model that nonetheless preserves important features of (5).
In summary of this discussion, (5) produces groups whose density is bounded as mass M increases but the model is difficult to analyze and simulate. There is a long mathematical history of approximating dynamics in the limit of slow variation to simplify governing equations [35,59]. As we will describe in detail in this paper, an expansion of the convolution in (5) the long-wave limit yields, after nondimensionalization, Retaining only the first term in this expansion yields a nonlinear diffusion equatioṅ which is ill-posed due to negative (backward) diffusion at small densities. However, retaining the second term in the expansion yields stabilizing fourth-order degenerate diffusion. We will concentrate on this degenerate fourth-order equation, It is curious that for this model, well-posedness at short wavelengths comes from the fourth-order term derived from nonlocal attraction. As attraction operates at longer wavelengths than repulsion, one must retain higher derivatives to model it faithfully. Mathematically, the fourth-order term describes attractive forces that are responsive not only to variations in density but to the convexity of the density. This yields a social force analogous to a surface tension, and has the effect of damping variation in the density. A moment expansion similar to the approximation (7) is used to derive a different model for the movement of ecological populations in [52]. A great deal is known about (9). As we will discuss below, it is related to both the Cahn-Hilliard equation [56] which arises as a model of phase separation and to thinfilm equations which model surface-tension driven wetting of a substrate by a viscous fluid [34,55,59]. While the existence theory for our model is more strongly related to previous work on thin film equations, the dynamics are akin to phase separation and therefore we refer to our model as a Cahn-Hilliard model. Our model (9) also may be considered a degenerate Sivashinsky equation [57,62].
After nondimensionalization, (9) is equivalent tȯ which reveals a little more of its structure. It is the composition of a negative-definite second-order self-adjoint operator (∇ · u∇) with the first variation of an energy, This equation is an example of the more general class of Cahn-Hilliard models [56] which take the form,u Here, M (u) is a non-negative mobility, so ∇·M (u)∇ is a negative-definite second-order self-adjoint operator. The equation minimizes an energy where F (u) = f (u) du is most commonly taken to be a double well potential in the Cahn-Hilliard literature. Eq. (12) was first developed in [25,28] to describe binary alloys. A hallmark of this model is that it exhibits phase separation, where the density segregates into regions in which u assumes one of the two minimum values of F , and these regions are separated by narrow transition layers [26,27,43].
Our model (10) borrows a feature more familiar in the thin film context. Biologically, the density u must remain nonnegative, much like a fluid film thickness. In our model, this constraint is enforced by the degenerate mobility M (u) ≡ u. This feature not only ensures that the density remains nonnegative, but it also allows regions of zero density, akin to dry spots in fluid films. Degenerate mobilities have been considered in both the Cahn-Hilliard context [38,56] and in the thin films context [6,7,11,14,19,20,42,45,50].
The work in [58] nicely summarizes and extends the analytical progress that has been made over the past two decades on a class of thin film equations including (10), primarily in one spatial dimension. For our model in one dimension, [58] shows that there exist nonnegative solutions with compact and potentially disjoint support. Where u is positive, these solutions are infinitely differentiable. Moreover, these solutions are continuously differentiable everywhere. Stated differently, the slope at the contact point, where u touches down to zero, is zero. Additionally, the touchdown is generically quadratic and the speed of propagation of the contact points is finite. To our knowledge, theory for uniqueness and existence in higher dimensions is lacking. Still, motivated by these one dimensional results, we will primarily consider solutions that are infinitely differentiable for positive values and touch down with zero contact slope.
The rest of this paper is organized as follows. In Section 2, we perform a longwave expansion of (5) leading to the degenerate Cahn-Hilliard equation (9), augmenting both models with an external potential describing the environment. Section 3 calculates the linear stability of the two models, builds a framework for examining them as energy minimization problems, and describes how numerical computation of minimizers will be performed. Section 4 studies in detail the minimizers for the case of no external potential. In brief, minimizers of both models share many properties, including compact support and a common peak density in the large mass limit. With the Cahn-Hilliard model established as a reasonable approximation of the nonlocal one, in Section 5 we consider an external potential and show that when a potential well is sufficiently steep, it can be replaced with an approximate obstacle potential consisting of regions where the potential is zero separated by inhibited regions where the potential is effectively infinite. Finally, Section 6 uses these obstacle potentials to investigate how peak population density is affected by resource distribution. A random distribution of resources tends to increase peak density above its intrinsic value while a periodic distribution can decrease it.
2. Derivation of the local model. We begin with the nonlocal model of [65] augmented with an external potential W (x) which models taxis due to environmental factors,u where u(x, t) is the population density and x ∈ R n . The parameter R > 0 controls the strength of short-range repulsion. For now, we assume the potential W to be continuously differentiable, although eventually we will consider a discontinuous limit. The kernel Q a in the convolution above describes social attraction over finite (rather than infinitesimal) distances. We place several assumptions on this kernel: 1. Q a is radial, that is, Q a ≡ Q a (r), r = |x|. This is the case for social interactions that are isotropic in space. 2. Q a (r) ≥ 0 so that Q a describes (pure) attraction.
3. As |r| → ∞, Q a approaches a constant since organisms at very long distances do not interact with each other. Without loss of generality, we assume Q a → 0 as |r| → ∞. 4. Q a ≤ 0 as a result of the previous two assumptions. This is a dimensioned model; by rescaling x, t, and u, we may write a dimensionless version of (14) asu where we have chosen R = 1. Also, we choose our nondimensionalization such that Q a has unit volume, that is, where M 0 denotes the zeroth moment. In addition, we specify the second moment of the potential in order that Q a has a length scale of order unity. This condition is This rescaling, which we take without loss of generality, will be convenient for our subsequent calculations. Presently, we will approximate (15) with a local partial differential equation by applying the convolution theorem to −∇Q a * u and performing a long-wave expansion in Fourier space [51,54]. Define the Fourier transform of a function f (x), Since our kernel Q a is radial in physical space, its Fourier transform is radial in the transformed variable k. So where k = |k| and r = |x|. Expanding the complex exponential, and substituting into (19) will yield a moment expansion for Q a . We will retain only the first three terms of this expansion. For the first term, due to our choice of scaling. For the second term, because the term linear in k vanishes due to symmetry; in fact, all terms odd in k vanish. Finally, for the third (quadratic in k) term, we have Here, the right hand side of (23a) expands the square of the dot product as a double sum, (23b) utilizes symmetry to recognize the vanishing non-diagonal moments, (23c) uses radial symmetry of Q a , (23d) substitutes the definition of M 2 , and (23e) uses the normalization of M 2 in (17).
Putting these results together, and therefore, we approximate which constitutes a long-wave approximation. With this approximation, the original governing equation (15) becomes the new, local equatioṅ This equation may be rearranged aṡ which is a Cahn-Hilliard equation of type (12) with degenerate mobility M (u) = u, f (u) = u 2 /2 − u, and an external potential W (x). We will explore the degree to which minimizers of this truncated model approximate those of the nonlocal model (15). While the former are easier to compute owing to the purely local interaction rules, it is important to realize what is lost in the local approximation. In (15), disjoint clumps of density interact with each other exponentially weakly owing to the nonlocality. However, in (27), disjoint clumps do not interact at all owing to the local behavior and degenerate diffusion.

Basic model characteristics.
To recapitulate, we are studying the nonlocal modelu and the degenerate Cahn-Hilliard equation which is a long-wave approximation to it, where u(x, t) is a non-negative density. To complete the formulation of the problem, one must specify the domain, the boundary conditions, the external potential W and initial conditions for the density u. For much of this paper we consider all of R n or periodic domains in one and two dimensions. Solutions to problems with degenerate diffusion often have compact support [2,42,69]. For the Cahn-Hilliard problem (29) we expect u(x, t) to be continuously differentiable everywhere and infinitely differentiable for all x in the support of u. For the nonlocal model (28), we expect u(x, t) may have jump discontinuities at the edge of its support, but will be infinitely differentiable for all x in the support of u [8,9]. Leveraging the compact support of the solution, when we wish to consider densities confined to a finite domain rather than specifying boundary conditions at the edge of a domain, we will either consider a large periodic domain or we will specify that W (x) is sufficiently large to drive the density to zero except in some subdomain of interest. Both models are conservation laws which are well known to conserve total mass (see, e.g., [65]). If where Ω is the domain (or the support of the solution) then For many problems with a reflection symmetry, including our models (28) and (29), the location of the center of mass will also be conserved when the domain is R n and the mass is finite. In this case, for both (28) and (29), if we define a straightforward calculation shows that from which we deduce that in the absence of an external potential, the center of mass will be stationary. In the presence of an external potential, mass tends to migrate towards minima of the potential.
For the nonlocal model (28), a result also stated in [65]. The Fourier transform of the attractive kernel Q a depends on spatial dimension. From our nondimensionalization in Section 2, This guarantees that the dispersion relations for (28) and (29) are identical to O(k 4 ).
To further analyze the model one must specify Q a . A commonly used example which we also adopt is the Laplace potential. In our nondimensionalization, this potential is where in which case which is a power of the Lorentzian function. We see that Q a (0) = −1 and Q a (k) is increasing in k, and tends to zero. Consequently, ifū > 1, σ(k) < 0 for all k > 0 and the homogeneous steady state is linearly stable. However, ifū < 1, there is a longwave instability. There exists a finite, continuous band of wave numbers extending from k = 0 to for which σ(k) > 0. For the Cahn-Hilliard model (29), the calculations are significantly simpler. We obtain The stability theory is analogous to the nonlinear model. Ifū > 1, σ(k) < 0 for all k > 0 and the homogeneous steady state is linearly stable. However, ifū < 1, there is a long-wave instability. There exists a finite, continuous band of wave numbers extending from k = 0 to k = √ 1 −ū for which σ(k) > 0. It is worthwhile to consider why a fourth-order truncation of (28) is superior to a second-order one. For a second-order truncation, the dispersion relation is simply The stability threshold is stillū = 1 as above. However, in the linearly unstable regime, all modes have positive growth rates which grow unboundedly with k, similar to the backwards heat equation. Thus, the linear problem for the second-order truncation is ill-posed in this case, shedding light on why an approximation to fourth order is desirable. Our fourth-order truncation is analogous to the strategy used in deriving amplitude equations for marginally long-wave unstable systems via modulation theory [35]. Retaining the destabilizing second-order and stabilizing fourth-order terms yields the most parsimonious truncation that is linearly well-posed.

Energy.
We will show that both the nonlocal model (28) and the Cahn-Hillard model (29) can be interpreted as gradient flows. Both may be written in the formu where L is a linear operator. For (28), Lu = Q a * u and for (29), For both models, f (u) = u 2 /2. In (42), ∇ · u∇ is a negative-definite self-adjoint second-order operator on non-negative functions u. We now define quantities useful for constructing an energy. First, define a symmetric quadratic form for which Lu is the first variation. For (28) define and for (29) define Next, choose whose first variation is f (u). The total energy for both models is The first variation of the energy is and both evolution equations can be written as It follows that the evolution is a gradient flow in a weighted metric [8,9,30,42,61,67]. The time derivative of the energy is Consequently, energy is decreasing except for stationary states where δE/δu is constant on the support of the solution.
3.3. Energy minimization and numerics. Because our governing equations are gradient flows, we expect almost all initial conditions to evolve to energy minimizers. We exploit energy minimization to explore the equilibria of both models, namely (28) with the Laplace potential (36), and (29). This strategy plays a prominent role for the remainder of this paper. In order for a solution u to be an equilibrium, it must minimize the energy (46), satisfy the appropriate mass constraint, and be nonnegative. In summary, the constrained minimization problem is subject to the constraints: We solve this problem using the fmincon minimization routine of Matlab's optimization toolbox, which accepts an objective function as well as constraints. We use periodic domains in one and two dimensions with equally spaced grid points. For the nonlocal model on a periodic domain, we replace the attractive potential Q a with its periodized analog and evaluate the convolution in Q via matrix multiplication. For the Cahn-Hilliard model, we evaluate derivatives using a centered difference on each interval in one dimension, or each cell in two dimensions. To aid convergence of the algorithm, we also supply the minimization algorithm with the gradient of our collocated approximation of the energy. We typically choose a random initial condition to seed the algorithm. Sometimes we solve the minimization problem on a coarse grid and then refine the grid by factors of two in order to obtain highly-resolved final states.
Our methodology neglects the dynamics of the approach to equilibria, and relatedly, the issue of whether those equilibria are accessible on relevant biological time scales. Nonetheless, we see the exploration of minimizers as an important first step in understanding (28) and (29). 4. Energy minimizers in the absence of an external potential.

Motivation and numerical computations.
Recalling that the nonlocal equation (28) and the degenerate Cahn-Hilliard equation (29) have steady, compactly supported solutions and motivated by the aggregation of biological organisms, we now study the simplest case of a single clump having total mass M on an effectively infinite domain with no external potential (W = 0). A phase plane analysis of (28) suggests that in one spatial dimension, for each mass M , there is a single energy-minimizing clump [65]. For a generalization of (28) to diffusion of power-law form, rigorous proofs show existence and uniqueness in higher dimensions, subject to appropriate assumptions [3,5]. Similarly, degenerate Cahn-Hilliard models such as (29) are gradient flows that minimize an energy; see [50,56] and references therein. Our primary strategy for studying (28) and (29) will be to exploit energy minimization.
To begin, we numerically compute minimizers as described in Section 3, varying total mass M as our parameter for both models, and in both one and two dimensions. When varying M , we choose values of domain length L such that the spatially homogeneous solution is linearly unstable to a single mode. Strictly speaking, the minimizers of the nonlocal model (28) are affected by their periodic images in our computation. In practice, with our choice of attractive potential (36), these effects are exponentially small and invisible on the plots shown. For the Cahn-Hilliard model (29), the compactly supported solution does not interact with its periodic images.
For both models in one dimension, example minimizers appear in  (28) and (29) differ at small masses, but coincide for large M , except in the transition region. Figure 2 scans two key properties of the clumps as a function of M , namely the peak density u ∞ and the size of the support supp(u) . The peak density increases but saturates to a value of 1.5 with increasing M ; we will later analyze this important property. The support increases with M , and the growth appears to be linear at large M (as we later confirm).
Example minimizers in a fully two dimensional geometry appear in Figure 3. The minimizers are radially symmetric. Panels (a) and (b) show clumps for small mass M for (28) and (29) respectively. To give a better sense of the mass distribution, Figure 4 To better understand these results, it is helpful to analyze the steady state problem. To find conditions for energy minimizers, we recall the discussion of [8,9]. Let Here,ū is an equilibrium solution of fixed mass M and δu is a small perturbation of zero mass, so that overall mass is conserved. Thus, where Ω is the support ofū. Then to leading order, the difference in energy between the equilibrium solution and the perturbed one is In order for the energy to be stationary, ∆E must vanish for all perturbations δu.
As δu has zero mass, a necessary and sufficient condition is that δE/δu in (53) is a constant, which we call λ. That is, in the support of u. The quantity λ has a physical interpretation: it is the energy per unit mass at each point in space [8,9]. The condition (54)  be a minimizer. For sufficient conditions, one must examine the continuation of λ outside of the support as well as the second variation of the energy. These ideas are discussed at length in [8,9]. While we have presented numerical results thus far, for the remainder of this section, we analyze the minimizers of (28) and (29).
From (47) and (54), and for W = 0, minimizers satisfy To make analytical progress, we consider the small and large mass limits in turn, with some ideas following [24,65].
When M is small, the characteristic length scale of the solution is small. Taylor expand the kernel as Radial profiles of the radially-symmetric two dimensional minimizers in Figure 3. Minimizers of the nonlocal model (28)  Under this approximation, (56) becomes where we have used the fact that 1 * u = M . Since (|x|) xx = 2δ(x), differentiate twice to obtain Without loss of generality, place the maximum of u at the origin; this implies u x (0) = 0 by symmetry. For convenience, rescale u by u(0) = u ∞ and x by u Multiply by (p 2 ) ξ , integrate once, and apply the conditions at ξ = 0 to obtain the constant of integration, leading to Take the square root of both sides and separate variables to find an implicit solution, As p → 0, ξ approaches the half-width in our rescaled variable, which implies supp(u) = 2 u 1/2 ∞ 3 2 1 0 p dp The mass constraint is By combining (63) and (64), write supp(u) and u ∞ as functions of mass M , supp(u) = 6 π Γ 5 6 Γ 2 3 3 8 For large mass, recall from Figures 1 and 2 that solutions approach wide plateaus with peak density u ∞ ≈ 1.5. To understand this behavior, we analyze the large M limit. When M is large, spatial scales are large. Thus Q a ≈ −δ(x) and (56) becomes from which it immediately follows that u is constant on the support of the solution. To determine what value that constant has, explicitly minimize the energy (46) subject to the mass constraint. The energy is However, the mass constraint dictates that for a this solution, supp(u) = M/u, which we substitute to obtain Thus which yields the critical point u = 3/2. Since the solution is constant, u | ∞ = 3/2, and supp(u) = 2M/3. A curious feature of the approximation we have used is that it is invariant under any area-preserving map of the spatial domain; put differently, the approximation suggests that any arrangement of disjoint clumps having constant density u = 3/2 has the same energy. However, we know there is an energetic cost for each clump (neglected in our analysis) proportional to the measure of its perimeter. Thus, we expect that the global minimizer is a rectangular profile of density u = 3/2.

Minimizers of the Cahn-Hilliard equation in one dimension. For the Cahn-Hilliard model (29), the energy (46) in one dimension becomes
As discussed in Section 1, (29) has compactly supported solutions with zero contact slope. Here, we establish that energy minimizers share this property. We derive this result in one dimension, though one may extend the calculation to higher dimensions.
Consider an energy-minimizing solutionū(x) with a support x ∈ [α, β]. We compute the first variation with respect to changes in the endpoints in addition to the solution itself. Let α =ᾱ + δα , β =β + δβ, u =ū + δu. Computing the first variation of the energy, we find where In deriving (71), we have used the fact that u vanishes at the endpoints x =ᾱ and x =β. For δE to vanish, it is necessary that u x (ᾱ) = u x (β) = 0. To see this, choose δu = 0 and suppose u x (ᾱ) and u x (β) are nonzero. In this case, choosing δα > 0 would reduce the energy, as would choosing δβ < 0. (The change in mass resulting from this shifting of the endpoints slightly into the support is a second-order effect.) However, this violates the assumption thatū is a minimizer. Therefore, u x (ᾱ) = u x (β) = 0. Throughout much of the rest of our analysis of (29), we make use of the result we have just shown, namely that energy-minimizing solutions have zero contact slope. From (47) and (54), and for W = 0, minimizers satisfy This differential equation has arisen previously in studies of the Cahn-Hilliard equation [44] and the thin film equation [50] and has solutions that can be expressed in terms of elliptic functions. Multiply through by u x , integrate once, and enforce that the energy minimizing solution has zero contact slope to obtain By symmetry, u x (0) = 0 from which it follows that at the center of the clump solution. Therefore, we know the relationship between u ∞ and λ, namely We derive expressions for the size of the support of the solution supp(u) and for the mass M in terms of u ∞ . For convenience, rearrange (74) as For supp(u) , we have where K represents the complete elliptic integral of the first kind. For mass M we have where E represents the complete elliptic integral of the second kind. As u ∞ increases from 0 to 3/2, the mass M increases from 0 to infinity, as seen in Figure 2(a). Ideally, we would invert these relationships to solve for u ∞ and supp(u) as a function of M , but this is cumbersome due to the elliptic integral functions. Therefore, as we did previously for the nonlocal model (28), we will analyze the small and large mass limits, which yield simple explicit expressions for u ∞ and supp(u) .
For small masses, neglect the nonlinear term in (73) to obtain As mentioned above, the solution that minimizes E has zero contact angle and must also respect the mass constraint. Up to translation, the unique solution is Thus, u ∞ = M/π and supp(u) = 2π. For large masses, the u xx term in (73) is negligible and the solution is identical to the large mass limit of the nonlocal case considered above. That is, as M → ∞, u ∞ → 3/2 and supp(u) → 2M/3.

Minimizers in two dimensions.
We briefly discuss the small and large mass limits in two dimensions. The small mass limit for the nonlocal model (28) is difficult to compute explicitly because of convolution with the kernel Q a ∝ e −|x| . An alternative approach is to use a different kernel such as the Bessel kernel used in [24,31]. For our kernel, numerical explorations and a scaling argument both suggest that in the limit M → 0, supp(u) and u ∞ both approach zero in a fashion analogous to the problem in one dimension. The small mass limit for the Cahn-Hilliard model (29) is also analogous to its corresponding one dimensional case, satisfying The solution is where the radius of the support, r * , is the first zero of J 1 (r), and J 0 and J 1 are Bessel functions of the first kind. The support and peak density are .
The large mass calculations for two dimensions are nearly identical to those for one dimension. For (28) and (29) as M → ∞, u ∞ → 3/2 and the minimizer is a disc of density u = 3/2, having radius We highlight two key results from this section. First, the Cahn-Hilliard model (29) approximates the nonlocal model (28) well, especially for large masses and away from an aggregation's boundaries. Second, in the limit of large mass, both models produce groups with an energetically preferred peak density of 3/2. This value is the benchmark against which we address a biologically motivated question: can an external potential modeling the environment be used to reduce the peak density of an aggregation? 5. Energy minimizers in an external potential well. With (29) established as a reasonable approximation of (28), at least for large enough masses and away from transition regions, we now use the degenerate Cahn-Hilliard model as a testbed to investigate the effects of the environment; that is, we now allow W = 0. More specifically, we focus on cases where W is a potential well, and show that when this well has sufficiently tall and steep sides, it can be replaced by an approximate obstacle potential that is either zero or a sufficiently large, and in fact effectively infinite, constant. In this scenario, the support of a minimizer is constrained to avoid locations where W is large. For the remainder of this paper, we simply refer to "obstacle potentials," but remind the reader that in our computational implementation, these are approximate obstacle potentials with large, rather than infinite, height.
The velocity in (29) only depends on the gradient of W . Furthermore, adding a constant to W affects the energy (70) by adding a multiple of the conserved mass. Therefore, W is arbitrary up to a constant. Without loss of generality, take the minimum value of W to be zero. As a prototypical example, consider a Gaussian barrier placed at each end of the interval [0, L], that is, where A measures the volume of the barriers and σ is the characteristic width. Figure 5 shows the ensuing minimizers of (29). Panel (a) varies the height A with fixed width σ = L/40 on a domain of length L = 10 with mass M = 15. As A increases, the density at the boundary decreases and is eventually driven to u = 0. For sufficiently large A, W drives the minimizer to zero over an interval, creating vacancies near the barriers. Panel (b) is similar, but varies σ with fixed A = 10, which is sufficiently large to ensure u = 0 at the boundary. As the barrier width σ decreases, the vacant interval shrinks to a point, and W pins the minimizer to u = 0 at the boundaries. For the smallest few values of σ, W is nonzero only at a the edges of the computational domain; that is to say, W approaches infinity at the domain endpoints and is zero within numerical tolerance elsewhere.
To restate: as we decrease σ, the narrow Gaussian, which we have numerically under-resolved, approaches an obstacle potential, and the minimizer approaches one which is pinned to u = 0 at the boundary. This limit suggests that we may replace our smooth potential W with an approximate obstacle potential that is large on some subset of the gridpoints and zero elsewhere. In this limit, continuity of W is lost. One consequence is that the analysis in Section 4.3 demonstrating zero contact slope for minimizers no longer holds. Indeed, it is apparent in Figure 5(b) that minimizers may have nonzero slope where they are pinned by the potential at the domain boundary.
The strategy of adopting an obstacle potential is effective in two dimensions as well. Figure 6 shows a minimizer of (29) with mass M = 9 on a domain with sides of length L = 3. In the absence of an external potential, the minimizer has constant density u(x) = 1. Panel (a) takes the potential W (x) = 10 3 at computational gridpoints that trace out the white heart-shaped curve, and W = 0 elsewhere. The potential drives the minimizer to u = 0 along this curve, apparent as the dark blue region, and the peak density is u ∞ ≈ 2.2. Panel (b) is similar, but takes W (x) = 10 3 on and exterior to the curve, creating a heart-shaped well. The potential confines the minimizer to the well, drives it to zero along the boundary, and pushes the peak density even higher, to u ∞ ≈ 4.7.
The next section focuses on cases where W is an obstacle potential, and examines strategies for reducing the peak density by engineering W .
6. The effect on peak density of random versus periodic potentials. Our biological motivation for this study is the dynamics of locust populations. Intuitively, one might think that the distribution of many organisms in the wild is affected by the distribution of resources on which they subsist, and indeed, this issue is part of the rich field of biogeography. The way that resources shape the distribution of locusts is intimately tied to locusts' phase polyphenism. Phase polyphenism refers to individuals of the same genotype manifesting different phenotypes [1,60]. For locusts, a pivotal phenotypic difference is the insect's social behavior, or lack thereof. A locust may exist in a solitarious phase in which it seeks isolation or a gregarious phase in which it seeks other locusts. Locusts in the gregarious phase may form large, migratory swarms that decimate crops. Therefore, dense social aggregations on the ground are of interest as they are precursors to dangerous flying groups.
Previous work [66] identifies a population density threshold above which a collective transition to a dangerous, gregarious group takes place. Thus, one way to prevent airborne swarms might be to suppress the density of locusts on the ground, so that the group never reaches the critical threshold. Because social locusts display taxis not only to each other, but to food, one can ask whether there exist external potentials W , modeling food sources, which serve to minimize the peak density u ∞ of the Here, A is sufficiently large to drive u to zero near the boundary. As σ decreases, the vacant segment shrinks. The smallest few values of σ to yield effectively nonzero W only at a the edges of the computational domain; that is to say, W is large at the domain endpoints and zero within numerical tolerance elsewhere. This choice pins the minimizer to be zero just at the endpoints, where a finite contact slope is observed.
population. We investigate this problem within the framework of (29). The locust model in [66] includes additional effects -crucially, the phase change from solitarious to gregarious locusts -but an understanding of the interaction between locusts and food sources even in the absence of phase change is an appropriate starting point for investigation. Figure 7 reveals some of the intricacies of the problem. Consider a one dimensional domain and obstacle potentials W (x) describing a sequence of square wells with dividers between them, intended to model cropland divided into plots. Each panel shows a minimizer of (29). At low and moderate masses M (top and middle rows), applying a large number of square wells in W (x) to the fixed domain reduces the  (29) for two obstacle potentials W (x). The domain has sides of length L = 3 with N = 160 computational grid points along each axis. Total mass is M = 9. (a) W (x) = 10 3 on the heart-shaped curve, and W = 0 elsewhere. The potential pins the minimizers to u = 0 along this curve. (b) Similar to (a), but the potential is now nonzero on and exterior to the curve, creating a heart-shaped well. In the absence of a potential, the minimizer is a constant, u(x) = 1. In (a), W drives the peak density higher, to u ∞ ≈ 2.2. In (b), W drives the peak density even higher, to u ∞ ≈ 4.7.
peak density u ∞ , whereas at high masses (bottom row) it has the opposite effect. Furthermore, at low and moderate densities, social aggregation can cause the mass to clump in a subset of the wells, as seen in panels (b) and (e), which may serve to either reduce or increase peak density.
In summary, the goal of this section is to choose external potentials W to minimize the peak density u ∞ using the Cahn-Hilliard model (29) as a testbed. For tractability, we restrict attention to the class of obstacle potentials. The next three subsections focus on a one dimensional domain with periodic W . In the the final subsection, we compare results for periodic potentials with those for random potentials in both one and two dimensions. Random potentials tend to increase peak density u ∞ while periodic potentials can decrease it.
6.1. Energy minimizers in a single square well. For the Cahn-Hilliard model (29), when supp(u) for a given mass on an infinite domain is less than the well width , the minimizer is identical to that found in Section 4. However, when the support is larger, the minimizer feels the boundaries of the square well. For sufficiently large W at the edge of the well (that is, sufficient contrast between desert and planted regions), u will be driven to zero.
From (47) and (54), minimizers satisfy subject to the mass constraint and the added restriction that u = 0 at the edges of the well because of the external potential. An implicit solution for u as a function of λ, mass M , and well size could be found in terms of elliptic integrals. However it is simpler to consider small and large mass limits, as in Section 4, and augment our calculations with numerical computations. First, we calculate the energy E of the minimizing solution for several cases depending on M and . Case IA: Small M , < 2π. In Section 4, we showed that for an infinite domain, the small M solutions have support of size 2π. We expect that for our present case, the solution will occupy the entire domain. To examine the small mass limit, linearize (93) to obtain the boundary value problem in addition to the mass constraint. The solution is Direct computation via substitution into (70) yields the energy for this solution, Crucially, the energy changes from positive to negative as increases through π.
Case IB: Small M , > 2π. The results are identical to the infinite interval result (so long as the support of the solution is less than ). Substituting (87) into (70), we find the energy Case II: Large M . Numerical simulation shows the density profile to be constant in the potential well except for a small boundary layer near the sides of the well of thickness O( /M ). Exploiting this fact, we have and thus the energy is 6.2. Energy per unit mass in a single well. Momentarily, we will make predictions about energy minimizers in multiple wells. To minimize the total energy we will show that a system with multiple wells distributes mass among those wells to minimize the energy per unit mass. Thus, as another preliminary step, we use the results of the previous subsection to calculate the energy per unit mass in a single well. Define the energy per unit mass, From the previous subsection, we can determine the behavior of E(M ) for the limiting cases of small and large M . We will see that E(M ) may have a minimum, which will be a crucial characteristic when we consider the case of multiple wells. We refer to the value of M that minimizes E(M ) as M * . Case I: < π. A sample computation of E(M ) is shown in Figure 8(a). It is monotonically increasing. This result is consistent with the small mass result (96), from which it follows that which is linearly increasing with M . Also, from the large mass result (99), which is again increasing. We conclude that M * = 0, that is, having the minimum possible mass in the well minimizes E(M ). Case II: π < < 2π. A sample calculation of E(M ) is shown in Figure 8(d). Relations (101) and (102) still hold, but crucially, the coefficient in the small mass limit is negative rather than positive, as it was in Case I. Thus, we expect M * to occur between the decreasing behavior at small mass and the increasing behavior at large mass, as evidenced in the figure.
Case III: > 2π (but finite). Recall that for an infinite domain and in the limit of small mass, the minimizer is given by (87) and the support is 2π. It is not surprising, (j-l) Case IV, chosen to be effectively infinite, much larger than supp(u). Here, E(M ) decreases monotonically to −3/8 and u ∞ increases monotonically to 3/2. Additionally, supp(u) widens without bound and, asymptotically, grows linearly. therefore, that the same small mass result holds for intervals wider than 2π. Plugging the minimizer into the energy and dividing by M yields which is decreasing. As the mass increases, the support widens as well, and eventually fills the domain, as shown in Figure 8(i). For sufficiently large mass, the density is nearly constant with u(x) ≈ M/ (except at the boundaries) and the approximation (102) still holds. As in Case II above, E(M ) is decreasing at small mass, increasing at large mass and reaches a minimum at M * between these two regimes. Case IV: large (effectively infinite). The small mass approximation (103) still holds. As M increases, the minimizer follows the branch of solutions found for the infinite interval. Moreover, by the energy argument in Section 4, the solution approaches a single rectangle of height u = 3/2, having, by direct calculation, Thus, E(M ) decreases monotonically but approaches this constant, and so M * is effectively infinite. It will be energetically favorable for the mass to form a clump of density u = 3/2 inside a single well as was seen previously.
6.3. Energy minimizers in multiple wells. We now take the external potential W to consist of a periodic sequence of square wells of width . The key question is whether the mass will distribute itself equally amongst the wells or whether it is energetically preferred to concentrate in a subset of the wells. We use the results from the previous subsections to answer this question.
More concretely, suppose that there are n equal wells in the domain and that for each well, E(m) is the energy associated with a minimizer of mass m located in that well. The total energy for a distribution of masses m i > 0 is because total mass is fixed. To proceed, consider an ansatz where the mass is equally distributed between k of the n wells. In this case, where µ is the mass contained in each occupied well. From this equation, we deduce that a necessary condition to minimize E(M ) is to minimize the energy per unit mass E(µ) over admitted values of µ.
For the sake of argument, suppose one can choose µ = M * . In this case, E(M ) is globally minimized. In practice, the values of µ are quantized since µ = M/k where k is a positive integer less than or equal to n. For large enough n, the system appears to choose µ that approximates M * well. Now apply these ideas to the numerical results in Figure 7. For panels (a,d,g) in the first column, there is one square well, that is, = L = 12π. For panels (b,e,h) in the second column, there are eight square wells, that is = L/8 = 3π/2, for which we have estimated M * ≈ 2.97 from the data in Figure 8(d). In panel (b), as M ≈ M * , the mass fills a single well. In panel (e), as 4M * < M < 5M * , the mass fills five of the wells. In panel (h), as M > 8M * , the mass equipartitions among all eight wells. Finally, for panels (c,f,i) in the third column, there are 24 square wells, that is, = L/24 = π/2 for which M * = 0. Therefore, for all values of M , mass equipartitions among all wells.
6.4. Comparing random and periodic potentials. We now demonstrate that a periodic pattern of resources in the environment can reduce peak density below what might be seen with a random distribution of resources. This effect is evident in Figure 9, which compares these two cases in one and two dimensions. For all four panels, the average density is u = 1/2, which is below the peak density of 3/2 seen in minimizers in the absence of an external potential (see Section 4). Additionally, we fix φ, the volume fraction of the computational domain on which the obstacle potential W is effectively infinite (here, 10 4 ) and we take W = 0 elsewhere. These choices of W model, respectively, regions barren and lush with vegetation. In Figure 9, we set φ ≈ 0.08.
First, consider the one dimensional case. In panel (a), mass equipartitions among all 40 available potential wells, and u ∞ ≈ 0.76. For the randomly distributed potential used in (b), some wells are left vacant and the peak density increases to u ∞ ≈ 1.95. Now, for two dimensions, panel (c) demonstrates that mass equipartitions among all 25 potential wells in the five-by-five grid, and u ∞ ≈ 1.2. However, for the randomly distributed potential used in (d), there exist hot spots with peak densities as high as u ∞ ≈ 1.7.
In short, periodic potentials can suppress the peak density below 3/2, while random potentials can concentrate density at values greater than 3/2. Figure 10 provides a systematic investigation for varying volume fraction φ in one and two dimensions. In both cases, at low φ, u ∞ ≈ 3/2 (dotted line) for both random (shown as squares) and periodic (shown as circles) potentials. As φ increases, the periodic potential suppresses peak density below 3/2 whereas the random potential increases it above. This result suggests the viability of carefully engineering resource layout as a strategy for reducing a population's peak density.

Conclusion.
Understanding how populations respond to social forces and environmental cues such as resource distribution is essential for modeling biological groups. For example, for desert locusts, spatially varying resource distributions can concentrate populations, which in the wild can lead to collective gregarization and dangerous outbreaks [36]. We have laid out a framework for studying these problems by exploiting energy minimization of the governing models.
More specifically, we began with aggregation model describing nonlocal social attraction, local, nonlinear repulsion, and an external potential modeling the environment. This model possesses an energy that is minimized by the dynamics. However, the model poses (at least) three challenges. First, time-dependent simulations are costly, and may converge slowly to attractors. Second, attractors typically have compact support. Standard simulation methods may produce oscillations emanating from contact points which typically violates nonnegativity of the solution. Third, the nonlocality couples all spatial locations, rendering numerical computations expensive, especially in two or more dimensions.
To address the first two challenges, we exploit the energy minimizing dynamics by focusing on the analysis and numerical computation of minimizers. This avoids costly computations of dynamics, and additionally, minimization algorithms allow us to enforce nonnegativity of the solution as a simple constraint in the optimization procedure. To address the third challenge, we perform a long-wave approximation of our model to obtain a fourth-order degenerate Cahn-Hilliard equation which is local in space while preserving the energy minimizing character of the original model. In summary, we have developed an approach to studying aggregations that replaces simulating a nonlocal equation with minimizing a local energy subject to a nonnegativity constraint.
In the absence of an external (that is, environmental) potential, minimizers of the nonlocal and Cahn-Hilliard equations agree well at large masses and away from the  Fig. 10. Peak density u ∞ of minimizers of the Cahn-Hilliard equation (29) for two classes of obstacle potentials W (x). We vary φ, the volume fraction of the computational domain in which the environmental potential W is effectively infinite (here, 10 4 ). Elsewhere, W = 0. These choices of W model, respectively, regions barren and lush with vegetation. Open circles correspond to choosing for W a periodic pattern of square wells, representing engineered cropland with barriers between crop beds. Squares correspond to choosing for W a uniformly random distributed pattern, producing an irregular, clumpy vegetation landscape. Error bars represent maximum and minimum values over ten realizations of W for the random case. Dotted lines represent u ∞ = 3/2, the energetically preferred peak density for large mass solutions in a large domain. (a) Minimizers in one dimension, with domain length L = 100, population mass M = 50, simulated with N = 480 gridpoints. (b) Minimizers in two dimensions, with domain length L = 36, population mass M = 648, and N = 120 gridpoints along each axis. For both (a) and (b), the average population density is 1/2. In both panels, a randomly generated W tends to increase peak density u ∞ above 3/2 while periodic square wells tend to decrease it below 3/2. contact points. In the limit of large domains and modest densities, these minimizers approach a compactly supported clump of population with peak density 3/2 and steep edges.
An external potential drives mass to accumulate near potential minima, often creating collections of compactly supported clumps. Many environmental landscapes can be partitioned into regions that are either inviting or inhospitable to biological organisms. This characterization arises naturally in numerical computations for sufficiently steep potential wells. Thus, we model these landscapes via obstacle potentials that are effectively infinite (inhospitable) on some subset of the domain and zero (inviting) elsewhere.
In this obstacle potential testbed, we ask how resource distribution effects peak population density. In the absence of an external potential and at modest average population densities, aggregations form with peak density 3/2, as mentioned above. A spatially random obstacle potential modeling the natural environment drives peak densities higher, even when the inhospitable portion of the domain is less than 10% of the total. However, a periodic pattern of square wells modeling divided crop beds can reduce the peak density below 3/2.
Anthropologists have chronicled the use of periodic planting patterns since ancient times, and typically attribute this strategy to better water management and thus increased crop yield. We are intrigued, but admit that it is pure speculation, to ask whether this strategy might confer any pest control advantages. Regardless, the biological literature notes that for locusts, "fractal dimension" of resources in the environment can drive large peak population densities [36]. While perhaps just a caricature, our simple model suggests that carefully distributing resources can inhibit populations from reaching their intrinsic peak density and plausibly could help control undesirable biological phenomena such as gregarious locust outbreaks that are triggered by surpassing a density threshold.