Global well-posedness, pattern formation and spiky stationary solutions in a Beddington–DeAngelis competition system

This paper investigates a reaction-advection-diffusion system that describes the evolution of population distributions of two competing species in an enclosed bounded habitat. Here the competition relationships are assumed to be of the Beddington–DeAngelis type. In particular, we consider a situation where first species disperses by a combination of random walk and directed movement along with the population distribution of the second species which disperse randomly within the habitat. We obtain a set of results regarding the qualitative properties of this advective competition system. First of all, we show that this system is globally well-posed and its solutions are classical and uniformly bounded in time. Then we study its steady states in a one-dimensional interval by examining the combined effects of diffusion and advection on the existence and stability of nonconstant positive steady states of the strongly coupled elliptic system. Our stability result of these nontrivial steady states provides a selection mechanism for stable wavemodes of the time-dependent system. Moreover, in the limit of diffusion rates, the steady states of this fully elliptic system can be approximated by nonconstant positive solutions of a shadow system that admits boundary spikes and layers. Furthermore, for the fully elliptic system, we construct solutions with a single boundary spike or an inverted boundary spike, i.e., the first species concentrates on a boundary point while the second species dominates the remaining habitat. These spatial structures model the spatial segregation phenomenon through interspecific competitions. Finally, we perform some numerical simulations to illustrate and support our theoretical findings.

