QUALITATIVE ANALYSIS FOR A DIFFUSIVE PREDATOR-PREY MODEL WITH HUNTING COOPERATION AND HOLLING TYPE III FUNCTIONAL RESPONSE

. The Spatio-temporal pattern induced by self-diﬀusion of a predator-prey model with Holling type III functional response that incorporates the hunting cooperation between predators has been investigated in this paper. For the local model without structure, the stability of non-negative equilibria with or without collaborative hunting in predators is studied. For the Spatio-temporal model, we analyze the eﬀect of hunting cooperation term on diﬀusion-driven Turing instability of the homogeneous positive equilibria. To get an idea about patterns formation near the Turing bifurcation, we derive and give a detailed study of the amplitude equation using the multiple-scale analysis. Our result shows that hunting cooperation plays a crucial role in determining the stability and the Turing bifurcation of the model, which is in sharp contrast to the case without cooperation in hunting. Furthermore, some numerical simulations are illustrated to visualize the complex dynamic behavior of the model.


Introduction
Modeling predator-prey systems have been one of the most important research themes in ecology. It is known that both prey and predator adjust their behaviors to gain the maximal benefits and to increase their biomass for each. Taking into account the fact that, in the real world, predator and prey are not completely static, they may move from one position to another, which indicates that the population density is not homogeneous. For this reason, spatial mathematical modeling is a convenient tool for describing Spatio-temporal population dynamics. An appropriate mathematical structure to explain the spatial aspect of population dynamics is specified by reaction-diffusion systems. Furthermore, reaction-diffusion models were initially applied to describe the ecological pattern formation by [30]. Moreover, Spatio-temporal patterns appear almost everywhere in nature, but their description and understanding still raise important and basic questions. Turing patterns in predator-prey models have been deeply investigated in the literature. To quote a few [1,5,6,14,15,20].
Keywords and phrases: Predator-prey system, Holling type III functional response, cooperative hunting, self-diffusion, amplitude equation, Turing instability.
More recently, there has been increasing interest in modeling predator-prey systems taking into account spatial heterogeneity and hunting cooperation mechanisms to describe the dynamics of the two populations. For example, [28] studied the Spatio-temporal predator-prey model with Holling type II functional response and investigated the effect of hunting cooperation in pattern formation. [7,29,33] studied the spatial predatorprey model with Holling type I functional response and concluded that the model with hunting cooperative may generate Turing instability. More and more recently, [34] studied a Spatio-temporal predator-prey model with Allee effect in prey and hunting cooperation in predator and taking into account the Holling type III functional response, they proved that the model shows both Turing and non-Turing patterns produced by the diffusion. Moreover, [17] studied a Spatio-temporal predator-prey model with hunting cooperation in the specialist predator population with three different types of parametrization of functional responses and found that the patterns evolve only when lower diffusivity of predators than prey has been considered. The purpose of this paper is not only analytically explore the existence and stability of a positive steady state but also to study the impact of self-diffusion as well as hunting cooperation on the formation of spatial patterns. The model reads as follows: where U and V denote the density of prey and predator at time t and position (u, v) respectively. The parameter r 1 > 0 is the per capita intrinsic growth rate of prey, K 1 > 0 is the carrying capacity of the environment for prey, α 1 is a parameter describing the strength of predators cooperation in hunting, h 1 > 0 represents the predator handling time, e 1 is the conversion rate of prey biomass to predator biomass, m 1 is the per capita mortality rate of predators, λ 1 is the successful attack rate in the absence of hunting cooperation, that is the number of preys consumed by a predator by a unit of searching time. While the term (λ 1 + α 1 V ) describes the attack rate in the case when predators cooperate. We consider a Holling type III functional response to represent the saturation effect and assume that it depends on both predators and prey densities φ(U, V ) = (λ 1 + α 1 V ) U 2 1 + h 1 (λ 1 + α 1 V ) U 2 . It is easy to see that hunting cooperation disappears when α 1 = 0. Moreover, ∇ 2 = ∂ 2 ∂u 2 + ∂ 2 ∂v 2 denotes the Laplacian operator in two-dimensional space, δ 1 and δ 2 are the self-diffusion coefficients of the prey and the predator, respectively.
In what follows, we non-dimensionalize the model (1.1) by introducing dimensionless variables: and dimensionless parameters: where the time is described by m 1 t instead of t. We get the following system where Ω is a bounded domain in R 2 with smooth boundary ∂Ω and represents the domain where these two species inhabit and ν is the outward unit normal on ∂Ω. We assume homogeneous Neumann boundary conditions, which means that model (1.4) is self-contained and has no population flux across the boundary ∂Ω, see [18].
The non-spatial model with hunting cooperation has been studied numerically by [2,26], and they concluded that the stable coexisting equilibrium is not destabilized by foraging facilitation.
The paper is organized as follows. In Section 2, we summarize the dynamics of the model without diffusion. Section 3 is devoted to studying the local stability of the uniform steady state of system (1.4). In Section 4, we carry out a nonlinear analysis using multiple-scale analysis to derive the amplitude equations. In Section 5, we perform some numerical simulations to illustrate our analytical results. Finally, we provide a discussion and conclusion in Section 6.

