Detecting minimum energy states and multi-stability in nonlocal advection–diffusion models for interacting species

Deriving emergent patterns from models of biological processes is a core concern of mathematical biology. In the context of partial differential equations, these emergent patterns sometimes appear as local minimisers of a corresponding energy functional. Here we give methods for determining the qualitative structure of local minimum energy states of a broad class of multi-species nonlocal advection–diffusion models, recently proposed for modelling the spatial structure of ecosystems. We show that when each pair of species respond to one another in a symmetric fashion (i.e. via mutual avoidance or mutual attraction, with equal strength), the system admits an energy functional that decreases in time and is bounded below. This suggests that the system will eventually reach a local minimum energy steady state, rather than fluctuating in perpetuity. We leverage this energy functional to develop tools, including a novel application of computational algebraic geometry, for making conjectures about the number and qualitative structure of local minimum energy solutions. These conjectures give a guide as to where to look for numerical steady state solutions, which we verify through numerical analysis. Our technique shows that even with two species, multi-stability with up to four classes of local minimum energy states can emerge. The associated dynamics include spatial sorting via aggregation and repulsion both within and between species. The emerging spatial patterns include a mixture of territory-like segregation as well as narrow spike-type solutions. Overall, our study reveals a general picture of rich multi-stability in systems of moving and interacting species.