1. Introduction. In this paper, we consider the following reaction-advectiondiffusion system x ∈ ∂Ω, t > 0, u(x, 0), v(x, 0) ≥, ≡ 0, x ∈ Ω, Here Ω is a bounded domain in R N , N ≥ 1 and ∂Ω is its piecewise boundary, with n denoting the unit outer normal on ∂Ω. D i , α i , β i , a i , b i and c i , i = 1, 2, are positive constants. χ is a constant and φ is a smooth function such that φ(v) > 0 for all v > 0. The initial conditions u 0 (x) and v 0 (x) are non-negative smooth functions which are assumed to be not identically zero. System (1.1) describes the evolution of population densities of two mutually interfering species over a bounded habitat Ω. u(x, t) and v(x, t) are population densities of the competing species at space-time location (x, t) ∈ Ω × R + , and diffusion rates D 1 and D 2 measure intensity of unbiased dispersal of the species. It is assumed that species u sense the population pressure from v and directs its dispersal accordingly. In particular, u moves up or down along the population gradient of interspecies v and χ is a constant that measures the intensity of such directed movement. One can choose χ < 0 if u invades the dwelling habitat of v and χ > 0 if u escapes the habitat of v. Therefore u takes active dispersal strategy to cope with population pressure from v, either to seek or avoid interspecific competition. φ reflects the variation of the directed movement strength to the population density of v.
For almost all mechanistic models that describe population dynamics, functional responses play key roles in the spatial-temporal and qualitative properties of the population distributions. One important feature of these competition dynamics is to take into account the environment's maximal load or carrying capacity. Among the most common types are Lotka-Volterra [45,7], Hollings [26,55], and Leslie-Gower [4], etc. To investigate mutual interference among intraspecies, Beddington [5] and DeAngelis et al. [14] introduced the following functional response of focal species u preying natural resources where R is the resource density, a is the attacking rate and h is the handling time, while b measures the interference rate. To account for both intra-and interspecific competition, de Villemereuil and Lopez-Sepulcre [46] generalized the B-D response into f 0 (R, u, v) = 1 1 + ahR + bu + cv , where b and c represent the intra-and interspecific competition rates, respectively. Moreover, they performed field experiments on the guppy-killifish system and collected data that fit the responses above in the absence and presence of interspecies. Functional response f 0 (R, u, v) models a competition relationship between u and v based on the idea that an increase in the population density of each species should suppress the growth rate of both individuals since they consume the same resources.
In this paper, we assume that species u and v interact with the Beddington-DeAngelis functional responses as in (1.2), where ecologically α 1 and α 2 represent intrinsic death rates of species u and v and β 1 and β 2 interpret their growth rates; a 1 and a 2 represent the compound effects of resource handling time and attack rate, while b 1 and c 2 account for intra-specific competition and b 2 and c 1 are the coefficients of interspecific competition. We are motivated to investigate the dynamics of population distributions to (1.1) and its stationary system due to the effect of the biased movement of species u. To manifest this effect, we assume that the resources are spatially homogeneous hence all the parameters are set to be positive constants.
To model the coexistence and segregation phenomenon through interspecific competition, various reaction-diffusion systems with advection or cross-diffusion have been proposed and studied. For example, [48,50,51] studied the global existence of (1.1) with Lotka-Volterra kinetics f (u, v) = a 1 − b 1 u − c 1 v u and g(u, v) = a 2 − b 2 u − c 2 v v. Existence and stability of nonconstant positive steady states have also been established through rigorous bifurcation analysis; moreover, they obtained transition layer steady states to the system. Another example of this type is the SKT model proposed by Shigesada, Kawasaki and Teramoto [41] in 1979 that studied the directed dispersals due to mutual interactions. See [30,31,36,38,41] etc. for works and recent developments on the SKT model. It is also worthwhile to mention that predator-prey models with Beddington-DeAngelis type functional response have been theoretically studied in [8,15,20] etc. and there are also several works on the prey-taxis models [22,47,49], and we refer the reader to [11] and the survey paper [10] for detailed discussions on reaction-advection-diffusion systems of population dynamics.
From the viewpoint of mathematical modeling, it is interesting and important to investigate the spatially inhomogeneous distribution of population densities such as the coexistence and segregation of mutually interfering species, in particular, due to the dispersal and species competition. Time-dependent solutions can be used to model the segregation phenomenon through the blow-up solutions in a finite or infinite time. That being said, a species population density converges to a δfunction or a linear combination of δ-functions. This approach has been taken for Keller-Segel chemotaxis system that models the directed cellular movements along the concentration gradient of stimulating chemicals in their environment [18,17,37]. An alternative approach is to show that the solutions exist globally and converge to bounded steady states. Then positive steady states with concentrating or aggregating structures such as boundary spikes, transition layers, etc. can be used to model the segregation phenomenon. This approach has been adopted by [31,38,52] etc. The blow-up solution or a δ-function is evidently connected to the species segregation phenomenon. However, it is not an optimal choice from the viewpoint of mathematical modeling since it challenges the rationality that population density can not be infinity. On the other hand, it brings challenges to numerical simulations and makes it impossible to analyze the solutions and models after the blow-up. Therefore, we are motivated to study the global well-posedness and stationary solutions of (1.1)-(1.2) in this paper.
The rest of this paper is organized as follows. In Section 2, we obtain the existence and uniform boundedness of positive classical solutions to (1.1) over multidimensional bounded domains-See Theorem 2.1. In Section 3, we study the existence and stability of nonconstant positive steady states of (1.1) over Ω = (0, L). Our method is based on the local theory of Crandall and Rabinowitz [12] and its new version recently developed by Shi and Wang in [40]. Our stability results give a selection mechanism for stable wavemodes of system (1.1)-See Proposition 2 and Theorem 3.2. Section 4 is devoted to the existence and asymptotic behaviors of nonconstant steady states with large amplitude. It is shown that (1.1)-(1.2) admits boundary spike and boundary layer solution if D 1 , χ are comparably large and D 2 is small-see Theorem 4.1. In Section 5, we present some numerical studies of the model and briefly discuss our results and their applications. Some interesting problems are also proposed for future studies. In the sequel, C and C i denote positive constants that may vary from line to line.

2.
Global existence and boundedness. In this section, we investigate the global existence of positive classical solutions to system (1.1)-(1.2). We will show that the L ∞ -norms of u and v are both uniformly bounded in time. The first set of our main results states as follows.
Theorem 2.1. Let Ω ⊂ R N , N ≥ 1, be a bounded domain with a piecewise smooth boundary ∂Ω. Assume that the function φ is smooth and positive for all v ≥ 0. Then for each pair (u 0 , v 0 ) ∈ C(Ω) × W 1,p , p > N , the IBVP (1.1)-(1.2) admits a unique and global solution (u(x, t), v(x, t)) for all t > 0; moreover, the solution is classical and uniformly bounded in time, and both u and v are positive onΩ for all t ∈ (0, ∞).
Many reaction-diffusion systems of population dynamics have maximum principles, which can be used to prove the global existence and boundedness of their classical solutions. However, the presence of advection term makes system (1.1) non-monotone and it inhibits the application of the comparison principle, at least in proving the boundedness of u. A review of the literature suggests that there are two well-established methods to prove the global existence for reaction-advectiondiffusion systems. One method is to derive the L ∞ through the L p -iterations of some compounded functions of the population densities; another method is to apply standard theory on semigroups generated by {e −∆t } t≥0 and L p -L q estimates on the abstract form of the system. We refer the reader to [16] for classical results and [19,53] for recent developments. Our proof of global existence of (1.1)-(1.2) involves both techniques.