Existence and stability of the uniform equilibria of the model without diffusion
Throughout the remaining part of this section, we present some results concerning the existence of uniform steady states and the local stability of each one of them. These equilibria correspond to spatially homogeneous steady states of the Spatio-temporal model (1.4). Clearly, we can easily show that the reaction-diffusion system (1.4) has trivial equilibrium E 0 = (0, 0) and axial equilibrium E 1 = (K, 0).
On the other hand, by calculating the variational matrix, we conclude that the trivial equilibrium E 0 = (0, 0) is always unstable and axial equilibrium The trivial equilibrium E 0 = (0, 0) is always unstable, which means that the two populations do not disappear simultaneously. On the other hand, the axial equilibrium that the predator population goes extinct, while the prey population reaches its carrying capacity. Otherwise, , which leads us to expect the coexistence of the two populations.
In the sequel, we study the existence and the number of uniform positive equilibria with respect to the model parameters. We separate our discussion into two cases.

Model without hunting cooperation (α 1 = 0)
In the absence of hunting cooperation, system (1.4) has the unique positive steady state E * = (x * , y * ), where The community matrix at the positive steady state E * = (x * , y * ) is given by On the other hand, In this case, the positive steady state is locally asymptotically stable. This situation may occur if, for example, a predator takes a lot of time to handle its prey over a domain where the carrying capacity of the prey is small enough. Furthermore, if K > 2h (2h − 1) √ 1 − h , we obtain tr (J 1 E * ) > 0, that is the positive steady state is unstable. When the parameter K passes through the critical value K 0 = 2h (2h − 1) √ 1 − h , the system around E * = (x * , y * ) enters a Hopf-bifurcation. In fact, the characteristic equation around E * = (x * , y * ) reads as It is clear that if then we have two imaginary eigenvalues and it is easy to obtain This confirms that the transversality condition is satisfied, i.e., a Hopf-bifurcation occurs around the uniform positive steady state. The positive steady state E * = (x * , y * ) is stable for K < K 0 and it is unstable for K > K 0 . The system passes from stability to instability when the parameter K increases and exhibits periodic solutions near the bifurcation value K 0 , which indicates that the parameter K plays an important role in controlling the stability of the unique positive steady state.
In the presence of cooperative behavior in predators, the existence, as well as the number of positive equilibria, are summarized in the following proposition. The variational matrix at the positive steady state E * = (x * , y * ) is given by Where (2.6) (2.7) Then So, if we assume that 0 < h < 1 2 and r > It means that if the predator consumes its prey more rapidly and the growth rate of the prey is large enough, thus the uniform positive steady state E * is locally asymptotically stable if the level of the prey at equilibrium is small.
We note that the non-dimensional parameter K = K 1 e 1 λ 1 m 1 , comprising the dimensional carrying capacity K 1 , the per capita mortality rate of predators m 1 , the attack rate λ 1 , and the conversion efficiency rate e 1 .
So, the parameter K is interpreted as the predator's maximum reproductive number. It is the average number of offspring an individual predator can reproduce over its lifetime when the predator is introduced into a prey population at its carrying capacity. From Section 2.1, when the predators do not cooperate collaboratively, we conclude that predators can survive and coexist with prey when the K values are greater than the amount With cooperative behavior, we show in Proposition 2.1 that model (1.1) may have either one or two positive steady states even if Therefore, cooperative hunting of the predator can lead to the extinction of the population or may promote the long time survival of species in the biological sense. So, it's interesting to study the dynamic behavior of the Spatio-temporal model (1.4) with hunting cooperation in the case when K < 1 √ 1 − h , to find out the situation where the hunting cooperation and the self-diffusion can be beneficial to ensure the coexistence of the population and avoid their extinction.

