Electronic Journal of Differential Equations, Vol. 2022 (2022), No. 67, pp. 1–19. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu STATIONARY AND OSCILLATORY DYNAMICS OF NICHOLSON’S BLOWFLIES EQUATION WITH ALLEE EFFECT CUIPING GUO, SHANGJIANG GUO Abstract. In this article we analyze the bistable dynamics of a Nicholson’s blowflies equation with Allee effect. Using Lyapunov-LaSalle invariance prin- ciple, we study the stability and basins of attraction of multiple equilibria. Also we study the existence, stability, and multiplicity of nontrivial steady- state solution and periodic solutions. These solutions generate long transient oscillatory patterns and asymptotic stable oscillatory patterns. 1. Introduction A biological species is said to exhibit an Allee effect in a habitat patch if its per capita rate of growth is decreasing at low densities. If the growth rate in the patch is negative when the density is below some threshold value, the species is said to exhibit a strong Allee effect. If the per capita growth rate remains positive at low densities, the species is said to exhibit a weak Allee effect. In a system with Allee effect, there usually exists a threshold under which the species will go to extinction. Namely, the species could not be expected to establish itself in the patch if introduced into the patch at a low enough density. In other words, it could not invade the patch. Moreover, it is typical that in such a system there are multiple stable equilibrium points. The study of Allee effect on the dy- namic behaviors of mathematical models provides enhanced insight into how Allee effects arise in the model outcomes, and is helpful for preventing the extinction of endangered species. In this article we consider the dynamical behavior of the delayed diffusive system ∂ ∂t u(t, x) = d∆u(t, x)− δu(t, x) + puγ(t− τ, x)e−au(t−τ,x), x ∈ Ω, ∂ ∂n u(t, x) = 0, x ∈ ∂Ω, (1.1) for t ≥ 0, where d, p, τ , a, and δ are positive constants, Ω is a connected bounded open domain in RN (N ≥ 1) with a smooth boundary ∂Ω. In biology, u(t, x) represents the population density of a species at time t and location x, δ is the per capita adult death rate, d describes the random movement of individuals, τ is a delay which usually represents the generation time. 2020 Mathematics Subject Classification. 34K18, 35B32, 35B35, 35K57, 35Q92, 92D40. Key words and phrases. Bistability; delay effect; Hopf bifurcation; stability. ©2022. This work is licensed under a CC BY 4.0 license. Submitted April 6, 2022. Published September 28, 2022. 1 2 C. GUO, S. GUO EJDE-2022/67 When γ = 1, it becomes the classical Nicholson’s blowflies equation which has been extensively studied in the literature [5, 12, 13]. For example, Guo [5] investi- gated the existence, stability, and multiplicity of nontrivial (spatially homogeneous or nonhomogeneous) steady-state solution and periodic solutions for a reaction- diffusion model with nonlocal delay effect and Dirichlet/Neumann boundary con- dition and then illustrate the general results by applications to the classical Nichol- son’s blowflies equation with one-dimensional spatial domain. Yi and Zou [12] studied the stability and existence of Hopf bifurcation of a delayed diffusive Nichol- son’s blowflies equation with Neumann boundary condition by applying maximun principal and some subtle inequalities. When γ > 1, the population growth rate is negative or decreasing function at low population size or density, which was found by Allee [1] and is termed as the strong Allee effect [2]. For the model (1.1) with d = 0 and Allee effect γ > 1, Terry [9] considered a special case of γ = 2 and found conditions of population extinction and persistence. For the case of 0 < γ < 1, Buedo-Fernández and Liz [3] established sharp global stability conditions for the positive equilibrium of equation (1.1) without diffusion. However, there are very few results for system (1.1) with general γ. In system (1.1), with the change of variables û = au, t̂ = δt, τ̂ = δτ, d̂ = d δ , p̂ = p δaγ−1 , and removing the hat, we obtain the equation ∂ ∂t u(x, t) = d∆u(t, x)− u(t, x) + puγ(t− τ, x)e−u(t−τ,x), x ∈ Ω, t ≥ 0, ∂ ∂n u(x, t) = 0, x ∈ ∂Ω t ≥ 0. (1.2) In this article, we analyze the dynamics of model (1.2) with γ > 1, including the existence, stability, Hopf bifurcation of steady states. We shall see that the change of p results in the presence of spatially homogeneous/nonhomogeneous steady-state solutions, while the presence of the time delay τ may cause some nonlinear oscilla- tions and hence can be regarded as a source of instability and oscillatory response of system. In a differential equation, either a steady state or a periodic orbit can be ob- tained by determining the zeros of an appropriate map and applying the Lyapunov- Schmidt procedure (for more details see [6]). In this article, we employ the Lyapunov- Schmidt reduction and the implicit function theorem to investigate the existence, uniqueness, bifurcation of spatially nonhomogeneous steady-state solutions (see Theorems 5.1 and 5.2). We shall see that the spatially nonhomogeneous steady- state solutions obtained by this means have an explicit algebraic form, which is helpful in the investigation of their stability and Hopf bifurcation. We set p∗ = eγ−1(γ − 1)1−γ , p∗ = eγγ1−γ , p∗n = (γ + 1 + σn)1−γ exp{γ + 1 + σn} for n ∈ N0 := N ∪ {0}. We shall see that the extinction steady state u∗0 is always locally stable, is globally asymptotically stable when p is small (that is, p < p∗), and attracts all small initial values regardless of τ when p > p∗. Therefore, the species goes to extinction if p is small (that is, p < p∗), and it persists if p is large and the initial population is appropriate. EJDE-2022/67 NICHOLSON’S BLOWFLIES EQUATION WITH ALLEE EFFECT 3 In addition to the extinction steady state u∗0 = 0, there exist two positive steady states u∗2(p) > u∗1(p) > 0 when the growth rate parameter p > p∗ and there is no positive steady state when p < p∗. The intermediate steady state u∗1(p) is always unstable; and the large steady state u∗2(p) is locally stable for p∗ < p < p∗0 and all τ > 0, but it becomes unstable for p > p∗0 and large τ . Moreover, u∗2(p) attracts all large initial values regardless of τ for a more restricted range of p. Therefore, our results confirm that the dynamics of system (1.2) is bistable due to the Allee effect structure. When p is large, time-delay-inducing oscillatory patterns occur near both positive equilibria u∗1(p) and u∗2(p), generating unstable and stable periodic solutions around respective steady states. For large τ , unstable periodic solutions around u∗1(p) cause long transient oscillatory patterns before solutions eventually converge to one of asymptotic stable states: the extinction state or the persistence state, which makes a bistable structure. The persistence state can be either the equilibrium u∗2(p) or a stable limit cycle around u∗2(p). Namely, for large p > 0, a large initial density always keeps the population persist in an oscillatory fashion. The remaining parts of this article are organized as follows. Some preliminaries are present in Section 2. Section 3 is devoted to the local/global stability and Hopf bifurcation of constant steady states. In Section 4 we first obtain the a priori bounds of nonnegative steady state solutions, which also identifies the regions of parameters of nonexistence of positive non-constant steady-state solutions, and then we use Leray-Schauder degree theory to investigate the existence of non-constant steady-state solutions. In section 5, we employ the Lyapunov-Schmidt reduction to obtain the existence, multiplicity, and concrete structural forms of spatially nonhomogeneous steady-state solution. Section 6 is devoted to the Hopf bifurcation at the nontrivial constant steady-state solutions. Finally, conclude our results and discuss some future work in Section 7. 2. Preliminaries Denote by Hk(Ω) (k ≥ 0) the Sobolev space of the L2-functions f(x) defined on Ω whose derivatives dn dxn f (n = 1, . . . , k) also belong to L2(Ω). Denote the spaces X = H2(Ω) ∩H1 0 (Ω) and Y = L2(Ω), where H1 0 (Ω) = {u ∈ H1(Ω)| ∂∂nu(x) = 0 for all x ∈ ∂Ω}. For any subspace Z of X or Y, we also define the complexification of Z to be ZC := Z ⊕ iZ = {x1 + ix2|x1, x2 ∈ Z}. For the complex-valued Hilbert space YC, we use the standard inner product 〈u, v〉 = ∫ Ω ū(x)v(x)dx. For convenience, let f(u) = puγe−u. (2.1) Note that f ′(u) = puγ−1e−u(γ − u), f ′′(u) = puγ−2e−u[γ(γ − 1)− 2γu+ u2], then we see that f ′(u) > 0 for all u ∈ (0, γ) and f ′(u) < 0 for all u ∈ (γ,+∞), f(u) ≤ fmax = f(γ) = pγγe−γ for all u > 0, and that f ′′(u) has two zero points ξ− and ξ+ satisfying 0 < ξ− < γ < ξ+ (see Figure 1). In fact, ξ± = γ ±√γ. The initial value of system (1.2) is φ ∈ C, where C = {φ ∈ C([−τ, 0],X) : φ(θ, x) ≥ 0 for all (θ, x) ∈ [−τ, 0]× Ω}. First, we have the following positivity and boundedness of solutions to system (1.2). 4 C. GUO, S. GUO EJDE-2022/67 - 6 u f(u) 0 u∗ 1(p) u∗ 2(p) u∗ 3(p) Figure 1. The graph of function f(x) = puγe−u with p > p∗ and γ > 1. Lemma 2.1. The solution of system (1.2) with initial value φ ∈ C is positive for all t > 0 and x ∈ Ω and is ultimately uniformly bounded. Proof. It follows from system (1.2) that ∂ ∂t u(x, t) ≥ d∆u(t, x)− u(t, x) for all t ∈ [0, τ ] and x ∈ Ω. By the strong maximum principle, we see that u(x, t) > 0 for all t ∈ [0, τ ] and x ∈ Ω. Repeating the above step for t ∈ [jτ, (j + 1)τ ] with j ∈ N, we obtain that u(t, x) is positive for all t > 0 and x ∈ Ω. In addition, note that ∂ ∂t u(x, t) ≤ d∆u(t, x)− u(t, x) + fmax for all t > 0 and x ∈ Ω. Then the strong maximum principle implies that lim sup t→+∞ u(t, x) ≤ fmax for all x ∈ Ω. Namely, u(t, x) is ultimately uniformly bounded. � Assume u is a constant steady state of system (1.2). Obviously, u = 0 is always a constant steady state. If u > 0, then it satisfies p = u1−γeu. Let g(u) = u1−γeu, then we have g′(u) = u−γeu(u− γ + 1). Obviously, limu→0+ g(u) = +∞, g(+∞) = +∞, g′(u) < 0 for u ∈ (0, γ − 1), and g′(u) > 0 for u ∈ (γ − 1,+∞). At u = γ − 1, g(u) attains its minimum value gmin = eγ−1(γ − 1)1−γ . Thus, the solution to p = u1−γeu undergoes a saddle-node bifurcation at p = gmin. Namely, we have the following results on the existence and multiplicity of nonnegative constant steady states of (1.2) (see Figure 2). Lemma 2.2. For system (1.2) with γ > 1, u∗0 = 0 is a steady state for all p > 0. Moreover, (i) If p < p∗ := eγ−1(γ − 1)1−γ , there is no positive constant steady state. (ii) If p = p∗, there exists a unique positive constant steady state u∗ = γ − 1. EJDE-2022/67 NICHOLSON’S BLOWFLIES EQUATION WITH ALLEE EFFECT 5 - 6 u g(u) 0 u∗ 1(p) γ − 1 u∗ 2(p) p∗ p Figure 2. The graph of function g(u) = u1−γeu with γ > 1. (iii) If p > p∗, there are exactly two distinctive positive constant steady states u∗1(p) and u∗2(p) satisfying 0 < u∗1(p) < γ − 1 < u∗2(p), f ′(u∗1(p)) > 1 > f ′(u∗2(p)), d dp u∗1(p) < 0 < d dp u∗2(p). (iv) If p∗ < p < p∗ := eγγ1−γ then 0 < u∗1(p) < γ − 1 < u∗2(p) < γ. (v) If p > p∗ then 0 < u∗1(p) < γ − 1 < γ < u∗2(p). The following lemma shows that when the initial function is less than the smaller positive constant steady state u∗1(p), the solution will converge to u∗0, which implies that the population will become extinct when the initial density of the population is small. This phenomenon is an important manifestation of Allee effect. Lemma 2.3. (i) Assume that p > p∗ and the initial function satisfies 0 < φ(θ, x) < u∗1(p) for all θ ∈ [−τ, 0] and x ∈ Ω. Then the solution of system (1.2) satisfies 0 < u(t, x) < u∗1(p) for all t ≥ 0 and x ∈ Ω. (ii) Assume that p > p∗ and the initial function satisfies u∗1(p) < φ(θ, x) < u∗3(p) for all θ ∈ [−τ, 0] and x ∈ Ω, where u∗3(p) > u∗2(p) satisfies f(u∗3(p)) = u∗1(p). Then the solution of system (1.2) satisfies u∗1(p) < u(t, x) < u∗3(p) for all t ≥ 0 and x ∈ Ω. (iii) Assume that p∗ < p < p∗ and the initial function satisfies u∗1(p) < φ(θ, x) < γ for all θ ∈ [−τ, 0] and x ∈ Ω. Then the solution of system (1.2) satisfies u∗1(p) < u(t, x) < γ for all t ≥ 0 and x ∈ Ω. Proof. (i) The conclusion that u(t, x) > 0 for all t > 0 follows from Lemma 2.1. It follows from 0 ≤ φ(θ, x) < u∗1(p) that f(φ(θ, x)) < f(u∗1(p)) = u∗1(p) for all θ ∈ [0, τ ] and x ∈ Ω. In view of system (1.2), we have ∂ ∂t u(x, t) < d∆u(t, x)− u(t, x) + u∗1(p) 6 C. GUO, S. GUO EJDE-2022/67 for all t ∈ [0, τ ] and x ∈ Ω. By the strong maximum principle, we see that u(x, t) < u∗1(p) for all t ∈ [0, τ ] and x ∈ Ω. Repeating the above step for t ∈ [jτ, (j + 1)τ ] with j ∈ N, we obtain that u(t, x) < u∗1(p) for all t > 0 and x ∈ Ω. (ii) It follows from u∗1(p) < φ(θ, x) < u∗3(p) that f(φ(θ, x)) > f(u∗1(p)) = u∗1(p) for all θ ∈ [0, τ ] and x ∈ Ω. In view of system (1.2), we have ∂ ∂t u(x, t) > d∆u(t, x)− u(t, x) + u∗1(p) for all t ∈ [0, τ ] and x ∈ Ω. By the strong maximum principle, we see that u∗1(p) < u(x, t) < fmax < u∗3(p) for all t ∈ [0, τ ] and x ∈ Ω. Repeating the above step for t ∈ [jτ, (j + 1)τ ] with j ∈ N, we obtain that u∗1(p) < u(t, x) < u∗3(p) for all t > 0 and x ∈ Ω. (iii) Note that p < p∗ if and only if fmax < γ, then using a similar argument as above, we can obtain the conclusion (iii). The proof is complete. � Lemma 2.3(iii) implies that C+ = {ϕ ∈ C|u∗1(p) < φ(θ, x) < γ for all (θ, x) ∈ [−τ, 0]× Ω} is a positively invariant set of system (1.2) with p∗ < p < p∗. In the following theorem, we see that u∗2 attracts all initial conditions in C+. Theorem 2.4. Assume that p∗ < p < p∗, then for every τ > 0 and initial value φ ∈ C+, the solution u(t, x) to (1.2) converges to u∗2(p) as t→∞. Proof. It is easy to see that function h(u) = u − 1 − lnu is strictly decreasing on (0, 1), is strictly increasing on (1,+∞), and has a global minimum 0 at u = 1. We define a Lyapunov functional V (u)(t) = ∫ Ω [ u∗2(p)h (u(t, x) u∗2(p) ) + u∗2(p) ∫ t t−τ h (f(u(s, x)) u∗2(p) ) ds ] dx. Obviously, V (u∗2(p)) = 0, and V is positive definite with respect to u ∈ (u∗1(p), γ). The derivative of V along the solutions of system (1.2) is d dt ∣∣ (1.2) V (u)(t) = ∫ Ω [( 1− u∗2(p) u(t, x) ) [−u(t, x) + f(u(t− τ, x))] ] dx + u∗2(p) ∫ Ω [ h (f(u(t, x)) u∗2(p) ) − h (f(u(t− τ, x)) u∗2(p) )] dx− du∗2(p) ∫ Ω ∆u(t, x) u(t, x) dx = ∫ Ω [( 1− u∗2(p) u(t, x) ) [u∗2(p)− u(t, x)] ] dx + ∫ Ω [( 1− u∗2(p) u(t, x) ) [f(u(t− τ, x))− u∗2(p)] ] dx + u∗2(p) ∫ Ω [ h (f(u(t, x)) u∗2(p) ) − h (f(u(t− τ, x)) u∗2(p) )] dx− du∗2(p) ∫ Ω |∇u(t, x)|2 u2(t, x) dx ≤ − ∫ Ω ( 1− u∗2(p) u(t, x) )2 u(t, x)dx+ ∫ Ω [f(u(t− τ, x))− u∗2(p)]dx − ∫ Ω u∗2(p) u(t, x) [f(u(t− τ, x))− u∗2(p)]dx EJDE-2022/67 NICHOLSON’S BLOWFLIES EQUATION WITH ALLEE EFFECT 7 + ∫ Ω [ u∗2(p)h (f(u(t, x)) u∗2(p) ) − f(u(t− τ, x)) + u∗2(p) + u∗2(p) ln f(u(t− τ, x)) u∗2(p) ] dx = − ∫ Ω ( 1− u∗2(p) u(t, x) )2 u(t, x)dx− ∫ Ω u∗2(p) u(t, x) [f(u(t− τ, x))− u∗2(p)]dx + u∗2(p) ∫ Ω [ h (f(u(t, x)) u∗2(p) ) + ln f(u(t− τ, x)) u∗2(p) ] dx = u∗2(p) ∫ Ω [ h (f(u(t, x)) u∗2(p) ) − h (f(u(t− τ, x)) u(t, x) ) − h (u(t, x) u∗2(p) )] dx ≤ u∗2(p) ∫ Ω [ h (f(u(t, x)) u∗2(p) ) − h (u(t, x) u∗2(p) )] dx. Note that f(u∗2(p)) = u∗2(p) and f is monotonic increasing in (0, γ). If u∗1(p) < u < u∗2(p) then u∗1(p) < u < f(u) < u∗2(p) and hence h ( f(u) u∗2(p) ) < h ( u u∗2(p) ) . If u∗2(p) < u < γ then u∗2(p) < f(u) < u < γ and hence h ( f(u) u∗2(p) ) < h ( u u∗2(p) ) . Thus, we have d dt ∣∣ (1.2) V (u)(t) ≤ 0 along an orbit u(t, x) of system (1.2) with any initial value φ ∈ C+. This implies that the solution u(t, x) to (1.2) converges to u∗2(p) as t→∞. � 3. Stability of constant steady states Let 0 = λ0 < λ1 < λ2 < · · · < λj < · · · satisfying limj→+∞ λj = +∞ be the eigenvalues of −∆ under the homogeneous Neumann boundary condition and φj be the normalized eigenfunction corresponding to λj . Let u∗ be the possible constant steady state of (1.2), then the linearized equation of (1.2) at u∗ is given by ∂ ∂t u(x, t) = d∆u(t, x)− u(t, x) + f ′(u∗)u(t− τ, x), x ∈ Ω, t ≥ 0, ∂ ∂n u(x, t) = 0, x ∈ ∂Ω t ≥ 0, where f(u) = puγe−u. Thus, u∗ is locally asymptotically stable in W 1,ς(Ω) with ς > N if all the solutions λ to the following characteristic equations admit negative real parts, Pn(λ, τ) := λ+ 1 + dλn − f ′(u∗)e−λτ = 0, n ∈ N0. (3.1) It follows from Pn(u + iv, τ) = 0 that u + 1 + dλn = f ′(u∗)e−τu cos(τv) and v = f ′(u∗)e−τu sin(τv), and hence (u+1+dλn)2+v2 = f ′2(u∗). This implies that all the zeros of Pn(·, τ) have negative real parts when 1 + dλn > |f ′(u∗)|. It follows from Pn(iω, τ) = 0 with ω > 0 that 1 + dλn = f ′(u∗) cos(τω) and ω = −f ′(u∗) sin(τω), and hence ω = ωn(u∗) := √ f ′2(u∗)− (1 + dλn)2. Thus, Pn(·, τ) has a pair of simple purely imaginary zeros ±iωn(u∗) when τ = τn,k(u∗) and |f ′(u∗)| > 1 + dλn, k ∈ N0, where τn,k(u∗) = 2kπ + θn(u∗) ωn(u∗) , 8 C. GUO, S. GUO EJDE-2022/67 θn(u∗) = [1 + sgn f ′(u∗)]π − sgn f ′(u∗) arccos δ + dλn f ′(u∗) . Thus, we have the following results. Lemma 3.1. (i) If 1 + dλn > |f ′(u∗)| for some n ∈ N0, then for every τ ≥ 0, all the zeros of Pn(·, τ) have negative real parts. (ii) If 1 + dλn < f ′(u∗) for some n ∈ N0, then for every τ ≥ 0, at least one zero of Pn(·, τ) has a positive real part. (iii) If f ′(u∗) < −δ − dλn for some n ∈ N0, then Pn(·, τ) has only zeros with negative real parts when 0 ≤ τ < τn,0(u∗), and exactly 2k zeros with positive real parts when τn,k−1(u∗) < τ ≤ τn,k(u∗), k ∈ N. Note that f ′(u∗0) = 0 < 1 + dλn for all n ∈ N0, then it follows from Lemma 3.1 that u∗0 is a locally asymptotically stable steady state of (1.2). Furthermore, we have the following result. Theorem 3.2. u∗0 is a locally asymptotically stable with respect to (1.2) for all p > 0 and τ ≥ 0. Moreover, (i) u∗0 is a globally asymptotically stable steady state of (1.2) for all τ ≥ 0 when p < p∗. (ii) Assume that p > p∗, then for every τ > 0 and initial value φ ∈ C0, the solution u(t, x) to (1.2) converges to u∗0 as t→∞, where C− = {ϕ ∈ C| 0 < φ(θ, x) < u∗1(p) for all (θ, x) ∈ [−τ, 0]× Ω}. Proof. (i) When p < p∗, we have f(u) u = p g(u) ≤ p∗ gmin = 1 and so f(u) < u for all u > 0. We shall establish the global stability of u∗0 by constructing a Lyapunov functional V (u)(t) = ∫ Ω u(t, x)dx+ ∫ 0 −τ ∫ Ω f(u(t+ s, x))dxds. By taking the time derivative of V along solutions of system (1.2), we have d dt ∣∣ (1.2) V (u)(t) = ∫ Ω [f(u(t, x))− u(t, x)] dx ≤ 0. where ≤ is actually = if and only if u = 0. Thus, by the Lyapunov-LaSalle invari- ance principle, u∗0 is globally asymptotically stable. (ii) Lemma 2.3(i) implies that C− is a positively invariant set of system (1.2) with p > p∗. When p > p∗, we have f(u) < u for 0 < u < u∗1(p). Thus, along an orbit u(t, x) of system (1.2) with any initial value φ ∈ C−, we also have d dt ∣∣ (1.2) V (u)(t) ≤ 0 This implies that the solution u(t, x) to (1.2) converges to u∗1(p) as t→∞. � It follows from p > p∗ and 0 < u∗1(p) < γ−1 < u∗2(p) that f ′(u∗1(p)) = γ−u∗1(p) > 1, which together with Lemma 3.1 implies that u∗1(p) is an unstable steady state of (1.2) for all p > 0 and τ ≥ 0. From Lemma 3.1 that we have the following result. EJDE-2022/67 NICHOLSON’S BLOWFLIES EQUATION WITH ALLEE EFFECT 9 Theorem 3.3. u∗1(p) is an unstable steady state of (1.2) for all p > p∗ and τ ≥ 0. Moreover, for each (n, k) ∈ N2 0 satisfying u∗1(p) < γ − 1 − dλn, system (1.2) with p > p∗ undergoes Hopf bifurcation near u = u∗1(p) and τ = τn,k(u∗1(p)). Finally, we consider the stability of the steady state u∗2(p) under the condition that p > p∗. Note that f ′(u∗2(p)) = γ − u∗2(p) < 1, then it follows from Lemma 3.1 that we see that (i) If u∗2(p) < γ+1, then for every τ ≥ 0, u∗2(p) is locally asymptotically stable; (ii) If γ + 1 < u∗2(p) < γ + 1 + dλ1, then u∗2(p) is locally asymptotically stable for 0 ≤ τ < τ2,0,0, and is unstable for τ > τ2,0,0, where τ2,0,0 = τ0,0(u∗2(p)); (iii) If γ+1+dλn < u∗2(p) < γ+1+dλn+1 for some n ∈ N, then u∗2(p) is locally asymptotically stable for 0 ≤ τ < max{τ2,j,0| 0 ≤ j ≤ n}, and is unstable for τ > max{τ2,j,0| 0 ≤ j ≤ n}, where τ2,j,0 = τj,0(u∗2(p)) for 0 ≤ j ≤ n. For each n ∈ N0, we set p∗n = (γ + 1 + dλn) 1−γ exp {γ + 1 + dλn} , then u∗2(p) = γ + 1 + dλn if and only if p = p∗n. Moreover, it is easy to see that p∗ < p∗n < p∗n+1 for all n ∈ N0. Thus, we have the following results. Theorem 3.4. (i) If p∗ < p < p∗0, then for every τ ≥ 0, u∗2(p) is locally asymptotically stable; (ii) If p∗0 < p < p∗1, then u∗2(p) is locally asymptotically stable for 0 ≤ τ < τ2,0,0, and is unstable for τ > τ2,0,0; (iii) If p∗n < p < p∗n+1 for some n ∈ N, then u∗2(p) is locally asymptotically stable for 0 ≤ τ < max{τ2,j,0| 0 ≤ j ≤ n}, and is unstable for τ > max{τ2,j,0| 0 ≤ j ≤ n}. Moreover, for each (j, k) ∈ {0, 1, 2, . . . , n} × N0, system (1.2) undergoes Hopf bifurcation near u = u∗2(p) and τ = τ2,j,k. 4. Nonconstant steady states Steady state solutions of (1.2) satisfy d∆u(x)− u(x) + f(u(x)) = 0, x ∈ Ω, ∂ ∂n u(x) = 0, x ∈ ∂Ω, (4.1) where, as in section 2, f(u) = puγe−u. In this section, we discuss the existence and nonexistence of nonconstant positive solutions of (4.1). Throughout the remaining part of this paper, the solutions refer to the classical solutions, by which we mean solutions in C2(Ω) ∩ C1(Ω). Similar to [11, Lemma 3.5], we have the following a priori estimate for nonnegative solutions for (1.1). Lemma 4.1. Assume that u(x) is a non-negative steady state solution of (1.2). (i) If p < p∗ then u(x) is exactly the constant solutions u∗0; (ii) If p ≥ p∗ then either u(x) is one of constant solutions u∗0, u∗1(p), and u∗2(p), or u(x) satisfies 0 < u(x) < u∗2(p) for all x ∈ Ω. Proof. Let x1 ∈ Ω be a maximum point of u; u(x1) = maxx∈Ω u(x). Then by using the maximum principle to (4.1), one can see f(u(x1))− u(x1) ≥ 0. If p < p∗ then f(u) < u for all u > 0 and hence u(x) = 0 for all x ∈ Ω. If p ≥ p∗ then we obtain 0 ≤ u ≤ u∗2(p) in Ω. If there exists x2 ∈ Ω such that u(x2) = 0, then u(x) ≡ 0 from strong maximum principle. Thus, if u is neither 0 nor u∗2(p), then from the 10 C. GUO, S. GUO EJDE-2022/67 strong maximum principle, we have 0 < u(x) < u∗2(p) for x ∈ Ω. The proof is complete. � Based on the above preparation, we are ready to derive a priori upper and lower bounds for all positive solutions to (4.1). Theorem 4.2. Assume that p ≥ p∗, then there exist two positive constants Ĉ and Č with Ĉ > Č such that any positive solution u(x) of (4.1) satisfies Č ≤ u(x) ≤ Ĉ for all x ∈ Ω. Proof. It follows from Lemma 4.1 that u(x) ≤ Ĉ := u∗2(p) for all x ∈ Ω. Let c(x) = puγ−1(x) exp{−u(x)} − 1. Then |c(x)| ≤ 1 + (γ − 1)γ−1 exp{1− γ}. From Harnack inequality (see [7, 8]), there exists a positive constant C such that sup Ω u(x) ≤ C inf Ω u(x). Hence it remains to prove that supΩ u(x) > c for some c > 0, which is independent of choice of solution. Similarly to the beginning of the proof of Lemma 4.1, we have f(u(x1))− u(x1) ≥ 0 and hence u(x1) ≥ u∗1(p), where x1 ∈ Ω is a maximum point of u and u(x1) = supx∈Ω u(x). Thus, we see that supΩ u(x) > 1 2u ∗ 1(p). The proof is complete. � Now we can show the nonexistence of positive nonconstant steady-state solutions when the diffusion coefficients d1 and d2 are large. Theorem 4.3. If p > p∗, then there exists d∗ > 0 such that the only nonnegative solutions to (4.1) with d > d∗ are u∗0, u∗1(p), and u∗2(p). Proof. Let u be a non-negative solution of (4.1) and denote ū = |Ω|−1 ∫ Ω u(x)dx. Then ∫ Ω [u(x)− ū]dx = 0. Multiplying (4.1) by u1 − ū1 and applying Lemma 4.1, we obtain d ∫ Ω |∇[u(x)− ū]|2dx = ∫ Ω [u(x)− ū][f(u(x))− u(x)]dx = ∫ Ω [u(x)− ū][f(u(x))− f(ū)]dx− ∫ Ω [u(x)− ū]2dx ≤ C1 ∫ Ω [u(x)− ū]2dx, where C1 depends on max{puγ−1e−u|γ − u| | 0 ≤ u ≤ u∗2(p)}. It follows from the Poincaré inequality that there exists a positive constant C2, independent of d, such that ∫ Ω [u(x)− ū]2dx ≤ C2 ∫ Ω |∇u(x)|2dx. Therefore, d ∫ Ω |∇[u(x)− ū]|2dx ≤ C3 ∫ Ω |∇[u(x)− ū]|2dx, where C3 is independent of d. Thus, if d > C3 then ∇[u(x)− ū] = 0, and hence u is a constant solution. The proof is complete. � EJDE-2022/67 NICHOLSON’S BLOWFLIES EQUATION WITH ALLEE EFFECT 11 In what follows, we use Leray-Schauder degree theory to show the existence of non-constant steady-state solutions when p > p∗. Recall the definition of Č and Ĉ from Theorem 4.2, we set X = {u ∈ C1(Ω) : ∂u ∂n = 0 on ∂Ω}, X0 = {u ∈ X : u ≥ 0}, X1 = {u ∈ X : Č/2 ≤ u ≤ 2Ĉ in Ω}, System (4.1) can be rewritten as −∆u = f(u)− u. (4.2) It is easy to see that u is a positive solution of (4.1) if and only if F(d, u) := u− (Id−∆)−1 [1 d f(u) + d− 1 d u ] = 0, u ∈ X0, where (Id−∆)−1 is the inverse of Id−∆ in X with the Neumann boundary condition. As F(d, ·) is a compact perturbation of the identity operator, the Leray-Schauder degree deg(F(d, ·),X1) is well defined from Theorem 4.2, and by the homotopy invariance, it is constant for all d when p > p∗. Thus, we have the following result. Lemma 4.4. Assume that p > p∗ and d > d∗. Then deg(F(d, ·),X1) = 2. Proof. In view of Theorem 4.3, we recall that if d > d∗ then F(d, ·) has exactly two zeros u∗1(p) and u∗2(p) in X1, and hence deg(F(d, ·),X1) = index(F(d, ·), u∗1(p)) + index(F(d, ·), u∗2(p)) Direct computation gives Fu(d, u∗j (p)) = Id− [1−Xj(d)] (Id−∆)−1, j = 1, 2. where Xj(d) = 1− f ′(u∗j (p)) d . It is easy to see that Fu(d, u∗2(p)) is invertible, i.e., 0 is not an eigenvalue of Fu(d, u∗j (p)). According to the Leray-Schauder index formula, we know index(F(d, ·), u∗j (p)) = (−1)γj(d), (4.3) where γj(d) is the number of negative eigenvalues (counting the algebraic multi- plicity) of the operator Fu(d, u∗j (p)). Then we shall count the negative eigenvalues of Fu(d, u∗j (p)). The linearized eigenvalue problem Fu(d, u∗j (p))ψ = σψ with ψ ∈ X can be represented by (1− σ)∆ψ = − [σ −Xj(d)]ψ. (4.4) Using a similar argument as that in Section 3, we see that σ is an eigenvalue of Fu(d, u∗j (p)) if and only if (σ − 1)λk + σ −Xj(d) = 0. for some k ∈ N0. Thus, eigenvalues of Fu(d, u∗j (p)) are given by σj,k(d) = λk +Xj(d) λk + 1 , where 0 = λ0 < λ1 < · · · < λn < · · · with limn→∞ λn = +∞ are the eigenvalues of the linear operator −∆ subject to the homogeneous boundary condition ∂ ∂nu = 0 on ∂Ω. 12 C. GUO, S. GUO EJDE-2022/67 Note that limd→∞X(d) = 0, then γj(d) = 0 for sufficiently large d and hence that index(F(d, ·), u∗j (p)) = 1 if d is large enough. This means that deg(F(d, ·),X1) = 2. The proof is complete. � Note that X1(0) = −∞, X2(0) = +∞, X1(d) < 0 < X2(d) and that X1(d) (respectively, X2(d)) is monotone increasing (respectively, decreasing) with respect to d. Obviously, index(F(d, ·), u∗2(p)) = 1 for all d > 0. Moreover, the monotone decreasing sequence {dj}∞j=j0 defined by dj := sup{d > 0 : X1(d) < −λj} for j ∈ N0 (4.5) satisfies limj→∞ dj = 0. Now, we have the following existence result for the non- constant steady state solutions: Theorem 4.5. Assume that p > p∗, then there exists at least one nonconstant solution of (4.1) if d ∈ (dj+1, dj) and j is even. Proof. If p > p∗, suppose for contradiction that there is no nonconstant solution of (1.1) under the assumptions of Lemma 4.4. Since u∗1(p) and u∗2(p) are solutions of (4.1) in X1, we have deg(F(d, ·),X1) = index(F(d, ·), u∗2(p)) + index(F(d, ·), u∗2(p)) = 2 (4.6) if d > d∗. It follows from the previous discussion that index(F(d, ·), u∗2(p)) = 1 for all d > 0. To derive index(F(d, ·), u∗1(p)), we count the number of negative eigenvalues of Fu(d, u∗1(p)) from the viewpoint of (4.4). If d ∈ (dj+1, dj), then (4.5) implies that all negative eigenvalues of Fu(d, u∗1(p)) consist of σ1,i(d), i = 0, 1, . . . , j. Therefore, the number γ1(d) of negative eigenvalues is γ1(d) = j + 1. If j is even, then (4.3) and (4.6) lead to deg(F(d, ·),X1) = 0. However, this contradicts Lemma 4.4. Then by the contradiction argument, we obtain at least one nonconstant solution if d ∈ (dj+1, dj) and j is even. � Remark 4.6. In Theorem 4.5, we investigate the existence of non-constant steady state solutions of (4.1) when p > p∗. However, from Theorem 4.5 we cannot draw any conclusion about the number and stability of these solutions. The main reason is that the explicit algebraic form of these solutions cannot be derived by means of topological methods. In the subsequent section, we will employ Lyapunov-Schmidt reduction to investigate the existence and multiplicity of non-constant steady state solutions. 5. Steady-state bifurcation It follows from the previous section that there is no steady-state bifurcation near u = u∗0 and u = u∗2(p) (in the later case, p > p∗). This section is devoted to the steady-state bifurcation near u = u∗1(p) when p > p∗. In view of Lemma 3.1, we see that system (1.2) with p > p∗ undergoes a steady-state bifurcation near u = u∗1(p) when u∗1(p) = γ − 1− dλn > 0 (equivalently, p = pn) for some n ∈ N, where pn = (γ − 1− dλn)1−γ exp{γ − 1− dλn}. Take p as a bifurcation parameter, and v(x) = u(x)− u∗1(p). Then system (4.1) is equivalent to the system d∆v(x)− (v(x) + u∗1(p)) + f(v(x) + u∗1(p)) = 0, x ∈ Ω, ∂ ∂n v(x) = 0, x ∈ ∂Ω. (5.1) EJDE-2022/67 NICHOLSON’S BLOWFLIES EQUATION WITH ALLEE EFFECT 13 We define F : X× R→ Y by F (v, p) = d∆v − (v + u∗1(p)) + f(v + u∗1(p)) for v ∈ X. Obviously, F (0, p) = 0 for all p > p∗. The Fréchet derivative of F with respect to v evaluated at (v, p) which is given by Lpv = d∆v + [γ − u∗1(p)− 1]v, v ∈ X. It is easy to see that Lp is invertible if and only if p = pn for some n ∈ N0. In the remaining part of this section, we focus on the following case: (H1) There exists n ∈ N such that pn > p∗ and that λn is a simple eigen- value of the linear operator −∆ subject to the homogeneous Neumann boundary condition on ∂Ω, with the associated eigenvector ϕn satisfying∫ Ω ϕ2 n(x)dx = 1. In this case, the kernel of Lpn is given by span{ϕn}, and that Lp is a self-adjoint operator, i.e., 〈u,Lpv〉 = 〈Lpu, v〉. Thus, we have the decompositions X = K ⊕ X1 and Y = K⊕ Y1, where K = span{ϕn} and X1 = {y ∈ X|〈v, y〉 = 0 for all v ∈ K}, Y1 = {y ∈ Y|〈v, y〉 = 0 for all v ∈ K}. Obviously, the operator Lpn : X→ Y is Fredholm with index zero. Lpn |X1 : X1 → Y1 is invertible and has a bounded inverse. Now, we can use Lyapunov-Schmidt reduction methods to reduce solving F (v, p) = 0 to finding solutions to F(v, p) = 0 in an open neighborhood U of (0, pn) in K×R, where F(v, p) = (I −Q)F (v + h(v, p), p), Q denotes the projection operators from Y onto Y1, and h : U→ X1 is a continuously differentiable map such that h(0, p) = 0 and QF (v + h(v, p), p) ≡ 0, (5.2) Substituting v = ϑϕ∗ ∈ K with ϑ ∈ R into F(v, p) = 0 and then calculating the inner product with ϕ∗ on Ω, we have g(ϑ, p) = 0, where G: R2 → R is explicitly given by G(ϑ, p) = ∫ Ω ϕn(x)F (ϑϕn(x) + h(ϑϕn(x), p), p)dx. (5.3) Notice that G(0, p) = 0, it follows that G: R2 → R takes the form G(ϑ, p) = ϑ[%(p− pn) +Kϑ+ χϑ2 + o(ϑ2)], where % = d dp ∣∣ p=pn u∗1(p) = −γ − 1− dλn pndλn < 0, K = f ′′(γ − 1− dλn) 2 ∫ Ω ϕ3 n(x)dx, χ = f ′′(γ − 1− dλn) 2 ∫ Ω hvv(ϕn(x), p)ϕ4 ∗(x)dx+ f ′′′(γ − 1− dλn) 6 ∫ Ω ϕ4 ∗(x)dx. It follows from (5.2) that Lpnhvv(ϕn, pn) +QFvv(ϕn, pn) = 0 and hence that hvv(ϕn, pn) = −f ′′(γ − 1− dλn)L−1 pn [ ϕ2 1 ] . Next we distinguish two cases to investigate the existence of nontrivial zero of G(·, p). We start with f ′′(γ − 1− dλn) 6= 0, that is, 1 + dλn 6= √ γ. Then we have K 6= 0. By using the implicit function theorem we see that there exist a constant ε > 0 and a continuously differentiable mapping ϑ: (pn − ε, pn + ε)→ R, such that G(ϑp, p) ≡ 0 for p ∈ (pn − ε, pn + ε). 14 C. GUO, S. GUO EJDE-2022/67 In fact, we have ϑp = %(pn − p) K + o(|p− pn|). (5.4) Theorem 5.1. If γ > 1 + dλn 6= √ γ for some n ∈ N satisfying (H1), then there exist a constant ε > 0 and a continuously differentiable mapping p → ϑp from (pn − ε, pn + ε) to R such that (4.1) has a nontrivial solution up = γ − 1− dλn + ϑpϕn + h(ϑpϕn, p), which exists for p ∈ (pn − ε, pn) ∪ (pn, pn + ε) and satisfies lim p→pn up = γ − 1− dλn. Next we consider the case where γ > 1 + dλn = √ γ and χ 6= 0, then the zeros of G(·, p) undergo a pitchfork bifurcation near p = pn. More precisely, if χ < 0 (respectively, χ > 0), then near the origin, one trivial zero ϑ = 0 and two nontrivial zeros ϑ = ϑ±p of G(·, p) exist for p < pn (respectively, p > pn), only one trivial zero ϑ = 0 exists for p > pn (respectively, p < pn). It follows from f ′′(γ − 1− dλn) = 0 that K = 0 and hence that χ = f ′′′(γ − 1− dλn) 6 ∫ Ω ϕ4(x)dx. In view of (5.3), we have the following result. Theorem 5.2. If γ > 1+dλn = √ γ for some n ∈ N satisfying (H1), and f ′′′(γ−1− dλn) < 0 (respectively, f ′′′(γ − 1 − dλn) > 0), then there exist a positive constant ε and two continuously differentiable mappings p → ϑ±p from (pn − ε, pn] to R (respectively, from [pn, pn + ε) to R) such that (4.1) has two nontrivial solutions u±p = γ − 1− dλn + ϑ±p ϕn + h(ϑ±p ϕn, p), which exists for p ∈ (pn − ε, pn] (respectively [pn, pn + ε) and satisfies lim p→pn u±p = γ − 1− dλn. Remark 5.3. In view of Theorem 3.3, u∗1(p) is unstable, and hence all the bifur- cated nonconstant steady states established in Theorems 5.1 and 5.2 are unstable. 6. Hopf bifurcation This section is devoted to the Hopf bifurcation at the nontrivial steady-state solution u∗2(p) of (1.2) under the following assumption (H2) There exists n ∈ N such that p > p∗n and that λn is a simple eigenvalue of the linear operator −∆ subject to the homogeneous Neumann boundary condi- tion on ∂Ω, with the associated eigenvector ϕn satisfying ∫ Ω ϕ2 n(x)dx = 1. Under this assumption, (3.1) with u∗ = u∗2(p) has a pair of simple purely imaginary solutions ±iωn(u∗2(p)) when τ = τn,k(u∗2(p)), k ∈ N0. Moreover, there exist ς > 0 and a continuously differentiable mapping λ: (τn,k(u∗2(p))− ς, τn,k(u∗(p)) + ς)→ C such that λ(τn,k(u∗2(p))) = iωn(u∗2(p)) and λ(τ) is a solution to (3.1) with u∗ = u∗2(p) for all τ ∈ (τn,k(u∗2(p))− ς, τn,k(u∗2(p)) + ς). In particular, Re{λ′(τn,k(u∗2(p)))} > 0. A Hopf bifurcation is said to be forward if there exist time-periodic solutions when parameter value τ > τn,k(u∗2(p)), and to be backward if τ < τn,k(u∗2(p)). For convenience, we set u∗ = u∗2(p), τ∗ = τn,k(u∗2(p)), ω∗ = ωn(u∗2(p)). EJDE-2022/67 NICHOLSON’S BLOWFLIES EQUATION WITH ALLEE EFFECT 15 Let T = 2π/ω∗, and CT (respectively, C1 T ) be the set of continuous (respectively, differentiable) T -periodic mappings from R into X2 C. If we denote ‖x‖0 = max t∈[0,T ] {‖u(t)‖XC} for u ∈ CT , and ‖u‖1 = max{‖u‖0, ‖u̇‖0} for u ∈ C1 T , then CT and C1 T are Banach spaces when they are endowed with the norms ‖ · ‖0 and ‖ · ‖1, respectively. It is easy to see that CT is a Banach representation of the group S1 with the action given by θ · u(t) = u(t+ θ) for θ ∈ S1. In view of complexity in analyzing (1.2), we introduce the inner product (·, ·): CT × CT → R defined by (v, u) = 1 T ∫ T 0 〈v(t), u(t)〉dt for u, v ∈ CT . Let ν ∈ (−1, 1), v(t) = u((1 + ν)t). Then equation (1.2) can be rewritten as (1 + ν) dv dt = [d4− 1]v + f(u∗ + v(t− (1 + ν)τ)). Define F: C1 T × R2 → CT by F(u, α, ν) = −(1 + ν) du dt + [d4− 1]v + f(u∗ + v(t− (1 + ν)τ)). (6.1) By varying the newly introduced small variable ν, one keeps track not only of solutions of (1.2) with period T but also of solutions with nearby period. In fact, solutions to F(u, α, ν) = 0 correspond to T 1+ν -periodic solutions of (1.2). It follows that F is S1-equivariant: θ · F(u, τ, ν) = F(θ · u, τ, ν), for all θ ∈ S1. Define Lτu = −du dt + [d4− 1]v + f ′(u∗)v(t− τ). The elements of kerLτ correspond to solutions of the linear system Lτu = 0 sat- isfying u(t) = u(t + T ). For convenience in computation we shall allow functions with range C instead of R. With respect to the inner product (·, ·): CT × CT → R, the adjoint operator of Lτ is L∗τu = d dt u(t) + [d4− 1]v + f ′(u∗)v(t+ τ). It follows that kerLτ∗ = span{ζ0, ζ0} and kerL∗τ∗ = span{ζ∗0 , ζ ∗ 0}, where ζ0, ζ ∗ 0 ∈ CT are defined as ζ0(t) = ϕneiω∗t, ζ∗0 (t) = w∗ϕneiω∗t, and w∗ = (1 + τ∗f ′(u∗)e−iω∗τ∗)−1. Obviously, spaces kerLτ∗ , rangeLτ∗ , and W = (kerL∗τ∗)⊥∩C1 T are S1-invariant subspaces of CT . Moreover, CT = kerLτ∗⊕rangeLτ∗ and C1 T = kerLτ∗ ⊕W. By Lyapunov-Schmidt reduction, we can reduce our Hopf bifurcation problem to the problem of finding zeros of the map B : kerLτ∗ × R2 → kerL∗τ∗ given by B(v, τ, ν) ≡ (I − P )F(v +W (v, τ, ν), τ, ν), (6.2) 16 C. GUO, S. GUO EJDE-2022/67 where P denotes the projection operator from CT onto rangeLτ∗ along kerL∗τ∗ , W : kerL× R2 →W is a continuously differentiable S1-equivariant map such that W (0, τ∗, 0) = 0 and PF(v +W (v, τ, ν), τ, ν) ≡ 0. (6.3) We refer to B as the bifurcation map of the system (1.2). It follows from the S1-equivariance of F and W that the bifurcation map B is also S1- equivariant. Moreover, B(0, τ∗, 0) = 0 and Bv(0, τ∗, 0) = 0. For each φ ∈ kerLτ∗ , φ = ςζ0 + ςζ0, where ς = 〈ζ∗0 , φ〉. Let G(ς, τ, ν) = (ζ∗0 ,B(ςζ0 + ςζ0, τ, ν)). Thus, we only need to consider the existence of nontriv- ial solutions to G(ς, τ, ν) = 0. It follows that Gς(0, τ∗, 0) = 0, Gς(0, τ∗, 0) = 0. (6.4) It is easy to see that G(·, τ, ν) is S1-equivariant. Using a similar arguments to that in [4], we can find two functions R, I : R3 → R such that G(ς, τ, ν) = R(|ς|2, τ, ν)ς + I(|ς|2, τ, ν)iς. (6.5) It follows from Gς(0, τ∗, 0) = 0 that R(0, τ∗, 0) = 0 and I(0, τ∗, 0) = 0. Let ς = reiθ. Then solving G is equivalent to either solve r = 0 or R(r2, τ, ν) = 0 and I(r2, τ, ν) = 0. Note that Gτ (ς, τ∗, 0) = (ζ∗0 ,Fτ (v, τ∗, 0)) = ςλ′(τ∗) +O(|ς|2), Gν(ς, τ∗, 0) = (ζ∗0 ,Fβ(v, τ∗, 0)) = −iω∗ς +O(|ς|2). Then the Jacobi determinant of functions R and I with respect to τ and ν is det [ Rτ (0, τ∗, 0) Rν(0, τ∗, 0) Iτ (0, τ∗, 0) Iν(0, τ∗, 0) ] = −ω∗Re{λ′(τ∗)} 6= 0. The implicit function theorem implies that there exists a unique function τ = τ(r2) and ν = ν(r2) satisfying τ(0) = τ∗ and ν(0) = 0 such that R(r2, τ(r2), ν(r2)) ≡ 0, I(r2, τ(r2), ν(r2)) ≡ 0 (6.6) for all sufficient small r. Therefore, g(ς, τ(|ς|2), ν(|ς|2)) ≡ 0 for ς sufficiently near 0. Therefore, system (1.2) has a bifurcation of periodic solutions. Namely, we have the following result. Theorem 6.1. In addition to assumption (H2), a Hopf bifurcation for (1.2) occurs at τ = τn,k(u∗2(p)). Namely, in a neighborhood of (u, τ) = (u∗2(p), τn,k(u∗2(p))) there is a branch of periodic solutions Uτ (x, t) satisfying Uτ (x, t) → u∗2(p) as τ → τn,k(u∗2(p)). The period Tτ of Uτ (x, t) satisfies that Tτ → 2π/ωn(u∗2(p)) as τ → τn,k(u∗2(p)). Note that dλn 6= 0, then the Hopf bifurcating periodic solutions Uτ (x, t) is spatially nonhomogeneous. In view of (6.3), we have PF(ςζ0 + ςζ0 + W (ςζ0 + ςζ0, τ, ν), τ, ν) ≡ 0. Write W (ςζ0 + ςζ0, τ∗, 0) and G(ς, τ∗, 0) as W (ςζ0 + ςζ0, τ∗, 0) = ∑ k+l≥2 1 k!l! Wklς kς l, G(ς, τ∗, 0) = ∑ k+l≥2 1 k!l! Gsklςkς l. It follows from (6.5) that G21 = R1(0, τ∗, 0) + iI1(0, τ∗, 0), where R1(u, τ, ν) = Ru(u, τ, ν) and I1(u, τ, ν) = Iu(u, τ, ν). Therefore, R1(0, τ∗, 0) = Re{G21} and EJDE-2022/67 NICHOLSON’S BLOWFLIES EQUATION WITH ALLEE EFFECT 17 I1(0, τ∗, 0) = Im{G21}. From (6.6), we can calculate the derivatives of τ(r2) and ν(r2) and evaluate at r = 0: τ ′(0) = −2 Re{G21}, ν′(0) = −2 Im{z′(τ∗)G21}. The bifurcation direction is determined by sgn τ ′(0), and the monotonicity of period of bifurcating closed invariant curve depends on sgn ν′(0). Note that G(ς, τ∗, 0) = iω∗ς + ( ζ∗0 ,F(ςζ0 + ςζ0 +W (zζ0 + zζ0, τ∗, 0), τ∗, 0) ) , then we have G21 = ( ζ∗0 , 2Sτ∗ [ζ0,W11] + Sτ∗ [ζ0,W20] + Eτ∗ [ζ0, ζ0, ζ0] ) , (6.7) where Sτ∗ [ϕ,ψ] = f ′′(u∗)ϕ(−τ∗)ψ(−τ∗), Eτ∗ [ϕ,ψ, φ] = f ′′′(u∗)ϕ(−τ∗)ψ(−τ∗)φ(−τ∗). We still need to compute W11 and W20. In fact, W20 = −L−1 τ∗ PSτ∗ [ζ0, ζ0], W11 = −L−1 τ∗ PEτ∗ [ζ0, ζ0]. Note that G20 = (ζ∗0 ,Sτ∗ [ζ0, ζ0]) = 0 and g11 = (ζ∗0 ,Sτ∗ [ζ0, ζ0]) = 0. Namely, Sτ∗ [ζ0, ζ0], Sτ∗ [ζ0, ζ0] ∈ rangeLτ∗ . Hence, the projection P on each of Sτ∗ [ζ0, ζ0] and Sτ∗ [ζ0, ζ0] acts as the identity. Thus, W20 = −f ′′(u∗)[d∆− 1 + f ′(u∗)e−2iω∗τ∗ − 2iω∗] −1ϕ2 ne2iω∗(·−τ∗), W11 = −f ′′(u∗)[d∆− 1 + f ′(u∗)]−1ϕ2 n, (6.8) and so G21 = f ′′′(u∗)w∗ exp{−iω∗τ∗} ∫ Ω ϕ4 n(x)dx − 2[f ′′(u∗)]2w∗ exp{−iω∗τ∗} ∫ Ω ϕ2 n(x)[d∆− 1 + f ′(u∗)]−1ϕ2 n(x)dx − [f ′′(u∗)]2w∗ exp{−3iω∗τ∗} ∫ Ω ϕ2 n(x)[d∆− 1 + f ′(u∗)e−2iω∗τ∗ − 2iω∗] −1 × ϕ2 n(x)dx. It is a standard result (see, for example, Guo and Wu [6] and Wu [10]) that for classical Hopf bifurcations, subcritical bifurcating periodic solutions are unstable, while supercritical bifurcating periodic solutions have the same stability as the trivial solution had before the bifurcation. We summarize the above discussion as follows. Theorem 6.2. Under the assumption (H2), for each fixed k ∈ N0, there exists a branch of periodic solutions, parameterized by τ , bifurcating from the point (u, τ) = (u∗(p), τn,k(u∗2(p))). Moreover, (i) Re{G21} determines the direction of the bifurcation and the stability of bifurcating periodic solutions: the bifurcating periodic solutions exist for τ > τn,k(u∗2(p)) (respectively, τ < τn,k(u∗2(p))), and have the same stability as the nontrivial steady state u∗2(p) had before the bifurcation (respectively, are unstable) if Re{G21} < 0 (respectively, > 0); 18 C. GUO, S. GUO EJDE-2022/67 (ii) Im{λ′(τn,k(u∗2(p)))G21} determines the period of the bifurcating periodic so- lutions along the branch: the period is greater than (respectively, smaller than) 2π/ωn(τn,k(u∗2(p))) if it is positive (respectively, negative). In Theorem 3.3, we know that u∗1(p) is an unstable steady state of (1.2) for all p > p∗ and τ ≥ 0. Moreover, for each (n, k) ∈ N2 0 satisfying u∗1(p) < γ − 1 − dλn, system (1.2) with p > p∗ undergoes Hopf bifurcation near u = u∗1(p) and τ = τ1,n,k(u∗1(p)). Using a similar argument, we can discuss the Hopf bifurcation at the nontrivial steady-state solution u∗1(p) of (1.2) and obtain the following result. Theorem 6.3. For each (n, k) ∈ N2 0 satisfying u∗1(p) < γ − 1− dλn, system (1.2) with p > p∗ undergoes Hopf bifurcation near u = u∗1(p) and τ = τn,k(u∗1(p)). Namely, in a neighborhood of (u, τ) = (u∗1(p), τn,k(u∗1(p))) there is a branch of unstable periodic solutions Uτ (x, t) satisfying Uτ (x, t)→ u∗1(p) as τ → τn,k(u∗1(p)). The period Tτ of Uτ (x, t) satisfies that Tτ → 2π/ωn(u∗1(p)) as τ → τn,k(u∗1(p)). 7. Discussion This article presents a delayed diffusive model, for which the growth rate not only depends on the present quantities but also its past quantities. The presence of time delay implies that the model depends on the solution at an earlier time in addition to derivatives. In fact, time delay occurs in many real life process, hence mathematical models with time delay may give much more reasonable insights into biological models. The present contribution presents a new mathematical model to understand how Allee effects arise in the model outcomes. Some substantial changes in the dynamics are observed in the model as the parameters vary. In this regard, steady-state bifurcation has been observed for the model. However, more complicated dynamics is observed with delay induced instability which enables the appearance of Hopf bifurcations in the system. In this article, we compute analytically primary branches of steady state solu- tions and periodic solutions. In particular, we present the bifurcation direction for each branch of steady state solutions and periodic solutions. Here, we should mention that the analytical results for the direction of Hopf bifurcation and sta- bility properties of the bifurcating solutions can be determined using the theory of normal form and central manifold reduction. In this concept, a second order approximation of the center manifold can be exploited and computation of the first Lyapunov ratio can be used to determine properties of the periodic branches aris- ing from Hopf points. In this article, however, we employ the Lyapunov-Schmidt reduction to obtain the existence of the bifurcating solutions and to see how the time delay significantly affects the direction of Hopf bifurcations. Motivated by the results given in this paper, another potential extension of the work would be to further investigate global continuation of nontrivial steady state solutions and nonlinear waves in system (1.2). In fact, we can use the global bifurcation theory to show that nontrivial steady state solutions exist for p not only near to but also far away from pn, and that these bifurcations of nonlinear waves exist for τ not only near to but also far away from the critical values τn,k(u∗2(p)). It is assumed in (H1) and (H2) that λn is a simple eigenvalue of the linear operator −∆ subject to the homogeneous Neumann boundary condition on ∂Ω. A natural question is what happens to system (1.2) when eigenvalue λn is m-multiple (m ≥ 2). In this case, the associated eigenspace and hence the kernel space K is EJDE-2022/67 NICHOLSON’S BLOWFLIES EQUATION WITH ALLEE EFFECT 19 m-dimensional, and then we can employ Lyapunov-Schmidt reduction to obtain a bifurcation mapping g : Rm × R → Rm. However, it is difficult to investigate the existence, multiplicity, and patterns of the bifurcation map g : Rm × R → Rm. Some of these results will be reported later. Acknowledgements. This work was partially supported by the National Natural Science Foundation of China (Grant No. 12071446), and by the Fundamental Re- search Funds for the Central Universities, China University of Geosciences (Wuhan) (Grant No. CUGST2). References [1] W. C. Allee; The Social Life of Animals, Heinemann, London, 1931. [2] F. Brauer, C. Castillo-Chávez; Mathematical Models in Population Biology and Epidemiol- ogy. Springer-Verlag, New York (2001). [3] S. Buedo-Fernández, E. Liz; On the stability properties of a delay differential neoclassical model of economic growth, Electron. J. Qual. Theory Differ. Equ., 2018 (2018), 1–14. [4] M. Golubitsky, I. Stewart, D. G. Schaeffer; Singularities and Groups in Bifurcation Theory, Springer-Verlag, New York, 1988. [5] S. Guo; Patterns in a nonlocal time-delayed reaction-diffusion equation, Z. Angew. Math. Phys. 69 (2018), 10 [6] S. Guo, J. Wu; Bifurcation Theory of Functional Differential Equations, Springer-Verlag, New York, 2013. [7] C. S. Lin, W. M. Ni, I. Takagi; Large amplitude stationary solutions to a chemotaxis system, J. Differential Equations, 72 (1988), 1–27. [8] R. Peng, J. P. Shi, M. X. Wang; On stationary patterns of a reaction-diffusion model with autocatalysis and saturation law, Nonlinearity 21 (2008), 1471–1488. [9] A. J. Terry; Impulsive adult culling of a tropical pest with a stage-structured life cycle, Nonlinear Anal. Real World Appl., 11 (2010), 645–664. [10] J. Wu; Theory and Applications of Partial Functional-Differential Equations, Springer- Verlag, New York, 1996. [11] F. Q. Yi, J. J. Wei, J. P. Shi; Bifurcation and spatiotemporal patterns in a homogeneous diffusive predator-prey system, J. Differential Equations, 246 (2009), 1944-1977. [12] T. Yi, X. Zou; Global attractivity of the diffusive Nicholson blowflies equation with Neumann boundary condition: a non-monotone case, Journal of Differential Equations, 245(2008), 3376–3388. [13] X.-Q. Zhao; Global attractivity in a class of nonmonotone reaction-diffusion equations with time delay, Canadian Applied Mathematics Quarterly, 17(2009), 271–281. Cuiping Guo School of Science, China University of Geosciences, Beijing 100083, China Email address: gcp@cugb.edu.cn Shangjiang Guo (corresponding author) School of Mathematics and Physics, China University of Geosciences, Wuhan, Hubei 430074, China. Center for Mathematical Sciences, China University of Geosciences, Wuhan, Hubei 430074, China Email address: guosj@cug.edu.cn 1. Introduction 2. Preliminaries 3. Stability of constant steady states 4. Nonconstant steady states 5. Steady-state bifurcation 6. Hopf bifurcation 7. Discussion Acknowledgements References