Preliminary results and Local existence.
First of all, we collect the local existence of classical solutions to (1.1)-(1.2) and their extension criterion based on Amann's theory in the following theorem.
Let Ω be a bounded domain in R N , N ≥ 1 with smooth boundary ∂Ω and assume that φ is a continuous function. Then for any initial data (u 0 , v 0 ) ∈ C(Ω)×W 1,p (Ω), p > N , satisfying u 0 , v 0 ≥, ≡ 0 onΩ, there exists a constant T max ∈ (0, ∞] and a unique solution Proof. Let w = (u, v). (1.1) can be rewritten as (2.1) is a triangular normally parabolic system since the eigenvalues of A are positive, then the existence part follows from Theorem 7.3 and Theorem 9.3 of [2] and the extension criterion follows from Theorem 5.2 of [3]. Moreover, one can apply the standard parabolic boundary L p estimates and Schauder estimates to see that u t , v t and all the spatial partial derivatives of u and v are bounded inΩ × (0, ∞) up to the second order, hence (u, v) has the regularities as stated in the Theorem. Finally, one can use parabolic strong maximum principle and Hopf's boundary point lemma to show that u > 0 and v > 0 onΩ×(0, T max ). This completes the proof of Theorem 2.2.

2.2.
A priori estimates. We collect some basic properties and a priori estimates of the local classical solutions obtained in Theorem 2.2.
moreover, for any p ∈ (1, ∞), there exists a positive constant C(p) dependent on Proof. To show (2.2), it is sufficient to show that Ω u(x, t)dx is uniformly bounded for t ∈ (0, ∞) since u(x, t) > 0 according to Theorem 2.2. Integrating the first equation in (1.1) over Ω yields To prove (2.3), we first see that v L 1 (Ω) is bounded for t ∈ (0, T max ) by the same arguments as above. Using the v-equation in (1.1), we have from the integration by parts that therefore, v(·, t) L p (Ω) is bounded over (0, T max ) for any p ∈ (1, ∞). To establish the L ∞ -bound on u, we need to estimate ∇v L p for p large. For this purpose, we convert the v-equation in (1.1) into the following abstract form v(·, t) = e −D2At v 0 + t 0 e −D2A(t−s) D 2 v(·, s) + g(u(·, s), v(·, s)) ds, (2.5) where g(u, v) is given in (1.2). After applying the L p -estimates in Newmann semigroups (e.g., [16,19,53]) on (2.5), we have the following result. and Ω such that By taking p > N in (2.6), we readily have the following result. Proof of Theorem 2.1. Thanks to (2.7) in Lemma 2.5, we only need to show that sup t∈(0,Tmax) u(·, t) L ∞ < ∞, whereas T max = ∞ readily follows; moreover, the regularities of (u, v) follow from Theorem 2.1. Without loss of generality, we assume that φ(v(·, t))∇v(·, t) L ∞ (Ω) ≤ 1 for all t ∈ (0, T max ) in light of Lemma 2.5. Then for any p > 1 we test the first equation of (1.1) by pu p−1 and integrate it over Ω by parts to have where C 1 (p) = 1 q β1 α1b1p q p |Ω| with q = p p−1 and we have applied the Young's inequality ab ≤ a p + b q q(p ) q p , ∀a, b, > 0. We recall Corollary 1 in [9] due to the Gagliardo-Nirenberg inequality: for any > 0, there exists a positive constant C 0 dependent on N and Ω such that (2.10) Finally, by applying the standard Moser-Alikakos iteration [1] to (2.10), we can show that u(·, t) L ∞ (Ω) is uniformly bounded for t ∈ (0, ∞). This completes the proof of Theorem 2.1.
3. Small-amplitude steady states and their stabilities. In this section, we investigate nonconstant positive steady states of (1.1)-(1.2) over one-dimensional finite domain in the following form where denotes the derivative taken with respect to x. We have assumed in (3.1) that α i = β i = 1, i = 1, 2 without loss of generality. Indeed, through the rescaling x ∈ Ω, the one-dimensional steady state of which recovers (3.1).
We will see that large advection rate χ drives the emergence of nonconstant positive solutions to (3.1). There are four constant solutions to system (3.1): (0, 0), is the unique positive solution provided that We assume this condition throughout the rest of our paper. The first inequality in (3.4) implies that b 2 c 1 < b 1 c 2 , therefore the interspecific competition is weak compared to the intraspecific competition. We call this condition the weak competition case. Similarly, we call the latter condition in (3.4) the strong competition case. The same weak and strong competition cases are proposed in the studies of SKT competition models in [41].