Local stability of the model with diffusion around the positive steady state
The Spatio-temporal model (1.4) has two trivial equilibria, the origin E 0 = (0, 0) and the axial equilibrium E 1 = (K, 0). The number of possibly uniform positive steady states are determined in Proposition 2.1. Note that a positive steady state E * of the reaction-diffusion system (1.4) without structure is spatially homogeneous, which is constant in space and time. We assume that the positive steady state E * = (x * , y * ) is asymptotically stable for model (1.4) without structure, which means that the spatially homogeneous steady state is asymptotically stable concerning spatially homogeneous small perturbations. Therefore, introducing a small heterogeneous perturbation around the homogeneous steady state may lead to guaranteeing the conditions for Turing instability. Which can be expressed as follows: where 0 < ε 1 , ε 2 1, λ k is the growth rate of perturbations in time t, "i" is the imaginary unit complex number, and k u (resp, k v ) is the wave number in the u (resp, v) direction, with k = (k u , k v ).
In the sequel, we study the local stability of the positive steady state E * for the Spatio-temporal system (1.4).

Stability of the model without hunting cooperation
Substituting (3.1) and (3.2) into system (1.4) and linearizing it around the positive steady state E * = (x * , y * ), we obtain the variational matrix as follows: The characteristic equation associated with E * is given by and Sufficient conditions to ensure that tr (A k ) < 0 are given by Then, the uniform positive steady state E * = (x * , y * ) is locally asymptotically stable, with taking into account that Consequently, we conclude that under some conditions on h and K, the self-diffusion cannot cause the Turing instability in the Spatio-temporal model (1.4) without hunting cooperation among predators. On the other hand, in the absence of cooperative behavior in the predator population, there is no spatial-temporal structure, i.e., the self-diffusion cannot make the species distribution uniform over the domain Ω as the time t increases.
Summarizing the above discussion, we announce the following proposition on the local stability of the Spatiotemporal system (1.4) without hunting cooperation around the uniform positive steady state E * = (x * , y * ).

Stability of the model with hunting cooperation
The variational matrix around the positive steady state E * = (x * , y * ) reads as where a 2 , b 2 , c 2 , and d 2 are given by (2.6) and (2.7). The corresponding characteristic equation is Note that and Therefore, the roots of the characteristic equation (3.6) are given by Sufficient stability conditions in the absence of self-diffusion are given in Section 2.2. Thus, Turing instability occurs only if one of the following conditions is violated: In fact, tr (B k ) < 0 for all k = 0 if we assume that tr (J 2 E * ) < 0. Therefore, the necessary and sufficient conditions for the occurrence of Turing instability is Det(B k ) < 0, for some k = 0. In other words, there exists some k which makes Det(B k ) < 0 then this confirms the occurrence of the Turing instability of the Spatio-temporal system (1.4) around its uniform positive steady state. As equation (3.8) is a quadratic polynomial in k 2 , we note that Thus, the critical condition for the occurrence of Turing instability is From inequality (3.12), we remark that the discriminant of the equation C(k 2 ) = 0 is positive, noted by Therefore, the above analysis gives the following results Lemma 3.2. If E * exist and the following conditions are hold

