Weakly nonlinear analysis of a two-species non-local advection-diffusion system

Nonlocal interactions are ubiquitous in nature and play a central role in many biological systems. In this paper, we perform a bifurcation analysis of a widely-applicable advection-diffusion model with nonlocal advection terms describing the species movements generated by inter-species interactions. We use linear analysis to assess the stability of the constant steady state, then weakly nonlinear analysis to recover the shape and stability of non-homogeneous solutions. Since the system arises from a conservation law, the resulting amplitude equations consist of a Ginzburg-Landau equation coupled with an equation for the zero mode. In particular, this means that supercritical branches from the Ginzburg-Landau equation need not be stable. Indeed, we find that, depending on the parameters, bifurcations can be subcritical (always unstable), stable supercritical, or unstable supercritical. We show numerically that, when small amplitude patterns are unstable, the system exhibits large amplitude patterns and hysteresis, even in supercritical regimes. Finally, we construct bifurcation diagrams by combining our analysis with a previous study of the minimisers of the associated energy functional. Through this approach we reveal parameter regions in which stable small amplitude patterns coexist with strongly modulated solutions.


Introduction
Spontaneous pattern formation occurs throughout nature [22], with examples ranging from animal coat patterns [35] to territory formation [27], cell sorting [6] and swarm aggregation [33].Therefore uncovering and analysing the mechanisms behind pattern formation is a central challenge in the life sciences where applied mathematics can play a role.Typically, research into pattern formation proceeds first by assessing which parameters may cause patterns to emerge spontaneously from a homogeneous steady state, using linear pattern formation analysis, sometimes called 'Turing pattern analysis' [35].This determines whether patterns may emerge at short times from arbitrarily small perturbations.However, it is also important biologically to show whether these patterns are stable.One approach to pattern stability is via weakly nonlinear analysis: a stable supercritical bifurcation branch suggest that asymptotic patterns will emerge continuously as the bifurcation parameter is changed, whereas an unstable subcritical branch suggests that large amplitude asymptotic patterns may appear abruptly as the bifurcation point is crossed, their amplitude being a discontinuous function of the bifurcation parameter.This discontinuity in amplitude with respect to parameter change indicates that a biological system might suddenly change its behaviour in a dramatic fashion with only a small change in the underlying mechanisms.
Many biological mechanisms generate attractive or repulsive forces governing phe-nomena such as chemotaxis ( [14,21]), bacterial orientation ( [2]), swarms of animals ( [29]), and motion of human crowds ( [20]).These mechanisms are driven by electrical, chemical or social interactions.These interactions arise from individual organisms collecting information from their environment, such as the presence of other individuals, food or chemicals.After gathering information, individuals move towards regions that contain important components for survival or move away from less favourable areas, thus creating spatially inhomogeneous distributions of individuals, which may have a certain degree of regularity in space and/or time (e.g.[33,28]).This process of acquiring information from the environment is generally nonlocal, as motile organisms are usually able to inspect a portion of their environment, either by prolonging their protrusions, as in the case of cells [8], or by using their sight, hearing or smell, as with animals [26].
In recent years there has been an increasing interest in the mathematical modelling of nonlocal advection as a movement model with nonlocal information [5,33,6,10,8].Recently, the following class of nonlocal advection-diffusion equations was proposed as a general model of interacting populations [28] γ ij ∇(K * u j ) , i = 1, . . ., N.
Here, u i (x, t) denotes the density of population i at position x and time t, for i ∈ {1, . . ., N } and D i > 0 is the diffusion rate of u i .Individuals can detect the presence of other individuals, whether conspecifics or not, over a spatial neighborhood described by spatial averaging kernel K, which is a symmetric, non-negative function modelling the sensing range.The term K * u j denotes the convolution between K and u j and describes the nonlocal interactions of u i with u j .The parameters γ ij are the inter/intra-species interaction parameters, giving the density-dependent rate at which species i advects towards (if γ ij < 0), or away from (if γ ij > 0), species j.
Model (1) implicitly focuses on time scales whereby birth and death processes are negligible.Nonetheless, it has a wide range of possible applications in that it generalizes a variety of existing models describing many different phenomena, such as animal home ranges [5], territory formation [15,27,30], and cell sorting [6].On the mathematical side, well-posedness of System (1) was analyzed in [17] and [23].
When the kernel K is sufficiently smooth, [17] shows that the system admits classical, positive and global solutions in 1D dimension, and local strong solutions in any higher dimension.When the kernel is non-smooth, in [23] it is proven that System (1) has weak solutions that exist globally in time.
From the perspective of pattern formation, numerical analysis shows that System (1) exhibits a great variety of spatio-temporal patterns, depending on the model parameters.These include segregated and aggregated stationary patterns, periodic time oscillating solutions, and aperiodic spatio-temporal behaviours [28], [17], [9].
In many cases the system admits an energy functional [18,9], which can be used to gain analytic insight into the steady asymptotic patterns that can form from this system.Although [18] focused on the N = 2 case, the methods are more generally applicable in principle.
Here, we perform a bifurcation analysis of one of the cases analyzed in [18], namely where N = 2, γ ij = γ ji and γ ii = 0.For simplicity, we also assume that D 1 = D 2 .
We use weakly nonlinear analysis to derive the equations governing the amplitude of the stationary solutions.Through analysis of the amplitude equations, we determine the nature of bifurcations generating branches of non-homogeneous solutions from a homogeneous state, then recover the shape of the non-homogeneous solutions and their stability.We validate our results through numerical analysis, setting K to be the top-hat distribution [18].Finally, we combine our results with results of [18] that were derived from an energy principle, to construct bifurcation diagrams that incorporate all the existing analysis of this system.
An interesting feature of our analysis is that the equation governing the modulation of small-amplitude patterns is not always the real Ginzburg-Landau (GL) equation.This contrasts with many examples of weakly nonlinear analysis, where the GL equation provides the amplitude of the stationary pattern and its stability: in subcritical regimes, the pattern solution is always unstable; in supercritical regimes, a periodic pattern is stable if its wavenumber lies within the Eckhaus band; [34,22,3,4,19,11].In our case, the real GL equation does not always provide a correct description of the pattern near the onset.This is because our system possesses a conservation law, i.e. mass is conserved for all time.This conservation law gives rise to a large-scale neutral mode (the zero mode) that can affect the stability of the pattern, so must be included into the analysis [12,24].Therefore, the resulting amplitude equations will consist of the GL equation coupled to an equation for the large-scale mode.
In [24] the authors used symmetry and scaling arguments to derive the amplitude equations governing systems with a conserved quantity.They proved that there exist stable stationary solutions in the form of strongly modulated patterns (i.e.patterns that consist of multiple Fourier modes), and these exist away from the branch that bifurcates from the constant steady state.The existence of strongly modulated patterns for System (1) has also been shown in [18] by analyzing the minimizers of an energy functional associated with the system.Here we build on this by investigating the existence and stability of small amplitude patterns, and showing that when these solutions are unstable, the system evolves towards either large amplitude or strongly modulated patterns.In addition, our analysis shows that, in some parameter regions, stable small amplitude patterns can coexist with stable strongly modulated solutions.
A similar two-species aggregation model was studied recently in [6].Their model differs from our model (2) in regard of the diffusion term.In [6] the terms D∂ xx u i for i = 1, 2 are replaced by density dependent diffusion terms D∂ x (u i ∂ x (u 1 + u 2 )).The pattern forming mechanism is similar to our model, however, the arising aggregations have compact support.
This paper is organised as follows.Linear stability analysis is given in Section 2 and a weakly nonlinear analysis in Section 3. In these two sections, the analysis is carried out with a generic kernel, in order to provide some general results that can be used for future works.Section 4 focuses on detailed analysis where K is the top-hat distribution.We analyse the amplitude equations, recover the bifurcation diagrams and compare analytical results with numerical solutions.We finally combine the analysis performed here with the results obtained in [18] to recover more exhaustive pictures of the bifurcation diagrams.In Section 5, we outline further extensions of this work and discuss possible applications of our results to natural systems.