3.1.
Global attractor in the purely diffusive system. We shall show that the emergence of nonconstant positive solutions to (3.1) is driven by large advection rate χ. To see this, we study the existence of nonconstant positive solutions to (3.2) with χ = 0, i.e., the following diffusive system x ∈ (0, L).
Our main result states as follows.
and it is the only positive steady state of (3.5) provided with either (i). b1 b2 > c1 c2 or and we rewrite (3.5) as Then the matrix . According to Theorem 3.1 of [30], (ū,v) is the only steady state of (3.5) if either (i) or (ii) holds. On the other hand, it follows from straightforward calculations that the linearized stability matrix corresponding to system (3.5) which has two negative eigenvalues if either (i) or (ii) occurs, then by the same analysis in Theorem 2.5 in [32] or [42], one can show that system (3.5) generates a strongly monotone semi-flow on is globally asymptotically stable. This completes the proof of Theorem 3.1.
Our results indicate that the global dynamics of the diffusion system (3.5) is dominated by the ODEs in the weak competition case and also in the strong case if one of the diffusion rates is large. However, if both D 1 and D 2 are small, (ū,v) is unstable in the strong competition case. We surmise that positive solutions with nontrivial patterns may arise when the system parameters and the domain geometry are properly balanced. Nonconstant positive solutions with spikes are investigated for the diffusion system (3.1) with Lotka-Volterra dynamics by various authors. See [23,33,34,35].
3.2. Advection-driven instability. Theorem 3.1 states that diffusion itself does not change the dynamics of the spatially homogeneous solution and no Turing's instability occurs for the diffusive system (3.5) in most cases. We proceed to investigate the effect of advection rate χ on the emergence of nonconstant positive solutions to (3.1). We shall show that this equilibrium loses its stability as the advection rate χ crosses a threshold value. To this end, we first study the linearized stability of the equilibrium (ū,v).
We have the following result on the linearized instability of (ū,v) to (3.1).
Proof. According to the standard linearized stability analysis, the stability of (ū,v) is determined by the eigenvalues of the following matrix In particular, (ū,v) is unstable if H k has an eigenvalue with positive real part for some k ∈ N + . It is easy to see that the characteristic polynomial of (3.7) takes the form p k (λ) has a positive root if and only if p k (0) = D k < 0, then (3.6) follows from simple calculations and the proof completes.
We have from Proposition 1 that (ū,v) loses its stability when χ surpasses χ 0 . In the weak competition case b1 b2 > 1−a1 1−a2 > c1 c2 , we see that χ 0 is always positive, hence (ū,v) remains locally stable for χ being small. In the strong competition case b1 b2 < 1−a1 1−a2 < c1 c2 , χ 0 < 0 if both D 1 and D 2 are sufficiently small. This corresponds to the fact that (ū,v) is unstable for χ = 0 in (3.1) in this case. By the same stability analysis above, we can show that the appearance of advection does not change the stability of the rest equilibrium points. It is also worthwhile to point out that Proposition 1 carries over to higher dimensions with ( kπ L ) 2 replaced by the k-eigenvalue of Neumann Laplacian.