Then, Turing instability occurs for model (1.4).
If we consider the parameter D 1 as the Turing bifurcation parameter of the Spatio-temporal system (1.4), then, there exist a threshold value D T 1 . Hence, when D 1 = D T 1 , we can determine the Turing region, and the explicit expression of D T 1 is obtained by solving the equation C(k 2 T ) = 0 as follows (3.14) From the condition C (k 2 T ) = 0, we get the following critical wave number k T (3.15) and the interval of wave number for which Tuning instability takes place is (k − , k + ) , where , then a 2 < 0. Thus, from the third necessary condition for the occurrence of Turing instability (see Lem. 3.2), we find that Therefore, from the ecological point of view, we conclude that prey disperses faster than predators (i.e., D 1 > D 2 ). On the other hand, the "activator" is the species promoting growth, and the "inhibitor" is the one repressing growth. It is well known that a necessary condition for the occurrence of Turing instability is that inhibitor has to diffuse faster than the activator [19]. Two situations with consequent qualitatively different reactions may occur: predators spread faster than prey or vice versa [32]. Although a large proportion of the bibliography falls in the first case, the other one is not an uncommon situation and is an essential point of our study. Specifically, to ensure the Turing instability conditions, we find that the prey diffuses faster than the predator. This situation may occur if, for example, the prey population is in a favorable environment which helps it to move faster than the predator species. Or another explanation for this phenomenon is the predators may use other hunting methods. The strategy they may adopt for hunting is to approach the prey by hiding until reaching some distance before they trigger the attack. While other predators prefer waiting in the same place until the passage or the arrival of their prey species. The lookout is a hunting technique that consists of a predator remaining hidden in silence and waiting for prey species to pass within range.
So, according to the above analysis, we announce the following theorem. Proof. Proof directly follows from the previous discussion, hence omitted.
In the above theorem, we have obtained the conditions for self-diffusion Turing instability. So, for a better understanding of this phenomenon, we have presented the plots of Det(B k ) and Re(λ + (k 2 )) (largest real part of the eigenvalue λ + ) for different values of bifurcation parameter D 1 .
In Figure 1 We still do not know which kind of Turing pattern will happen. For this reason, we shall derive the amplitude equations of Turing patterns close to the onset D 1 = D T 1 , which interprets the stability of various forms of Turing patterns of the Spatio-temporal system (1.4).

Amplitude equations
The amplitude equation is an important tool to understand the dynamics of reaction-diffusion systems for parameter values near the Turing bifurcation curve. Furthermore, the stability analysis of the amplitude equations give an idea of various patterns that will evolve near the Turing bifurcation boundary, for example, spot, labyrinth, and a mixture of spots and stripes. Hence, multiple-scale perturbation analysis is useful to derive the amplitude equations to study the dynamics of the active slow modes. For more details about multiple-scale perturbation analysis, we refer to [22,23,35].
where W 0 defines the direction of the eigenmodes in concentration space and A j = A x j , A y j T ,Ā j = Āx j ,Ā y j T , respectively denote the amplitudes associated with modes k j , −k j . Turing patterns are well described by a system of three active resonant pairs of modes (k j , −k j ) , j = 1, 2, 3 two of which make an angle of 2π 3 and | k j |= k T , j = 1, 2, 3. Further, i is the imaginary unit, and r = (u, v) is the spatial vector in two dimensions. Setting X = x − x * and Y = y − y * . Then, to obtain the amplitude equations, we should first write the Taylor expansion of model (1.4) at the equilibrium point E * as follows: where a 2 , b 2 , c 2 , and d 2 are the same as defined in (2.6), (2.7) and Then, system (4.2) can be converted to the following form: where L is the linear operator given by 4) and N (X, Y ) is the nonlinear vector as follows Near the Turing bifurcation threshold D 1 = D T 1 , we perturb the bifurcation parameter D 1 along with X, Y and t according to the following forms with a small perturbation ε as follows: where |ε| 1. Then the Taylor expansion of L with respect to ε is where and H = ∇ 2 0 0 0 (4.9) Remark 4.1. We take the amplitude A j (j = 1, 2, 3), to be a variable that changes slowly with respect to time.
Next, we substitute the equations corresponding to X and Y defined in system (4.6) in the nonlinear vector N (X, Y ), and we get where 11) and g 20 x 1 x 2 + g 11 (x 1 y 2 + x 2 y 1 ) + g 02 y 1 y 2 + 1 6 g 30 x 3 1 + 1 2 g 21 x 2 1 y 1 + Substituting (4.6) into (4.3), then we get the different orders of ε as follows: First order of ε: where (x 1 , y 1 ) t is the linear combination of the eigenvectors corresponding to the eigenvalue zero. The solution of the linear problem (4.13) satisfying the Neumann boundary conditions is given by where Similarly, for the order ε 2 , we obtain In view of Fredhom solvability conditions, see [21], the vector function of the right-hand side of equation (4.15) must be orthogonal to the zero eigenvectors of operator L T * which is the adjoint operator of the operator L T . Then, the zero eigenvectors of L T * are , and c.c represents the conjugate term of exp(−ik j .r). Based on orthogonal conditions, we obtain where F j X , F j Y denote the coefficients corresponding to exp(ik j .r) in F X , F Y respectively. From equation (4.18), we get the following equations: where m 2 = 1 2 f 20 φ 2 + f 11 φ + 1 2 f 02 and n 2 = 1 2 g 20 φ 2 + g 11 φ + 1 2 g 02 .