Linear stability analysis
We consider System (1) with two interacting populations, u 1 and u 2 , that either mutually avoid or attract with the same strength (i.e.γ 12 = γ 21 ).We set γ := γ 12 = γ 21 and fix D 1 = D 2 =: D, and γ 11 = γ 22 = 0. Therefore, System (1) reads as ( We work on the one dimensional spatial domain Ω = − l 2 , l 2 and impose periodic boundary conditions We consider an even and non-negative kernel K such that where the constant α denotes the sensitivity radius.We assume that α < l/2.Due to the periodic boundary conditions, we also assume that K(x) is wrapped around periodically over the domain.
The periodic boundary conditions (Equation ( 3)) ensure that in System (2) the total mass of each population u i is conserved in time.Indeed the following identities Hence where the constant p i denotes the size of population u i , for i = 1, 2.
Equation (6) implies that system (2) has a unique equilibrium point given by

Nondimensionalization
We start our analysis by rescaling the original system (2) using the following non-dimensional coordinates and variables Note that, instead of α, one could have rescaled using any other constant that is proportional to the standard deviation of K(x) instead, which may be useful if K(x) does not have compact support, for example.
In the non-dimensional spatial domain, we define the following kernel By Equation ( 9), we see that Supp( K) = [−1, 1] and that By ( 8) and ( 9), it follows that the convolution product becomes where * ∼ denotes the convolution operator in the rescaled spatial coordinate.
By substituting Equations ( 8), ( 9) and (11) in Equations ( 2), we obtain the following non-dimensional system where x ∈ − l 2α , l 2α .By the relations in Equation ( 8), the boundary conditions now read as: The boundary conditions (Equation ( 13)) imply that the total mass of each population ũi is conserved in time.Therefore, for i = 1, 2 and all t ≥ 0, the following identities hold where the second equality uses the identities in Equation ( 8) and the third equality uses Equation (6).By Equation ( 14) it follows that the non-dimensional system in (12) has a unique equilibrium point given by To simplify the notation, we define γ := γ lD and L := l α , and by dropping the tildes, the non-dimensional system (12) reads as where x ∈ − L 2 , L 2 .The boundary conditions for System (16) read as:

Linear stability analysis
We now perform a linear stability analysis of system (16) about the equilibrium point (see Equation ( 15)).To this end, we consider a perturbation of the homogeneous solution (18) of the following form subject to boundary conditions (17), where u (0) is a constant vector, λ ∈ R is the growth rate and q is the wavenumber of the perturbation.By substituting Equation (19) into Equation ( 16) and neglecting nonlinear terms, we obtain the following eigenvalue problem where and where the second equality uses the fact that K(x) is an even function and then K(x) sin(qx) is an odd function.
The wavenumbers q must be chosen in such a way that the periodic boundary conditions in Equation ( 17) are satisfied, and thus we have a discrete set of admissible wavenumbers given by The equilibrium ū (Equation ( 18)) is unstable when λ ± (q m ) > 0 for some m ∈ Z ≥0 .Note that λ ± (q 0 ) = 0 so the system never becomes unstable at wavenumber q 0 .For m > 0, if K(q m ) ̸ = 0, we denote by γ ± m the instability thresholds of the wavenumber q m , which are defined as Therefore the equilibrium ū (Equation ( 18)) is unstable when In the following section, we will perform a weakly nonlinear analysis to study the evolution of the perturbation w when the equilibrium ū becomes linearly unstable.We will adopt γ as bifurcation parameter and denote by q c the first admissible wavenumber that is destabilized as |γ| is increased.By Equation ( 25), we note the critical wavenumber q c is defined as where the set I is defined in (24).We also underline that q c depends on the choice of kernel K and may not be unique.We will denote by γ ± c the corresponding bifurcation thresholds, that is

Amplitude equations
In this section we perform a weakly nonlinear analysis based on the method of multiple scales.Close to the threshold of instability, that is in the weakly non-linear regime, we will use an expansion technique to recover an approximated solution, characterized by a slowly varying amplitude, and the equations governing the amplitude of the solution.Through the analysis of these equations (usually referred to as amplitude equations), we recover the amplitude and stability of the stationary solutions.
The idea behind the multiple scale method comes from the observation that, just above an instability threshold, a nonlinear state is given by a superposition of modes whose wavenumbers q lie in a narrow band q − ≤ q ≤ q + (see [13] Cap 6).The resulting nonlinear state is a solution governed by one or more unstable modes and characterized by an amplitude that varies slowly in space, due to the superposition of modes with almost identical wavenumbers.Also, the amplitude evolves slowly in time because, close to the onset of instability, all growth rates are small.
Generally just beyond a bifurcation threshold, if the band of unstable wavenumbers [q − , q + ] around q c has width O(ε), where ε ≪ 1, the positive growth rates are O(ε 2 ).Therefore, the solution evolves as where X = εx is a long spatial scale, T = ε 2 t is a slow temporal scale, Ã(X, T ) is a complex function and denotes the slow modulation of the critical mode e iqcx , and Ã * is the complex conjugate of Ã.Also, in the limit of ε → 0, this solution must satisfy the boundary conditions in Equation (17).
However, in systems with a conservation law, so that λ(0) = 0, long-scale modes evolve on long timescales, and must be included in the analysis (see also [24]).Therefore solutions to System ( 16)-( 17) evolve as where B(X, T ) is a real function and denotes the slow modulation of the mode corresponding to the zero wavenumber, q = 0.
Recall that the homogeneous steady state is linearly stable for γ − c < γ < γ + c , and becomes unstable for γ < γ − c or γ > γ + c .In the following Theorem, we derive an approximation of the solutions close to the instability thresholds (γ ≈ γ + c or γ ≈ γ − c ) and the equations governing the amplitude of the solutions.Since the analysis is broadly the same, we do not distinguish between γ + c and γ − c and use γ c to denote both the thresholds.This Theorem also shows that the ansatz in Equation (30) correctly describes solutions in the weakly nonlinear regime.
Proof.Recall the definition of w from Equation (19).Separating the linear part from the non linear part, System ( 16) can be rewritten as where the actions of linear operator L γ and the non-linear operator Q γ on the vectors r = (r 1 , r 2 ) T and s = (s 1 , s 2 ) T are defined as Choosing γ such that γ − γ c ∼ ε 2 , we write the following expansion From the definition of ε, it follows that either γ (2) = γ c or γ (2) = −γ c .In particular, ).We then employ the method of multiple scales and adopt a long spatial scale X = εx and multiple temporal scales T 1 , T 2 , . . .such that As ε → 0, temporal and spatial derivatives decouple as We employ a regular asymptotic expansion of w in terms of ε where and must satisfy the boundary conditions in Equations (17).
By Equations ( 38) and ( 41), we see that the operators L γ and Q γ in (37) decouple in orders of ε as By substituting Equations ( 41), ( 38), ( 40) and (43) into Equation (36), we obtain Next we collect the terms at each order of ε and obtain a sequence of equations for each w i .At order ε, we obtain the homogeneous linear problem where the function w 1 , has the form as in (42).Therefore, we have: where the second equality uses with K defined in (22).The fourth equality in Equation ( 45) is satisfied if and only if Non-trivial solutions to Equation (47) exist when either the determinant of the matrix is zero or q m = 0. Recalling the definition of q m (24) and γ c (28), we see that nontrivial solutions exist only for q m = q 0 and q m = q c .Therefore, the function w 1 that satisfies this linear problem where the complex conjugate of A, and where K is defined in (22).Since γ c K(q c ) √ ū1 ū2 = ±1 (see Equation ( 28)), ρ can be defined up to a constant.We shall choose the following normalization At this stage, the amplitudes A(X, T 1 , T 2 ) and A 0 (X, T 1 , T 2 ), and the vector ρ 0 are still unknown.
At order ε 2 we obtain the following problem with where the second equality uses Equation (46), the third equality is true because, by Equation (50), the term on the second line is equal to zero, and the fourth equality uses the definition of ρ (Equation (51)).
Notice that any a ̸ = 0 satisfying the condition in (54) is a constant multiple of Therefore Equation ( 52) only has a solution when ρ 01 = ρ 02 = 0 and ∂ T 1 A = 0, that is the amplitude A does not depend on T 1 .From now on, we will denote T 2 by T for simplicity and write A(X, T ) instead of A(X, T 2 ).
Therefore, the linear problem in Equation (52) reduces to Finally, by Equation (56) it follows that the function w 2 , having the form as in (42), is given by where B 0 (X, T ) is a real function and are constant vectors.Notice that ∂ xx L γc [ψ 0 B 0 (X, T )] = 0, for any ψ 0 and B 0 (X, T ).
At order ϵ 3 , we find the following problem where By Equation (50), it follows that the third term of the second equality of Equation (60) is the null vector.In order to simplify the notation, we rewrite Equation (61) as: The linear problem in Equation (60) admits a solution if and only the Fredholm condition ⟨G, a⟩ = 0 is satisfied, where a is defined in Equation (54).Note that the terms G 2 (A 2 ) X e 2iqcx +G 2 (A * 2 ) X e −2iqcx and G 3 A 3 e 3iqcx +G 3 A * 3 e −3iqcx are hortogonal to a. Therefore, the Fredholm condition ⟨G, a⟩ = 0 for Equation (60) gives the following amplitude equation where At order ϵ 4 , we have the following problem Since the function w 4 is as in (42), in Equation (65) all terms independent of x must be equal to zero, that is When ū1 = ū2 , we can choose ψ 01 = ψ 02 and, by setting B := ψ 01 B 0 , we obtain the following amplitude equations where and σ and Λ are given in Equation (64).Notice that ν = δ/ψ 01 (see Equation ( 64)), with ψ 01 = ψ 02 and ū1 = ū2 .On the other hand, if ū1 ̸ = ū2 , Equation (66) is satisfied when ψ 01 = ψ 02 = 0 and (|A| 2 ) XX = 0. □
In the supercritical regime, as the homogeneous steady state becomes unstable, stationary small amplitude patterns emerge and correspond to solutions of Equation (33) with A = a 0 e iϕ , where ϕ ∈ R is the phase of the pattern and the amplitude a 0 is real and must satisfy a 2 0 = σ/Λ.These small amplitude solutions are always stable ( [34]).
Analogously, stationary small amplitude patterns correspond to solutions of Equation (34) with A = a 0 e iϕ and B = 0, where ϕ ∈ R and a 2 0 = σ/Λ.However, in this case the stationary patterns might be destabilized by large-scale modes ( [12]).In the following Proposition we will derive a stability condition for these stationary solutions.
Proof.By Theorem 3.1, if ū1 = ū2 , the amplitude of the stationary solutions to System ( 16) is governed by Equation (34).When σ > 0 and Λ > 0, stationary small amplitude patterns exist and correspond to solutions of (34) with A = a 0 e iϕ and B = 0, where ϕ ∈ R and a 2 0 = σ/Λ.To study the stability of this stationary solution, we consider the following perturbation We substitute the perturbation (70) in Equations ( 34), and by linearizing in a and b we obtain: We consider a perturbation of the form a(X, T ) = e λT (V e iQX + W * e −iQX ) and b(X, T ) = e λT (U e iQX + U * e −iQX ), (72 where λ is the growth rate of the perturbation, U, V, W ∈ C and Q ≥ 0 denotes a spatial mode.Notice that a is a complex perturbation, while b is real.Upon substituting Equations (72) in Equations (71), we obtain the following eigenvalue problem from which we recover the growth rates Recalling that a 2 0 = σ/Λ, a simple calculation shows that λ+ The analysis so far is valid for any non-negative, symmetric kernel K satisfying Equation (4).In the following section, we adopt the top-hat distribution and use the results obtained so far to recover the instability thresholds and to predict the shape of the emerging patterns.
For readers convenience, we conclude this section with Table (1

The top hat distribution
In this section we analyze System (2) with The parameter α, modelling the sensing radius of an organism, is such that α < l/2, where l is the length of the domain.As in Section 2, we will work in dimensionless coordinates, so that our study system is given by Equations ( 16) and the dimensionless averaging kernel is (76)

Linear stability analysis
Linear stability analysis of System ( 16) around the equilibrium point ū = (p 1 , p 2 ) (Equation ( 18)), gives the following eigenvalues (see Equation ( 23)) where Recall that the admissible wavenumbers are q m = 2πm/L, with m ∈ N.

Analysis of the amplitude equations and bifurcations
By Theorem 3.1, when ε = | γ−γc γc | ≪ 1 (where γ c = γ ± c ), the solutions to System (16) have the following form Recall from (32) that the constants ρ 1 , ρ 2 are defined as Note that in the mutual avoidance case (γ > 0), γ c = γ + c > 0 and then ρ 2 < 0, which implies that u 1 and u 2 show a spatial oscillation that is out of phase.On the other hand, in the mutual attraction regime (γ < 0), γ c = γ − c < 0 and then ρ 2 > 0, which means that the spatial pattern for u 1 and u 2 are in phase.Theorem 3.1 also says that A(X, T ) and B(X, T ) are governed by the following equations 2. If ū1 = ū2 , where the coefficients σ, Λ, ν, µ and η are defined in Equation (35) As discussed in Section 3, the sign of Λ determines the type of bifurcation: for Λ > 0 the system exhibits a supercritical bifurcation, while for Λ < 0 the system undergoes a subcritical bifurcation (see also Table (1)).The sign of Λ depends on ū1 , ū2 and on the length of the domain, L (see the definition of Λ in Equation ( 35)).
For γ > 0, if ū1 = ū2 then the qualitative behaviour of Λ(L) remains unchanged as ū1 = ū2 are varied.In fact, Figure 2 As shown in Section 3, if Λ(L) is positive then small amplitude patterns emerge from the homogeneous steady state beyond the bifurcation threshold.These solutions are always stable when ū1 ̸ = ū2 but can be unstable when ū1 = ū2 .Proposition 3.1 shows that when ū1 = ū2 the stability of small amplitude patterns is determined by the coefficients of the amplitude equations in (85) and that, in particular, these solutions are unstable if Γ = Λµ ην − 1 < 0. By using the definitions of Λ, ν, µ and η in Equation ( 35), we recover Γ = (1 + K1 (q 1 ))(2 K1 (2q 1 ) + K1 (q 1 )) 2 K1 (q 1 )( K1 (2q 1 ) + K1 (q 1 )) Note that Γ does not depend on ū1 .Indeed, since q 1 = 2π/L, it follows that Γ depends only on L. In Figure 3 we show the graphs of Γ versus L for γ > 0 in (a), and γ < 0 in (b).We also recall that we are analyzing the sign of Γ in supercritical regimes (Λ > 0), for this reason we plot the curve Γ(L) only in those intervals in which Λ > 0. The graph in Figure 3(a) shows that in the mutual avoidance case (γ > 0), small amplitude patterns exist and are unstable for 3 < L < 3.5, and that they become stable as L > 3.5.Figure 3(b) shows that in the mutual attraction scenario (γ < 0), Γ(L) is always negative and therefore small amplitude patterns are always unstable.These results are summarized in Figure 4.
In summary, our analysis shows that the nature of the transition and the stability of the bifurcation patterns depend mainly on L. These results can be read and reinterpreted in terms of the parameters of the original system (2), recalling that L = α/l, where α is the sensing radius and l is the length of the dimensional spatial domain.Therefore, the qualitative behaviour of the system under study strongly depends on the measure of the sensing radius compared on the length of the domain.28)) versus the domain length L. When the magnitude of γ is small, the homogeneous steady state is linearly stable.As the magnitude of γ increases, the system undergoes a bifurcation and the homogeneous steady state becomes unstable as γ crosses γ ± c .For γ > 0 (a), when L is small the system undergoes a subcritical bifurcation.As L increases, the bifurcation becomes supercritical, and the emerging patterns will be unstable.As L increases further, the system undergoes a supercritical bifurcation leading to the emergence of stable patterns.
For γ < 0 (b), when L is small the system undergoes a supercritical bifurcation generating unstable small amplitude patterns.As L increases, the bifurcation becomes subcritical