Steady state bifurcation.
To establish the existence of nonconstant positive solutions to (3.1), we first use the bifurcation theory of Crandall-Rabinowitz [12] by taking χ as the bifurcation parameter. To this end, we rewrite (3.1) into the following abstract form By the same arguments that lead to (iv) of Lemma 5.1 in [9] or Lemma 2.3 of [52], one can show that For bifurcation to occur at (ū,v, χ), we need the implicit function theorem to fail on F at this point thus require the condition N D (u,v) F(ū,v, χ) = 0. Let (u, v) be a nontrivial solution in this null space, then it satisfies the following system   (3.9) Expanding u and v into the following series and substituting them into (3.9) yields k = 0 is ruled out in (3.9) thanks to (3.4). For k ∈ N + , (3.10) has nonzero solutions if and only if its coefficient matrix is singular, which implies that local bifurcation can occur only at Moreover, the null space is one-dimensional and is spanned by Having the candidates for bifurcation values χ k , we now show that local bifurcation does occur at (ū,v, χ k ) in the following theorem, which establishes nonconstant positive solutions to (3.1).
Assume that (3.4) holds, and for any two integers k = j ∈ N + , where (ū,v) is the positive equilibrium of (3.1) given in (3.3). Then for each k ∈ N + , there exists δ > 0 small such that (3.1) admits nonconstant solutions The solutions are continuous functions of s in the topology of X × X × R and have the following expansions for |s| being small, with (ū k ,v k ) and Q k defined in (3.12); moreover, all nontrivial solutions of (3.1) Proof. Thanks to the arguments above, our results follow from Theorem 1.7 of Crandall and Rabinowitz [12] once we prove the following transversality condition where (ū k ,v k ) is given in (3.12) and R(·) denotes the range of the operator. We argue by contradiction and assume that (3.15) fails. Therefore there exists a nontrivial pair (u, v) to the following problem Testing the first two equations in (3.16) by cos kπx L over (0, L) yields (3.17) The coefficient matrix is singular due to (3.11), then we reach a contradiction in (3.17) and this proves (3.15). On the other hand, we need condition (3.13) such that χ k = χ j for all integers k = j, which is also required for the application of the local theory in [12] for single eigenvalue. The proof of Proposition 2 is completed.  (s, x), v k (s, x), χ k (s)) obtained in Proposition 2. Here the stability or instability means that of the nonconstant solution viewed as an equilibrium of the time-dependent system of (3.1). F is C 4 -smooth in s if φ is C 4 , therefore, according to Theorem 1.18 in [12], (u k (s, x), v k (s, x), χ k (s)) is C 3 -smooth and we have the following expansions (3.18) where (ψ i , ϕ i ) ∈ Z as defined in (3.14) and K i is a constant for i = 1, 2.
Note that Taylor's expansion gives where the o(s 3 )-terms in (3.18)  Theorem 3.2. Assume that all the conditions in Proposition 2 hold. Suppose that χ k0 = min k∈N + χ k of (3.8). Then for all positive integers k = k 0 and small In Table 1, we list the wavemode numbers k 0 and the corresponding minimum bifurcation values for different interval lengths. Figure 1 gives numerical simulations on the evolutions of stable patterns over different intervals. In Figure 3, we plot the spatial-temporal solutions to illustrate the formation of a single boundary spike to u and a boundary layer to v. The boundary spike and layer here correspond to those obtained in Theorem 3.1. if χ is chosen to be slightly larger than χ k 0 . We can also see that larger domains support higher wave modes. Figure 1 is given to illustrate the wavemode selection mechanism. Remark 1. Theorem 3.2 provides a rigourous selection mechanism for stable nonconstant positive solutions to system (3.1). If (u k (s, x), v k (s, x)) is stable, then k must be the integer that minimizes the bifurcation value χ k over N + , that being said, if a bifurcation branch is stable, then it must be the left-most branch.  Table 1. χ is chosen to be slightly larger than χ0 and the rest system parameters are chosen to be the same as in Table 1. Initial data are small perturbations of (ū,v).
Multiplying (3.21) by cos kπx L and integrating it over (0, L) by parts give rise to where the eigenvalue λ satisfies