Dynamical analysis of amplitude equations
In this subsection, we investigate the dynamics of the amplitude equations and determine which kind of patterns will be obtained in our predator-prey model. For this reason, we suppose that each amplitude in (4.31) can be decomposed to where ρ j and θ j represent mode and the corresponding phase angle, respectively. Substituting (4.34) into system of equations (4.31) and separating the real and imaginary parts, we obtain four differential equations as follows where θ = θ 1 + θ 2 + θ 3 and τ 0 , µ, g 1 , g 2 , and χ are defined in (4.32) and (4.33).
The dynamical system (4.35) possesses five kinds of solutions represented in the following proposition Proposition 4.2.
Then, the spatial patterns can be determined by the values of µ, g 1 , g 2 , and χ. From Figure 2, we show that in region (µ 1 , µ 2 ) exists a bistable state of system (1.4). In other words, when the control parameter µ lies in this region, the hexagon patterns H + π and the stationary steady state are all stable. While, the hexagon patterns H − π are unstable. Moreover, the hexagon patterns with θ = 0 (i.e., H 0 ) are always unstable when µ > µ 2 . On the other hand, when µ lies in the region (µ 2 , µ 3 ), the stripe patterns S p are unstable and the hexagon H + π patterns are stable. In region (µ 3 , µ 4 ), the system admits another bistable state between the hexagon patterns H 0 and stripe patterns S p . When the parameter µ increases to the critical point µ 4 , we observe that the hexagon patterns H + π are unstable.

Numerical simulations of pattern formation
In this section, numerical simulations are given to illustrate our analytic results obtained in previous sections. Furthermore, to depict the effect of self-diffusion and hunting cooperation on the dynamics of the Spatiotemporal model (1.4). We consider the 2D spatial space [0, 70] × [0, 70], the domain size has been chosen large enough to give enough room for the species to propagate in space. Zero flux boundary conditions have been chosen for all simulations, i.e., the case of immigration or emigration of either species is not taken into account. Moreover, the numerical simulations were performed by using the finite-difference method with a suitably chosen small time step ∆t = 0.001 and space step size ∆h = 0.5. Also, the initial condition is all taken as a random perturbation around the positive steady state E * . More preciously, we have employed statistically uncorrelated Gaussian white noise perturbation in space, which is mathematically noted in the two-dimensional case as follows where γ 1 and γ 2 are very small real numbers and ε ij , η ij are statistically uncorrelated Gaussian white noise perturbations with zero mean and fixed variance in two-dimensional space. Throughout this section, in every snapshot, the blue color corresponds to the low density, and the red color corresponds to the high density of species.

Pattern selection
In terms of results in Sections 3 and 4 and by using the ode45 method in Matlab, we obtain the following numerical results on the spatial patterns.
According to the above analysis, there should be only hexagon patterns H + π for both prey and predator populations, which confirms that the numerical simulations in Figure 4 correspond to the theoretical results. More precisely, blue spots on the red background shown in Figure 4c-e represent that the prey population lies in isolated circular regions with low density, called "cold spots". Similarly, in Figure 4d-f, red spots on the blue background are obtained and called "hot spots". An interesting feature is that the regions where cold spots occur for the prey population coincide with the hot spots for the predator population, which means areas of low prey density correspond to areas of high predator density. Hence, which can be explained by an inverse correlation in the patterns produced by prey and predator populations when the self-diffusion coefficients are taken to be equal up to a certain level. On the other hand, this phenomenon can be caused by the impact of several ecological factors. For example, the fear effect due to predation risk enhances the survival probability of the prey population, and it can greatly reduce the reproduction of the prey population. Also, the fear of the predator population may influence the physiological situation of the prey species, and it may cause a long-term loss for prey species. For example, in the Greater Yellowstone Ecosystem, wolves (Canis lupus) influence the reproductive physiology of elks (Cervus elaphus) [10]. Moreover, the frightened prey species naturally forage less, for which their growth rate decreases and embraces some survival mechanisms like starvation [9,11]. A higher level of acute risk for predation can create prey species to quit habitats or foraging sites momentarily, returning when the acute risk has passed and the prey species are relatively safe [11].
On the other hand, using another set of parameter values, which can lead to other numerical results that illustrate the stripes-spots patterns, as follows. In Figure 5, we observe the emergence of stripes patterns for both populations, which indicates that the numerical simulations are in agreement with theoretical results.