Numerical Simulations
In this Section, we perform a numerical investigation of system (16).To solve numerically System (16), we use the spectral method and numerical schemes presented in [18].By employing a continuation technique, we recover numerical bifurcation diagrams which are compared with the bifurcation diagrams obtained via the weakly nonlinear analysis.We show that our weakly nonlinear analysis provides accurate approximations of stable steady-state solutions in supercritical stable regimes, as long as we stay close to the bifurcation threshold.We also analyse those bifurcations that generate unstable small amplitude patterns.In these cases, we numerically detect the existence of stable large amplitude solutions, which are not predicted by the weakly nonlinear analysis, but which were predicted by an energy method in [18].First, we analyze the scenarios depicted in Figures 2(b) (γ > 0) and (d) (γ < 0), in which ū1 ̸ = ū2 .These figures show subcritical bifurcations for sufficiently small values of L, then a shift to a supercritical regime, as L increases, and again a subcritical regime, as L increases further.Recall that if ū1 ̸ = ū2 then supercritical bifurcations always give rise to stable small amplitude solutions.
Figure 5 shows bifurcation diagrams obtained by fixing ū1 = 0.1 and ū2 = 10 and by changing L, in the mutual avoidance regime (γ > 0).This case corresponds to the scenario shown in Figure 2(b) (center).Dashed and solid lines represent unstable and stable branches, respectively, computed analytically, while the dots are computed numerically.For L = 2.7, the weakly nonlinear analysis predicts a subcritical bifurcation, and the numerical simulations confirm this result.In fact, just beyond the instability threshold (γ > γ c ≈ 3.20), we find stable large amplitude solutions, which persist when we decrease the control parameter γ below the instability threshold (Figure 5(a)).For L = 5, the analysis predicts a supercritical bifurcation and, again, the numerical simulations confirm this result.In Figure 5(b) we see, indeed, a good matching between the analytical branch and the numerical solutions, as long as γ is sufficiently close to the bifurcation threshold γ c ≈ 1.32.Finally, for L = 15 the subcritical bifurcation predicted by our analysis is also detected numerically (see Figure 5(c)).Here, we observe bistability between the homogeneous steady state and non-homogeneous solutions below the instability threshold γ c ≈ 1.03.γ > 0, ū1 ̸ = ū2    (for Λ > 0 and Γ > 0) (see Figure 4).In particular, for γ > 0, system (16) undergoes subcritical bifurcations for 2 < L < 3, unstable supercritical bifurcations for 3 < L < 3.5, and stable supercritical bifurcations for L > 3.5 (see Figure 4(a)).
In Figure 7 we analyze System (16) with γ > 0 and ū1 = ū2 = 10, for L = 3.1 in (a), and L = 4 in (b).In Figure 7(a) (left) we show the spatio-temporal evolution of a numerical solution whose initial condition is a small perturbation of the weakly nonlinear solution with L = 3.1.We observe that the numerical solution moves away from the initial condition and evolves toward a large amplitude pattern.The initial condition and the final stationary state are reported in Figure 7(a) (center).Therefore, when the supercritical branch is unstable, the system supports large amplitude patterns.These solutions exist even below the bifurcation threshold, as shown by the bifurcation diagram in Figure 7(a) (right).These large amplitude solutions are not predicted by the weakly nonlinear analysis.However we conjecture that they might be obtained analytically by expanding the weakly nonlinear analysis to higher orders.
In Figure 7(b) (left) we show the spatio-temporal evolution of a numerical solution whose initial condition is a small perturbation of the weakly nonlinear solution with L = 4.In this case, the analysis predicts that the small amplitude pattern is stable.
In the numerical simulation we observe that the solution moves towards a small amplitude pattern, which is well approximated by the weakly nonlinear analysis.
This result confirms the stability predicted by our analysis.The initial condition and the final stationary state are reported in Figure 7(b) (center).Finally, a comparison between the analytical and numerical bifurcation diagrams is shown in 7(b) (right).
γ < 0, ū1 = ū2  and stable branches, respectively, which are computed analytically, while the dots are computed numerically.As the length of the domain increases, the system changes its qualitative behaviour.
In (a): L = 5 and the system exhibits a supercritical bifurcation at γ = γ c ≈ −13.2, giving rise to a branch of unstable small amplitude solutions.In (b), L = 10 and at γ = γ c ≈ −10.7 the system exhibits a subcritical bifurcation