Introduction
A central purpose of mathematical biology is to provide a way of linking biological processes to emergent patterns (Levin 1992;Murray 2001). In cell biology, such insights can illuminate the mechanisms behind the growth of cancerous tumours, and inform the development of interventions to slow or halt that growth (Altrock et al. 2015;Byrne 2010;Painter and Hillen 2013). In ecology, the insights on mechanisms behind animal space use can be valuable for species conservation (Bellis et al. 2004;Macdonald and Rushton 2003;Zeale et al. 2012), ensuring maintenance of biodiversity (Hirt et al. 2021;Jeltsch et al. 2013), and controlling biological invasions (Hastings et al. 2005;Lewis et al. 2016;Shigesada and Kawasaki 1997).
For partial differential equation (PDE) models of biological systems, one useful method to link process to pattern is to construct an energy functional for a system, if it exists. Then the local minima of this energy functional give possible final configurations of the system. Our focus here is to develop techniques for finding such local energy minima in a particular system of PDEs describing symmetric nonlocal multi-species interactions, with the parallel biological aim of being able to detect and describe the possible long-term patterns that may emerge from underlying processes.
The PDE system we focus on is a multi-species system of nonlocal advection diffusion equations recently introduced (Potts and Lewis 2019) and slightly generalised by Giunta et al. (2021a). This system models the spatial structure of ecosystems over timescales where births and deaths are negligible and has the following functional form for i ∈ {1, . . . , N }, where D i and γ i j are constants, and u i (x, t) is the density of a species of moving organisms in location x at time t. Individuals detect the presence of others over a spatial neighborhood described by spatial averaging kernel K , which is a symmetric, non-negative function with K L 1 = 1. The magnitude of γ i j gives the rate at which species i advects towards (resp. away) from species j if γ i j < 0 (resp. γ i j > 0). Whilst the detection of individuals may be direct, e.g. through sight smell or sound, Potts and Lewis (2019) showed that the above formalism can also be used when interactions are mediated by marks in the environment or memory of past interactions. Note that, as well as modelling different species of organism, Eq.
(1) can also be used to model N different groups within a species, or to describe more complex situations where organisms may be spatially delineated by something other than species, e.g. mixed-species territorial flocks of birds (Mokross et al. 2018). However, we use the term 'species' for simplicity. Equation (1) generalises a variety of existing models. In the case N = 1 and γ 11 < 0, Eq. (1) is an aggregation-diffusion equation (Carrillo et al. 2018(Carrillo et al. , 2019 and also arises in model of animal home ranges (Briscoe et al. 2002). For N = 2 and γ 12 , γ 21 > 0, Eq. (1) can be related to models of territory formation (Ellison et al. 2020;Potts and Lewis 2016b;Rodríguez and Hu 2020) and cell sorting (Burger et al. 2018) (the latter also includes γ 12 , γ 21 < 0). The case of arbitrary N with γ i j = 1 has also been recently studied in the context of territories (Ellefsen and Rodríguez 2021). Finally, the N = 2 case with γ 12 and γ 21 having different signs has been studied in the context of predator-prey dynamics (Di Francesco and Fagioli 2016). So there is a wide range of possible applications arising from Eq. (1).
Whilst our approach is quite general in potential applicability, there are various specific biological questions that might be addressed by classifying minimum energy solutions. A simple example is that of animal territory formation. How much avoidance is necessary for segregated territories to form? Is the emergence of territories history dependent? Do symmetric avoidance mechanisms always lead to symmetric territories? As another example, in the case of mutualistic species, we can ask similar questions. How much attraction is necessary for aggregation? Is it history dependent? All of these questions can benefit from the insight provided by classifying minimum energy solutions to Eq. (1), as well as more complex questions regarding multi-species questions that may exhibit a mixture of attraction and avoidance mechanisms.
The model given by Eq. (1) has been shown to exhibit rich pattern formation properties, including aggregation, segregation, oscillatory patterns and non-periodic spatio-temporal solutions suggestive of strange attractors (Potts and Lewis 2019). In Potts and Lewis (2019), for the simple case where N = 2, γ ii = 0, and γ 12 = γ 21 , an energy functional was constructed that is decreasing in time, bounded below, and becomes a steady state of Eq. (1) as t → ∞. Furthermore, numerical experiments suggest that only stationary patterns emerge in this case (Potts and Lewis 2019). Here, our first task is to generalise this N = 2 energy functional to arbitrary N , but where γ i j = γ ji for all i, j, ∈ {1, . . . , N }. Related work by Jüngel et al. (2022) found two more energy functionals which are based on the Shannon entropy on the one hand and a Rao-like entropy on the other. However, our focus here is on the generalization of the energy function from Potts and Lewis (2019).
Once this energy functional has been constructed, our second task is to minimise it to ascertain the functional form of the local minimum energy solutions. For this, we work in the local limit, i.e. where K tends towards a Dirac-δ function. We give a numerical technique for showing that, if we start with a class of stable steady state solutions for different K , then take the local limit, we return a piecewise constant function. This technique makes use of the theory of Gröbner bases and associated methods from computational algebraic geometry. It is a generalisation of a method first used in Potts and Lewis (2016b). In situations where the local limit is piecewise constant, local minima of the energy functional can be found by searching through the space of piecewise constant functions. We show that this can sometimes be done analytically, using some basic examples in one spatial dimension to illustrate the methods. Even in case N = 2, this process reveals a range of situations where there are multiple local energy minima, all of which we verify via numerics away from the local limit. Overall, the methods presented here enable users to detect local minimum energy states of Eq. (1), including multiple minima, in any situation where γ i j = γ ji . This paper is organized as follows. We begin with linear stability analysis, in Sect. 2. This sets the stage by showing that the γ i j = γ ji case (for all i, j) leads to stationary pattern formation at small times (from perturbations of the homogeneous steady state) as long as the species have the same-sized populations. In Sect. 3, we construct an energy functional associated with Eq. (1) in the case γ i j = γ ji (for all i, j) and analyze its properties, particularly that it decreases in time and is bounded below. Noteably, unlike the linear analysis, this does not require the species to have the same-sized populations. This section ends with a conjecture about the structure of the attractor, which is somewhat stronger than what we are able to show in this paper, but for which we have numerical evidence to suggest it might be true. In Sect. 4, we describe our technique for finding stable steady states, assuming that the local limit of stable steady states is piecewise constant, generalising a method used in Potts and Lewis (2016a). In Sect. 5, we give a method for proving that this local limit is piecewise constant, demonstrating our proof for N = 2 and arbitrary γ i j , then for N = 3 with specific examples of γ i j .

Notation and assumptions
We use the following notation conventions throughout. Let S ⊂ R n be a measurable set. Then we denote the measure of S by |S|, so that where 1 : R n → R is the constant function 1(x) = 1. Let ⊂ R n and f : L p ( ) → R. We use the following norms Let M ∈ N and g = (g 1 , g 2 , . . . , g M ) : (L p ( )) M → R. Then we define To ease notation, we usually write g L p ( ) instead of g (L p ( )) M , if the meaning is clear from the context. We also may drop explicit dependence on .

Linear stability analysis
Inhomogeneous solutions of PDEs can emerge when a change in a parameter causes the loss of stability of a homogeneous steady state, leading to the formation of inhomogeneous solutions (sometimes referred to as Turing patterns after Turing (1952)), which can be either stationary or periodically oscillating in time. In this section, we will analyze the linear patterns supported by Eq. (1).
In Eq. (1), the total mass of each species i is conserved in time, indeed on the periodic domain T, on which conditions (3) hold, the following identities are satisfied Hence, for all i = 1, . . . , N , where the constant p i is the population size of species i. Therefore, Eq.
(1) has an homogeneous steady statē unique for each value of p i (determined by the initial condition). To study the stability ofū, we introduce the vector where u (0) is a constant vector, λ ∈ R is the growth rate of the perturbation, x = (x 1 , . . . , x n ) ∈ T and κ = (κ 1 , . . . , κ n ) is the wave vector, whose components are the wave numbers of the perturbation and must satisfy the boundary conditions (3). We thus have Substituting Eq. (7) into Eq. (1) and neglecting nonlinear terms, we obtain the following eigenvalue problem where and whereK For each κ, the eigenvalue with greatest real part (called the dominant eigenvalue) determines whether or not non-constant perturbations of the constant steady state at wavenumber κ will grow or shrink at short times. If the dominant eigenvalue has positive real part and non-zero imaginary part, then these perturbations oscillate in time as they emerge. If the dominant eigenvalue is real, such oscillations will not occur at short times. Now, ifū i =ū j and γ i j = γ ji for all i, j = 1, 2, . . . , N then L is symmetric, so all its eigenvalues are real (Artin 2011). Therefore non-constant perturbations of the constant steady state will not oscillate at short times. In practice, situations where the dominant eigenvalue is real and positive are often accompanied by non-constant stable steady states. Although this does not follow by necessity (Giunta et al. 2021b), this observation nonetheless suggests that this case provides a good starting point in searching for non-constant stationary patterns.
In the following sections, we will study the γ i j = γ ji case through an energy functional analysis, showing how this can give us insights into the structure of nonconstant stable steady states. It turns out that for this analysis, we do not need the additional assumptionū i =ū j .
We conclude this section by analysing the N = 2 case in detail, to provide some results required in later sections. In this case, the characteristic polynomial of the matrix L is whose roots are giving the eigenvalues of L. The condition γ 12 = γ 21 ensures that the argument of the square root is always positive and therefore the eigenvalues λ ± are real. As a concrete example, if p 1 = p 2 = 1, L 1 = · · · = L N = 1, D 1 = D 2 , γ 12 = γ 21 and γ 11 = γ 22 then the system admits a linear instability if there exists at least one κ > 0 such that 3 Energy functional In this section, we will define an energy functional associated to Eq. (1) with γ 12 = γ 21 , and show that it is continuous, bounded below, decreases in time, and that its stationary points coincide with those of Eq. (1). This gives evidence to suggest that Eq. (1) with γ 12 = γ 21 will tend towards a steady state, which will be inhomogeneous in space if the constant steady stateū is linearly unstable. During this section, we will assume a positivity result, namely that u i (x, 0) > 0 implies u i (x, t) > 0, for all i = 1, . . . , N , for all t > 0. This result has been already proved in one spatial dimension (Giunta et al. 2021a). This proof relies on a Sobolev embedding theorem only valid in one dimension, so other tools will be needed to give a proof in arbitrary dimensions. Indeed, at the time of writing, this positivity result has not yet been established in arbitrary dimensions.
First, we re-write Eq. (1) as follows Then we define the following energy functional where x = (x 1 , x 2 , . . . , x n ).
The first term D i u i ln u i is the entropy of each of the populations on their own and the second term γ i j (K * u j )u i denotes the interaction energy between the populations (Carrillo et al. 2020). The factor 1 2 before the sum is required so that we can leverage the γ i j = γ ji symmetry later on.

Proposition 1
The energy functional E, defined in Eq. (16), is a continuous function of the variables u 1 , u 2 , . . . , u N .
Proof First we show that the following functions are continuous as long as u i is positive across space and time Equation (17) is continuous since it is the product of continuous functions. For Eq. (18), we first observe that if K ∈ L 1 and u ∈ L p , with 1 ≤ p ≤ ∞, then by Young's convolution inequality. Moreover, since K where the last equality uses K L 1 = 1. Equation (20) shows that u → K * u is a Lipschitz function and thus a continuous function. Therefore Eq. (18) is continuous because it is the product of continuous functions. This shows that the integrand in Eq. (16) is continuous. Now let 1 ≤ p ≤ ∞ and g : L p ( ) → L p ( ) be a continuous function. Define a function G : L p ( ) → R by It remains to show that G is continuous. To this end, let > 0 and u ∈ L p ( ). Then since g is continuous, there exists δ > 0 such that for any Remark 1 Note that whilst we have used K L 1 = 1, the previous proposition also holds for any K ∈ L 1 .
Proposition 2 Suppose γ i j = γ ji , for all i, j = 1, . . . , N . For any positive (for each component) initial data (u 1,0 , . . . , u N ,0 ), the energy functional E [u 1 (x, t), u 2 (x, t), . . . , u N (x, t)] is non-increasing over time, where (u 1 , u 2 , . . . , u N ) is the trajectory of Eq. (1) starting from (u 1,0 , . . . , u N ,0 Here, the second equality uses that T g(K * h)dx = T h(K * g)dx as long as K (x) = K (−x) for x ∈ R n . The fourth equality uses γ i j = γ ji and the sixth uses Eq. (15). Before continuing the computations in Eq. (22), we simplify notation by setting Observing that we continue the previous computation to give The final inequality uses the assumption that u i > 0. The second equality uses integration by parts. The third equality follows from the following equalities and we observe that each term in Eq. (26) is equal to zero due to the periodic boundary conditions in Eq. (3).
Equation (25) shows that E is decreasing over time unless which is a steady state of Eq. (15), or equivalently of Eq. (1).

Remark 2
Proposition 2 rules out the existence of non-stationary, time-periodic solutions. Indeed, as E is monotonic decreasing, if there exist t, τ > 0 such that so Eq. (27) holds and (u 1 (x,t),. . . ,u N (x,t)) is a stationary solution.
Proof We first observe that for all γ ∈ R, the following inequalities hold The first inequality uses the fact that γ ≥ − |γ |, for all γ ∈ R, the second uses Hölder's inequality and the third uses Young's convolution inequality. Moreover, since u i > 0, condition (5) ensures that u i (x, t) L 1 = p i for all t ≥ 0 and thus the right-hand side of Eq. (28) is finite. Finally, by observing that inf u i >0 {u i ln(u i )} = −e −1 and also by using Inequality (28), we obtain the following estimates where the last equality uses the integral condition (5). Thus E is bounded below.
Proof Since K L ∞ < ∞ , Prop. 3 ensures that the following set is bounded below. Due to the Completeness Axiom of the real numbers, the set in (31) has an infimum l u 0 , which is determined by the initial condition u 0 . Moreover, by Proposition 2, E is a non-increasing monotonic function of time, so tends to its infimum l u 0 as t → ∞.
Proposition 4 shows that for any initial data u 0 ∈ L 1 (T) N the trajectory starting from u 0 evolves over time towards a configuration that is a local minimiser of E, with energy E = l u 0 . We also observe that if E reaches the minimum value l u 0 at a finite time T , then the trajectory becomes stationary. Indeed, if E(u(T )) = l u 0 then E(u(t)) ≡ l u 0 for all t ≥ T . Hence, the minimum at E = l u 0 corresponds to a steady state that is Lyapunov stable (i.e. any solution that starts arbitrarily close to the steady state will remain arbitrarily close). However, it does not guarantee asymptotic stability (i.e. any solution that starts arbitrarily close to the steady state tend toward the steady state).
In the next Section, we will propose a method to determine the structure of these minimum energy states of Eq. (1).
Finally, we note that the convergence of E towards a finite minimum value does not guarantee that every solution converges towards a steady state when γ i j = γ ji , as opposed to fluctuating in perpetuity. Nevertheless, this is something we would like to establish. Indeed, in all our numerical investigations, both here (in Sect. 4) and in previous works (Potts and Lewis 2019;Giunta et al. 2021a), we have only every observed (numerically) stable steady state solutions emerging, and have never observed perpetually fluctuating solutions. Therefore, we conclude this section formulating the following conjecture. This is left as an open problem, but one possible means of attack might be the via the S 1 -equivariant theory of Buttenschön and Hillen (2021), applied there to a single-species system with a similar (but not identical) non-local advection term.

A method to find minimum energy states
In this section, we will propose a method to gain insight into the possible structures of minimum energy to Eq. (1). We build on methods first proposed in Potts and Lewis (2016a, Sect. 3.4) and recent existence results of Jüngel et al. (2022). We work in one spatial dimension and assume the assumptions of Sect. 1.1.
As shown in the previous section, the energy will always tend towards a local minimum, leading to a minimum energy state for the system, which is also a steady state.
When solving Eq. (1) for the top-hat kernel As α tends to zero, the solution appears to tend towards a piece-wise constant function (Panel a and c) or the limit of arbitrarily narrow, arbitrarily high piece-wise constant functions (Panel b and d).
The parameter values used in the simulations are D 1 = D 2 = 1, p 1 = p 2 = 1, γ 11 = γ 22 = 0, γ 12 = 1.05 in Panel a and c, γ 12 = −1.05 in Panel b and d numerically, we find that for decreasing α, the asymptotic steady state solutions look increasingly like piece-wise constant functions, or the limit of arbitrarily narrow, arbitrarily high piece-wise constant functions, with single or multiple peaks. These structures become more singular as α → 0. In Fig. 1, we see this for some simple examples. Note that as α → 0, the top-hat kernel in Eq. (32) becomes a Dirac delta measure, and the model (1) becomes a local cross-diffusion model. Hence we call this limit α → 0 as the local limit. Jüngel et al. (2022) derived a solution theory for non-smooth interaction kernels K , which includes the case of a top-hat kernel as in Eq. (32). They consider Eq. (1) for the case where there are constants π i such that the matrix (π i γ i j ) i j is positive definite. For that case they showed global existence of weak solutions in Sobolev spaces. They also show a local-limit result. As α → 0 there exists a subsequence of solutions of Eq. (1), with K as in Eq. (32), that converge to a solution of the local version of Eq. (1). The norm of this convergence varies depending on the space dimension. In n = 1 we can use any L p -norm and in dimensions n ≥ 2 we use the L n n−1 -norm. These limits are piece-wise constant solutions, and spike solutions, depending on the sign of γ i j . They arise as minimizers of the local version of the energy functional (Eq. (16)), which is Hence in the following we consider piece-wise constant energy minimizers, assuming that they are close to the minimizers of the non-local problem and we confirm this relation numerically. We also focus on the n = 1 case and write L = L 1 for simplicity.
We now explain our method in detail. First, Eq. (25) in one dimension tells us that any minimum energy solution, u i (x), occurs when for each i ∈ {1, . . . , N }. Next we take the local limit of Eq. (34), which in the case K = K α is the limit α → 0. In this limit, Eq. (34) becomes Therefore, either u i (x) = 0, or, for any subinterval on which u i (x) = 0, there exists a constant c i ∈ R such that In principle, there might exist infinitely many subintervals on which u i (x) = 0, and c i may vary between these different subintervals. However, for each set of constants c 1 , . . . , c N , Eq. (36) will typically have a finite number of common solutions (indeed, Sect. 5 shows how to determine whether we are in this 'typical' situation). Therefore, on each subinterval I in which u i (x) = 0, there exists a finite set of values u c i1 , . . . , u c ih , with h ∈ N, satisfying Eq. (36), such that where I il , for i = 1, . . . , N and l = 1, . . . , h, are disjoint subsets of I such that ∪ l I il = I for each i. By considering all such subintervals I together, Eq. (37) defines a class of piece-wise constant functions on [0, L]. The aim here is to examine which of these functions is a local minimum of the energy and satisfies all model assumptions.
The general case is too complicated to deal with in one go, so we demonstrate our method on some simple examples for the case of two species, N = 2. We start by studying the case γ 11 = γ 22 = 0, so there is neither self-attraction nor self-repulsion. We split this analysis further into the cases of mutual avoidance (γ 12 > 0) and mutual attraction (γ 12 < 0). Then we analyze the case where γ 11 , γ 22 = 0.

Analytic results in the local limit
Minimising the energy over the full class of functions given by Eq. (37) turns out to be too complicated. However, our numerics (see Fig. 1) suggest that the local limit (i.e. α → 0 in the case K = K α ) of any solution to Eq. (1) is a function of the following form where u c i ∈ R + and S i are subsets of [0, L], for i ∈ {1, 2}. Therefore we restrict our search by looking for the minimisers of the energy (Eq. (16)) in the class of piece-wise constant functions defined as in Eq. (38).
By Eq. (5), in Eq. (38) we require the following constraint recalling from Eq. (2) that |S| denotes the measure of a set S, not the cardinality, and p i denotes the total population size of species i. We wish to find the solutions of the form in Eq. (38), subject to Eq. (39), that are local minimisers of the energy, Eq. (16). Placing Eq. (38) into Eq. (16), and taking the spatially-local limit (i.e. α → 0 in the case K = K α ), gives where the first equality uses γ 12 = γ 21 , the second equality uses Eq. (38) and the third equality uses Eq. (39). In Eq. (40), notice that if we keep |S 1 | and |S 2 | fixed whilst lowering |S 1 ∩ S 2 | then the energy decreases. Thus, if |S 1 |+|S 2 |≤ L, we can construct disjoint sets S 1 and S 2 , and these will correspond to lower energy solutions than any pair of non-disjoint sets of equal measure. Furthermore, if |S 1 | + |S 2 | > L, we can construct sets S 1 and S 2 , such that |S 1 ∩ S 2 | = |S 1 | + |S 2 | − L and these will correspond to lower energy solutions than any other pair of sets of equal measure. Therefore henceforth, when |S 1 | + |S 2 | ≤ L, we will assume that S 1 ∩ S 2 = ∅, and when |S 1 | + |S 2 | > L, we will assume that |S 1 ∩ S 2 | = |S 1 | + |S 2 | − L.
To search for the local minimizers of the energy in Eq. (40), we thus define To constrain our search, notice that Eq. (39) and |S i |≤ L imply that The region of the (u c 1 , u c 2 )-plane defined by Eq. (42) is shown as white region in Fig. 2. Our strategy will be as follows. First we will look for the local minima of Eq. (41), subject to Eq. (42), in the case where |S 1 | + |S 2 | ≤ L. Then we will look in the region |S 1 | + |S 2 | > L. Combining these results will then give us a complete picture of the local minima of E(u c 1 , u c 2 ). Starting with |S 1 | + |S 2 | ≤ L, Eq. (39) shows that this case is equivalent to the following condition By analysing the partial derivatives of E(u c 1 , u c 2 ) in the region of the (u c 1 , u c 2 )-plane defined by Eq. (43), we see that there are no critical points in this region. Furthermore, Therefore minima in this region must lie on the boundary, p 1 /u c 1 + p 2 /u c 2 = L, which is shown as solid black line in Fig. 2. Analysis of the partial derivative of E(u c 1 , u c 2 ) on this boundary shows that E(u c 1 , u c 2 ) has a unique minimum point, given by This is also a local minimum of the region defined by Eq. (43). This can be shown by performing a Taylor expansion of E(u c 1 , u c 2 ) about the point M S in the region given by Since the slope of the tangent line to the curve , we choose two arbitrarily small constants, and δ, such Since M S lies on the boundary curve |S 1 |+|S 2 |= L (Fig. 2), we have so far only established that it is a minimum of the region where |S 1 |+|S 2 |≤ L. We now need to find out whether it is a minimum for the whole admissible region (the white region in Fig. 2).
To this end, we perform a Taylor expansion of E(u c 1 , u c 2 ) in a neighbourhood of M S within the region |S 1 |+|S 2 |≥ L, which is also the region where p 1 /u c 1 + p 2 /u c 2 ≥ L, by Eq. (39). Since the slope of the tangent line to the curve , we choose two arbitrary constants, and δ, such that Then the Taylor expansion of if γ 12 > D 1 D 2 L p 1 D 1 + p 2 D 2 , where the inequality uses D 2 1 p 1 + D 2 2 p 2 δ ≤ 0. We now examine whether there are any other minima of E(u c 1 , u c 2 ) in the region where |S 1 |+|S 2 | > L. By Eq. (42), the condition |S 1 | + |S 2 | > L is equivalent to p 1 /u c 1 + p 2 /u c 2 > L. Therefore we have the following constraints A direct calculation using partial derivatives shows that there are no local minima of E(u c 1 , u c 2 ) (Eq. (45)) in the interior of the region of the plane (u c 1 , u c 2 ) defined by Eq. (47). Therefore any local minimum must occur on the boundary. On the part of the boundary given by u c i = p i /L, for i = 1, 2, there is a unique minimum at This is also a local minimum of the region defined by Eq. (47). This can be shown by performing a Taylor expansion of where the inequality uses ≥ 0, δ ≥ 0, so that we remain in the u i ≥ p i /L region in Fig. 2. In summary, if 0 < γ 12 < (41)) has a unique minimum, given by M H . However, if γ 12 > D 1 D 2 L p 1 D 1 + p 2 D 2 then E(u c 1 , u c 2 ) has two local minima, given by M H and M S (see Fig. 2). Now, we recover the local minimizer u i (x) (Eq. (38)) of the energy (Eq. (33)). To give a concrete example, we use the parameter values p 1 = p 2 = D 1 = D 2 = L = 1. If (u c 1 , u c 2 ) = M H then u 1 (x) = u 2 (x) = 1, the homogeneous steady state, which we denote by S H . If (u c 1 , u c 2 ) = M S then with |S i |= 1/2, for i = 1, 2, and |S 1 ∩ S 2 |= 0. This is a class of solutions we denote by S 2,2 S , where the subscript S stands for segregation and the superscript 2, 2 denotes the finite positive value that functions u 1 (x) and u 2 (x) take, respectively. To avoid any confusion, we want to stress that the points M H (Eq. (48)) and M S (Eq. (44)) are local minima of E(u c 1 , u c 2 ) (Eq. (41)), while the functions S H and S 2,2 S are minimizers of the energy E[u 1 , u 2 ] (Eq. (40)).
In our example, if 0 < γ 12 < 1/2, E(u c 1 , u c 2 ) (Eq. (33)) has a unique minimum, given by S H . If γ 12 > 1/2 the energy has two local minima, given by S H and S 2,2 S . However, recall that S H and S 2,2 S are derived by minimizing the energy (Eq. (33)) in a particular class of piece-wise constant functions given by Eq. (38). Therefore, the steady states S H and S 2,2 S may not be minima of the full function space where solutions  (41)). The point M H , corresponding to the homogeneous steady state, is always a local minimum. Whether the point M S is a local minimum depends upon the value of γ 12 might live. However, the linear stability analysis performed in Sect. 2, and particularly Eq. (14), suggests that in the limit as α tends to zero, S H is stable if γ 12 < 1. This gives rise to the diagram of analytically-predicted steady states given by the red and black lines in Fig. 3.

Numerical verification
The analysis of Sect. 4.1.1 suggests that for p 1 = p 2 = D 1 = D 2 = L = 1, when 1/2 < γ 12 < 1 and the averaging kernel K is arbitrarity small, Eq. (1) should exhibit bistability between the homogeneous solution, S H , and an inhomogeneous solution arbitrarily close to S 2,2 S . Here, we verify this numerically. Figure 3 summarises our results. To produce this figure, we start with K = K α and γ 12 = 1.2, so that the homogeneous steady state is unstable. The initial condition is a small perturbation of the solution given in Eq. (50) which we run to numerical steady state. We then reduce the magnitude of γ 12 by γ 12 = 0.05 and solve the system again using a small random perturbation of the previous simulation as initial condition. We then repeat this process of reducing γ 12 and re-running to steady state until the system returns to the homogeneous steady state. This process of slowly changing one parameter and re-running to steady state is a type of numerical bifurcation analysis used in many previous studies, e.g. Painter and Hillen (2011). The numerical scheme we use for solving our particular system is detailed in Giunta et al. (2021a).
We examine three different values of α in Fig. 3. For each of these, we observe that the inhomogeneous solution persists below γ 12 = 1 and above γ 12 = 1/2 and, as predicted by our calculations of Sect 4.1.1, the system shows bistabilty and hysteresis. Furthermore, as α decreases (towards the local limit), the numerical branches appear to tend towards the branch calculated in Sect. 4.1.1.
Finally, in Fig. 4, we show some numerical stationary solutions for different values of α, as γ 12 varies in the range [0.55, 1.05]. We observe that, as α decreases, the numerical solution appears to tend to a piece-wise constant function of the class given in Eq. (50) and predicted by the analysis of Sect. 4.1.1.

Analytic results in the local limit
As in Sect. 4.1.1, here we will look for the minimizers of the local version of the energy (Eq. (33)) in the class of piece-wise constant functions defined as where u c i ∈ R + and S i are subsets of [0, L], for i ∈ {1, 2}. Placing Eq. (51) into Eq. (33), and repeating the same argument of Sect. 4.1.1, we obtain In this case, to minimize Eq. (52) we note that, since γ 12 < 0, E[u 1 , u 2 ] can be lowered by increasing |S 1 ∩ S 2 |, whilst keeping everything else the same. Therefore if we keep |S 1 | and |S 2 | unchanged, then |S 1 ∩ S 2 | is maximised when either S 1 ⊆ S 2 or S 2 ⊆ S 1 , so that |S 1 ∩ S 2 |= min i |S i |. Thus and therefore we have that E[u 1 , u 2 ] → −∞ as min{ p 1 u c 2 , p 2 u c 1 } → ∞. As we approach this limit, u c 1 , u c 2 become arbitrarily large, so u 1 and u 2 (Eq. (51)) become arbitrarily high, arbitrarily narrow functions with overlapping support. We will denote the limit of this solution by S ∞ A , in which the subscript A stands for aggregation and the ∞ superscript denotes that the solution becomes unbounded in the local limit. Thus E[u 1 , u 2 ] is minimized by S ∞ A whenever γ 12 is negative, regardless of its magnitude. One can also show, using a very similar argument to Sect. 4.1.1 (details omitted), that the homogeneous steady state, S H , is the only other possible local minimiser of the energy that satisfies Eq. (42), and this is only a local minimum when γ 12 > −L( p 1 D 1 + p 2 D 2 )/( p 1 p 2 ). However, linear stability analysis (Eq. (13)) suggests that, in the limit as α tends to zero, the homogeneous steady state is linearly stable only if γ 12 > −L √ D 1 D 2 /( p 1 p 2 ). Since Young's inequality for products implies that L √ D 1 D 2 /( p 1 p 2 ) < L( p 1 D 1 + p 2 D 2 )/( p 1 p 2 ), any time S H is linearly stable it is also a local energy minimiser within the set of functions given by Eq. (51). The red and black lines in Fig. 5a are the conclusion from combining all the results from Sect. 4.2.1, both energy functional and linear stability analysis, in the case where p 1 = p 2 = D 1 = D 2 = L = 1.

Numerical verification
The analysis of Sect. 4.2.1 suggests that when γ 12 > −L √ D 1 D 2 /( p 1 p 2 ), γ 12 < 0, and α is arbitrarily small, Eq. (1) should display bistability between the homogeneous solution and an inhomogeneous solution, whose structure tends towards S ∞ A as α → 0. Here we verify this conjecture numerically, with results shown in Fig. 5a and b.
To construct these figures, we perform a similar analysis to Sect. 4.1.2. We simulate Eq. (1) with K = K α (Eq. (32)) for small values of α. We use the parameter values p 1 = p 2 = D 1 = D 2 = L = 1, as in Sect. 4.1.1. For these values, the constant steadystate is stable to perturbations at all wavenumbers for −1 < γ 12 < 0. Therefore, we begin our analysis by setting γ 12 = −1.2, reducing the magnitude of γ 12 by a small amount ( γ 12 = 0.05) at each iteration of the analysis, as in Sect. 4.1.2.
Our results show that patterns persist beyond γ 12 = −1, and the extent of this persistence depends on α (Fig. 5a). As α is decreased, the numerical stationary states become higher, narrower functions with qualitatively similar shapes, as predicted by the previous analysis (Fig. 5b).

The case 11 , 22 = 0
The case γ 11 , γ 22 = 0 uses similar arguments to those in Sect. 4.1. We therefore just summarise the results here, leaving details of the calculations for "Appendix A".
In our computations, we consider the case γ 22 = γ 11 and fix the other parameter values as p 1 = p 2 = D 1 = D 2 = L = 1. The analysis of this case reveals five distinct classes of qualitatively-different stable solutions (Fig. 6a), each of which we have verified through numerical analysis (where throughout this section we use 'stable' to mean 'Lyapunov stable'). These are (i) territory-like segregation patterns, S 2,2 S , (b) Fig. 6 Panel a shows the five qualitatively-different local minimum energy states revealed by the analysis in "Appendix A". Note that the S 1,∞ S solution also allows for u 1 and u 2 to be swapped. These plots were produced by setting K = K α , α = 0.025 and by fixing the following parameter values: p 1 = p 2 = D 1 = D 2 = L = 1. For each graph of Panel a, we fixed different values of the parameters γ 11 and γ 12 , in particular we used: γ 11 = 0.2 and γ 12 = 1.05, for S 2,2 S ; γ 11 = −0.15 and γ 12 = 0.4, for S ∞,∞ S and S 1,∞ S ; γ 11 = 0.2 and γ 12 = −1.05, for S ∞ A ; γ 11 = 0.2 and γ 12 = 0.2, for S H . Panel b shows the minimum energy solutions to Eq. (1) in different subregions of the plane (γ 12 , γ 11 ), for N = 2, γ 22 = γ 11 and γ 12 = γ 21 , predicted by the analysis in "Appendix A". This graph is obtained by fixing the following parameter values: p 1 = p 2 = D 1 = D 2 = L = 1 the height of which remains finite as K becomes arbitrarily narrow, (ii) segregation patterns where the height of both species becomes arbitrarily high as K becomes arbitrarily narrow, denoted by S ∞,∞ S , (iii) segregation patterns where the height of just one species becomes arbitrarily high as K becomes arbitrarily narrow but the other remains at finite height, denoted by S 1,∞ S , (iv) aggregation patterns, S ∞ A , where the height of both species becomes arbitrarily high as K becomes arbitrarily narrow, and (v) the spatially homogeneous solution S H . Figure 6b shows the parameter regions in which the analysis from "Appendix A" predicts we should see these various solutions. Notice that there are regions in which we have two-, three-, and even four-fold stability. These calculations are verified numerically in Figs. 7 and 8. In particular, Figs. 7 and 8 show that, as α becomes smaller, so the numerical results become closer to our analytic predictions.
As shown in Fig. 6b, when species exhibit mutual attraction (γ 12 < 0), our analysis predicts two stationary states: the homogeneous distribution S H and the aggregation pattern S ∞ A . In particular, if γ 12 < 0 and species show mutual avoidance, i.e. γ 11 > 0, there always exists a region in the parameter space in which both stationary states, S H and S ∞ A , are stable. However, if the magnitude of self-avoidance γ 11 is relatively weaker than the rate of mutual-attraction γ 12 , aggregation is more favored than the homogeneous distribution. In this case, S ∞ A is the only stable steady state, while the S H solution is unstable.
On the other hand, in the mutual-and self-attraction case (γ 12 < 0, γ 11 < 0), bistability between the homogeneous distribution S H and the aggregation pattern S ∞ A is observed as long as the magnitudes of γ 12 and γ 11 are sufficiently small. However, if the rates of mutual and self-attraction become stronger, aggregation is favoured over the homogeneous distribution. Consequently, as the magnitudes of γ 11 and γ 12 increase, the homogeneous solution S H loses stability.
The scenario becomes even richer when γ 12 > 0. In particular, if the species exhibit mutual avoidance (γ 12 > 0) and self-avoidance (γ 11 > 0), the stable steady states predicted by our analysis are the homogeneous solution S H and segregation pattern S 2,2 S . When the strength of self-repulsion (γ 11 ) is relatively stronger than the tendency to avoid individuals from the other species (γ 12 ), the homogeneous distribution is favoured over aggregation with conspecifics, so that S H is the only stable steady state. However, if the rate of mutual avoidance γ 12 increases, the tendency to avoid individuals from the foreign species promotes the formation of spatial distributions in which the two species are segregated into distinct sub-regions of space. Indeed, Fig. 6b shows that as γ 12 increases, the segregation pattern S 2,2 S acquires stability. However, as long as the magnitude of self-avoidance is sufficiently strong, the homogeneous distribution remains stable. Indeed, we observe that there is a parameter region in which the system shows bistability between S H and S 2,2 S . Finally, if the strength of mutual avoidance γ 12 becomes sufficiently stronger than the propensity to avoid conspecifics, segregation becomes more favored over the homogeneous distribution. Indeed, as γ 12 increases, S H loses its stability.
In the mutual avoidance (γ 12 > 0) and self-attraction (γ 11 < 0) scenario, the stable states predicted by our analysis include S H (homogeneous) and S 2,2 S (territory-like segregation) as before, but also S ∞,∞ S (self-aggregated species that are segregated from one another) and S 1,∞ S (segregated species where only one population is selfaggregated). If the magnitudes of self-attraction γ 11 and mutual avoidance γ 12 are sufficiently small, the homogeneous distribution, S H , is also stable. However, for small values of γ 11 , as the rate of mutual avoidance γ 12 increases, we observe the same scenario discussed above: S 2,2 S gains stability and there exists a region in the parameter space in which both S 2,2 S and S H are stable. Finally S H loses stability as γ 12 increases further. We also observe that high rates of self-attraction γ 11 favour the formation of sub-regions with high densities of individuals. Therefore, when the magnitude of γ 11 is strong, S ∞ S and S ∞ H solutions are favored over the homogeneous distribution S H and the inhomogeneous distribution S 2,2 S , which become unstable. Finally, we verify this multi-stability numerically for small α, with results shown in Figs. 7 and 8. As in the γ 11 = γ 22 = 0 cases, the numerics follow our analytic predictions well, giving better approximations for smaller α.

The steady states in the local limit
In the previous section, we found piecewise constant energy minimisers of the local limit of Eq. (1). These can attain only a discrete set of values. Here, we confirm this observation by showing that, on each subinterval where the solution is differentiable, it must be constant.
For N = 2 we prove that the image of any minimum energy solution must lie in a finite set. This proof works for any parameter values D i and γ i j . We were not, however, able to prove this result in full generality for arbitrary N . Nonetheless, we do provide a method for constructing a proof for any particular set of parameter values, and put these ideas into practice in some example cases where N = 3.

The general setup
Let K (x) = δ(x), the Dirac delta function with mass concentrated at x = 0. Then in one spatial dimension Eq. (1) becomes Any local minimum energy solution to Eq. (54) is given by a set of functions We therefore require that, on any subinterval where u i (x) = 0, which implies that Equation (56) can be written in matrix form as and u = (u 1 , . . . , u N ) T . Equation (57) holds on each subinterval where u i (x) = 0. We wish to show that differentiable solutions are necessarily constant. Equation (57) only has a nontrivial solution if either det(A 1 ) = 0 or ∂u ∂ x = 0. The latter means that u is constant, so we need to investigate the condition det(A 1 ) = 0.

The case N = 2
To make things simple, we begin by focusing on the case N = 2. We use the notation A (2) 1 to mean the matrix A 1 (Eq. (57)) for N = 2, so that The condition det(A 1 ) = 0 then implies If u is differentiable then we can differentiate Eq. (59) with respect to x, leading to the following Combining Eq. (60) with the first row of the vector equation A Then det(A (2) 1 ) = 0, det(A (2) 2 ) = 0 is a system of two simultaneous equations in two unknowns. These have at most three solutions, as we show in "Appendix B".
The exact form of these solutions is rather cumbersome, so we omit writing them down explicitly. However, it is instructive to give a simple example, which we do in the case γ 11 = γ 22 = 0. Here, there is a single solution to det(A 2 ) = 0 of the following form Regardless of whether or not we impose the condition γ 11 = γ 22 = 0, the solution set of (u 1 , u 2 ) is a finite set. Therefore each differentiable part of a solution of Eq. (56) is constant.

The case N = 3
We now show how to extend the arguments of Sect. 5.2 to the N = 3 case. The expressions become too complicated in N = 3 to give a complete analysis, so we instead give some examples to demonstrate how one can ascertain whether or not image of u(x) is contained in a finite set. Similar to the strategy for N = 2, the aim is to construct a system of equations that constrain the possible solutions for u(x). For N = 3, this involves constructing three equations, which each take the form det(A , 2, 3}), and showing that this set of simultaneous equations has a finite number of solutions. Whilst for N = 2, we were able to calculate the number of solutions exactly by solving polynomial equations, this is not possible for N = 3 as the polynomials are usually of order 5 or more (Stewart 2015). Instead, we use the theory of Gröbner bases to prove the solution set is finite.