The effect of varied D 1 on the pattern formation
In this subsection, we show how different values of prey diffusion D 1 can lead to different patterns. To do this, we slightly modify the prey diffusion to obtain a configuration that meets the criteria established in Section 4 for spots, stripes, or mixed patterns.  µ range of µ Figure 6 Type of pattern for both prey and predator 4 0.2477 µ 1 < µ 2 < µ < µ 3 < µ 4 (a), (b) Cold spots and Hot spots 8 1.4943 µ 1 < µ 2 < µ < µ 3 < µ 4 (c), (d) Cold spots and Hot spots 32 8.9815 µ 1 < µ 2 < µ 3 < µ < µ 4 (e), (f) Cold spots and Hot spots From Figure 6a-d, we observe the emergence of spot patterns for both populations. Also, we have µ ∈ (µ 2 , µ 3 ), which confirms that the numerical simulations are in agreement with theoretical results. In the last snapshots (see Fig. 6e-f), we observe the appearance of spot patterns. Therefore, the numerical simulations cannot correspond to the theoretical analysis. This phenomenon cannot be explained by the amplitude equation, because the controlling parameter D 1 is far away from the critical value D T 1 = 3.2059.

The effect of hunting cooperation on Turing space
To explain the impact of cooperative behavior in predators on the dynamics of the Spatio-temporal system (1.4), we varied the intensity of hunting cooperation α and fixed the other parameter values.
In Figure 7a-c, we draw the different curves of the polynomial function Det(B k ) against the square of wave number k 2 for different values of prey diffusion parameter D 1 . In Figure 7d-f, we plot the real part of the eigenvalue Re(λ + (k 2 )) against the square of the wave number k 2 for different values of prey diffusion D 1 . When we increase the intensity of hunting cooperation between predators, we observe that the Turing space becomes large enough, (see Fig. 7). From an ecological point of view, the surface of this region increases when α increases,    meaning that it is difficult to find Turing patterns when the prey diffuses much faster than the predator group.
In other words, we deduce that the main effects of higher hunting cooperation are to accelerate the insurgence of patterns and to increase the instability of the system when the chosen D 1 is quite far from its critical value D T 1 .

The effect of varied α on the pattern formation
To illustrate the impact of the strength of predator-hunting cooperation on the distribution of patterns associated with the Spatio-temporal system (1.4), we varied the value of the parameter α and kept the other parameter values fixed. Recall that α = α 1 λ 1 e 1 m 1 λ 1 is proportional to the parameter describing the strength of the predator's cooperation in hunting.
Therefore, a simple calculation leads to the following results: In Figure 8, we show the evolution of spatial patterns of prey and predator at different values of controlling parameter of hunting cooperation between predators α with t = 500. In Figure 8a-b, we observe the emergence of stripes patterns for both populations with α = 13. Otherwise, we have µ = 0.3989 ∈ (µ 3 , µ 4 ) which confirm the existence of stripes patterns. Therefore, the numerical simulations are compatible with the results in Section 4. Similarly, Figure 8c-d presents the distribution of prey and predator with α = 7, and we observe the appearance of spot-stripes patterns for the prey population and the spot patterns for the predator population. Also, we have µ = 0.3099 ∈ (µ 2 , µ 3 ). Then, in this situation, the numerical simulations are not consistent with the theoretical results for the prey population.
In Figure 8e-f, we decrease the intensity of hunting cooperation between predators. We show that the spatial patterns present the existence of spot patterns, and we have µ = 0.2080 ∈ (µ 2 , µ 3 ), which also confirms that the theoretical results are compatible with numerical analysis for the prey population. From Figure 8, we conclude that predation rate has a significant influence on Turing selection and show our model dynamics exhibit a transition from stripes to spot-stripe and to spot replication. Biologically, it means that the strong hunting cooperation between predators gives some different pattern types (Spots-Stripes and spots). This confirms that the system's instability (1.4) increased with a higher level of hunting cooperation. More precisely, In an area where prey species are very abundant, despite an important hunting collaboration between predators, we notice the appearance of very intense stripes. Thus, ecologically, we can explain this phenomenon as if the prey species were in a favorable place to diffuse faster than the predator species, or that the fear due to the risk of predation imposes to the latter to hide until security is back. On the other hand, when there is a decrease in the intensity of hunting cooperation between predators, weaker spot patterns appear for both populations, which can be explained by the fact that the prey leaves in a safe place. Therefore, this leads to a decrease in the number of predator species, in the reason for searching for another foraging source.