Bistability between small amplitude patterns and strongly modulated solutions
The existence of non-constant solutions to system (16), far away from any bifurcation of the constant solution, was already detected and analyzed in [18] using an energy method.By minimising an energy functional associated with the system, nontrivial stationary solutions were revealed which, as L increases, tend to look increasingly like piecewise constant functions, when γ > 0, or spike solutions, when γ < 0. We call such solutions strongly modulated because they are given by the superposition of more than one unstable Fourier mode.In this section, we will combine numerical and analytic solutions inferred by both the weakly nonlinear analysis here and the results presented in [18] to construct more comprehensive bifurcation diagrams.
For this, we focus on the case γ > 0 and ū1 = ū2 .Here, the system exhibits supercritical bifurcations for large values of L (see Figure 2 (a)).Also, as shown in Figure 3 (a), these supercritical bifurcations generate stable small amplitude patterns.In [18] we showed that under the same conditions (that is L ≫ 1, γ > 0 and ū1 = ū2 ), the system supports strongly modulated patterns.Therefore we expect that for L sufficiently large, there exist parameter regions in which small amplitude patterns and strongly modulated solutions coexist and are stable.
We have verified this numerically and the results are shown in Figure 9.When L is not too large, the system admits small amplitude solutions that bifurcate supercritically from the homogeneous steady state and remains stable as the control parameter γ increases (see Figure 9 (a)).In this case, we do not find strongly modulated solutions.As L increases, the supercritical branch of patterns predicted by the weakly nonlinear analysis still exists and is stable as long as γ is sufficiently close to the bifurcation threshold (see Figure 9 (b)).However, a second branch appears higher up, representing the strongly modulated solutions predicted by [18].As L increases further, the branch of stable small amplitude solutions becomes smaller and smaller (Figure 9(c)), and the solutions transition to strongly modulated for values of γ closer to the bifurcation threshold.L becomes sufficiently large, the system support strongly modulated patterns which coexist with stable small amplitude patterns.