Example 1
For this example, we let D i = 1, γ ii = 0, γ 12 = γ 21 = γ 23 = γ 32 = 2, and γ 13 = γ 31 = 1. Then Since det(A 1 ) = 0, we have Again, assuming u is differentiable, we can differentiate Eq. (64) with respect to x which leads to the following Combining Eq. (65) with the first two rows of A Once again, we have that det (A (3) 2 ) = 0, leading to the following polynomial equation Differentiating Eq. (67) with respect to x gives where Combining Eq. (68) with the first two rows of A We now have a set of three polynomials such that the image of u(x) must lie on the common zeros of this set. In the N = 2 case (Sect. 5.2), we had just two polynomials, both of which were cubics, thus it is possible to find formulae for the common zeros. Here, however, we have a polynomial of degree six (det(A 3 )). Since there is no general solution to a sixth degree polynomial (Stewart 2015), we cannot solve the system det(A 2 ) = 0, det(A 3 ) = 0 directly.
Instead, we use a classical result from algebraic geometry, which says that the number of common zeros of S is finite iff for each i ∈ {1, 2, 3}, the Gröbner basis of the ideal I (S) generated by S contains a polynomial whose leading monomial is a power of u i (Adams and Loustaunau 1994). Computation of the Gröbner basis of an ideal generated by a set of polynomials is an algorithmic procedure that is encoded into various mathematical packages, such as Mathematica (Wolfram et al. 1999) or Macauley2 (Eisenbud et al. 2013).
We use Mathematica to calculate the Gröbner basis of I (S). The result is a set of five polynomials whose leading monomials are β 1 u 19 3 , β 2 u 2 u 2 3 , β 3 u 2 2 u 3 , β 4 u 4 2 and β 5 u 1 , where β 1 , . . . , β 5 are constants (some of which are of the order 10 26 so we refrain from writing down their exact numerical values). For each i, there is a polynomial in the Gröbner basis whose leading monomial is a power of u i . Therefore, the common zeros of S are finite and the image of u(x) is contained in a finite set. Since we have assumed u(x) is differentiable, it must also be constant.