Discussion and conclusion
In recent years, there has been considerable interest in spatial and temporal behaviors of interacting species in ecosystems. From the point of view of biology, it is meaningful to consider several significant factors that affect the qualitative behavior of each ecological system, such as refuge, time delay, the Allee effect, hunting cooperation between predators, etc. So, in the present paper, we have considered a Spatio-temporal predatorprey model with self-diffusion terms and Holling type III functional response incorporating hunting cooperation among predators. In fact, in the real world, there are many species that our proposed model (1.4) can fit. Notably, species that live in groups we mention, for example, the population of Orca and the Blue Whale. Where they hunt in packs and cooperate in achieving a successful hunt. For this specific predator, hunting cooperation is very efficient in catching bigger prey. Thus, to avoid hunting by a group of Orcas, the acceleration of the Blue Whale can be considered very high, see https://www.nationalgeographic.com/animals/article/ orcas-can-kill-blue-whales-the-biggest-animal-on-earth.
More specifically, to determine the qualitative dynamics of our proposed model, we have investigated the analysis step by step. For the model without structure, first, in the absence of hunting cooperation between predators, we find that system (1.4) admits a unique positive steady state if and only if 0 < h < 1 and K > 1 √ 1 − h and it is asymptotically stable if 0 < h < 1 2 (or 1 2 < h < 1 and K < 2h (2h − 1) √ 1 − h ). Ecologically, it means an area where the predator captures and consumes its prey quickly. Also, its prey is associated with a high carrying capacity, which is then sufficient to ensure the coexistence of both species. Second, when predators integrate into the hunting cooperation, the dynamics become more complex. So, from Proposition (2.1), we have observed that the hunting cooperation affects the existence of the positive steady state and it can lead to the occurrence of two positive steady states if αrK > 2, P (x 2 ) > 0, and K < 1 √ 1 − h i.e., P (K) < 0. In general, from an ecological point of view, in an environment where the prey population is abundant, then even if the predators cooperate intensively and handle their prey rapidly, there exists at least one positive steady state.
For the Spatio-temporal model, to find the different types of pattern formation, we will have derived the amplitude equations using the linear stability analysis. Specifically, stability conditions for different patterns formation are derived and condensed in a diagram spanned by control parameters, which allows us to predict which type of patterns can emerge.
In particular, through the analysis and numerical simulations in Section 3, we have obtained the expression of the critical value of prey diffusion D T 1 and found that the homogeneous steady state becomes unstable when the parameter D 1 is greater than the critical value D T 1 . Afterward, in Section 4, we found that this paper establishes the amplitude for the active modes, which determine the stability of the amplitudes. More particularly, we found that when the parameter µ is close to the critical points µ 2 and µ 3 , the numerical results correspond perfectly to the theoretical analysis, as shown in Figure 4. Therefore, To show the impact of prey diffusion D 1 , we have varied the parameter D 1 and obtained that when µ ∈ (µ 2 , µ 3 ), the numerical simulations are in agreement with theoretical results. But, when µ ∈ (µ 3 , µ 4 ), we have observed that the numerical illustrations don't correspond perfectly to the theoretical analysis, as shown in Figure 6e-f. This phenomenon cannot be explained by the amplitude equations, because the control parameter D 1 is far away from the critical value D T 1 . Further, we have illustrated the impact of a varied parameter describing the strength of predators in hunting and found a transition from stripes to the spots-stripes to spots patterns. Even if, we have µ ∈ (µ 2 , µ 3 ) which are not consistent with the theoretical results. Ecologically, it means that strong intensity of hunting cooperation between predators can accelerate the insurgence of patterns and also increase the Turing space (see Fig. 7).
In this paper, we only consider the pattern selection of the two-dimensional predator-prey model with selfdiffusion. It is questionable whether the system with cross-diffusion and adding some time delay represents the delay in forming a group and preparation for the attack. We will discuss this system in detail in future work.