Discussion
We have analysed bifurcations for a nonlocal advection diffusion system with two interacting populations that either mutually avoid or mutually attract.First, we analysed the linear stability of the homogeneous steady state and recovered the instability thresholds.Beyond these thresholds, the homogeneous steady state becomes unstable and the system is expected to form spatially inhomogeneous patterns.To predict the evolution of the system in the unstable regime, we used weakly nonlinear analysis to recover the equations governing the amplitude of the pattern and approximations of the inhomogeneous solutions.We found that the amplitude equations consist of a Ginzburg-Landau equation coupled with an equation for the zero mode.
Indeed, we obtained a sequence of linear problems whose general solutions must be a linear combination of the critical mode and the zero mode.This follows from the fact that the system under study obeys a conservation law.An equivalent result was shown in [25], where similar amplitude equations were derived using symmetry and scaling arguments.By means of the amplitude equations, we recovered the condition that ensures the stability of the patterns bifurcating from the homogeneous steady state.
To obtain concrete numerical results, we analysed the case where the spatialaveraging kernel, K, is a top-hat distribution.By combining analysis of the amplitude equation with numerical solutions, we showed that the system exhibits a variety of different types of bifurcations and bistability regimes, strongly depending on the ratio l/α.In particular, we found stable small amplitude patterns bifurcating supercritically from the homogeneous steady state at the onset of the instability.We also found subcritical regimes generating unstable small amplitude patterns, which coexist with both the stable homogeneous solution and stable large amplitude patterns.
In this case, numerics revealed an hysteresis effect due to the bistability between two stationary states.Finally, we also found supercritical bifurcations generating unstable small amplitude patterns.Beyond the instability threshold, we numerically detected stable large amplitude patterns that persist even when decreasing the bifurcation parameter below the instability threshold, revealing again a hysteresis effect similar to that found in the subcritical regime.
By combining weakly nonlinear analysis, numerical simulations and the energy functional analysis from [18], we obtained a comprehensive bifurcation picture.We found parameter regions exhibiting bistability between small amplitude patterns and strongly modulated solutions, when l/α ≫ 1.The range of bistability becomes smaller and smaller as l/α increases, because the small amplitude patterns lose their stability for values of the control parameter increasingly closer to the bifurcation threshold (Figure 9).Overall, our analysis reveals that our system may display discontinuous phase transitions either when α ≈ l or when the sensing range α is very small compared to the length l of the domain.
Our study provides an example of how to combine different and complementary approaches to recover more comprehensive pictures of the bifurcation diagrams.To extend these results further, it would be interesting to expand the weakly nonlinear analysis up to higher orders.Such an approach could reveal analytically some of the large amplitude branches here found numerically, as well as the branches of solutions connecting small and large amplitude patterns.Numerical continuation software, such as pde2path [36], gives another way of approaching this problem [31,11].Our analysis revealed parameter regions with bistability between two extended states, a scenario in which systems often exhibit snaking branches of localized solutions [7,37].Extending our weakly nonlinear analysis to higher orders may help locate the codimension-two point where the nascence of localised structures may take place, which would be an interesting subject for future work.
Our focus here has been on a particular example of Equation (1) [28], with just two populations and no self-interaction terms (N = 2, γ ii = 0).However, even in this relatively-simple system, we found an unexpectedly rich variety of patterning scenarios.Therefore, we conjecture that analysis of the system with N ≥ 3 populations and/or γ ii ̸ = 0 would reveal even more complex patterning and bifurcation structure.
Our next goal, indeed, is to analyse the more general scenarios (N ≥ 3, γ ii ̸ = 0).A possible way forward might be to analyse phase transitions by combining the tools used here with those from [9].In [9] the authors studied the phase transitions of System (1) has several applications to natural systems and, in particular, to ecological systems.Therefore the analysis presented in this paper, as well as possible future extensions, might help to address some important ecological questions regarding the emergence of territories, as well as their sizes and stability [26].Indeed, variations in territory size and shape can strongly affect population structure and dynamics [1], therefore understanding the mechanisms and consequences of these changes is crucial for informing the design of efficient conservation strategies.Our results support the hypothesis that the formation of territorial patterns is not just a consequence of a heterogeneity in resources distribution, but that they can emerge as a consequence of animal behaviour and mutual interactions [1,16,26].Our analysis also predicts that a small sensing range relative to the length of the domain can facilitate a territory instability, in agreement with other theoretical studies suggesting that poor sensory information can promote the range size instability ( [32]).In summary, the analysis of the class of models (1) with the techniques here presented and discussed can help to resolve biological and ecological questions that may be inaccessible to experimental investigation.