Example 2
In the previous example, we were able to show that the image of u(x) is contained in a finite set by showing it lies on the intersection of three polynomials, which is the minimum number of polynomials required in the case N = 3. However, sometimes three polynomials is not enough. Here, we detail an example which requires the construction of five polynomials to ensure the intersection of their zeros is a finite set.
Since det(A (3) 1 ) = 0, we have Differentiating Eq. (75) with respect to x leads to the following Combining Eq. (76) with the first two rows of A Once again, we have that det (A (3) 2 ) = 0, leading to the following polynomial equation Differentiating Eq. (78) with respect to x gives where Combining Eq. (79) with the first two rows of A We now have a set of three polynomials S = det(A 3 ) , such that the image of u(x) must lie on the common zeros of this set. The Gröbner basis of I (S) contains eight polynomials whose leading terms are β 1 u 2 u 9 3 , β 2 u 2 u 8 3 , β 3 u 2 u 8 3 , β 4 u 2 u 8 3 , β 5 u 2 u 8 3 , β 6 u 2 u 8 3 , β 7 u 2 u 8 3 , β 8 u 2 u 8 3 for constants β 1 , . . . , β 8 . Here, the Gröbner basis of I (S) does not contain a polynomial a with leading monomial that is a power of u i for any i = 1, 2, 3, so the common zeros of S do not necessarily form a finite set. Therefore we need to search for further polynomials on which the solution lies, to see if we can constrain the solutions into a finite set.
To this end, we combine Eq. (76) with the first and the third row of A Since det(A Differentiating Eq. (85) with respect to x gives where Combining Eq. (86) with the second and third row of A We now have a set of five polynomials S = det(A 32 ) , such that the image of u(x) must lie on the common zeros of this set. The Gröbner basis of I (S) consists of seven polynomials whose leading monomials are 32768u 9 3 , 12u 2 u 2 3 , 6u 2 2 u 3 , 96u 3 2 , −18u 1 , 18u 1 u 2 , 12u 2 1 . Since, for each i ∈ {1, 2, 3}, this set contains a power of u i , the common zeros of S are finite, and therefore the image of u(x) is contained in a finite set. Hence if u(x) is differentiable, it must be constant.

Discussion
A central aim of mathematical biology is to predict emergent features of biological systems, using dynamical systems models. Stable steady states provide an important class of emergent features, so identification of these is a key task of mathematical biology. However, for nonlinear PDEs, this is not usually an easy task (Robinson and Pierre 2003). Indeed, often this is replaced by the more tractable task of examining a system's behaviour close to the constant steady state, which enables linear or weakly nonlinear approximations. But it is the behaviour far away from the constant solution that is interesting biologically, as that is where the patterns exist that we perceive in biological systems.
Here, we have detailed a novel method to help find local minimum energy states, which are Lyapunov stable, in a system of nonlocal advection-diffusion equations for modelling N species (or groups) of mobile organisms, each of which move in response to the presence of others. Our study system is closely related to (and often directly generalises) a wide variety of previous models, including those for cell aggregation (Carrillo et al. 2018) and sorting (Burger et al. 2018), animal territoriality (Potts and Lewis 2016a) and home ranges (Briscoe et al. 2002), the co-movements of predators and prey (Di Francesco and Fagioli 2016), and the spatial arrangement of human criminal gangs (Alsenafi and Barbaro 2018). Therefore our results have wide applicability across various areas of the biological sciences.
Whilst analytic determination of stable steady states in PDEs remains a difficult task in general, numerical analysis always leaves the question open of whether one has found all possible steady states or whether there are more that the researcher has simply not stumbled upon. To help guide numerical investigations, we have constructed a method, combining heuristic and analytic features, that gives clues as to where stable steady states might be found in multi-species nonlocal advection-diffusion systems. We have demonstrated in a few examples that numerical investigations agree with the predictions of our method. Whilst our method does not give an analytic solution, it should be a valuable tool for finding stable steady states in biological models that can be modelled by nonlocal advection-diffusion systems.
Our method relies on constructing an energy functional for the PDE system. We were only able to do this in the case γ i j = γ ji for all i, j ∈ {1, . . . , N } and assuming that the kernel K is identical for all species. These constraints mean that each pair of species (or populations or groups) respond to one another in a symmetric fashion, either mutually avoiding or mutually attracting with identical strengths of avoidance or attraction, respectively. This generalises a recent result of Ellefsen and Rodríguez (2021) who construct an energy functional for the case where γ i j = 1 for all i, j ∈ {1, . . . , N }. We conjecture that this energy functional could be used to prove that the attractor of our study system is an unstable manifold of fixed points. However, we were unable to prove this here, so encourage readers to take on this challenge.
Whilst it may be possible to construct energy functionals in some example situations where γ i j = γ ji for some i, j, or where the kernel is not identical for all species (we leave this as an open question), we expect that it is not possible in general, since there are situations where the numerical analysis suggests the attractors do not consist of stable steady states, but patterns that fluctuate in perpetuity (Potts and Lewis 2019).
Perhaps the simplest situation where this has been observed is for N = 2, γ 11 , γ 22 < 0, and γ 12 < 0 < γ 21 (Giunta et al. 2021a), whereby both populations aggregate and one 'chases' the other across the terrain without either ever settling to a fixed location. Furthermore, to keep our analysis as simple as possible, we only applied the techniques of Sect. 4 to some concrete examples in n = 1 spatial dimension. Nonetheless, there is no a priori reason why the techniques in Sect. 4 could not be extended to higher dimensions in the future.
Whilst our method is designed for application to models of nonlocal advection, for which there are existence and regularity results (Giunta et al. 2021a), it works by examining the local limit of stable solutions. The reason for this is that these solutions are piecewise constant, so we can constrain our search for the minimum energy, enabling minimisers to be found analytically. The disadvantage is that the local limit of stable solutions is not itself the steady state solution of a well-posed system of PDEs: in the local limit, Eq. (1) becomes ill-posed. More precisely, it is unstable to arbitrarily high wavenumbers whenever the pattern formation matrix has eigenvalues with positive real part. Nonetheless, we have shown that the local limit of minimum energy solutions to the nonlocal problem is a useful object to study, even if it may not itself be the steady state solution of a system of PDEs.
It would be cleaner, however, if we were able to develop theory that did not require taking this local limit. For N = 1, Potts and Painter (2021) developed techniques that are analogous to the ones proposed here but in discrete space. In this case, the actual stable steady states of the discrete space system become amenable to analysis via an energy functional approach similar to the one proposed here. However, generalisations of this technique to N > 1 do not appear to be trivial from our initial investigations.
Another possible way forward is to use perturbation analysis, starting with the minimum energy solutions from the local limit, studied here, and perturbing them to give solutions to the full nonlocal system. One could then minimise the energy across this class of perturbed solutions (which would no longer be piecewise constant) to find stable steady states of the nonlocal system in Eq. (1). This is quite a nontrivial extension of the present methods, which we hope to pursue in future work. One possible avenue might be to use a kernel that allows the non-local model to be transformed into a higher-order local model (Bennett and Sherratt 2019;Ellefsen and Rodríguez 2021).
Figures 3, 5, 7 and 8 show numerical bifurcation analysis of our system in certain examples. This naturally leads to questions about the nature of these bifurcations. In particular, the discontinuity in amplitude that occurs as the constant steady state loses stability is something that is also seen with subcritical pitchfork bifurcations. In this case, the stable branches may be joined to one another by an unstable branch, or some more complicated structure. It would be valuable to investigate analytically whether this is the case. Standard tools include weakly non-linear analysis and Crandall-Rabinowitz bifurcation theory, both of which have been used successfully for nonlocal advection-diffusion equations (Buttenschön and Hillen 2021;Eftimie et al. 2009).
The system we study assumes that species advect in response to the population density of other species. However, it is agnostic as to the precise mechanisms underlying this advection. Previous studies show that Eq. (1) can be framed as a quasi-equilibrium limit of various biologically-relevant processes, such as scent marking or memory Lewis 2016a, b, 2019). This quasi-equilibrium assumption says, in effect, that the scent marks or memory map stabilise quickly compared to the probability density of animal locations. However, it would be valuable to examine the extent to which these processes might affect the emergent patterns away from this quasi-equilibrium limit. Along similar lines, it would also be valuable to examine the extent to which our results translate to the situation where we model each individual as a separate entity, as in an individual based model (IBM), rather than using a population density function, which is a continuum approximation of an IBM. We have recently begun developing tools for translating PDE analysis to the situation of individual based models, which could be useful for such analysis (Potts et al. 2022).
In summary, we have developed novel methods for finding nontrivial steady states in a class of nonlinear, nonlocal PDEs with a range of biological applications. As well as revealing complex multi-stable structures in examples of these systems, our study opens the door to various questions regarding the bifurcation structure, the effect of nonlocality, and the structure of the attractor. We believe these will lead to yet more significant, but highly fruitful, future work.
We will look for the local minimizers of the following energy functional, where K = δ, in the class of piece-wise constant functions defined as where u c i ∈ R + and S i are subsets of [0, L], for i ∈ {1, 2}. Recall that, by Eq. (5), in Eq. (A2) we require the following constraint where the first equality uses γ 12 = γ 21 and the third equality uses Eq. (A3).
Since the general analysis of this case is not straightforward, we instead set γ 11 = γ 22 and fix the other parameter values as p 1 = p 2 = D 1 = D 2 = L = 1. Therefore Eq. (A4) becomes In the following, we will look for the minimizers of Eq. (A5) and examine different cases demarcated by the signs of γ 11 and γ 12 .
To search for the local minimizers of the energy in Eq. (A5), we then define To constrain our search, notice that Eq. (A3), p i = 1 and |S i | ≤ L = 1 imply that We analyse E(u c 1 , u c 2 ) (Eq. (A6)) under the constraint in Eq. (A7), first in the region where |S 1 | + |S 2 | ≤ 1 and then in the region where |S 1 | + |S 2 | > 1. By combining these results we will have a complete picture of the local minima of E(u c 1 , u c 2 ). Note that by Eq. (A3), the case |S 1 | + |S 2 | ≤ 1 is equivalent to By analysing the partial derivatives of E(u c 1 , u c 2 ) in the region of the (u c 1 , u c 2 )-plane defined by Eq. (A8), one can check that there are no local minima in this region. Furthermore, E(u c 1 , u c 2 ) → ∞ as either u c 1 → ∞ or u c 2 → ∞. Therefore any minima in this region must lie on the boundary, 1/u c 1 + 1/u c 2 = 1 (solid black line in Fig. 2). Analysis of the partial derivative of E(u c 1 , u c 2 ) on this boundary shows that E(u c 1 , u c 2 ) has a unique minimum point, given by This is also a local minimum of the region defined by Eq. (A8). This can be shown by performing a Taylor expansion of E(u c 1 , u c 2 ) about the point M S . Since the slope of the line tangent to the curve 1/u c 1 + 1/u c 2 = 1 in M S is −1, we choose two constant, and δ, such that + δ ≥ 0 and the Taylor expansion gives where the inequality uses γ 11 > 0, + δ ≥ 0. However, since the point M S lies on the boundary curve |S 1 | + |S 2 | = 1, we do not yet know whether it is a minimum for the whole admissible region defined by Eq. (A7) (white region in Fig. 2). To this end, we examine whether M S is a minimum of E(u c 1 , u c 2 ) (Eq. (A6)) in the region where |S 1 |+|S 2 | > 1. By Eq. (A7), the condition |S 1 |+|S 2 | > 1 is equivalent to 1/u c 1 + 1/u c 2 > 1. Therefore we have the following constraints Since |S 1 ∩ S 2 | = |S 1 | + |S 2 | − 1, when |S 1 | + |S 2 | > 1 the function E(u c 1 , u c 2 ) (Eq. (A6)) can be rewritten as E(u c 1 , u c 2 ) = ln(u c 1 ) + ln(u c 2 ) + 1 2 γ 11 (u c 1 + u c 2 ) + γ 12 u c 1 u c 2 |S 1 ∩ S 2 | = ln(u c 1 ) + ln(u c 2 ) + 1 2 γ 11 (u c 1 + u c 2 ) + γ 12 u c 1 u c 2 (|S 1 |+|S 2 |−1), = ln(u c 1 ) + ln(u c 2 ) + 1 2 γ 11 (u c 1 + u c 2 ) + γ 12 u c 1 u c where the third equality uses |S i |= 1 u c i . To verify whether M S is also a minimum on the part of the domain given by Eq. (A10), we perform a Taylor expansion of E(u c 1 , u c 2 ) in a neighbourhood of M S within the region 1/u c 1 + 1/u c 2 < 1. Since the slope of the tangent line to the curve 1/u c 1 + 1/u c 2 = 1 at the point M S is −1, we choose two arbitrary constants, and δ, such that + δ ≤ 0. Then Taylor expansion of E(u c 1 , u c 2 ) is E(2 + , 2 + δ) ≈ E(2, 2) + ∂ u c 1 E(2, 2) + ∂ u c 2 E(2, 2)δ = E(2, 2) + 1 2 (1 + γ 11 − 2γ 12 )( + δ) if γ 11 < 2γ 12 − 1, where the inequality uses + δ ≤ 0. Next, we look for any other minima in the region defined by Eq. (A10). By analysing first partial derivatives, one can show that there are no local minima of E(u c 1 , u c 2 ) (Eq. (A11)) in the interior of this region. Therefore any local minima must occur on the boundaries. On the part of the boundary given by u c i = 1, for i = 1, 2, there is a unique minimum at This is also a local minimum of the region defined by Eq. (A10). This can be shown by performing a Taylor expansion of E(u c 1 , u c 2 ) about the point M H , to give E(1 + , 1 + δ) ≈ E(1, 1) + ∂ u c 1 E(1, 1) + ∂ u c 2 E(1, 1)δ = E(1, 1) + 1 + 1 2 γ 11 ( + δ) where the inequality uses γ 11 > 0, ≥ 0 and δ ≥ 0. Here, and δ are chosen to be non-negative so that we remain in the u c i ≥ 1 region (Fig. 2). Therefore, if γ 11 > 2γ 12 − 1, E(u c 1 , u c 2 ) (Eq. (A5)) has a unique minimum, given by M H . Whilst if 0 < γ 11 < 2γ 12 − 1, then E(u c 1 , u c 2 ) has two local minima, given by M H and M S . Finally, we write down the functions u i (x) (Eq. (A2)) which locally minimize the energy E[u 1 , u 2 ] (Eq. (A4)). If (u c 1 , u c 2 ) = M H then u 1 (x) = u 2 (x) = 1, the homogeneous steady state, which we denote by S H . If (u c 1 , u c 2 ) = M S then with |S i | = 1/2, for i = 1, 2, and |S 1 ∩ S 2 | = 0, denoted by S 2,2 S . In conclusion, if γ 11 > 2γ 12 − 1, the energy E(u c 1 , u c 2 ) (Eq. (A4)) has a unique minimum, given by S H . However, if 0 < γ 11 < 2γ 12 − 1 the energy has two local minima, given by S H and S 2,2 S . Furthermore, linear stability analysis (Eq. (14)) suggests that when α tends to zero, the homogeneous steady state is stable if γ 11 > γ 12 − 1. This gives rise to the diagram of analytically-predicted steady states given by the red and black lines in Fig. 7a.