QI WANG, LING JIN AND ZENGYAN ZHANG
If χ k = χ k0 = min k∈N + χ k , thenD k < 0 thanks to (3.11) hencep k (λ) always has a positive root λ(0) for all k = k 0 . From the standard eigenvalue perturbation theory in [24], (3.20) always has a positive root λ(s) for small s if k = k 0 . This finishes the proof of the instability part.
Differentiating (3.23) with respect to χ and then taking χ = χ k0 , we have that The coefficient matrix is singular thanks to (3.11), then we must have that which, in light of (3.12), implies thaṫ According to Theorem 1.16 in [13], the functions λ(s) and −sχ k0 (s)μ(χ k0 ) have the same zeros and signs near s = 0, and lim s→0 −sχ k0 (s)μ(χ k0 ) λ(s) = 1, forλ(s) = 0, then we conclude that sgn(λ(0)) = sgn(K 2 ) in light of K 1 = 0. Following the standard perturbation theory, one can show that sgn(λ(s)) = sgn(K 2 ) for s ∈ (−δ, δ), δ being small. This finishes the proof of Theorem 3.2. According to Theorem 3.2, if k 0 is the positive integer that minimizes χ k over N + , the bifurcation branch Γ k0 (s) around (ū,v, χ k0 ) is stable if it turns to the left and it is unstable if it turns to the right. However, for all k = k 0 , Γ k (s) is always unstable for s being small. It is unknown about the global behavior of the continuum of Γ k (s). We discuss this in details in Section 5.
Proposition 1 indicates that (ū,v) loses its stability as χ surpasses χ k0 . Theorem 3.2 shows that the stability is lost to the nonconstant steady state (u k0 (s, x), v k0 (s, x)) in the form (Q k0 , 1) cos k0πx L and we call this the stable wavemode of (1.1). If the initial value (u 0 , v 0 ) are small perturbations from (ū,v), then spatially inhomogeneous patterns can emerge through this mode, at least when χ is around χ k0 , therefore stable patterns with interesting structures, such as interior spikes, transition layers, etc. can develop as time evolves. Rigorous analysis of the boundary spike solutions will be analyzed in the coming section.
If the interval length L is small, we see that therefore χ 1 = χ k0 = min k∈N + χ k and Theorem 3.2 shows that the only stable pattern is (u 1 (s, x), v 1 (s, x)) which is spatially monotone. It is easy to see that k 0 increases if L increases. This indicates that small domains only support monotone stable solutions, while large domains support non-monotone stable solutions. Indeed, one can construct non-monotone solutions to (3.1) by reflecting and periodically extending the monotone ones at the boundary points 0, ±L, ±2L, ...

4.
Boundary spikes of the 1D systems. This section is devoted to positive solutions of (3.1) with large amplitude compared to the small amplitude bifurcating solutions. In particular, we study the boundary spike solution as D 1 and χ approach infinity with χ/D 1 ∈ (0, ∞) being fixed. Since the smallness of D 2 gives rise to boundary spikes for system (3.1) and our approach works for a wide class of the sensitivity functions, we put = √ D 2 and choose φ(v) ≡ 1 in (3.1) without loss of generality. Therefore, it is the goal of this section to investigate positive solutions with large amplitude of the following system   Figure 2 illustrates the formation of such spatial pattern.

4.1.
Convergence to a shadow system. Theorem 4.1 is an immediate consequence of several results and we first study (4.1) in the large limit of D 1 . To this end, we need the following a priori estimates.
Proof. We have from the maximum principles that therefore v is bounded in C 2 ([0, L]) thanks to the standard elliptic Schauder estimates.
On the other hand, we integrate the u-equation over (0, L) to see that Testing the u-equation by u and integrating it over (0, L) give rise to
We now study the asymptotic behaviors of positive solutions (u, v) of (4.1) by passing the advection rate D 1 and advection rate χ to infinity. We assume that χ D1 remains bounded in this process, hence both u and v are bounded as in Lemma 4.2.
uniformly for all χ i and D 1,i > 0, by the compact embedding one finds that, v i converges to some v ∞ in C 1 ([0, L]) as i → ∞, after passing to a subsequence if necessary. On the other hand, we integrate the u-equation in (4.1) over (0, x) to have Denoting w i = u i e rivi , we have that Sending i to ∞ in (4.3), we conclude that w i → 0 uniformly on [0, L], therefore w i = u i e rvi converges to a nonnegative constant λ . Moreover we can use standard elliptic regularity theory to show that v is C ∞ -smooth and it satisfies the shadow system (4.2).
Proposition 3 implies that when both D 1 and χ are sufficiently large, the steady state (u(x), v(x)) can be approximated by the structures of the shadow system (4.2). We want to remark that the arguments in this Proposition carry over to the case of higher-dimensional bounded domains. To find boundary spikes of (1.1)-(1.2), we study the asymptotic behavior of the shadow system with small diffusion rate and then investigate the original system with large diffusion and advection rates. This approach originates from the idea in [25], applied by [43,44] in reaction-diffusion systems and developed by [30,31] for reaction-diffusion system with cross-diffusions.

