Electronic Journal of Differential Equations, Vol. 2025 (2025), No. 38, pp. 1–19. ISSN: 1072-6691. URL: https://ejde.math.txstate.edu, https://ejde.math.unt.edu DOI: 10.58997/ejde.2025.38 TRAVELING WAVES OF A DIFFUSIVE MODIFIED LESLIE-GOWER MODEL WITH CHEMOTAXIS SHUNA WANG, JIANG LIU, JUN FANG, XIAOJIE LIN Abstract. In this article, we study a diffusive modified Leslie-Gower model with chemotaxis and large wave speed. By applying traveling wave transformation and changing the time scale, this modified Leslie-Gower model can be transformed into a singularly perturbed system. We establish the existence of heteroclinic orbits connecting different equilibria for the system with- out perturbation by constructing invariant regions and using the Poincaré-Bendixson theorem. Then the existence of traveling wave solutions for the diffusive modified Leslie-Gower system is demonstrated via the geometric singular perturbation theory and Fredholm theory. 1. Introduction In ecosystems there exist a large number of predation relationships among species. Numerous scholars have conducted a series of investigations into interactions between predators and prey [9, 11, 19, 30], which play an important role in revealing the dynamic evolution of populations. To describe that the population density of predator is restricted by preys, Leslie and Gower [15, 16] improved the classical Lotka-Volterra model [20, 28] and proposed the Leslie-Gower model du dt = u(α1 − βu− θ1v), dv dt = v(α2 − θ2v u ), (1.1) where u and v represent the density of the prey and predator, respectively. The increase of the predator population follows the logistic growth model. The carrying capacity of environment is positively correlated with the number of prey population. Based on the Leslie-Gower model (1.1), many modified versions have been considered. As mentioned in [1], the predator can seek alternative food for survival instinct if its preferred food is extremely scarce. Due to the scarcity of preferred food, the growth rate of predator can decline. To model this ecological phenomenon, authors added a constant k to the system (1.1), and yielded du dt = u(α1 − βu− θ1v), dv dt = v(α2 − θ2v u+ k ), (1.2) where k represents the maximum decline rate to indicate the environmental protection. Models (1.1) and (1.2) are both ordinary differential equations, which are constructed based on the assumption that the spatial distribution of populations is uniform. However, the uneven dis- tribution of resources leads to an unbalanced spatial distribution of organisms in nature. In order to survive, the population may migrate and diffuse from high density areas to low density areas. 2020 Mathematics Subject Classification. 35K57, 92D25, 35C07, 34D15. Key words and phrases. Modified Leslie-Gower model; chemotaxis; traveling waves; geometric singular perturbation. ©2025. This work is licensed under a CC BY 4.0 license. Submitted February 21, 2025. Published April 10, 2025. 1 2 S. WANG, J. LIU, J. FANG, X. LIN EJDE-2025/38 Mathematically, this phenomenon can be described by the self-diffusion of species. Therefore, many researchers focused on the following reaction-diffusion system with Leslie-Gower term ∂u ∂t = d1∆u+ u(α1 − βu− θ1v), ∂v ∂t = d2∆v + v(α2 − θ2v u ). (1.3) Several investigations have been carried out with respect to model (1.3). Du and Hsu [4] proved the existence of postively steady-state solutions with prescribed spatial patterns. Moreover, they compared with the classical Lotka-Volterra model, some crucial differences in the dynamic behavior were observed. Hamizah M. Safuan et al. [26] found the existence of traveling waves and derived the minimum wave speed. Furthermore, they also carried out the stability of traveling waves. With the use of the Lyapunov function and the transformation technique, Zhou and Wei [33] obtained a new global stability result of the positive equilibrium and generalized results to the more general reaction-diffusion system. In addition to the self-diffusion of predator and prey, the diffusion of predator may be influenced by prey. This influence gives rise to chemotactic movement, making the predator population migrate towards areas of high prey density. In 1987, based on the data of prey gliding in ladybugs, Kareiva and Odell [13] firstly proposed a predator-prey system with prey-taxis to explain the predation between ladybugs and aphids. In the ecological system, chemotaxis plays an significant role in maintaining ecological balance and biodiversity, such as regulating prey species (pest) to prevent outbreaks. Furthermore, models incorporating chemotaxis may produce different spatial patterns and bring new dynamical behaviors. Recently, Li [17] studied a two-species system with chemotaxis ∂u ∂t = ∆u+ u(α1 − βu− θ1v), ∂v ∂t = ∆v − χ∇ · (v∇u) + ϑuv − ϱv, (1.4) where the chemotaxis term χ∇ · (v∇u) expresses the tendency of the predator to move upward or downward along the gradient direction of prey, and χ denotes the prey-taxis coefficient. When χ > 0 (χ < 0), the chemotaxis is characterized as attractive (repulsive). From a biological perspective, χ > 0 implies that predator tends to diffuse towards areas where the density of prey is higher to enhance capture rate; Conversely, χ < 0 represents the predator tends to diffuse towards areas where the density of prey is lower to avoid prey group defense. That is, the spatiotemporal dynamics of predators is influenced by the prey density. Under Neumann boundary conditions, Li proved that system (1.4) admits a unique global bounded classical solution. Currently, various biological models with chemotaxis have been proposed and extensively studied [6, 18, 22]. Zhao and Hu in [31] investigated the classical solutions of a prey-taxis model with Sigmoid function response under homogeneous Neumann boundary conditions. Moreover, they demonstrated that the prey-taxis sensitivity coefficient destabilizes the stability of the homogeneous steady state when the prey defend. Qiu and Guo [24] focused on a Leslie-Gower model with chemotaxis. Based on the asymptotic analysis and bifurcation theory, they obtained the local and global bifurcation of nonconstant steady-states by taking the chemotaxis coefficient as a bifurcation parameter. Inspired by the above mentioned papers, we investigate the diffusion model with chemotaxis ∂u ∂t = d1∆u+ a1u(1− b1u− r1v), ∂v ∂t = d2∆v + χ∇ · (v∇u) + a2v(1− r2v u+ k ), (1.5) where u and v respectively represent the density of the prey and predator at spatial position x and time t. d1 and d2 are diffusion coefficients of populations. ∇ · (v∇u) accounts for the chemotaxis phenomenon. χ > 0 stands for the sensitivity of chemotaxis. a1 and a2 response the intrinsic growth rate of prey and predator. b1 reflects the intensity of competition among members of prey species. r1 and r2 are the maximum values of the average decrease rate for organisms. k indicates the environmental protection. All the parameters mentioned above are positive. EJDE-2025/38 TRAVELING WAVES OF A DIFFUSIVE LESLIE-GOWER MODEL 3 For biological population systems, traveling waves can describe the behaviors of species, such as propagation, migration, and invasion. Therefore, studying the traveling waves of population systems has great significance. In recent years, the qualitative properties of the traveling waves of the Leslie-Gower model have also been widely studied [2, 8, 25, 27]. For example, Tian and Zhang [27] explored the traveling waves for a modified Leslie-Gower model with continuous diffusion or discrete diffusion based on the squeeze method and a Lyapunov function. Guo and Cheng [8] considered a prey-predator system of Leslie-Gower type under shifting environment, and the existence of traveling waves was proved by utilizing the method of constructing the suitable upper and lower solution with a monotone iteration technique. On the other hand, the interactions of populations occur on different time scales in nature. The dynamical behaviors observed in these systems also reflect these multi-scale characteristics [3, 10]. The main theoretical basis for studying such problems is the geometric singular perturbation theory [7, 12]. In recent years, the geometric singular perturbation theory has also been widely applied to many other equations [5, 21, 29, 32]. In this article, we focus on the existence of traveling waves for the modified Leslie-Gower model with chemotaxis (1.5) and large wave speed. The highlights of our paper mainly lie in the following aspects. • Several Leslie-Gower models investigated in previous works can be considered as special instances of our model. This includes the models presented in references [26, 23]. • We transform the modified Leslie-Gower model with chemotaxis (1.5) into a fast-slow system. And the existence of traveling waves that connect different equilibrium points can be obtained by applying the geometric singular perturbation theory. • We obtain two types of traveling waves when system (1.5) has no co-existence equilibrium and three types of traveling waves when system (1.5) admits the co-existence equilibrium. Then we generalize the results obtained in [18], which proved the existence of traveling waves between the co-existence equilibrium and (0, 0). The rest of this article is structured as follows. In Section 2, we provide the geometric singular perturbation theory and present some preliminaries on system (1.5). In Section 3, the existence of heteroclinic solutions for the slow subsystem (2.12) is proved, based on the phase plane analysis along with Poincaré-Bendixson theorem. In Section 4, the existence of the traveling waves for system (1.5) is established via the application of geometric singular perturbation theory and Fredholm theory. In Section 5, we present a main summary of our findings. 2. preliminaries In this section, we outline the geometric singular perturbation theory. Moreover, we investigate the existence and stability of boundary equilibria and the co-existence equilibrium. 2.1. Geometric singular perturbation theory. Firstly, we introduce the definitions and con- clusions related to geometric singular perturbation [12, 14]. Consider the differential equations x′ = f(x, y, ε), y′ = εg(x, y, ε), (2.1) where ′ = d dt denotes the differentiation with respect to the fast time scale t. Fast variable x ∈ Rp, slow variable y ∈ Rq with p, q ≥ 1, f : Rp × Rq ×R → Rp, g : Rp × Rq ×R → Rq. The functions f and g are sufficiently smooth. And the positive parameter ε is sufficiently small, which can be regarded as separation of time scales. Set t = τ ε , (2.1) has the equivalent form εẋ = f(x, y, ε), ẏ = g(x, y, ε), (2.2) where · = d dτ denotes the differentiation with respect to the slow time scale τ . Based on the relationship between t and τ , system (2.1) and system (2.2) are referred to as the fast system and slow system. 4 S. WANG, J. LIU, J. FANG, X. LIN EJDE-2025/38 Setting ε = 0 in (2.1) and (2.2), one has the fast subsystem (layer system) x′ = f(x, y, 0), y′ = 0, (2.3) and the slow subsystem (reduced system) 0 = f(x, y, 0), ẏ = g(x, y, 0). (2.4) Definition 2.1. The set defined by the algebraic equations in the slow subsystem (2.4) M = {(x, y) ∈ Rp × Rq : f(x, y, 0) = 0} is called the critical set, which is also the set of equilibria for the fast subsystem (2.3). If M is a submanifold of Rp × Rq, then M is referred to as a critical manifold. Definition 2.2. The critical manifold M is normally hyperbolic if the linearization of (2.3) at all point on M has exactly p eigenvalues off the imaginary axis. Lemma 2.3 ([12, 14]). Consider the system (2.1), if submanifold M0 of critical manifold M is compact and normally hyperbolic, then for 0 < r < +∞ and 0 < ε≪ 1, the following conclusions hold, (i) there exists a slow manifold Mε that lies within O(ε) of M0 is diffeomorphic to M0, and Mε is locally invariant under the flow of (2.1). This means that the flow of system (2.1) only can enter or leave through the boundary of Mε; (ii) there exist Cr-smooth unstable manifolds Wu loc(Mε) and stable manifolds W s loc(Mε), which are O(ε) close to and diffeomorphic to Wu loc(M0) and W s loc(M0), respectively. In addition, Wu loc(Mε) and W s loc(Mε) are locally invariant with respect to the flow of system (2.1). Definition 2.4. A traveling wave solution of system (1.5) can be expressed as z = x− ct, (u, v)(x, t) = (U, V )(z), (2.5) where c > 0 represents the wave speed. Applying the traveling wave transformation (2.5) to system (1.5), we obtain the wave profile system −cUz = d1Uzz + a1U(1− b1U − r1V ), −cVz = d2Vzz + χ(UzVz + UzzV ) + a2V (1− r2V U + k ), (2.6) where Uz = dU dz and Uzz = d2U dz2 . Setting ξ = − z c , system (2.6) is transformed into εÜ = U̇ − a1U(1− b1U − r1V ), εV̈ = 1 ρ V̇ − εγ(U̇ V̇ + ÜV )− a2 ρ V (1− r2V U + k ), (2.7) where ˙ = d dξ , ρ = d2 d1 , γ = χ d1ρ and ε = d1 c2 satisfying 0 < ε ≪ 1. System (2.7) is equivalent to the slow syetem U̇ = U1, V̇ = V1, εU̇1 = U1 − a1U(1− b1U − r1V ), εV̇1 = −εγU1V1 − γV [U1 − a1U(1− b1U − r1V )] + 1 ρ V1 − a2 ρ V (1− r2V U + k ), (2.8) EJDE-2025/38 TRAVELING WAVES OF A DIFFUSIVE LESLIE-GOWER MODEL 5 where ξ indicates the slow time scale. Setting ζ = ξ ε in system (2.8) results in the fast system U ′ = εU1, V ′ = εV1, U1 ′ = U1 − a1U(1− b1U − r1V ), V1 ′ = −εγU1V1 − γV [U1 − a1U(1− b1U − r1V )] + 1 ρ V1 − a2 ρ V (1− r2V U + k ), (2.9) where ′ = d dζ , and ζ indicates the fast time scale. Setting ε = 0 in slow system (2.8) and fast system (2.9), we respectively obtain the slow subsystem U̇ = U1, V̇ = V1, 0 = U1 − a1U(1− b1U − r1V ), 0 = 1 ρ V1 − a2 ρ V (1− r2V U + k ), (2.10) and the fast subsystem U ′ = 0, V ′ = 0, U1 ′ = U1 − a1U(1− b1U − r1V ), V1 ′ = −γV [U1 − a1U(1− b1U − r1V )] + 1 ρ V1 − a2 ρ V (1− r2V U + k ). (2.11) Based on Definition 2.1, the critical manifold of fast system (2.9) is M0 = { (U,U1, V, V1) : U1 = a1U(1− b1U − r1V ), V1 = a2V (1− r2V U + k ) } , which is two-dimensional and consists of equilibria of the fast subsystem (2.11). The Jacobian matrix of system (2.11) at any point of M0 reads 0 0 0 0 0 0 0 0 −a1(1− 2b1U − r1V ) a1r1U 1 0 a1γV (1− 2b1U − r1V )− a2r2V 2 ρ(U+k)2 W (U, V, U1) −γV 1 ρ  , where W (U, V, U1) = −γU1 + a1γU(1− b1U − 2r1V ) + a2 ρ ( 2r2VU+k − 1). After computations, we obtain that the matrix has two positive eigenvalues 1, 1 ρ , and two zero eigenvalues. Therefore, the dimension of the fast variables is equal to the number of eigenvalues off the imaginary axis. By Definition 2.2, we deduce that the critical manifold M0 is normally hyperbolic. According to Lemma 2.3, for 0 < ε ≪ 1, there exists a perturbed locally invariant manifold Mε which is O(ε)-close to M0. 2.2. Equilibria of the slow subsystem (2.12). The dynamics behavior restricted on M0 is determined by the slow subsystem U̇ = a1U(1− b1U − r1V ), V̇ = a2V (1− r2V U + k ). (2.12) The existence and topological properties of boundary equilibria for the slow subsystem (2.12) can be characterized by the following lemma. Lemma 2.5. The slow subsystem (2.12) admits three boundary equilibria E0(0, 0), E1( 1 b1 , 0) and E2(0, k r2 ). Among them, (i) E0 is always an unstable node; 6 S. WANG, J. LIU, J. FANG, X. LIN EJDE-2025/38 (ii) E1 is always a saddle; (iii) E2 is a stable node for 0 < r2 < r1k, E2 is a saddle for r2 > r1k, and E2 is a saddle-node for r2 = r1k. Proof. It is evident that the slow subsystem (2.12) admits three boundary equilibria E0(0, 0), E1( 1 b1 , 0) and E2(0, k r2 ). (i) The linearized matrix of the slow subsystem (2.12) at E0(0, 0) is J(E0) = ( a1 0 0 a2 ) , which has two positive real eigenvalues a1 and a2, thus the equilibrium E0 is an unstable node. (ii) The linearized matrix of the slow subsystem (2.12) restricted to E1( 1 b1 , 0) is J(E1) = ( −a1 −a1r1 b1 0 a2 ) , which has one negative eigenvalue −a1 and one positive eigenvalue a2, thus the equilibrium E1 is a saddle. (iii) The Jacobian matrix of the slow subsystem (2.12) at E2(0, k r2 ) is J(E2) = ( a1(1− r1k r2 ) 0 a2 r2 −a2 ) . If 0 < r2 < r1k, then 1 − r1k r2 < 0. It follows that the matrix J(E2) has two negative real eigenvalues a1(1− r1k r2 ) and −a2. Thus the equilibrium E2 is a stable node. If r2 > r1k, then 1 − r1k r2 > 0. One obtains that the matrix J(E2) has one positive eigenvalue a1(1− r1k r2 ) and one negative eigenvalue −a2, that is the equilibrium E2 is a saddle. If r2 = r1k, then 1− r1k r2 = 0. And the Jacobian matrix J(E2) has one negative real eigenvalue −a2 and one zero eigenvalue, which indicates that E2 is a high-order equilibrium. To conduct a detailed analysis regarding the type of this high-order equilibrium E2, we apply the linear translation x = U , y = V − k r2 , which moves the equilibrium E2(0, k r2 ) to the origin (0, 0). Then by rescaling time via ξ = r2(x+ k)τ , the slow subsystem (2.12) becomes dx dτ = −a1r2x(x+ k)(b1x+ r1y), dy dτ = a2(r2y + k)(x− r2y). (2.13) Let X = 1 r2 x and Y = − 1 r2 x+ y, the system (2.13) can be rewritten as dX dτ = −a1r2k(b1r2 + r1)X 2 − a1r1r2kXY − a1r2 2(b1r2 + r1)X 3 − a1r1r2 2X2Y, dY dτ = −a2r2kY − a2r2 2Y 2 + a2r2k(b1r2 + r1)X 2 + a2r2 2k(b1r2 + r1)X 3 + a2r1r2 2X2Y. (2.14) Using dη = −a2r2kdτ in the system (2.14) yields dX dη = a1(b1r2 + r1) a2 X2 + a1r1 a2 XY + a1r2(b1r2 + r1) a2k X3 + a1r1r2 a2k X2Y, dY dη = Y + r2 k Y 2 − (b1r2 + r1)X 2 − r2(b1r2 + r1) k X3 − r1r2 k X2Y. (2.15) Since the coefficient of X2 is a1(b1r2+r1) a2 > 0, E2 is a saddle-node. Moreover, the parabolic sector is in the right halfplane and two hyperbolic sectors are in the left halfplane. In order to study the EJDE-2025/38 TRAVELING WAVES OF A DIFFUSIVE LESLIE-GOWER MODEL 7 dynamics near (0, 0), we introduce the polar coordinates X = Rcosθ and Y = Rsinθ. Then the system (2.15) can be given by dR dη = −a2r2kRsin2θ − a1r2k(b1r2 + r1)R 2cos3θ − a2r2 2R2 sin3 θ + [a2r2k(b1r2 + r1)− a1r1r2k]R 2cos2θ sin θ +O(R3), dθ dη = −a2r2ksinθ cos θ +O(R). (2.16) The polar coordinate transformation replaces dynamical behaviors near (0, 0) with the unit circumference S1. It is obvious that there are four singularities (0, θi) on invariant set {0} × S1, where θi = 0, π2 , π, 3π 2 , i = 1, 2, 3, 4. For 0 < R ≪ 1, we can obtain that dR dη < 0 holds along θ = θ1, θ2, θ4, and dR dη > 0 holds along θ = θ3. In conclusion, the trajectory distribution nearby the unit circumference is illustrated in Figure 1 (a). Moreover, by reducing the unit circumference to (0, 0), we can get the distribution of trajectories near the origin (0, 0) of the XOY plane, which is shown in Figure 1 (b). Therefore, any trajectory originating from (X0, Y0) near (0, 0) (X0 ≥ 0) will tend to the origin (0, 0). It follows that the equilibrium E2(0, k r2 ) of the slow subsystem (2.12) is a attractive saddle-node. The proof is complete. Figure 1 (c) illustrates the dynamics of the system (2.12) around E2 when a1 = a2 = 2, b1 = r1 = 1, r2 = 5 and k = 5. □ Now we show the existence and linear stability of the co-existence equilibrium. Lemma 2.6. (i) If r2 ≤ r1k, the slow subsystem (2.12) has no co-existence equilibrium. (ii) If r2 > r1k, the slow subsystem (2.12) has one co-existence equilibrium E∗( r2−r1k r1+b1r2 , 1+b1k r1+b1r2 ), which is stable. Proof. The co-existence equilibrium E∗(U∗, V∗) is a point that the non-trivial prey and predator isoclines intersect. It must satisfy V = − b1 r1 U + 1 r1 , V = 1 r2 U + k r2 . By direct calculations, we can obtain that U∗ = r2−r1k r1+b1r2 and V∗ = 1+b1k r1+b1r2 . Since the co-existence equilibrium E∗(U∗, V∗) lies in the first quadrant, for the case r2 ≤ r1k, the slow subsystem has no co-existence equilibrium; and for the case r2 > r1k, the slow subsystem (2.12) admits one unique co-existence equilibrium E∗( r2−r1k r1+b1r2 , 1+b1k r1+b1r2 ). We analyze linear stability of the co-existence equilibrium when r2 > r1k. The Jacobian matrix of the slow subsystem (2.12) restricted to E∗( r2−r1k r1+b1r2 , 1+b1k r1+b1r2 ) is as follows J(E∗) = ( −a1b1(r2−r1k) r1+b1r2 −a1r1(r2−r1k) r1+b1r2 a2 r2 −a2 ) . Then the determinant and trace of J(E∗) are det(J(E∗)) = a1a2(b1r2 + r1)(r2 − r1k) r2(r1 + b1r2) > 0, tr(J(E∗)) = −a1b1(r2 − r1k) r1 + b1r2 − a2 < 0. Therefore, the co-existence equilibrium E∗( r2−r1k r1+b1r2 , 1+b1k r1+b1r2 ) is stable. □ Remark 2.7. Since the equilibrium E2 is a high-order equilibrium when r2 = r1k, for the slow subsystem (2.12), the investigation of the dynamics in a small neighbourhood of E2 is challenging. We leave the case for future study. 8 S. WANG, J. LIU, J. FANG, X. LIN EJDE-2025/38 (a) (b) (c) Figure 1. (a) Diagram depicting the dynamics after the blow-up at (0, 0), (b) trajectory distribution near (0, 0), (c) phase portrait of the system (2.12) around E2 with a1 = a2 = 2, b1 = r1 = 1, r2 = 5 and k = 5. 3. Heteroclinic orbits of the slow subsystem (2.12) In this section, we explore the existence of heteroclinic orbits for the slow subsystem (2.12). Based on Lemma 2.6, the existence results are divided into two cases. Case 1: When r2 < r1k, the slow subsystem (2.12) has no co-existence equilibrium. In this case, we respectively demonstrate the existence of heteroclinic orbits from E0 to E2 and from E1 to E2 in Section 3.1, as shown in Theorem 3.3. Case 2: When r2 > r1k, the slow subsystem (2.12) possesses a unique co-existence equilibrium E∗. In this case, we have the existence of heteroclinic orbits between E0, E1, E2 and E∗ in Section 3.2, as shown in Theorem 3.5. 3.1. Existence of heteroclinic orbits for the slow subsystem (2.12) without co-existence equilibrium. From Lemma 2.5 and Lemma 2.6, we determine that if r2 < r1k, the slow subsystem (2.12) only possesses three boundary equilibria E0(0, 0), E1( 1 b1 , 0) and E2(0, k r2 ). Moreover, E0 is an unstable node, E1 is a saddle, and E2 is a stable node. From the phase portrait shown in Figure 2, we find that there are two different types of orbits connecting E0(0, 0) to E2(0, k r2 ) and E1( 1 b1 , 0) to E2(0, k r2 ). EJDE-2025/38 TRAVELING WAVES OF A DIFFUSIVE LESLIE-GOWER MODEL 9 Figure 2. Phase portrait of the system (2.12) with a1 = a2 = 2, b1 = r1 = 1, r2 = 3 and k = 5. We define the region Ω1 ∈ R2 as Ω1 = {(U, V ) : 0 ≤ U ≤ 1 b1 , 0 ≤ V ≤ b1k + 1 b1r2 }, which is displayed in Figure 3 (a). Lemma 3.1. If r2 < r1k, then Ω1 is a positively invariant region of the slow subsystem (2.12). Proof. The region Ω1 is surrounded by four lines Γ1 = {(U, V )|U = 0, 0 ≤ V ≤ b1k + 1 b1r2 }, Γ2 = {(U, V )|V = 0, 0 ≤ U ≤ 1 b1 }, Γ3 = {(U, V )|U = 1 b1 , 0 ≤ V ≤ b1k + 1 b1r2 }, Γ4 = {(U, V )|V = b1k + 1 b1r2 , 0 ≤ U ≤ 1 b1 }. (a) (b) Figure 3. Schematic of invariant regions for system (2.12) with r2 < r1k, where Ω1 in (a), and Ω2 in (b). We shall demonstrate that no trajectory of the slow subsystem (2.12) initiated within the region Ω1 can leave Ω1 through the boundaries Γi(i = 1, 2, 3, 4). This indicates that for any point 10 S. WANG, J. LIU, J. FANG, X. LIN EJDE-2025/38 (U, V ) ∈ Γi(i = 1, 2, 3, 4), it satisfies the condition ni · (U̇ , V̇ ) ∣∣ Γi ≤ 0, where ni indicates the outward normal vector of Γi at the point (U, V ). It is obvious that Γ1 and Γ2 are both invariant sets of system (2.12) and no trajectories of (2.12) can intersect the boundaries Γ1 or Γ2. The boundary Γ3 has the outer normal vector n3 = (1, 0). Then for all points (U, V ) ∈ Γ3, we find that n3 · (U̇ , V̇ ) |Γ3= −a1r1 b1 V ≤ 0. The boundary Γ4 has the outer normal vector written as n4 = (0, 1). For any point (U, V ) ∈ Γ4, we have n4 · (U̇ , V̇ ) |Γ4= a2(b1k + 1)(b1U − 1) b1 2r2(U + k) ≤ 0. Thus, the direction of the vector field along the boundaries Γ3 and Γ4 points into Ω1. It follows that no orbits of the slow subsystem (2.12) can leave the boundaries Γi, i = 1, 2, 3, 4. The proof is complete. The invariant region Ω1, which is surrounded by four red lines, with arrows indicating the flow direction along the boundaries is shown in Figure 3 (a). □ Subsequently, we prove that the unstable manifold of E1( 1 b1 , 0) starting from E1 will enter the region Ω1. Lemma 3.2. The unstable manifold emanating from E1( 1 b1 , 0) of the slow subsystem (2.12) will intersect the region Ω1, that is Wu(E1) ∩ Ω1 ̸= ∅. Proof. According to Lemma 2.5, the linearization matrix at E1 has two different eigenvalues −a1 < 0 and a2 > 0. Consequently, the boundary equilibrium E1 possesses a one-dimensional stable manifold and a one-dimensional unstable manifold. The stable manifold approaches E1 along the U-axis. The local unstable manifold of E1 is tangential to ς, which is the corresponding eigenvector of the positive eigenvalue a2 and reads as ς := (−1, (a1 + a2)b1 a1r1 )T . Thus Wu(E1) ∩ Ω1 ̸= ∅, for the slope of ς is − (a1+a2)b1 a1r1 < 0. □ Using Lemmas 3.1 and 3.2, we obtain the existence of heteroclinic orbits between E0, E1 and E2 for the slow subsystem (2.12). Theorem 3.3. If r2 < r1k, then the following statements hold. (i) The slow subsystem (2.12) has a unique heteroclinic orbit (U, V ) connecting the equilibrium E1( 1 b1 , 0) to E2(0, k r2 ). (ii) The slow subsystem (2.12) has infinite number of heteroclinic orbits (U, V ) that connects the equilibrium E0(0, 0) to E2(0, k r2 ). Proof. (i) We denote P (U, V ) = a1U(1− b1U − r1V ), Q(U, V ) = a2V (1− r2V U + k ). Assume that there exists a closed orbit L within the region Ω1. The rotation number of the closed orbit L with respect to the vector field (P,Q) is 1. In fact, the rotation number of L with respect to the vector field (P,Q) must be 0 due to no equilibrium in Ω1, which leads to a contradiction. This implies that if there exists a closed orbit, its interior must contain at least one equilibrium. Therefore, it is evident that the slow subsystem (2.12) has no closed orbit in Ω1. Our discussion for the type of equilibria in Lemma 2.5 reveals that E1 is a saddle, characterized by a one-dimensional unstable manifold. Furthermore, Lemma 3.2 indicates that the unstable manifold will enter Ω1. Therefore, according to the Poincaré-Bendixson theorem, the ω limit of the unstable manifold of E1 in Ω1 is identified as E2. Consequently, the slow subsystem (2.12) admits a unique heteroclinic orbit connecting E1 to E2, which is denoted by Γ̃. EJDE-2025/38 TRAVELING WAVES OF A DIFFUSIVE LESLIE-GOWER MODEL 11 (ii) We define a region Ω2, which is bounded by two segments Γ2, E0E2 and the orbit Γ̃. Figure 3 (b) displays Ω2 with red curves, with arrows indicating the flow direction on the boundary. Since Γ2, Γ̃ and E0E2 are all orbits of the slow subsystem (2.12), it is clear that Ω2 is invariant as orbits cannot intersect. All equilibria in this region are located on the boundaries, hence the slow subsystem (2.12) does not have any closed orbit in Ω2. Based on Lemma 2.5, we know that E0 is an unstable node which is repelling, while E2 is a stable node which is attracting. Therefore, by applying the Poincaré-Bendixson theorem, the slow subsystem (2.12) possesses infinitely many heteroclinic orbits connecting the equilibria E0 and E2. □ A numerical simulation to illustrate Theorem 3.3 is presented in Figure 4. Figure 4. Numerical simulation diagrams of two types of heteroclinic orbits, as described in Theorem 3.3, generated from system (2.12) with different initial values, where a1 = a2 = 2, b1 = r1 = 1, r2 = 3 and k = 5. 3.2. Existence of heteroclinic orbits for the slow subsystem (2.12) with a unique co- existence equilibrium E∗. Based on Lemma 2.5 and Lemma 2.6, when r2 > r1k, the slow subsystem (2.12) has a unique co-existence equilibrium E∗( r2−r1k r1+b1r2 , 1+b1k r1+b1r2 ) and three boundary equilibria E0(0, 0), E1( 1 b1 , 0), E2(0, k r2 ). Among them, E0 is an unstable node, E1 and E2 are two saddles, and co-existence equilibrium E∗ is stable. As depicted in Figure 5, the phase portrait of system (2.12) reveals that there exist three distinct types of heteroclinic orbits. Subsequently, we shall present a detailed demonstration. From P (U, V ) = a1U(1− b1U − r1V ), Q(U, V ) = a2V (1− r2V U + k ), we obtain (PU +QV )|U∗ = a1(1− 2b1U∗ − r1V∗) + a2(1− 2r2V∗ U∗ + k ) = a1(1− 2b1 r2 − r1k r1 + b1r2 − r1 1 + b1k r1 + b1r2 ) + a2[1− 2r2(r1 + b1r2) r2 − r1k + k(r1 + b1r2) 1 + b1k r1 + b1r2 ] = −a1b1(r2 − r1k) r1 + b1r2 − a2 < 0. Therefore, the sign of PU +QV remains unchanged within a sufficiently small neighborhood of E∗. According to Dulac theorem, there is no closed orbits in this neighborhood. Moreover, since E∗ is locally stable, the trajectory originating from the point ( r2−r1k r1+b1r2 + ϵ, 1+b1k r1+b1r2 + ϵ) must converge to E∗, where the constant ϵ is sufficiently small. We denote this trajectory by L̂. 12 S. WANG, J. LIU, J. FANG, X. LIN EJDE-2025/38 Figure 5. Phase portrait of the system (2.12) with a1 = a2 = 2, b1 = r1 = 1, r2 = 6 and k = 4. Furthermore, for b1U + r1V − 1 ≥ 0, if U ≥ U∗, we obtain b1U̇ + r1V̇ = r1a2V (1− r2 U + k 1− b1U r1 ) = a2V U + k [(r1 + r2b1)U + r1k − r2] ≥ a2V U + k [(r1 + r2b1) r2 − r1k r1 + b1r2 + r1k − r2] = a2V U + k (r2 − r1k + r1k − r2) = 0. Thus, for the trajectory L̂, we have U̇ ≤ 0, that is the U -component of L̂ is monotonically decreasing. Similarly, for r2V − U − k ≥ 0, if U ≥ U∗, we obtain r2V̇ − U̇ = −a1U(1− b1U − r1V ) = −a1U(1− b1U − r1 U + k r2 ) = a1U r2 [(r1 + b1r2)U + (r1k − r2)] ≥ a1U r2 [(r1 + b1r2) r2 − r1k r1 + b1r2 + (r1k − r2)] = a1U r2 (r2 − r1k + r1k − r2) = 0. Then, for the trajectory L̂, we have V̇ ≤ 0, that is the V -component of L̂ is monotonically decreasing. Consequently, there exists ξ̂ ∈ R such that the trajectory L̂ intersects the line U = 1 b1 at a point ϕξ̂( r2−r1k r1+b1r2 + ϵ, 1+b1k r1+b1r2 + ϵ) = P̂ ( 1 b1 , â) as ξ → −∞. Then, we have the following lemma. Lemma 3.4. If r2 > r1k, then the region Σ1 is a positively invariant region for the slow subsystem (2.12). EJDE-2025/38 TRAVELING WAVES OF A DIFFUSIVE LESLIE-GOWER MODEL 13 Proof. We define a region Σ1 ∈ R2, which is bounded by three segments and a trajectory Li (i = 1, 2, 3, 4), where L1 := {(U, V )|U = r2 − r1k r1 + b1r2 , 0 ≤ V ≤ b1k + 1 r1 + b1r2 }, L2 := {(U, V )|V = 0, r2 − r1k r1 + b1r2 ≤ U ≤ 1 b1 }, L3 := {(U, V )|U = 1 b1 , 0 ≤ V ≤ â}, L4 := {(U, V )|ϕξ( r2 − r1k r1 + b1r2 + ϵ, 1 + b1k r1 + b1r2 + ϵ), ξ ∈ [ξ̂,+∞)} . The outer normal vector of L1 is N1 = (−1, 0). For any point (U, V ) ∈ L1, we have N1 · (U̇ , V̇ )|L1 = −a1 r2 − r1k r1 + b1r2 (1− b1 r2 − r1k r1 + b1r2 − r1V ) ≤ r2 − r1k r1 + b1r2 (1− b1 r2 − r1k r1 + b1r2 − r1 1 + b1k r1 + b1r2 ) = 0. Based on Lemma 3.1, no trajectory of the slow subsystem (2.12) can leave through the bound- aries L2 and L3. Therefore, any trajectory of the slow subsystem (2.12) starting inside Σ1 cannot leave Σ1 through its boundaries Li (i = 1, 2, 3). Since L4 is a trajectory of the slow subsystem (2.12) and trajectories cannot intersect, Σ1 is a positively invariant region. The proof is complete. A schematic of Σ1 is shown in Figure 6 (a), where the flow directions along boundaries are marked by arrows. □ (a) (b) Figure 6. Schematic of invariant regions for system (2.12) with r2 > r1k, where Σ1 in (a), and Σ2 in (b). The slow subsystem (2.12) admits heteroclinic orbits connecting between the equilibria E0(0, 0), E1( 1 b1 , 0), E2(0, k r2 ), and E∗( r2−r1k r1+b1r2 , 1+b1k r1+b1r2 ), which can be stated as follows. Theorem 3.5. If r2 > r1k, then the following statements hold. (i) the slow subsystem (2.12) has a single heteroclinic orbit that connects the equilibrium E1( 1 b1 , 0) to E∗( r2−r1k r1+b1r2 , 1+b1k r1+b1r2 ). (ii) the slow subsystem (2.12) has infinitely many heteroclinic orbits connecting the equilibrium E0(0, 0) to E∗( r2−r1k r1+b1r2 , 1+b1k r1+b1r2 ). (iii) the slow subsystem (2.12) has a unique heteroclinic orbit connecting the equilibrium E2(0, k r2 ) to E∗( r2−r1k r1+b1r2 , 1+b1k r1+b1r2 ). Proof. (i) Since there is no equilibrium inside Σ1, it is evident that the slow subsystem (2.12) cannot have a closed orbit in Σ1. By applying the Poincaré-Bendixson theorem, any trajec- tory starting from the unstable manifold of E1 must eventually approach the stable co-existence 14 S. WANG, J. LIU, J. FANG, X. LIN EJDE-2025/38 equilibrium E∗. Hence, the slow subsystem (2.12) features a heteroclinic orbit connecting E1 to E∗. Additionally, it follows that E1 possesses a one-dimensional unstable manifold, therefore this heteroclinic orbit is unique, and we denote it by L̃. (ii) Define the region Σ2 as the area enclosed by three segments Γ2, E0A, E∗A and one orbit L̃, which is shown in Figure 6 (b). According to Lemma 3.1, the trajectories can not leave Σ2 from the segments Γ2 and E0A. Additionally, the outwards vector of E∗A is written by NE∗A = (0, 1). For all points (U, V ) ∈ NE∗A, one has NE∗A · (U̇ , V̇ ) = a2 b1k + 1 r1 + b1r2 [ 1− r2(b1k + 1) (U + k)(r1 + b1r2) ] ≤ a2 b1k + 1 r1 + b1r2 [1− r2(b1k + 1)(r1 + b1r2) (r2 + b1r2k)(r1 + b1r2) ] = 0. Since orbits cannot intersect, it is evident that Σ2 is a positively invariant region for the slow subsystem (2.12). All equilibria lie on the boundaries. Hence, the slow subsystem (2.12) does not possess any closed orbits within Σ2. Drawing from Lemma 2.5 and Lemma 2.6, we recognize that E0 is repelling, while E∗ is attracting. Consequently, applying the the Poincaré-Bendixson theorem, the slow subsystem (2.12) contains an infinite number of heteroclinic orbits that connect the equilibria E0 and E∗. (iii) Since E2 is a saddle when r2 > r1k, the equilibrium E2 of system (2.12) has a one- dimensional unstable manifold. Through a straightforward calculation, we know that the eigen- vector ν of the Jacobian matrix J(E2) is (1, a2 a1(r2−r1k)+a2r2 )T , which is directed towards the interior of the region Σ2. Therefore, the unstable manifold of E2 tangent to ν will enter Σ2. Applying Poincaré-Bendixson theorem, there exists only one heteroclinic orbit from E2 to E∗. The proof is complete. A numerical simulation demonstrating the heteroclinic orbits of the slow subsystem (2.12) is shown in Figure 7. □ Figure 7. Numerical simulation diagrams of three types of heteroclinic orbits, as described in Theorem 3.5, where a1 = a2 = 2, b1 = r1 = 1, r2 = 6 and k = 2. EJDE-2025/38 TRAVELING WAVES OF A DIFFUSIVE LESLIE-GOWER MODEL 15 4. Traveling wave solutions of system (1.5) In this section, we focus on investigating the presence of traveling wave solution for system (1.5) with the large wave speed. This problem is equivalent to studying the heteroclinic orbits for the slow system (2.8). The main tools we use are the geometric singular perturbation theory and Fredholm theory. Referring to Lemma 2.3, we have that there admits a slow manifold Mε, which is O(ε)-close to M0 for sufficiently small ε > 0. Furthermore, Mε and M0 are diffeomorphic. The slow manifold Mε can be expressed as Mε = { (U,U1, V, V1) : U1 = P (U, V ) + ϕ̄(U, V, ε), V1 = Q(U, V ) + ψ̄(U, V, ε) } , where P (U, V ) = a1U(1− b1U − r1V ), Q(U, V ) = a2V (1− r2V U + k ), ϕ̄ and ψ̄ are both smooth functions of ε that satisfy ϕ̄(U, V, 0) = 0 and ψ̄(U, V, 0) = 0. Thus these functions ϕ̄ and ψ̄ can be expanded in a Taylor series with respect to the ε as follows ϕ̄(U, V, ε) = εϕ1(U, V, ε) +O(ε2), ψ̄(U, V, ε) = εψ1(U, V, ε) +O(ε2). By substituting U1 = P (U, V ) + ϕ̄(U, V, ε), V1 = Q(U, V ) + ψ̄(U, V, ε), into the slow system (2.8), we obtain εU̇1 = ε( ∂U1 ∂U U̇ + ∂U1 ∂V V̇ ) = ε [ (PU + εϕ1U )U̇ + (PV + εϕ1V )V̇ ] +O(ε2) = ε(PUU1 + PV V1) +O(ε2), and εV̇1 = ε( ∂V1 ∂U V̇ + ∂V1 ∂V V̇ ) = ε [ (QU + εψ1U )U̇ + (QV + εψ1V )V̇ ] +O(ε2) = ε(QUU1 +QV V1) +O(ε2), where PU = a1(1− 2b1U − r1V ), PV = −a1r1U, QU = a2r2V 2 (U + k)2 , QV = a2(1− 2r2V U + k ). By examining the equations of the slow system (2.8), we find that εU̇1 = U1 − a1U(1− b1U − r1V ) = a1U(1− b1U − r1V ) + ϕ1(U, V )ε+O(ε2)− a1(1− b1U − r1V ) = ϕ1(U, V )ε+O(ε2), and εV̇1 = −εγU1V1 − γV [U1 − a1U(1− b1U − r1V )] + 1 ρ V1 − a2 ρ V (1− r2V U + k ) = −εγ [ a1U(1− b1U − r1V ) + ϕ1(U, V )ε+O(ε2) ][ a2V (1− r2V U + k ) 16 S. WANG, J. LIU, J. FANG, X. LIN EJDE-2025/38 + ψ1(U, V )ε+O(ε2) ] − γV [ a1U(1− b1U − r1V ) + ϕ1(U, V )ε +O(ε2)− a1U(1− b1U − r1V ) ] + 1 ρ [a2V (1− r2V U + k ) + ψ1(U, V )ε+O(ε2)] − a2 ρ V (1− r2V U + k ) = −εγ [ a1a2UV (1− b1U − r1V )(1− r2V U + k ) ] − γV ϕ1(U, V )ε + 1 ρ ψ1(U, V )ε+O(ε2) = [ − γa1a2UV (1− b1U − r1V )(1− r2V U + k )− γV ϕ1(U, V ) + 1 ρ ψ1(U, V ) ] ε+O(ε2). Comparing the coefficients of ε in the above equations, we obtain ϕ1(U, V ) = PUU1 + PV V1, ψ1(U, V ) = γρa1a2UV (1− b1U − r1V )(1− r2V U + k ) + γρV (PUU1 + PV V1) + ρ(QUU1 +QV V1), and Mε = { (U,U1, V, V1) : U1 = P (U, V ) + ε(PUU1 + PV V1) +O(ε2), V1 = Q(U, V ) + ε [ γρa1a2UV (1− b1U − r1V )(1− r2V U + k ) + γρV (PUU1 + PV V1) + ρ(QUU1 +QV V1) ] +O(ε2) } . Then the slow system (2.8) restricted to Mε is described by dU dξ = a1U(1− b1U − r1V ) + ε(PUU1 + PV V1) +O(ε2), dV dξ = a2V (1− r2V U + k ) + ε { γρa1a2UV (1− b1U − r1V )(1− r2V U + k ) + γρV (PUU1 + PV V1) + ρ(QUU1 +QV V1) } +O(ε2). (4.1) It is evident that the system (4.1) is reduced to the slow subsystem (2.12) as ε→ 0. Next, we demonstrate that for 0 < ε≪ 1, the system (4.1) admits a heteroclinic orbit connecting E1( 1 b1 , 0) to E2(0, k r2 ) under the condition r2 < r1k. For ε = 0, by Theorem 3.3, it follows that the system (4.1) has a heteroclinic orbit connecting E1 to E2, written as (U0, V0). For 0 < ε≪ 1, we set U = U0 + εÛ1 + . . . , V = V0 + εV̂1 + . . . . (4.2) Substituting this expansion in the system (4.1) and equating the coefficients of ε, we derive the differential system that determines Û1 and V̂1 dφ(ξ) dξ −H(ξ)φ(ξ) = C(ξ), (4.3) where φ(ξ) = ( Û1 V̂1 ) , H(ξ) = ( a1(1− 2b1U0 − r1V0) −a1r1U0 a2r2V0 2 (U0+k)2 a2(1− 2r2V0 U0+k ) ) , C(ξ) = ( c1 c2 ) , EJDE-2025/38 TRAVELING WAVES OF A DIFFUSIVE LESLIE-GOWER MODEL 17 c1 = a1 2U0(1− 2b1U0 − r1V0)(1− b1U0 − r1V0)− a1a2r1U0V0(1− r2V0 U0 + k ), c2 = γρa1a2U0V0(1− b1U0 − r1V0)(1− r2V0 U0 + k ) + γρV0 [ a1 2U0(1− 2b1U0 − r1V0)(1− b1U0 − r1V0) − a1a2r1U0V0(1− r2V0 U0 + k ) ] + ρa1a2r2 U0V0 2 (U0 + k)2 (1− b1U0 − r1V0) + ρa2 2V0(1− 2r2V0 U0 + k )(1− r2V0 U0 + k ). Subsequently, we demonstrate that system (4.3) admits a solution that satisfies the boundary condition Û1(±∞) = 0, V̂1(±∞) = 0. We define the operator L = d dξ −H(ξ), and let L2 represent the space of square-integrable functions, with inner production given by∫ +∞ −∞ (X(ξ), N(ξ))dξ, where the notation (·, ·) denotes the Euclidean inner product on R2. According to the Fredholm theory, the system (4.3) has a solution if and only if∫ +∞ −∞ (X(ξ), C(ξ))dξ = 0, for all functions X(ξ) ∈ R2 that lie in the kernel of the adjoint operator L. It can be easily confirmed that the adjoint operator L∗ is given by L∗ = − d dξ −HT (ξ), where HT (ξ) = ( a1(1− 2b1U0 − r1V0) a2r2V0 2 (U0+k)2 −a1r1U0 a2(1− 2r2V0 U0+k ) ) . To determine the kernel of L∗, we should seek all X(ξ) satisfying L∗X(ξ) = 0, that is, dX(ξ) dξ = −HT (ξ)X(ξ). (4.4) Then the question of persistence reduces to the solvability of the equation (4.4). It is clear that the zero solution is a solution to equation (4.4). Since the matrix HT (ξ) is nonconstant, it is challenging to find the general solution to equation (4.4). However, we are specifically interested in solutions that meet the boundary condition X(±∞) = 0, and indeed, the sole such solution is the zero solution. Recall that (U0(ξ), V0(ξ)) is the solution of the system (2.12) derived in Theorem 3.3. Although we do not have an explicit expression for it, it follows that (U0(ξ), V0(ξ)) tends to E2(0, k r2 ) as ξ → +∞. As ξ → +∞, the matrix −HT (ξ) finally becomes a constant matrix JT = ( −a1(1− r1k r2 ) 0 −a2 r2 a2 ) . (4.5) Through a direct computation, we find the matrix J admits two positive real eigenvalues λ1 = −a1(1 − r1k r2 ) and λ2 = a2. As a result, the only solution that meets X(+∞) = 0 is the zero solution. This implies that the Fredholm orthogonality condition holds trivially∫ +∞ −∞ (X(ξ), C(ξ))dξ = ∫ +∞ −∞ (0, C(ξ))dξ = 0, 18 S. WANG, J. LIU, J. FANG, X. LIN EJDE-2025/38 and the system (4.4) exists a solution satisfying the boundary condition Û1(±∞) = 0, V̂1(±∞) = 0. Therefore, for sufficiently small ε > 0, the system (4.1) has a heteroclinic orbit that connects E1( 1 b1 , 0) to E2(0, k r2 ). For system (4.1) with 0 < ε ≪ 1, the existence of the heteroclinic orbits connecting from E0(0, 0) to E2(0, k r2 ), as well as the existence of the heteroclinic orbit linking E0(0, 0), E1( 1 b1 , 0) and E2(0, k r2 ) to E∗( r2−r1k r1+b1r2 , 1+b1k r1+b1r2 ), can be established in a similar way. Consequently, we get the following theorem. Theorem 4.1. There exists a large wave speed c∗ > 0. Under the assumption c > c∗, with regard to the traveling waves (U, V )(z) = (u, v)(x, t) of system (1.5), where z = x − ct, the following results hold. (i) If r2 > r1k, then the system (1.5) has two types of traveling waves. The first type connects the equilibria E1( 1 b1 , 0) and E2(0, k r2 ). The second type connects E0(0, 0) and E2(0, k r2 ). (ii) If r2 < r1k, then the system (1.5) has three types of traveling waves. The first type connects E0(0, 0) and E∗( r2−r1k r1+b1r2 , 1+b1k r1+b1r2 ). The second type connects E1( 1 b1 , 0) and E∗( r2−r1k r1+b1r2 , 1+b1k r1+b1r2 ). The third type connects E2(0, k r2 ) and E∗( r2−r1k r1+b1r2 , 1+b1k r1+b1r2 ). 5. Conclusion In this article, we studied a diffusive modified Leslie-Gower model with chemotaxis (1.5), in which the diffusion of the predator is also influenced by prey. Through our research, it is revealed that for a large sufficiently wave speed c > 0, there exist five kinds of traveling waves for the system (1.5). These traveling waves of the diffusive modified Leslie-Gower model with chemotaxis (1.5) are associated with the heteroclinic orbits of the slow subsystem (2.12). In Section 3, for the slow subsystem without co-existence equilibrium, we established the existence of heteroclinic orbits connecting E0(0, 0) and E2(0, k r2 ), E1( 1 b1 , 0) and E2(0, k r2 ), as detailed in Theorem 3.3. For the slow subsystem with one unique co-existence equilibrium, we obtained that the existence of heteroclinic orbits connecting E0(0, 0) and E∗( r2−r1k r1+b1r2 , 1+b1k r1+b1r2 ), E1( 1 b1 , 0) and E∗( r2−r1k r1+b1r2 , 1+b1k r1+b1r2 ), E2(0, k r2 ) and E∗( r2−r1k r1+b1r2 , 1+b1k r1+b1r2 ), as shown in Theorem 3.5. To prove the existence of homoclinic orbits, we adopted the method of invariant regions and Poincaré-Bendixson theorem. In Section 4, based on the geometric singular perturbation theory and Fredholm theorem, we proved that the diffusive modified Leslie-Gower model with chemotaxis (1.5) admits the different traveling waves for 0 < ε≪ 1. The results are demonstrated in Theorem 4.1. Acknowledgements. The authors wish to thank to the anonymous referees for their valuable comments and helpful suggestions that improved the quality of this article. This work was sup- ported by the Natural Science Foundation of China (Grant No. 12271220) and by the Postgraduate Research and Practice Innovation Program of Jiangsu Province (Grant No. KYCX24-3010). References [1] M. Aziz-Alaoui, M. Okiye; Boundedness and global stability for a predator-prey model with modified Leslie- Gower and Holling-type II schemes, Appl. Math. Lett., 16 (2003) 1069-1075. [2] S. Biswas, L. Bhutia, T. Kar, et al.; patiotemporal analysis of a modified Leslie-Gower model with cross- diffusion and harvesting, Phys. D, 470 (2024) 134381. [3] S. Chen, J. Li, Singular perturbations of generalized Holling type III predator-prey models with two canard points, J. Differential Equations, 371 (2024) 116-150. [4] Y. Du, S. Hsu; A diffusive predator-prey model in heterogeneous environment, J. Differential Equations, 203 (2004) 331-364. [5] Z. Du, J. Li, X. Li; The existence of solitary wave solutions of delayed Camassa-Holm equation via a geometric approach, J. Funct. Anal., 306 (2018) 988-1007. [6] Z. Du, J. Liu, Y. Ren; Traveling pulse solutions of a generalized Keller-Segel system with small cell diffusion via a geometric approach, J. Differential Equations, 270 (2021) 1019-1042. [7] N. Fenichel; Geometric singular perturbation theory for ordinary differential equations, J. Differential Equa- tions, 31 (1979) 53-98. EJDE-2025/38 TRAVELING WAVES OF A DIFFUSIVE LESLIE-GOWER MODEL 19 [8] Q. Guo, H. Cheng; Existence of forced waves and their asymptotics for Leslie-Gower prey-predator model with some shifting environments, Discrete Contin. Dyn. Syst. Ser. B, 29 (2024) 2419-2434. [9] Y. Hua, Z. Du, J. Liu; Dynamics of the epidemiological Predator-Prey system in advective environments, J. Math. Biol., 89 (2024) 28-59. [10] Y. Hua, X. Lin, J. Liu, et al; Dynamics of traveling waves for predator-prey systems with Allee effect and time delay, Electron. J. Differential Equations, 2024 No. 33 (2024) 1-19. [11] J. Huang, M. Lu, C. Xiang, et al; Bifurcations of codimension 4 in a Leslie-type predator-prey model with Allee effects, J. Differential Equations, 414 (2025) 201-241. [12] C. K. R. T. Jones; Geometric singular perturbation theory, Springer, Berlin, 1995. [13] P. Kareiva, G. Odell; Swarms of Predators Exhibit ”Preytaxis” if Individual Predators Use Area-Restricted Search, Am. Nat., 130 (1987) 233-270. [14] C. Kuehn; Multiple time scale dynamics, Appl. Math. Sci., Springer, Switzerland 191 (2015). [15] P. Leslie; Some further notes on the use of matrices in population mathematics, Biometrika, 35 (1948) 213-245. [16] P. Leslie, J. Gower; The properties of a stochastic model for the predator-prey type of interaction between two species, Biometrika, 47 (1960) 219-234. [17] D. Li; Global stability in a multi-dimensional predator-prey system with prey-taxis, Discrete Contin. Dyn. Syst., 41 (2021) 1681-1705. [18] D. Li, N. Tan, H. Qiu; Traveling wave solutions in a modified Leslie-Gower model with diffusion and chemotaxis, Z. Angew. Math. Phys., 75 (2024) 165-188. [19] J. Liu, J. Wu, X. Lin, et al; Traveling wave solutions of a diffusive predator-prey system with Holling II type functional response, Proc. Amer. Math. Soc., 153 (2025) 577-589. [20] A. Lotka; Elements of Physical Biology, Williams and Wilkins Company, Philadelphia, 1925. [21] Q. Qiao; Traveling waves and their spectral instability in volume-filling chemotaxis model, J. Differential Equations, 382 (2024) 77-96. [22] Q. Qiao, X. Zhang; Traveling waves to a chemotaxis-growth model with Allee effect, J. Differential Equations, 416 (2024) 1747-1770. [23] H. Qiu, S. Guo; Steady-states of a Leslie-Gower model with diffusion and advectio, Appl. Math. Comput., 346 (2019) 695-709. [24] H. Qiu, S. Guo; Bifurcation structures of a Leslie-Gower model with diffusion and advection, Appl. Math. Lett., 135 (2023) 108391. [25] J. Roy, S. Dey, B. Kooi, et al.; Prey group defense and hunting cooperation among generalist-predators induce complex dynamics: a mathematical study, J. Math. Biol., 89 (2024) 22-64. [26] H. Safuan, I. Towers, Z. Jovanoski, et al.; On travelling wave solutions of the diffusive Leslie-Gower model, Appl. Math. Comput., 274 (2016) 362-371. [27] Z. Tian, L. Zhang; Traveling Wave Solutions for a Continuous and Discrete Diffusive Modified Leslie-Gower Predator-Prey Model, Qual. Theory Dyn. Syst., 23 (2024) 253-301. [28] V. Volterra; Fluctuations in the Abundance of a Species considered Mathematically, Nature, 47 (1926) 558-560. [29] K. Wang, Z. Du, J. Liu; Traveling pulses of coupled FitzHugh-Nagumo equations with doubly-diffusive effect, J. Differential Equations, 374 (2023) 316-338. [30] Y. Yao, J. He; Bifurcations of codimension three in a Leslie-Gower type predator-prey system with herd behavior and predator harvesting, J. Appl. Anal. Comput., 15 (2025) 734-761. [31] Z. Zhao, H. Hu; Boundedness, stability and pattern formation for a predator-prey model with Sigmoid func- tional response and prey-taxis, Electron. J. Differential Equations, 37 (2023) 1-20. [32] H. Zhang, Y. Xia; Persistence of kink and anti-kink wave solutions for the perturbed double sine-Gordon equation, Appl. Math. Lett., 141 (2023) 108616. [33] W. Zhou, X. Wei; Global stability in a diffusive predator-prey model of Leslie-Gower type, Partial Differ. Equ. Appl. Math., 7 (2023) 100472. Shuna Wang School of Mathematics and Statistics, Jiangsu Normal University, Xuzhou, Jiangsu 221116, China Email address: 2020211141@jsnu.edu.cn Jiang Liu School of Mathematics and Statistics, Jiangsu Normal University, Xuzhou, Jiangsu 221116, China Email address: jiangliu@jsnu.edu.cn Jun Fang School of Mathematics and Statistics, Jiangsu Normal University, Xuzhou, Jiangsu 221116, China Email address: 2020231222@jsnu.edu.cn Xiaojie Lin (corresponding author) School of Mathematics and Statistics, Jiangsu Normal University, Xuzhou, Jiangsu 221116, China Email address: linxiaojie@jsnu.edu.cn 1. Introduction 2. preliminaries 2.1. Geometric singular perturbation theory 2.2. Equilibria of the slow subsystem (2.12) 3. Heteroclinic orbits of the slow subsystem (2.12) 3.1. Existence of heteroclinic orbits for the slow subsystem (2.12) without co-existence equilibrium 3.2. Existence of heteroclinic orbits for the slow subsystem (2.12) with a unique co-existence equilibrium E 4. Traveling wave solutions of system (1.5) 5. Conclusion Acknowledgements References