Figure 5 :
Figure 5: Comparison between analytical and numerical bifurcation diagrams of system (16) with density-dependent advection strength γ > 0, and nonlocal kernel K = K 1 (see Equation (76)), ū1 = 0.1 and ū2 = 10, for different values of the length of the domain L.These scenarios correspond to Figure 2 (b) (center).Dashed and solid lines represent unstable and stable branches, respectively, which are computed analytically, while the dots are computed numerically.As the length of the domain increases, the system changes its qualitative behaviour.In (a): L = 2.7 and the system exhibits a subcritical bifurcation at γ = γ c = 3.19933.In (b), L = 5 and at γ = γ c = 1.32131 a branch of stable solutions bifurcates from the homogeneous state .In (c), L = 15 and the system exhibits a subcritical bifurcation at γ = γ c = 1.02985.

Figure 6
Figure 6 shows bifurcation diagrams obtained by fixing ū1 = 0.1 and ū2 = 10, for three different values of L, in the mutual avoidance regime (γ < 0).This case corresponds to the scenario shown in Figure 2(d) (center).The numerical simulations, again, confirm the results of the weakly nonlinear analysis: we have detected subcritical transitions for L = 2 and L = 10, and a stable branch bifurcating supercritically for L = 5, whose amplitude is well approximated by the weakly nonlinear analysis.