4.2.
Boundary spike layers of the shadow system. We proceed to construct boundary spikes of the shadow system (4.2). Our results state as follows.
To prove Proposition 4, we first choose λ > 0 to be a predetermined fixed but arbitrary constant and establish nonconstant positive solutions v (λ; x) to the following problem Then we find λ = λ > 0 such that (λ , v (λ ; x)) satisfies the integral constraint

then (4.4) becomes
It is equivalent to study the existence and asymptotic behaviors of (4.6) in order to prove Proposition 4. We first collect some facts about f (λ; s) introduced in (4.7) . We denote Lemma 4.3. Let r ∈ (0, ∞). Then for each λ ∈ (0, c2 b2r e rv * − r c 2 −1 ), g(λ; s) has two positive roots s 1 < s 2 and 1 + g(λ; s) has two positive roots s 3 < s 4 such that s 1 < s 3 < s 4 < s 2 ; moreover as λ → 0, (s 1 , s 3 ) → (0, 1 c2 ) and s 2 , s 4 → ∞. Proof. It is easy to see that the roots of g(s) or 1+g(s) must be positive if they exist; moreover, if 1 + g(s) has positive roots, so does g(s). Hence we only need to show that 1 + g(s) has positive roots and it is equivalent to show that its minimum over R + is negative. Let s * be the critical point of 1+g(s), we have from straightforward calculations that therefore 1 + g(s) is convex and its minimum value is 1 + g(s * ) = 1 + c2 r − c2 r ln( c2e rv * b2λr ) < 0 as desired. Since s i is continuous in λ, putting λ = 0 in 1 + g(s) gives s 1 = 0 and s 3 = 1 c2 . Moreover, s 2 > s 4 > s * → ∞ as λ → 0. This finishes the proof.
In light of the unique solution W 0 to (4.11), we now construct a boundary spike to (4.9) hence a boundary layer to (4.4). To this end, we choose a smooth cut-off function ρ(x) such that ρ(x) ≡ 1 for |x| ≤ L we want to show that (4.10) has a solution in the formw (λ, x) = W ,λ (x)+ ψ(λ, x). Then ψ satisfies L ψ + P + Q = 0, (4.12) where L = 2 d 2 dx 2 +f w (λ; W ,λ (x)), (4.13) and According to (4.12)-(4.15), P and Q measure how W ,λ (x) approximates solution w (λ, x). Our existence result is a consequence of several lemmas. Set We want to point out that Lemma 4.6 generalizes Lemma 5.3 in [48] which holds for L p , p > 1. Assuming Lemmas 4.6-4.8, we prove the following results of positive solutions to (4.10).
Proposition 5. Let r ∈ (0, ∞) and a 2 ∈ (0, 1). Suppose that (4.8) is satisfied. There exists a small 4 > 0 such that for all ∈ (0, 4 ), hence S is a contraction mapping on B(R 0 ) for small positive . We conclude from the Banach fixed-point theorem that S has a fixed point ψ in B, which is a smooth solution of (4.10). It is easy to show thatw satisfies (4.16) and this finishes the proof of Proposition 5.
3 ), since W 0 decays exponentially at ∞, we can also see that P is bounded. Therefore sup x∈[0,L] |P | is bounded.
where we have used the fact that |f ww (Ψ i )|, |f www (Ψ i )| are bounded for Ψ i ∈ B(R 0 ).