A.2 Mutual attraction ( 12 < 0)
In this section, we analyze the local minimizers of the energy (Eq. (A4)) for γ 12 < 0, γ 11 ∈ R and γ 12 = γ 21 . We observe that the energy in Eq. (A4) decreases as |S 1 ∩ S 2 | increases, whilst keeping everything else constant. Therefore if we keep |S 1 | and |S 2 | unchanged, then |S 1 ∩ S 2 | is maximised when either S 1 ⊆ S 2 or S 2 ⊆ S 1 , so that |S 1 ∩ S 2 |= min i |S i |. Thus by repeating the same argument presented in Sect. 4.2.1 for γ 12 < 0 and γ 11 = 0, we see that E[u 1 , u 2 ] → −∞ as min{u c 1 , u c 2 } → ∞. As we approach this limit, u c 1 , u c 2 become arbitrarily large, so u 1 and u 2 (Eq. (A2)) become arbitrarily high, arbitrarily narrow functions with overlapping support. We will denote the limit of this solution by S ∞ A . One can also show, using a very similar argument to "Appendix A.1" (details omitted), that the homogeneous steady state, S H , is the only other possible local minimiser of the energy that satisfies u c i ≥ 1, for i = 1, 2, and this is only a local minimum when γ 12 > −γ 11 − 2. However, linear stability analysis (Eq. (13)) suggests that, in the limit as α tends to zero, the homogeneous steady state is linearly stable only if γ 12 > −γ 11 − 1. Therefore, any time S H is linearly stable, it is also a local energy minimiser within the set of functions given by Eq. (A2). These results give rise to the diagram of analytically-predicted steady states given by the red and black lines in Fig. 7b-c.