Figure 6 :
Figure 6: Comparison between analytical and numerical bifurcation diagrams of system (16) with γ < 0, K = K 1 (see Equation (76)), ū1 = 0.1 and ū2 = 10, for different values of the length of the domain L.These scenarios correspond to Figure 2 (d) (center).Dashed and solid lines represent unstable and stable branches, respectively, which are computed analytically, while the dots are computed numerically.As the length of the domain increases, the system changes its qualitative behaviour.In (a): L = 2.5 and the system exhibits a subcritical bifurcation at γ = γ c = −4.2758.In (b), L = 5 and at γ = γ c = −1.32131 a branch of stable solutions bifurcates from the homogeneous state.In (c), L = 10 and the system exhibits a subcritical bifurcation at γ = γ c = −1.06895.

(a) L = 3 . 1 (b) L = 4 Figure 7 :
Figure 7: Numerical investigation of system (16) in the mutual avoidance regime (γ > 0) with ū1 = ū2 = 10, for two different values of L. In (a): L = 3.1 and the analysis predict an unstable supercritical bifurcation at γ = γ c = 0.225754.On the left, numerical simulation showing that the system moves away from the unstable solution and evolves toward a large amplitude pattern.In the center, initial condition and the final stationary state.On the right, comparison between analytical and numerical bifurcation diagrams.In (b): L = 4 and the analysis predict a stable supercritical bifurcation at γ = γ c = 0.15708.On the left, numerical simulation showing that the system moves towards the stable small amplitude solution.In the center, initial condition and the final stationary state.On the right, comparison between analytical and numerical bifurcation diagrams.

Figure 8 :
Figure 8: Comparison between analytical and numerical bifurcation diagrams of system (16) with γ < 0, K = K 1 (see Equation (76)), ū1 = ū2 = 10, for different values of the length of the domain L.These scenarios correspond to Figure 2 (c) (right).Dashed and solid lines represent unstable

the
Mckean-Vlasov equation by analysing the minimizers of the energy associated to the problem.Combining this with weakly nonlinear analysis might shed light on the number of steady states at the onset of an instability, and consequently on type of phase transition occurring when the bifurcation parameter crosses the instability threshold.

Table 1 :
List and description of main parameters involved in the study of stability and bifurcations.