4.3.
Boundary spike and boundary layer to the full system. Now we prove Theorem 4.1 by showing that the solutions to (4.1) perturb from its shadow system (4.4) when D 1 is sufficiently large. First of all, we let s = 1 Define the projection operator P :
It is an interesting and also important mathematical question to study the stability of the single boundary spike solutions to the shadow system and the full system. To this end, one needs more detailed asymptotic expansions on the perturbed solutions than what obtained here, however this has be to postponed to future studies. We refer to [21,29,27,28,54] and the references therein for stability analysis of spiky solutions of some closely related reaction-diffusion systems and/or the shadow systems.

5.
Conclusions and discussions. In this paper, we investigate the interspecific competition through a reaction-advection-diffusion system with Beddington-DeAngelis response functionals. Species segregation phenomenon is modeled by the global existence of time-dependent solutions and the formation of stable boundary spike/layers over one-dimensional intervals.
There are several findings in our theoretical results. First of all, we prove that system (1.1)-(1.2) admits global-in-time classical positive solutions which are uniformly bounded globally for any space dimensions. For Ω = (0, L), we show that large directed dispersal rate χ of the species destabilizes the homogeneous solution (ū,v) and nonconstant positive steady states bifurcate from the homogeneous solution. Moreover, we study the linearized stability of the bifurcating solutions which shows that the only stable bifurcating solutions must have a wavemode number k 0 , which is a positive integer that minimizes bifurcating value χ k in (3.11)-Proposition 2 and Theorem 3.2. That being said, the constant solution loses its stability to the wavemode cos k0πx L over (0, L). Our result indicates that if the domain size L is sufficiently small, then only the first wavemode cos πx L can be stable, which is spatially monotone; moreover, large domain supports stable patterns with more aggregates than small domains. We prove that the steady states of (1.1)-(1.2) over (0, L) admits boundary spike in u and boundary layer in v if D 1 and χ are comparably large and D 2 is small. Therefore, we can construct multi-spike or multi-layer solutions to the system by reflecting and periodically extending the single spike over ±L, ±2L,... See Figure 4. These nontrivial structures can be used to the segregation phenomenon through interspecific competition. Compared to the transition layer in [48], our results suggests the boundary spike solutions to (1.1) is due to the Beddington-DeAngelis dynamics since a similar approach was adapted in [48].
We propose some problems for future studies. System (1.1) describes a situation that u directs its disperse strategy over the habitat to deal with the population pressures from v, which moves randomly over the domain. From the viewpoint of mathematical modeling, it is interesting to study the situation when both species takes active dispersals over the habitat. Then system (1.1) becomes a doubleadvection system with nonlinear kinetics. Analysis of the new system becomes much more complicated since no maximum principle is available for the second equation. We also want to point out that it is also interesting to study the effect of environment heterogeneity on the spatial-temporal behaviors of the solutions.
Our bifurcation analysis is carried out around bifurcation point (ū,v, χ k ). It is interesting to investigate the global behavior of the continuum of Γ k (s), denoted by C. According to the global bifurcation theory in [39] and its developed version in [40], there are the following alternatives for C: (i) it is either not compact in X × X × R + , or it contains a point (ū,v, χ j ), j = k. Actually, by the standard  elliptic embeddings, we can show that C is bounded in the axis of X × X if χ is finite. Therefore, it either extends to infinity in the χ-axis or it stops at (ū,v, χ j ). It is unclear to us which is the case for χ k . Considering Keller-Segel chemotaxis models over (0, L) without cellular growth, some topology arguments are applied in [9,52] etc. to show that all solutions on the first branch must be either monotone increasing or deceasing. Therefore χ 1 can not intersect at (ū,v, χ j ) for j = 1 since all the solutions around it must be non-monotone. This shows that Γ 1 extends to infinity in the χ direction. However, the appearance of complex structures of the kinetics in our model inhabits this methodology.
Our numerical simulations suggest that these boundary spikes and layers are globally stable for a wide range of system parameters. Rigorous stability analysis of these spikes is needed to verify our theoretical results. It is also interesting to investigate the large-time behavior of (1.1)-(1.2), which requires totally different approaches from what we apply in our paper.