Lionel Roques

Mickaël D Chekroun

On population resilience to external perturbations

Keywords: reaction-diffusion, heterogeneous media, harvesting models, stochastic environments, periodic environments

We study a spatially-explicit harvesting model in periodic or bounded environments. The model is governed by a parabolic equation with a space-dependent nonlinearity of KPP type, and a negative external forcing term. The domain is either the whole space R N , with periodic coefficients, or a bounded domain. Analyzing the stationary states, we define two main types of solutions: the "significant" solutions, which always stay above a certain small threshold value, and the "remnant" solutions, which are always below this value. Using sub-and supersolution methods, the characterization of the first eigenvalue and first eigenfunction of some linear elliptic operators, we obtain existence and nonexistence results, as well as results on the number of stationary solutions. We also characterize the asymptotic behavior of the evolution equation as a function of the forcing term amplitude. In particular, we define critical thresholds on the forcing term below which the population density converges to a significant state, while it converges to a remnant state whenever the forcing term lies above the highest threshold. These bounds were shown to be useful in studying the influence of environmental fragmentation on the long-time behavior of the population density, in terms of the forcing term amplitude. We also present numerical results in the case of stochastic environments.

Introduction

Overexploitation has led to the extinction of many species [START_REF]IUCN Red List of Threatened Species. A Global Species Assessment[END_REF]. Traditionally, models of ordinary differential equations (ODEs) or difference equations have been used to estimate the maximum sustainable yields from populations and perform quantitative analysis of harvesting policies and management strategies [START_REF] Getz | Population Harvesting: Demographic Models of Fish, Forests and Animal Resources[END_REF]. Occulting age or stage structures as well as delay mechanisms, which will not be treated by the present paper, the ODEs models are generally of the type:

dU dt = F (U ) -Y (U ), (1.1)
where U is the population biomass at time t, F (U ) is the growth function and Y (U ) corresponds to the harvest function. In these models, the most commonly used growth function is logistic, with F (U ) = U (µ -νU ) ( [START_REF] Beddington | Harvesting natural populations in a randomly fluctuating environment[END_REF], [START_REF] Murray | Minimum domains for spatial patterns in a class of reaction-diffusion equations[END_REF], [START_REF] Schaefer | Some considerations of population dynamics and economics in relation to the management of the commercial marine fisheries[END_REF]), where µ > 0 is the intrinsic growth rate of the population and ν > 0 models its susceptibility to crowding effects. Different harvesting strategies Y (U ) have been considered in the literature, and are used in practical resource management. A very common one is the constant-yield harvesting strategy, where a constant number of individuals are removed per unit of time: Y (U ) = δ, with δ a positive constant. This harvesting function naturally appears when a quota is set on the harvesters ( [START_REF] Robinson | Towards wildlife management in tropical forests[END_REF], [START_REF] Robinson | Sustainable harvest of neo-tropical mammals[END_REF], [START_REF] Stephens | Sustainable exploitation of social species: a test and comparison of models[END_REF]). Another frequently used harvesting strategy is the proportional harvesting strategy (also called constant-effort harvesting), where a constant proportion of the population is removed. It leads to a harvesting function of the type Y (U ) = δU .

Much less has been done in this field using reaction-diffusion models (but see [START_REF] Kurata | Optimal Spatial Harvesting Strategy and Symmetry-Breaking[END_REF], [START_REF] Neubert | Marine reserves and optimal harvesting[END_REF], [START_REF] Oruganti | Diffusive Logistic equation with constant effort harvesting, I: Steady States[END_REF]). The aim of this paper is to perform an analysis of some harvesting models, within the framework of reaction-diffusion equations.

One of the most celebrated reaction-diffusion model has been introduced by Fisher [START_REF] Fisher | The advance of advantageous genes[END_REF] and Kolmogorov, Petrovsky and Piskunov [START_REF] Kolmogorov | Etude de l'équation de la diffusion avec croissance de la quantité de matière et 29 son application à un problème biologique[END_REF] in 1937. Since then, it has been widely used to model spatial propagation or spreading of biological species into homogeneous environments (see books [START_REF] Murray | Minimum domains for spatial patterns in a class of reaction-diffusion equations[END_REF], [START_REF] Okubo | Diffusion and Ecological Problems -Modern Perspectives[END_REF] and [START_REF] Turchin | Quantitative Analysis of Movement: Measuring and Modeling Population Redistribution in Animals and Plants[END_REF] for review). The corresponding equation is

u t = D∇ 2 u + u(µ -νu), (1.2)
where u = u(t, x) is the population density at time t and space position x, D is the diffusion coefficient, and µ and ν still correspond to the constant intrinsic growth rate and susceptibility to crowding effects. In the 80's, this model has been extended to heterogeneous environments by Shigesada et al. [START_REF] Shigesada | Traveling periodic waves in heterogeneous environments[END_REF]. The corresponding model (which we denote by SKT-model in this paper) is of the type:

u t = D∇ 2 u + u(µ(x) -ν(x)u). (1.3)
The coefficients µ(x) and ν(x) now depend on the space variable x, and can therefore include some effects of environmental heterogeneity. More recently, this model revealed that the heterogeneous character of the environment played an essential role on species persistence, in the sense that for different spatial configurations of the environment, a population can survive or become extinct, depending on the habitat spatial structure ( [START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF], [START_REF] Cantrell | Spatial Ecology via Reaction-Diffusion Equations[END_REF], [START_REF] Roques | Species persistence decreases with habitat fragmentation: an analysis in periodic stochastic environments[END_REF], [START_REF] Shigesada | Biological invasions: theory and practice[END_REF]).

As mentioned above, the combination of a harvesting model with a Fisher-KPP population dynamics model, leading to an equation of the form u t = D∇ 2 u+u(µ-νu)-Y (x, u) has been considered in recent papers, either using a spatially dependent proportional harvesting term Y (x, u) = q(x)u in [START_REF] Neubert | Marine reserves and optimal harvesting[END_REF] [START_REF] Oruganti | Diffusive Logistic equation with constant effort harvesting, I: Steady States[END_REF], or a spatially dependent and time-constant harvesting term Y (x) = h(x) in [START_REF] Kurata | Optimal Spatial Harvesting Strategy and Symmetry-Breaking[END_REF]. In these papers, the models were considered in bounded domains with Dirichlet (lethal) boundary conditions.

Here we study a population dynamics model of the SKT-type, with a spatially dependent harvesting term Y (x, u):

u t = D∇ 2 u + u(µ(x) -ν(x)u) -Y (x, u). (1.4)
We mainly focus on a "quasi-constant-yield" case, where the harvesting term only depends on u for very low population densities (ensuring the nonnegativity of u). We consider two types of domains and boundary conditions. In the first case, the domain is bounded with Neumann (reflective) boundary conditions; this framework is often more realistic for modelling species which cannot cross the domain boundary. In the second case, we consider the model (1.4) in the whole space R N with periodic coefficients. This last situation, though technically more complex, is useful, for instance, for studying spreading phenomena ( [START_REF] Berestycki | Front propagation in periodic excitable media[END_REF], [START_REF] Berestycki | Analysis of the periodically fragmented environment model : II -Biological invasions and pulsating travelling fronts[END_REF]), and for studying the effects of environmental fragmentation, independently of the boundary effects. Lastly, note that the effects of variability in time of the harvesting function will be investigated in a forthcoming publication [START_REF] Chekroun | Spatially-explicit harvesting models. The influence of seasonal variations[END_REF].

In § 2, we define a quasi-constant-yield harvesting reaction-diffusion model. We prove, on a firm mathematical basis, existence and nonexistence results for the equilibrium equations, as well as results on the number of possible stationary states. We also characterize the asymptotic behavior of the solutions of (1.4). In § 3, we illustrate the practical usefulness of the results of § 2, by studying the effects of the amplitude of the harvesting term on the population density in terms of environmental fragmentation. Lastly, in § 4, we give new results for the proportional harvesting case Y (x, u) = q(x)u.

Mathematical analysis of a quasi-constant-yield harvesting reaction-diffusion model

For the sake of readability, the proofs of the results of § 2 are postponed to § 2.5.

Formulation of the model

In this paper, we consider the model:

u t = D∇ 2 u + u(µ(x) -ν(x)u) -δh(x)ρ ε (u), (t, x) ∈ R + × Ω.
(2.5)

The function u = u(t, x) denotes the population density at time t and space position x. The coefficient D, assumed to be positive, denotes the diffusion coefficient. The functions µ(x) and ν(x) respectively stand for the spatially dependent intrinsic growth rate of the population, and for its susceptibility to crowding effects. Two different types of domains Ω are considered: either Ω = R N or Ω is a smooth bounded and connected domain of R N (N ≥ 1). We qualify the first case as the periodic case, and the second one as the bounded case. In the periodic case, we assume that the functions µ(x), ν(x) and h(x) depend on the space variables in a periodic fashion. For that, let L = (L 1 , . . . , L N ) ∈ (0, +∞) N . We recall the definition:

Definition 2.1 A function g is said to be L-periodic if g(x + k) = g(x) for all x = (x 1 , • • • , x N ) ∈ R N and k ∈ L 1 Z × • • • × L N Z.
Thus, in the periodic case, we assume that µ, ν and h are L-periodic. In the bounded case we assume that Neumann boundary conditions hold: ∂u ∂n = 0 on ∂Ω, where n is the outward unit normal to ∂Ω. The period cell C is defined by

C := (0, L 1 ) × • • • × (0, L N ),
in the periodic case, and in the bounded case, we set

C := Ω,
for the sake of simplicity of some forthcoming statements. We furthermore assume that the functions µ and ν satisfy

µ, ν ∈ L ∞ (Ω) and ∃ ν , ν ∈ R s.t. 0 < ν < ν(x) < ν, ∀ x ∈ Ω. (2.6)
Regions with higher values of µ(x) and lower values of ν(x) will be qualified as being more favorable, while, on the other hand, regions with lower µ(x) and higher ν(x) values will be considered as being less favorable or equivalently, more hostile.

The last term in (2.5), δh(x)ρ ε (u), corresponds to a quasi-constant-yield harvesting term. Indeed, the function ρ ε satisfies:

ρ ε ∈ C 1 (R), ρ ′ ε ≥ 0, ρ ε (s)
= 0 for all s ≤ 0, and ρ ε (s) = 1 for all s ≥ ε, (2.7) where ε is a non-negative parameter. With such a harvesting function, the yield is constant in time whenever u ≥ ε, while it depends on the population density when u < ε. In the sequel, the parameter ε is taken to be very small. As we prove in the next sections, there are many situations where the solutions of the model always remain larger than ε. For these reasons, we qualify our model as quasi-constant-yield harvesting SKT-model, the "dominant" regime being the constant-yield one. Note that the function ρ ε ensures the non-negativity of the solutions of (2.5). From a biological point of view, ε can correspond to a threshold below which harvesting is progressively abandoned. Considering constant-yield harvesting functions without this threshold value would be unrealistic since it would lead to harvest on zero-populations.

Finally, we specify that δ ≥ 0 and that h is a function in L ∞ (Ω) such that there exists α > 0 with α ≤ h(x) ≤ 1 for all x ∈ Ω.

(2.8)

We call h the harvesting scalar field, and δ designates by this way the amplitude of this field. Before starting our analysis of this model, we consider the no-harvesting case, i.e. when δ = 0. We recall the main known results in this case. These results will indeed be necessary for the analysis of the quasi-constant-yield harvesting SKT-model.

The no-harvesting case

When δ = 0 in equation (2.5), our model reduces to the SKT-model described by equation (1.3). The behavior of the solutions of this model has been extensively studied in [START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF] and [START_REF] Berestycki | Analysis of the periodically fragmented environment model : II -Biological invasions and pulsating travelling fronts[END_REF].

Results are formulated in terms of first (smallest) eigenvalue λ 1 of the Schrödinger operator L µ defined by

L µ φ := -D∇ 2 -µ(x)I,
with either periodic boundary conditions (on the period cell C) in the periodic case or Neumann boundary conditions in the bounded case. This operator is the linearized one of the full model around the trivial solution. Recall that λ 1 is defined as the unique real number such that there exists a function φ > 0, the first eigenfunction, which satisfies

-D∇ 2 φ -µ(x)φ = λ 1 φ in C, φ > 0 in C, φ ∞ = 1, (2.9)
with either periodic or Neumann boundary conditions, depending on Ω. The function φ is uniquely defined by (2.9) ( [START_REF] Berestycki | Front propagation in periodic excitable media[END_REF]), and belongs to W 2,τ (C) for all 1 ≤ τ < ∞ (see [START_REF] Adams | Sobolev Spaces[END_REF] and [START_REF] Amann | Fixed point equations and nonlinear eigenvalue problems in ordered Banach spaces[END_REF] for further details). We set

φ := min x∈C φ(x).
We recall that a stationary state p of equation (1.3) satisfies the equation,

-D∇ 2 p = p(µ(x) -ν(x)p).
(2.10)

The following result on the stationary states of (2.10) is proved in [START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF].

Theorem 2.2 (i) If λ 1 < 0, the equation (2.10) admits a unique nonnegative, nontrivial and bounded solution, p 0 .

(ii) If λ 1 ≥ 0, the only nonnegative and bounded solution of (2.10) is 0.

Moreover, in the periodic case, the solution p 0 is L-periodic. Throughout this paper, p 0 always denotes the stationary solution given by Theorem 2.2, Part (i).

In order to emphasize that this solution can be "far" from 0 (see Definition 2.5, and the commentary following (2.14)), we give a lower bound for p 0 .

Proposition 2.3 Assume that λ 1 < 0, then p 0 ≥ -λ 1 φ ν in Ω.
The asymptotic behavior of the solutions of (1.3) is also detailed in [START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF]. It is proved that λ 1 < 0 is a necessary and sufficient condition for species persistence, whatever the initial population u 0 is: Theorem 2.4 Let u 0 be an arbitrary bounded and continuous function in Ω such that u 0 ≥ 0, u 0 ≡ 0. Let u(t, x) be the solution of (1.3), with initial datum u(0, x) = u 0 (x).

(i) If λ 1 < 0, then u(t, x) → p 0 (x) in W 2,τ loc (Ω), for all 1 ≤ τ < ∞, as t → +∞ (uniformly in the bounded case).

(ii) If λ 1 ≥ 0, then u(t, x) → 0 uniformly in Ω as t → +∞.

The situation (i) corresponds to persistence, while in the case (ii) the population tends to extinction. In the sequel, unless otherwise specified, we therefore always assume that λ 1 < 0, so that the population survives, at least when there is no harvesting. We are now in position to start our main analysis of steady states and related asymptotic behavior of the solutions of (2.5).

Stationary states analysis

As it is classically demonstrated in finite dimensional dynamical system theory and many problems in the infinite dimensional setting (see e.g. [START_REF] Temam | Infinite-Dimensional Dynamical Systems in Mechanics and Physics[END_REF]), the asymptotic behavior of the solutions of (2.5) is governed in part by the steady states and their relative stability properties. In that respect, we study in this section the positive stationary solutions of (2.5), namely the solutions of

-D∇ 2 p δ = p δ (µ(x) -ν(x)p δ ) -δh(x)ρ ε (p δ ), x ∈ Ω, (2.11)
in the periodic and bounded cases. When needed, we may note (2.11,δ) instead of (2.11). Note that, provided p δ ≥ ε in Ω, p δ is equivalently a solution of the simpler equation

-D∇ 2 p δ = p δ (µ(x) -ν(x)p δ ) -δh(x), x ∈ Ω.
(2.12)

This last equation has been analyzed in the case of Dirichlet boundary conditions in [START_REF] Oruganti | Diffusive Logistic equation with constant effort harvesting, I: Steady States[END_REF], in the particular case of constant coefficients µ and ν.

Because of the type of harvesting function considered here, we are lead to introduce the definition:

Definition 2.5 Set ε 0 := ε φ ≥ ε.
We say that a nonnegative function σ is

remnant whenever max C σ < ε 0 , whereas it is significant if it is a bounded function satisfying min C σ ≥ ε 0 .
Remark 1: The concepts of remnant and significant solutions, as well as the harvesting term δh(x)ρ ε (u), are not classical. In order to clarify these notions, we present in Fig. 1 a short graphical study of the non-spatial model: The abscissae of the points of intersection of the solid and dashed lines correspond, respectively, to remnant (if smaller than ε) and significant (if strictly larger than ε) stationary solutions of (2.13). We observe that the number of significant solutions is: one

dU dt = U (µ -νU ) -δρ ε (U ) =: k(U ), t ∈ R + , (2.13
, if δ < k(ε) (case δ = δ a ); two, if k(ε) ≤ δ < µ 2 /(4ν) (case δ = δ b ); one, if δ = µ 2 /(4ν) (case δ = δ * ); zero, if δ > µ 2 /(4ν). The number of nonzero remnant solutions is zero or more, if δ ≤ k(ε) (depending on the shape of ρ ε ); one or more, if δ > k(ε), since, from (2.7), ρ ′ ε (0) = 0.
We assumed here that ε 0 = ε.

with constant coefficients µ, ν > 0.

Since ε 0 is assumed to be small in our model, the remnant solutions of (2.11) correspond to very low population densities. On the other hand, significant solutions are everywhere above ε 0 . In particular, a constant-yield is ensured in that case. Contrarily to the ODE case, stationary solutions which are neither remnant nor significant may exist, as outlined in the next theorems. However, as we will see while studying the long-time behavior of the solutions of the model (2.5), they are of less importance (see Theorem 2.11 and § 3). The threshold ε 0 is different from ε in general. We had to define remnant and significant functions using ε 0 for technical reasons (see the proof of Theorem 2.10, part (ii), equation (2.31)). Since ε is assumed to be very small, it has no implication on the biological interpretation of our results. Moreover, most of our results still work when ε 0 is replaced by ε.

Let us now start our analysis of equation (2.11). In the sequel, we always assume that

ε 0 < -λ 1 φ 4ν , (2.14)
so that, in particular, from Proposition 2.3, the solution p 0 of (2.10) is significant.

We begin by proving that there exists a threshold δ * such that, if the amplitude δ is below δ * , the equation (2.11) admits significant solutions, while it does not in the other case.

Theorem 2.6 Assume that λ 1 < 0, then there exists δ * ≥ 0 such that (i) if δ ≤ δ * there exists at least a positive significant solution p δ ≤ p 0 of (2.11).

(ii) if δ > δ * , there is no positive significant solution of (2.

11).

Remark 2: There is no positive bounded solution of (2.11) whenever λ 1 ≥ 0.

Under stronger hypotheses, we are able to prove that (2.11) admits at most two significant solutions. In order to state this result, we need some definitions. Let G be the space defined by

G := H 1 (C), (2.15)
in the bounded case, and by

G := H 1 per = ψ ∈ H 1 loc (R N ) such that ψ is L-periodic , (2.16)
in the periodic case. Let us define the standard Rayleigh quotient, for all ψ ∈ G, ψ ≡ 0, and for all σ ∈ L ∞ (C),

R σ (ψ) := C D|∇ψ| 2 -σ(x)ψ 2 C ψ 2 .
(2.17)

According to the Courant-Fischer Theorem (see e.g. [START_REF] Belhachmi | Shape optimization problems for eigenvalues of elliptic operators[END_REF]), the second smallest eigenvalue λ 2 of the operator L µ can be characterized by:

λ 2 = min E k ⊂G,dim(E k )=2 max ψ∈E k , ψ ≡0 R µ (ψ). (2.18)
This characterization is equivalent to the classical one given in [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF].

We are now in position to state the theorem:

Theorem 2.7 Assume that λ 1 < 0 ≤ λ 2 , then, in the bounded case, the equation (2.11) admits at most two significant solutions. In the periodic case, (2.11) admits at most two L ′ -periodic significant solutions for all L ′ ∈ (0, +∞) N . Moreover, under these hypotheses, if two solutions p 1,δ and p 2,δ exist, they are ordered in the sense that, for instance, p 1,δ < p 2,δ in Ω.

Remark 3: Similar methods also allow us to assess a result on the number of solutions of equation (2.12). Indeed, if λ 1 < 0 ≤ λ 2 , then we obtain that (2.12) admits at most two non-negative bounded (and periodic in the periodic case) solutions. If these solutions exist, they are ordered.

In the periodic case, Theorem 2.7 also gives some information on the periodicity of the significant solutions of (2.11), which are actually found to have the same periodicity as the coefficients of the equation (2.11):

Corollary 2.8 Assume that λ 1 < 0 ≤ λ 2 .
Then, in the periodic case, the significant periodic solutions of (2.11) are L-periodic.

The fact that λ 1 < 0 is directly related to the instability of the trivial solution in the SKT-model. The additional condition λ 2 ≥ 0 in this theorem is linked to the existence of a stable manifold or center manifold of the steady state 0 of the SKT-model, in some appropriate functional spaces (see [START_REF] Temam | Infinite-Dimensional Dynamical Systems in Mechanics and Physics[END_REF]). Therefore, the assumptions of Theorem 2.7, and the Krein Rutmann theory, allow us to conclude that under these assumptions, the unstable manifold of 0 is of dimension equal to one or equivalently the stable manifold is of codimension one. Such results on multiplicity of solutions of elliptic nonlinear equations with a source or sink term have been investigated in the past, and are known nowadays as being of Ambrosetti-problem type. These results also involve manifolds of codimension 1 (in the functional space of forcing), and first and second eigenvalues (for the Laplace operator only) (see [START_REF] Nirenberg | Topics in Nonlinear Functional Analysis[END_REF] for a survey of these results).

In any event, Theorem 2.7 relies on the assumption that λ 2 ≥ 0. In the next proposition, we give conditions under which λ 2 may become positive. Proposition 2.9 (i) In the bounded case, if C is a (smooth) domain with

diameter d := max x,y∈C x -y R N , λ 2 (C) ≥ D π d 2 -max C µ. (ii) In the periodic case, λ 2 (C) ≥ D π L d 2 -max C µ, where L d denotes
the length of the longest diagonal of the period cell C.

For instance, when C = [0, 1] × [0, 1], we have d = L d = √ 2; thus, for D = 1 and max C µ = 4, we get λ 2 > 0.9. However, this lower bound is far from being optimal. Indeed, in all our computations of § 3, and under the same hypothesis on C and D, we always had λ 2 > 0, while max C µ = 10.

Sharper lower bounds for λ 2 can be found in [START_REF] Bramble | Bounds in the Neumann problem for second order uniformly elliptic operators[END_REF]; however, these bounds are also more sensitive to the geometry of the domain, and thus less general. They are therefore not detailed here.

We now introduce a result which is important for more applied ecological questions. Indeed, one of the main drawbacks of Theorem 2.6 is that it gives no computable bound for δ * . Obtaining information on the value of δ * is precious for more inclined ecological questions such as the study of the relationships between δ * and the environmental heterogeneities. The next theorem states some computable estimates of δ * .

Let us define

δ 1 := λ 2 1 φ ν(1 + φ) 2 and δ 2 := λ 2 1 4αν . (2.19)
Note that neither δ 1 nor δ 2 depend on δ and ε.

Theorem 2.10 (i) If λ 1 < 0 and δ ≤ δ 1 , then there exists a positive significant (and L-periodic in the periodic case) solution p δ of (2.11) such that

p δ ≥ - λ 1 φ ν(1 + φ)
.

(ii) If λ 1 < 0 and δ > δ 2 , the only possible positive bounded solutions of (2.11) are remnant.

The lower bound of Part (i), for p δ , does not depend on ε. Thus, there is a clear distinction between the remnant and significant solutions. Note that, of course,

δ 1 ≤ δ 2 .
The formulae (2.19) allow numerical evaluations. An important quantity to compute is the size of the gap δ 2 -δ 1 , and its fluctuations in terms of environmental configurations. This question is addressed in § 3 through a numerical study.

Asymptotic behavior

In this section, we prove that the quantity δ * in fact corresponds to a maximum sustainable yield, in the sense that when δ is smaller than δ * , the population density u(t, x) converges to a significant stationary state of (2.5) as t → ∞, whereas when δ is larger than δ * , the population density converges to a stationary state which is not significant. In fact, when δ is larger than the quantity δ 2 defined by (2.19) we even prove that the population converges to a remnant stationary state of (2.5).

We assume here that the harvesting starts on a stabilized population governed by the standard SKT-model with δ = 0. From Theorem 2.4, this means that we study the behavior of the solutions u(t, x) of our model (2.5), starting with the initial datum u(0, x) = p 0 (x). Since we have assumed that λ 1 < 0, it follows from Theorem 2.2, Proposition 2.3 and (2.14) that p 0 is well defined and significant.

Let us describe, with the next theorem, the long-time behavior of the population density.

Theorem 2.11 Let u(t, x) be the solution of (2.5) with initial datum u(0, x) = p 0 (x). Then u is non-increasing in t and, (i) if δ ≤ δ * , u(t, x) → p δ (x) uniformly in Ω as t → +∞, where p δ is the maximal significant solution of (2.11). Moreover p δ is L-periodic in the periodic case;

(ii) if δ > δ * , then the function u(t, •) converges uniformly in Ω to a solution of (2.11) which is not significant;

(iii) if δ > δ 2 , the function u(t, •) converges uniformly in Ω to a remnant solution of (2.11).

Remark 4: If, in addition, we assume that λ 2 ≥ 0, then Theorem 2.7 says that, whenever δ ≤ δ * , the equation (2.5) admits at most two significant stationary states (which are periodic stationary states in the periodic case). In that case, the stationary state p δ selected at large times is the higher one. If we do not assume that λ 2 ≥ 0, this stationary state can still be defined as "the maximal one" that can be constructed by a sub-and super-solution method (see [START_REF] Amann | Supersolution, monotone iteration and stability[END_REF]).

From the above theorem, we observe that, whenever δ ≤ δ * , the solution u(t, x) of (2.5), with initial datum p 0 , remains significant for all times t ≥ 0. This ensures a constant yield in time, and justifies the name of the model.

Similar results could be obtained for a wider class of initial data. Indeed, with similar methods, the convergence of u(t, x) to a significant solution of (2.11) can be obtained whenever δ ≤ δ * for all bounded and continuous initial data u(0, x) which are larger than the smallest significant solution of (2.11). In particular, when u(0, x) is larger than the maximal significant solution of (2.11), u(t, x) converges to this maximal significant solution as t → +∞. A more detailed analysis of the basin of attraction related to the maximal significant solution will be further investigated in the forthcoming paper [START_REF] Chekroun | Spatially-explicit harvesting models. The influence of seasonal variations[END_REF].

Theorem 2.11 shows that the practical determination of δ * is directly linked to the size of the gap δ 2 -δ 1 . As we will see in § 3, this gap (δ 1 , δ 2 ) can be very narrow in certain situations. In those cases, the numerical computation of δ 1 and δ 2 therefore gives a sharp localization of the maximum sustainable quota δ * ∈ [δ 1 , δ 2 ], that can be of non-negligible ecological interest.

Proofs of the results of § 2

Proof of Proposition 2.3: Let φ be defined by (2.9), with the appropriate boundary conditions. Set κ 0 := -λ 1 ν . Then the function κ 0 φ satisfies

-D∇ 2 (κ 0 φ) -µ(x)κ 0 φ + ν(x)(κ 0 φ) 2 = λ 1 κ 0 φ + ν(x)(κ 0 φ) 2 , = κ 0 φ(λ 1 + ν(x)κ 0 φ) ≤ 0.
Thus κ 0 φ is a subsolution of the equation (2.10) satisfied by p 0 . Since for M ∈ R large enough, M is a supersolution of (2.10), it follows from the uniqueness of the positive bounded solution p 0 of (2.10) that p 0 ≥ κ 0 φ ≥ -λ 1 φ ν .

Before proving Theorem 2.6, we begin with the following lemma.

Lemma 2.12 For all δ > 0, if p δ is a nonnegative bounded solution of (2.11), then p δ ≤ p 0 .

Proof of Lemma 2.12: Assume that there exists x 0 ∈ Ω such that p δ (x 0 ) > p 0 (x 0 ). The function p δ satisfies

-D∇ 2 p δ -p δ (µ(x) -ν(x)p δ ) = -δh(x)ρ ε (p δ ) ≤ 0,
thus p δ is a subsolution of the equation (2.10) satisfied by p 0 . Since for M ∈ R large enough, M is a supersolution of (2.10) we can apply a classic iterative method to infer the existence of a solution p ′ 0 of (2.10) (with Neumann boundary conditions in the bounded case since both p δ and M satisfy Neumann boundary conditions) such that p δ ≤ p ′ 0 ≤ M . In particular p ′ 0 (x 0 ) > p 0 (x 0 ), which is in contradiction with the uniqueness of the positive bounded solution of (2.10).

Proof of Theorem 2.6: Let us define δ * := sup{δ ≥ 0, (2.11) admits a significant solution}.

For δ = 0, we know from Proposition 2.3 that p 0 is a significant solution of (2.11). Moreover, for δ large enough, the nonexistence of significant solutions of (2.11) is a direct consequence of the maximum principle (it is also a consequence of the proof of Theorem 2.10, Part (ii)). Thus δ * is well defined and bounded.

Assume that δ * > 0, and let us prove that equation (2.11,δ * ) admits a significant solution. By definition of δ * , there exists a sequence (p δ k ) k∈N of solutions of (2.11,δ k ) with 0 < δ k ≤ δ * and δ k → δ * as k → +∞. Moreover, from Lemma 2.12, ε 0 ≤ p δ k ≤ p 0 for all k ≥ 0. Thus, from standard elliptic estimates and Sobolev injections, the sequence (p δ k ) k∈N converges (up to the extraction of some subsequence) in W 2,τ loc , for all 1 ≤ τ < ∞, to a significant solution p δ * of (2.11,δ * ). Now, let 0 ≤ δ < δ * . Then

-D∇ 2 p δ * -p δ * (µ(x) -ν(x)p δ * ) + δh(x) = (δ -δ * )h(x) < 0,
thus p δ * is a subsolution of (2.11,δ). Since p 0 is a supersolution of (2.11,δ), and p δ * ≤ p 0 , a classical iterative method gives the existence of a significant solution p δ of (2.11,δ) (with Neumann boundary conditions in the bounded case since both p 0 and p δ satisfy Neumann boundary conditions). This concludes the proof of Theorem 2.6.

Proof of Theorem 2.7: As a preliminary, we prove that if two solution exist, then they cannot intersect. Let p 1,δ and p 2,δ be two significant solutions of (2.11). In the bounded case, we assume that p 1,δ and p 2,δ satisfy Neumann boundary conditions. In the periodic case, we assume that there exists L ′ ∈ (0, +∞) N such that p 1,δ and p 2,δ are L ′ -periodic, we then denote the period cell by C ′ . Let us set q δ := p 2,δ -p 1,δ . Then q δ verifies

-D∇ 2 q δ -[µ(x) -ν(x)(p 1,δ + p 2,δ )]q δ = 0, (2.20)
thus, setting ρ(x) := µ(x) -ν(x)(p 1,δ + p 2,δ ), we obtain

-D∇ 2 q δ -ρ(x)q δ = 0, (2.21)
with the same boundary conditions that were satisfied by p 1,δ and p 2,δ .

Let λ 1 and λ 2 be respectively the first and second eigenvalues of the operator L ρ := -D∇ 2 -ρI. Let R σ (φ), be defined by equation (2.17). Since ρ(x) < µ(x) -2νε 0 for all x ∈ Ω, we get

R ρ (ϕ) ≥ R µ (ϕ) + 2νε 0 ,
for all ϕ ∈ G ′ , where G ′ := H 1 (C) in the bounded case and

G ′ := H 1 per = ϕ ∈ H 1 loc (R N ) such that ϕ is L ′ -periodic ,
in the periodic case. Thus, by the classical min-max formula (2.18), it follows that

λ 2 ≥ λ 2 + 2νε 0 > 0. (2.22)
Furthermore, from (2.21), 0 is an eigenvalue of the operator L ρ . Thus, (2.22) implies that λ 1 = 0. As a consequence, q δ is a principal eigenfunction of the operator L ρ . The principal eigenfunction characterization thus implies that q δ has a constant sign. Finally, we get that p 1,δ and p 2,δ do not intersect each other.

Let us now prove that equation (2.11) admits at most two significant solutions. Arguing by contradiction, we assume that there exist three significant (L ′ -periodic in the periodic case, for some L ′ ∈ (0, +∞) N ) solutions p 1,δ , p 2,δ , and p 3,δ of (2.11). From the above result, we may assume, without loss of generality, that p 3,δ > p 2,δ > p 1,δ > ε 0 . Set q 2,1 := p 2,δ -p 1,δ and q 3,2 := p 3,δ -p 2,δ , then these functions satisfy the equations -D∇ 2 q 2,1 -ρ 2,1 (x)q 2,1 = 0, (2.23) and -D∇ 2 q 3,2 -ρ 3,2 (x)q 3,2 = 0, (

with ρ 2,1 := µ(x) -ν(x)(p 1,δ + p 2,δ ) and ρ 3,2 := µ(x) -ν(x)(p 2,δ + p 3,δ ). Moreover, q 2,1 > 0 and q 3,2 > 0. Thus 0 is the first eigenvalue of the operators L ρ 2,1 := -D∇ 2 -ρ 2,1 I and L ρ 3,2 := -D∇ 2 -ρ 3,2 I with either Neumann or L ′ -periodic boundary conditions. From the strong maximum principle (see e.g. [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF]) (together with Hopf's Lemma in the bounded case, and using the L ′ -periodicity of q 3,2 in the periodic case), we obtain the existence of θ > 0 such that q 3,2 > θ. Since the operator L ρ 3,2 is self-adjoint, we have the following formula for its first eigenvalue λ 1 3,2 ,

λ 1 3,2 = min ϕ∈G ′ R ρ 3,2 (ϕ), thus λ 1 3,2 = min ϕ∈G ′ R ρ 2,1 (ϕ) + C ν(p 3,δ -p 1,δ )ϕ 2 C ϕ 2 ≥ min ϕ∈G ′ R ρ 2,1 (ϕ) + νθ, ≥ λ 1 2,1 + νθ,
where λ 1 2,1 is the first eigenvalue of the operator L ρ 2,1 . Since the first eigenvalues of the operators L ρ 2,1 and L ρ 3,2 are both 0, we deduce that 0 ≥ 0 + νθ > 0, hence a contradiction.

Proof of Corollary 2.8: Let p δ be a significant L ′ -periodic solution of (2.11), and let k ∈ N i=1 L i Z. From the L-periodicity of the equation (2.11), p δ (• + k) is also a solution of (2.11). By periodicity of p δ , the functions p δ and p δ (• + k) intersect each other. Thus, from Theorem 2.7, since δ and p δ (• + k) are both L ′ -periodic, p δ ≡ p δ (• + k). Therefore, p δ is a L-periodic function.

Proof of Proposition 2.9: In the bounded case, let C be the convex hull of the set C. It was proved in [START_REF] Payne | An optimal poincare inequality for convex domains[END_REF] This proves part (i) of Proposition 2.9.

In the periodic case, since H 1 per can be seen as a subset of H 1 (C), it follows from (2.18) that,

λ 2 ≥ min E k ⊂H 1 (C),dim(E k )=2 max ψ∈E k , ψ ≡0 R µ (ψ). (2.25)
The period cell C is convex but not smooth enough to assert that the righthand side of (2.25) is equal to the second eigenvalue in the bounded case. Let L d be the longest diagonal of C. Then C is included in a ball B L d of diameter L d . Thus, from formula (2.18), the right-hand side of (2.25) is larger than the second eigenvalue of L µ on B L d . From (i), the conclusion of (ii) follows.

Proof of Theorem 2.10, Part (i): Let λ 1 and φ be defined by (2.9), and let κ be a nonnegative real number such that κ > ε 0 . Then we have

-D∇ 2 (κφ) -κφ(µ(x) -κφν(x)) + δh(x)ρ ε (κφ) ≤ λ 1 κφ + κ 2 φ 2 ν(x) + δ ≤ κφ(λ 1 + κφν(x)) + δ ≤ max τ ∈I {τ (λ 1 + τ ν)} + δ,
(2.26) where I = {κφ(x), x ∈ C}. Setting g(τ ) := τ (λ 1 + τ ν), since φ ∞ = 1, and since g is a convex function, it follows from (2.26) that

-D∇ 2 (κφ) -κφ(µ(x) -κφν(x)) + δh(x)ρ ε (κφ) ≤ max{g(κ), g(κφ)} + δ.
(2.27)

Let us take κ 0 be such that g(κ 0 ) = g(κ 0 φ), namely κ 0 = - (2.28) from the hypothesis on δ of Theorem 2.10, Part (i). Therefore, κ 0 φ is a subsolution of (2.11) with either L-periodic or Neumann boundary conditions. Moreover, if M is a large enough constant, M is a supersolution of (2.11) with L-periodic or Neumann boundary conditions. Thus, it follows from a classical iterative method that there exists a solution p δ of (2.11), with the required boundary conditions, and which satisfies κ 0 φ ≤ p δ ≤ M in Ω. Moreover, in the periodic case, since κ 0 φ and M are L-periodic and since the equation (2.11) is also L-periodic, it follows that p δ is L-periodic. Theorem 2.10, Part (i) is proved.

λ 1 ν(1 + φ) (note that κ 0 φ > ε). We get -D∇ 2 (κ 0 φ) -κ 0 φ(µ(x) -κ 0 φν(x)) + δh(x) ≤ - λ 2 1 φ ν(1 + φ) 2 + δ ≤ 0,
Proof of Theorem 2.10, Part (ii): Assume that λ 1 < 0, δ > δ 2 and that there exists a positive bounded solution p δ of (2.11) which is not remnant, i.e.

there exists x 0 with p δ (x 0 ) ≥ ε 0 .

(2.29)

Since φ is bounded from below away from 0 and p δ is bounded, we can define

γ * = inf {γ > 0, γφ > p δ in Ω} > 0. (2.30)
It follows from the definition of γ * that γ * φ ≥ p δ in Ω, and in particular,

γ * φ(x 0 ) ≥ p δ (x 0 ) ≥ ε 0 . Since φ ∞ =1, we get γ * ≥ ε 0 . Thus, γ * φ ≥ ε 0 φ = ε, (2.31)
which implies ρ ε (γ * φ) = 1. Thus, h(x)ρ ε (γ * φ) ≥ α, and we get,

-D∇ 2 (γ * φ)-γ * φ(µ(x)-γ * φν(x))+δh(x)ρ ε (γ * φ) ≥ γ * φ(λ 1 +γ * φν(x))+δα,
on Ω. Moreover, since γ * φ > 0 and ν ≥ ν, we have γ * φ(λ

1 + γ * φν(x)) ≥ - λ 2 1 4ν . Using the fact that δ > δ 2 = λ 2 1 4αν
, we thus get:

-D∇ 2 (γ * φ)-γ * φ(µ(x)-γ * φν(x))+δh(x)ρ ε (γ * φ) ≥ - λ 2 1 4ν + δα > 0, (2.32)
on Ω. Therefore, γ * φ is a supersolution of (2.11). Set z := γ * φ -p δ . From the definition of γ * , we know that z ≥ 0, and that there exists a sequence (x n ) n∈N in Ω such that z(x n ) → 0 as n → +∞.

In the bounded case, up to the extraction of some subsequence, x n → x ∈ Ω as n → +∞. By continuity, z(x) = 0. Moreover, subtracting (2.11) to (2.32), we get

-D∇ 2 z + [ν(x)(γ * φ + p δ ) + χ(x) -µ(x)]z > 0 in Ω, (2.33)
where the function χ is defined by χ(x) = δh(x) ρε(γ * φ(x))-ρε(p δ (x))

γ * φ(x)-p δ (x) whenever γ * φ(x) -p δ (x) = 0, and χ(x) = ρ ′ ε (p δ (x)) otherwise. Since ρ ε is C 1 , χ is bounded. Thus b(x) := ν(x)(γ * φ + p δ ) + χ(x) -µ(x) is a bounded function.
Using the strong elliptic maximum principle, we deduce from (2.33) that z ≡ 0. Thus γ * φ ≡ p δ is a positive solution of (2.11). It is in contradiction with (2.32).

In the periodic case, we must also consider the situation where the sequence (x n ) n∈N is not bounded. Let (x n ) ∈ C be such that

x n -x n ∈ N i=1 L i Z.
Up to the extraction of some subsequence, we can assume that there exists x ∞ ∈ C such that x n → x ∞ as n → +∞. Set φ n (x) = φ(x + x n ) and p δ,n (x) = p δ (x + x n ). From standard elliptic estimates and Sobolev injections, it follows that (up to the extraction of some subsequence)

p δ,n converge in W 2,τ loc , for all 1 ≤ τ < ∞, to a function p δ,∞ satisfying -∇ 2 (Dp δ,∞ ) -p δ,∞ (µ(x + x ∞ ) -p δ,∞ ν(x + x ∞ )) + δh(x + x ∞ )ρ ε (p δ,∞ ) = 0, in R N , while γ * φ n converges to γ * φ ∞ := γ * φ(• + x ∞ ), and -∇ 2 (Dγ * φ ∞ )-γ * φ ∞ (µ(x+x ∞ )-γ * φ ∞ ν(x+x ∞ ))+δh(x+x ∞ )ρ ε (γ * φ ∞ ) > 0, in R N . Let us set z ∞ (x) := γ * φ ∞ (x)-p δ,∞ (x). Then z ∞ (x) = lim n→+∞ z(x+ x n ), therefore z ∞ ≥ 0 and z ∞ (0) = 0. Moreover, there exists a bounded function b ∞ such that -D∇ 2 z ∞ + b ∞ z ∞ > 0 in R N .
(2.34)

It then follows from the strong maximum principle that z ∞ ≡ 0 and we again obtain a contradiction. Finally, we necessarily have p δ ≤ ε 0 , and the proof of Theorem 2.10, Part (ii) is complete.

Proof of Theorem 2.11, Part (i): Assume that δ ≤ δ * . Let p δ be the unique maximal significant solution defined in the proof of Theorem 2.10, Part (i).

Then, from Lemma 2.12,

p δ (x) ≤ p 0 (x) = u(0, x), ∀ x ∈ Ω, (2.35)
which implies

p δ (x) ≤ u(t, x) in R + × Ω, (2.36)
since p δ is a stationary solution of (2.5). Moreover, since p 0 is a supersolution of (2.11), u is nonincreasing in time t, and standard parabolic estimates imply that u converges in W 2,τ loc (Ω), for all 1 ≤ τ < ∞, to a bounded stationary solution u ∞ of (2.5). Furthermore, from (2.36) we deduce that p δ ≤ u ∞ ≤ p 0 . Since p δ is the maximal positive solution of (2.11), it follows that u ∞ ≡ p δ . Moreover, in the periodic case, since p 0 and the equation (2.5) are L-periodic, u(t, x) is also L-periodic in x. Therefore the convergence is uniform in Ω. Part (i) of Theorem 2.11 is proved.

Proof of Theorem 2.11, Parts (ii) and (iii): Assume that δ > δ * . Since 0 is a stationary solution of (2.5) and u(0, x) = p 0 > 0, we obtain that u(t, x) > 0 in R + × Ω, and again, from standard parabolic estimates, we know that u converges in W 2,τ loc (Ω) (for all 1 ≤ τ < ∞) to a bounded stationary solution u ∞ ≥ 0 of (2.5) as t → +∞. Moreover, in the periodic case, from the L-periodicity of the initial data and of the equation (2.5), we know that u(t, •) and u ∞ are L-periodic. Therefore the convergence is uniform in Ω. It follows from Theorem 2.6, Part (ii) that u ∞ cannot be a significant solution of (2.11). Moreover, if δ > δ 2 , Theorem 2.10, Part (ii) ensures that u ∞ is a remnant solution of (2.11).

Numerical investigation of the effects of environmental fragmentation

We propose here to apply the results of §2, on the estimation of the maximum sustainable yield, to the study of the effects of environmental fragmentation. A theoretical investigation of the relationships between maximum sustainable yield and fragmentation is difficult to achieve (see Remark 5).

To overcome this difficulty, we propose a numerical study in the case of stochastic environments. Firstly, we show that the gap δ 2 -δ 1 , obtained from (2.19) and Theorem 2.10, remains small whatever the degree of fragmentation is. This gap corresponds to the numerical values of the harvesting quota δ for which we do not know whether the population density will converge to a significant or a remnant solution of the stationary equation (2.11).

Secondly, we show that there is a monotone increasing relationship between the maximal sustainable yield δ * and the habitat aggregation.

Remark 5: In a periodic environment, a simple way of changing the degree of fragmentation without changing the relative spatial pattern (favorable area/unfavourable area ratio), is to modify the size of the period cell C. Assume that µ(x) = η x L , for some 1-periodic function η with positive integral, and for some L > 0. This means that the environment consists in square cells of side L. Setting λ 1,L := λ 1 , and φ L := φ, we then have

-D∆φ L -η x L φ L = λ 1,L φ L on [0, L] N . The function ψ L (x) := φ L (Lx) thus satisfies -D∆ψ L -L 2 η(x)ψ L = L 2 λ 1,L ψ L in [0, 1] N , with 1-periodicity.
From the Rayleigh formula we thus obtain:

λ 1,L = min ψ∈H 1 per D L 2 [0,1] N |∇ψ| 2 [0,1] N ψ 2 - [0,1] N ηψ 2 [0,1] N ψ 2 ,
therefore λ 1,L < 0 (since ψ ≡ 1 ∈ H 1 per ), and λ 1,L decreases with L. It implies that δ 2 increases with L. The relationship between δ 1 and L is less clear since, φ L = min C φ L may not always be an increasing function of L.

In order to lessen the boundary effects, and to focus on fragmentation, we place ourselves in the periodic case. For our numerical computations, we assume that the environment is made of two components, favorable and unfavorable regions. It is expressed in the model (2.5) through the coefficient µ(x), which takes two values µ + or µ -, depending on the space variable x. We also assume that

µ + > µ -, ν(x) ≡ 1, h(x) ≡ 1 and D = 1.
Using a stochastic model for landscape generation [START_REF] Roques | Species persistence decreases with habitat fragmentation: an analysis in periodic stochastic environments[END_REF], we built 2000 samples of binary environments, on the 2-dimensional period cell C = [0, 1] 2 , with different degrees of fragmentation. In all these environments, the favorable region, where µ(x) = µ + , occupies 20% of the period cell. The environmental fragmentation is defined as follows. We discretize the cell C into n C = 50 × 50 equal squares C i . The lattice made of the cells C i is equipped with a 4-neighborhood system V (C i ) (see Fig. 2), with toric conditions. On each cell C i , we assume that the function µ either takes the value µ + or µ -, while the number n + = card{i, µ ≡ µ + on C i } is fixed to n C × 20 100 = 500. For each landscape sample ω = (µ(C i )) i=1...n C , we set s(ω) = 1 2

C i C i
C i ⊂C C j ∈V (C i ) 1 1{µ(C j ) = µ(C i )}
, the number of pairs of neighbors (C i , C j ) such that µ takes the same value on C i and C j (1 1{•} is the indicator function). The number s(ω) is directly linked to the environmental fragmentation: a landscape pattern is all the more aggregated as s(ω) is high, and all the more fragmented as s(ω) is small (Fig. 3). Thus, we shall refer to s as the "habitat aggregation index".

Remark 6: There exist several ways of obtaining hypothetical landscape distributions. The commonest are neutral landscape models, originally introduced by Gardner et al. [START_REF] Gardner | Neutral Models for the Analysis of Broad-Scale Landscape Pattern[END_REF]. They can include parameters which regulate the fragmentation [START_REF] Keitt | Spectral representation of neutral landscapes[END_REF]. We preferred to use a stochastic landscape model presented in [START_REF] Roques | Species persistence decreases with habitat fragmentation: an analysis in periodic stochastic environments[END_REF], since it allows an exact control of the favorable and unfavorable surfaces, and is therefore well adapted for analyzing the effects of fragmentation per se. This model is inspired from statistical physics. The number of pairs of similar neighbors s is controlled during the process of landscape generation. This quantity can be measured a posteriori on the landscape samples. Other measures of fragmentation could have been used, such as fractal dimension (see [START_REF] Mandelbrot | The Fractal Geometry of Nature[END_REF]). For a discussion on the different ways of measuring habitat fragmentation in real-world situations the interested reader can refer to [START_REF] Fahrig | Effects of habitat fragmentation on biodiversity[END_REF].

For our computations, we took µ + = 10 and µ -= 0, and we computed the corresponding values of λ i 1 , δ i 1 and δ i 2 on each landscape sample ω i of aggregation index s i , for i = 1 . . . 2000. The eigenvalues λ i 1 were computed with a finite elements method. We fitted the data sets {(s i , δ i 1 )} i=1...2000 and {(s i , δ i 2 )} i=1...2000 using ninth degree polynomials (it is enough to assess if the relations between s and δ 1 , δ 2 tend to be monotonic or not). The resulting fitted curves δ 1,f and δ 2,f are presented in Fig. 4. Under the assumption of normally distributed values of δ 1 and δ 2 for fixed s values, we computed a lower prediction bound (δ 1,lo ) for new observation of δ 1 and an upper prediction bound for δ 2 (δ 2,up ), with a level of certainty of 99%. Thus, given a configuration ω, with a fixed value of s, when δ is smaller than δ 1,lo , we take a 0.5% chance of being above δ 1 , while when δ is larger than δ 2,up , we take a 0.5% chance of being below δ 2 . The small thickness of the intervals (δ 1,lo , δ 2,up ) emphasizes the quality of the relationship between the habitat aggregation index s and the maximum sustainable yield δ * ∈ [δ 1 , δ 2 ]. This also indicates that the criteria of Theorems 2.10 and 2.11 are close to be optimal, at least in some situations.

Furthermore, as we can observe, the values of δ 1 and δ 2 tend to increase as s increases, and thus as the environment aggregates. Since δ * ∈ [δ 1 , δ 2 ], we deduce from the computations presented in Fig. 4 that δ * tends to increase with environmental aggregation.

These tests were performed for particular values of µ + and µ -. However, the thickness of the interval (δ 1 , δ 2 ) can be determined for all values of µ + , 3500 4000 4500 5000

1 1.1 1.2 1.3 s: habitat aggregation index δ 1,f δ 2,f δ 2,up δ 1,lo

Remnant solutions

Significant solutions

Figure 4: Solid lines: δ 1,f and δ 2,f correspond respectively to the data sets {(s i , δ i 1 )} i=1...2000 and {(s i , δ i 2 )} i=1...2000 , fitted with ninth degree polynomials. Dashed lines: δ 1,lo is a lower prediction bound for new observations of δ 1 and δ 2,up an upper prediction bound for new observations of δ 2 , with in both cases a level a certainty of 99%. µ -without further numerical computations, provided that µ + -µ -= 10. Indeed, let us set B := µ + -µ -. For a fixed value of B, let µ 0 (x) be a given L-periodic function in L ∞ (R N ) taking only the two values µ + 0 = B and µ - 0 = 0. Let λ 1,0 be the first eigenvalue of the operator -∇ 2 -µ 0 I on C, with L-periodicity conditions, φ 0 the associated eigenfunction with minimal value φ 0 and

δ 1,0 := λ 2 1,0 φ 0 (1 + φ 0 ) 2 and δ 2,0 := λ 2 1,0 4
.

We have the following proposition:

Proposition 3.1 Assume that µ(x) = µ 0 (x) + µ -, with µ -> λ 1,0 . Let δ 1
and δ 2 be defined by (2.19). Then we have δ

2 -δ 1 = 1 - µ - λ 1,0 2 (δ 2,0 -δ 1,0 ).
This result also indicates that the information on δ * is all the more precise as the growth rate function takes low values. However, the "relative thickness" of the interval (δ 1 , δ 2 ), compared to δ 1 , δ 2 -δ 1 δ 1 , does not depend on µ -, as it can be easily noticed.

Proof of Proposition 3.1: The relation λ 1 [µ(x)] = λ 1,0 -µ -is a direct consequence of the uniqueness of the first eigenvalue λ 1 . We assume that µ -> λ 1,0 , so that λ 1 [µ(x)] < 0. From the uniqueness of the eigenfunction φ associated to λ 1 , φ does not depend on µ -. Therefore, δ 1 and δ 2 satisfy

δ 1 = (λ 1,0 -µ -) 2 φ 0 (1 + φ 0 ) 2 and δ 2 = (λ 1,0 -µ -) 2 4
. The result immediately follows.

A few comments on the proportional harvesting model

In this model, the population density u is governed by the equation

u t = D∇ 2 u + u(µ(x) -ν(x)u) -q(x)u, x ∈ Ω, (4.37)
with L-periodicity of the functions µ(x), ν(x) and q(x) in the periodic case, and with Neumann or Dirichlet boundary conditions in the bounded case. Setting τ (x) := µ(x) -q(x), this model becomes equivalent to the SKT-model (1.3). Hence, many properties of the solutions of this model are described in the existing literature.

In particular the existence, nonexistence and uniqueness results of Theorems 2.2 and 2.4 apply. The condition λ 1 [µ(x) -q(x)] < 0 is therefore necessary and sufficient for species persistence. Furthermore, the theoretical results of [START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF], [START_REF] Cantrell | Spatial Ecology via Reaction-Diffusion Equations[END_REF], [START_REF] Roques | Mathematical analysis of the optimal habitat configurations for species persistence[END_REF], [START_REF] Roques | Species persistence decreases with habitat fragmentation: an analysis in periodic stochastic environments[END_REF] on the effects of habitat arrangement on species persistence are also true for this model. For instance, when the function µ(x) is constant, with µ(x) ≡ µ 1 > 0, and if the domain Ω is convex and symmetric with respect to each axis {x 1 = 0}, ..., {x N = 0}, the next result is a straightforward consequence of the paper [START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF],

Theorem 4.1 (i) In the periodic case, λ 1 [µ 1 -q * k (x)] ≤ λ 1 [µ 1 -q(x)]. (ii) In the bounded Dirichlet case, λ 1 [µ 1 -q * k (x)] ≤ λ 1 [µ 1 -q(x)], (iii) In the bounded Neumann case, if Ω is a rectangle, λ 1 [µ 1 -q ♯ k (x)] ≤ λ 1 [µ 1 -q(x)].
Here q * k denotes the symmetric decreasing Steiner rearrangement of the function q with respect to the variable x k , and q ♯ k denotes the monotone rearrangement of q with respect to x k (see [START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF] and [START_REF] Berestycki | On the monotone rearrangement in cylinders and applications[END_REF] for the definition of these rearrangements). These rearrangements of a function q not only preserve its mean value, but also its distribution function. This means that if, for instance, q corresponds to a "patch" function taking the values q 1 , q 2 and q 3 in some regions A 1 , A 2 and A 3 respectively, with A 1 + A 2 + A 3 = |C|, then the areas of the regions where the rearranged functions q * and q ♯ take the values q 1 , q 2 and q 3 remain equal to A 1 , A 2 and A 3 respectively. Theorem 4.1, combined with Theorem 2.4 say that the spatially rearranged harvesting strategies are better for species survival. This result can be helpful from a resource management point of view. Indeed, the authorities can rearrange the position of the harvested areas in order to improve the chances of population persistence. The result of Theorem 4.1 shows that, in the framework of these models, the creation of a large reserve gives persistence more chances than the creation of several small reserves, and is in accordance with the former results of [START_REF] Kurata | Optimal Spatial Harvesting Strategy and Symmetry-Breaking[END_REF] and [START_REF] Neubert | Marine reserves and optimal harvesting[END_REF] in the Dirichlet case. See Fig. 5 for some illustrations in the bounded case with Dirichlet and Neumann boundary conditions. The initial effort function q(x) takes two values, q + > 0 in the white area and q -= 0 in the shadowed regions, which correspond to reserves. (b): Position of the reserves after a symmetric decreasing Steiner rearrangement along the ∆ 1 and ∆ 2 axes, successively. The rearranged configuration (b) always give more chances of species persistence. In the figures (c) and (d), the boundary Γ is divided into two parts: Γ = Γ 1 ∪ Γ 2 . Γ 1 is represented with a solid line and can correspond to a coast, while Γ 2 is represented with a dashed line, and can correspond to a non-physical limit that the species cannot cross (Neumann boundary conditions). (c): The effort function q(x) again takes two values, q + > 0 in the white area and q -= 0 in the reserves. (d): Position of the reserves after monotone rearrangement along the horizontal and vertical axes, successively. The chances of persistence are better in the rearranged configuration (d).

∆ 1 ∆ 2 Γ (a) ∆ 1 ∆ 2 Γ (b) Γ 1 Γ 2 (c) Γ 1 Γ 2 (d)
This forcing term could be regarded as a "quasi-constant-yield" harvesting, depending only on the population density u when u is below a certain small threshold ε. The introduction of such a threshold ε was necessary for ensuring the nonnegativity of the solutions of our model, and therefore its actuality.

We carried out new mathematical results on the elliptic equation satisfied by the stationary states of the model, and on the associated parabolic equation. Both qualitative and quantitative results were obtained.

From the qualitative point of view, we described the behavior of the model solutions in terms of the harvesting amplitude δ. Two main types of stationary solutions were found: the remnant solutions, always below a small threshold ε 0 and therefore close to 0, and the significant solutions, always above this threshold, thus ensuring a time-constant yield. We discussed the maximum number of significant stationary solutions, which we found equal to 2, under an hypothesis of positivity of the second eigenvalue λ 2 of a linear operator. We further investigated the long-time behavior of the solution of our model, starting from a non-harvested population at equilibrium. We found a critical value δ * of the harvesting term amplitude, below which the population density tends over time to a significant stationary solution, and above which it converges to a stationary solution which is not significant. We also established quantitative formulae for some lower and upper bounds for δ * : δ 1 and δ 2 respectively. The threshold δ 2 has the additional property that, whenever the amplitude δ is above δ 2 , the population density decreases to a remnant stationary solution.

The quantitative aspects of our study mainly consisted in discussing the effect of environmental fragmentation on these thresholds δ 1 and δ 2 , and therefore on the interactions between environmental fragmentation and maximum sustainable yield. Namely, when computing the values of δ 1 and δ 2 on 2000 samples of stochastically obtained patchy environments, with different levels of fragmentation, we found an increasing relationship between these two coefficients and an environmental aggregation index s. This indicates that, for given areas of favorable and unfavorable regions, the harvesting quota that a species can sustain, while ensuring and time-constant yield, is higher when the favorable regions are aggregated.

The reader may note that, in our model, the species mobility was not affected by the environmental heterogeneity. Such a dependence could be modelled by using a more general dispersion term, of the form ∇ • (A(x)∇u), instead of D∇ 2 u, where A(x) stands for the diffusion matrix (see [START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF], [START_REF] Shigesada | Biological invasions: theory and practice[END_REF]). In fact, most of our results still work when the matrix A is of class C 1,α (with α > 0) and uniformly elliptic; i.e. when it exists τ > 0 such that A(x) ≥ τ I N for all x ∈ Ω. Indeed, Theorems 2.2, 2.4, 2.7, 2.10, 2.11 remain true under this more general assumption. However, the effects of environmental heterogeneity may differ, depending on the way A(x) and µ(x) are correlated (see [START_REF] Kinezaki | Spatial dynamics of invasion in sinusoidally varying environments[END_REF]). In the proportional harvesting case, the results of § 4 on the effects of the arrangements of the harvested regions may also not be valid with this dispersion term. However, in situations where A(x) takes low values (slow motion) when q(x) is low ("reserves", see § 4), as underlined in [START_REF] Roques | Mathematical analysis of the optimal habitat configurations for species persistence[END_REF], a simultaneous rearrangement of the functions A(x) and q(x) would lead to lower λ 1 values and therefore to higher chances of species survival.

Figure 1 :

 1 Figure 1: The logistic growth function U → U (µ -νU ) (solid line), and the harvesting function U → δρ ε (U ) for three values of δ (dashed lines).The abscissae of the points of intersection of the solid and dashed lines correspond, respectively, to remnant (if smaller than ε) and significant (if strictly larger than ε) stationary solutions of (2.13). We observe that the number of significant solutions is: one, if δ < k(ε) (case δ = δ a ); two, if k(ε) ≤ δ < µ 2 /(4ν) (case δ = δ b ); one, if δ = µ 2 /(4ν) (case δ = δ * ); zero, if δ > µ 2 /(4ν).The number of nonzero remnant solutions is zero or more, if δ ≤ k(ε) (depending on the shape of ρ ε ); one or more, if δ > k(ε), since, from (2.7), ρ ′ ε (0) = 0. We assumed here that ε 0 = ε.

2 .

 2 that the second Neumann eigenvalue of the Laplace operator -D∇ 2 on C was larger than D π d Since C ⊂ C, we have H 1 (C) ⊂ H 1 ( C). Using formula (2.18), we thus obtain that the second eigenvalue of L µ in the bounded case satisfies λ 2

Figure 2 :

 2 Figure 2: The 4-neighborhood system: an element C i of C and its four neighbors.

Figure 3 :

 3 Figure 3: Some samples of the landscapes used for the computations of δ 1 and δ 2 , with different values of the habitat aggregation index s. The black areas correspond to more favorable environment, where µ(x) = µ + .

Figure 5 :

 5 Figure 5: Examples of applications of Theorem 4.1, Parts (ii) and (iii) to reserves management. In the figures (a) and (b), the boundary Γ of Ω is lethal (Dirichlet boundary conditions). (a):The initial effort function q(x) takes two values, q + > 0 in the white area and q -= 0 in the shadowed regions, which correspond to reserves. (b): Position of the reserves after a symmetric decreasing Steiner rearrangement along the ∆ 1 and ∆ 2 axes, successively. The rearranged configuration (b) always give more chances of species persistence. In the figures (c) and (d), the boundary Γ is divided into two parts: Γ = Γ 1 ∪ Γ 2 . Γ 1 is represented with a solid line and can correspond to a coast, while Γ 2 is represented with a dashed line, and can correspond to a non-physical limit that the species cannot cross (Neumann boundary conditions). (c): The effort function q(x) again takes two values, q + > 0 in the white area and q -= 0 in the reserves. (d): Position of the reserves after monotone rearrangement along the horizontal and vertical axes, successively. The chances of persistence are better in the rearranged configuration (d).

DiscussionWe have proposed a model for the study of populations in heterogeneous environments, for populations submitted to an external negative forcing term.


D Bitar

N Kacem

N Bouhaddi

Nonlinear dynamics of two-dimensional granular periodic lattices

We study wave propagation in a damped driven discrete nonlinear two-dimensional (2D) periodic structure. The considered system is composed of identical spherical particles and subjected to harmonic horizontal base excitation, where its collective nonlinear dynamics is caused by the Hertzian contact law and the applied compressive load. The coupled equations governing the nonlinear vibrations of the considered system have been solved using an analytical-numerical solving procedure, which allows a detailed study of complex model interactions of the periodic structure. The multiple scales method coupled with standing wave decomposition were used to transform this nonlinear system into a set of complex differential equations, which has been numerically solved. Numerical simulations have been performed in order to analyze the stability, the modal interactions and the bifurcation topologies resulting from the collective dynamics of particles in contact, in a 2D configuration.

Introduction

The effect of structural periodicity on wave propagation has been extensively studied in a wide array of fields, such as vibrations of spring-mass systems, arrays of coupled pendulums, electrostatically coupled microbeams, electromagnetic waves in photonic crystals and sound waves in phononic crystals. Several researches were based on studying wave propagation in two-dimensional linear periodic structures. For instance, Leamy [START_REF] Leamy | Exact wave-based bloch analysis procedure for investigating wave propagation in twodimensional periodic lattices[END_REF] et al. detailed an exact wave-based approach for characterizing wave propagation in two-dimensional periodic lattices, while Ruzzene et al. [START_REF] Soranna | Wave beaming effects in two-dimensional cellular structures[END_REF] evaluated the dynamic behavior of 2D cellular structures with the focus on the effect of the dynamics of the propagation of elastic waves within the structure. Moreover, Zou et al. [START_REF] Liang | Control of the elastic wave bandgaps in two-dimensional piezoelectric periodic structures[END_REF] presented a detailed calculations of the dispersion relations of bulk waves propagating in two-dimensional (2D) piezoelectric composite structures and showed that the first bandgaps could be controlled according to need. In the field of structural mechanics, Langley [START_REF] Langley | The response of two-dimensional periodic structures to point harmonic forcing[END_REF] et al. considered a beam grillage constructed from strips of aluminum with bolted joints, while in mechanical vibrations Zhou et al. [START_REF] Ichchou | Multi-scale modelling for two-dimensional periodic structures using a combined mode/wave based approach[END_REF] extended the study of wave finite element method to the two-dimensional periodic beam grillage to calculate the wave dispersion characteristics.

Wave propagation in nonlinear periodic media is spare and as most of the structures are inherently nonlinear it is important to investigate the effect of such nonlinearities on the wave propagation characteristics. Nonlinear periodic structures may exhibit a variety of rich dynamic properties such as multi-stability, modal interactions, bifurcation topology transfer, existence of highly stable localized solutions, solitary waves and variation in wave speeds, propagation direction related to wave amplitude and nonlinearity and so on. For instance, in the field of acoustics, Narisetti et al. [START_REF] Leamy | A perturbation approach for analyzing dispersion and group velocities in two-dimensional nonlinear periodic lattices[END_REF] developed a perturbation approach for two-dimensional mono-atomic lattice of identical masses connected by nonlinear spring stiffness of type Duffing between two adjacent masses to evaluate the influence of nonlinearities on the location of band gaps, group velocity magnitudes and the direction of energy propagation. In the field of photonics and phononics, Maldovan et al. [START_REF] Maldovan | Simultaneous localization of photons and phonons in two-dimensional periodic structures[END_REF] demonstrated theoretically the simultaneous localization of photons and phonons in the same spatial region by introducing lattice defects in a periodic array of dielectric/elastic material that exhibits gaps for both electromagnetic and elastic waves. In Biomedical, Sheldon et al. [START_REF] Sheldon | Spatial heterogeneity enhances and modulates excitability in a mathematical model of the myometrium[END_REF] modeled the muscular layer of the uterus as a 2D periodic structures of cells coupled through resistors that represents gaps junctions, where they demonstrated that the spatial heterogeneity enhances and modulates excitability in a mathematical model of the muscular layer of the uterus. Moreover, in optics, Kantner et al. [START_REF] Yanchuk | Delay-induced patterns in a two-dimensional lattice of coupled oscillators[END_REF] showed that arbitrary stable spatiotemporal periodic patterns can be created in 2D lattices of coupled oscillators with inhomogeneous coupling delay.

Hertzian chains of spherical particles represent a famous example of nonlinear discrete structure, where the nonlinearity is caused by the contact law and the applied compressive load. The dynamic behavior onedimensional (1D) chains of spheres have been widely studied [START_REF] Nesterenko | Propagation of nonlinear compression pulses in granular media[END_REF]. Heterogeneous, 1D strongly compressed spheres chain has been used to generate solitary wave [START_REF] Nesterenko | Dynamics of Heterogeneous Materials, Dynamics of Heterogeneous Materials[END_REF][START_REF] Yang | Interaction of highly nonlinear solitary waves with linear elastic media[END_REF]. It was observed that the wave speed changes with increase in amplitude at certain frequency leading to small changes in the band edges [START_REF] Boechler | Analytical and experimental analysis of bandgaps in nonlinear on dimensional periodic structures[END_REF]. In addition, Theocharis et al. [START_REF] Theocharis | Nonlinear Periodic Phononic Structures and Granular Crystals[END_REF] describe the dynamic behavior of nonlinear periodic phononic structures, along with how such structures can be utilized to affect the propagation of mechanical waves. However, the study of 2D granular lattices has not been widely studied. Recently, Leonard et al. [START_REF] Leonard | Directional wave propagation in a highly nonlinear square packing of spheres[END_REF] studied the stress wave properties in a basic 2D granular crystals. Aware about the richness and the complexity of the collective dynamics occurring between localized modes in periodic nonlinear arrays [START_REF] Bitar | Collective dynamics of periodic nonlinear oscillators under simultaneous parametric and external excitations[END_REF], we study the nonlinear dynamics of 2D spherical particles in contact under compression.

In this context, we investigate wave propagation in a discrete nonlinear two-dimensional damped granular periodic lattice subject to an external harmonic excitation with dynamics governed by a Hertzian contact model. The coupled equations governing the nonlinear vibrations of the considered system have been solved using an analytical-numerical solving procedure, which allows a detailed study of complex model interactions of the periodic structure. The multiple scales method coupled with standing wave decomposition were used to transform this nonlinear system into a set of complex algebraic equations, which has been numerically solved. Numerical simulations have been performed in order to analyze the stability, the modal interactions and the bifurcation topologies resulting from the collective dynamics of particles in contact, in a 2D configuration.

Equivalent Mass-Spring system

Consider a 2D periodic structure of identical unit cells subjected to harmonic horizontal base excitation. Each unit cell consist of spherical particle coupled with two uniform beams and subjected to compressive forces applied to all ends (Figures 1(a) 1(b))). Hence, an equivalent spring-mass system can be constructed as shown in Figure 1(c). The displacement of the i, j th sphere from its equilibrium position in the initially compressed chain is defined by u i,j . The periodic system is subjected to harmonic horizontal base excitation ÿ(τ ). For dynamical displacements with small amplitudes relative to those due to static load, where a power series expansion was calculated in order to obtain the following i, j th particle equivalent equation of motion. m eq üi,j + k eq u + c x (-ui-1,j + 2 ui,j -ui+1,j ) + c y (-ui,j-1 + 2 ui,j -ui,j+1 ) where m eq is the effective mass of the two-beams-sphere system, k eq is the equivalent structural beams linear stiffness, c x = ξ x m eq ω 0 and c y = ξ y m eq ω 0 are the damping coefficients in the x and y directions respectively. ω 2 0 = keq+2kx+2ky meq is the natural frequency where k x and k y represent the linear coupling spring constants. In addition, k 2x , k 2y and k 3x , k 3y are respectively the coupling quadratic and cubic parameters in the x and y directions and Y 0 = m eq Ω 2 y 0 the harmonic based excitation amplitude (y(τ ) = y 0 cos(Ωτ )).

-k x (u i-1,j -2u i,j + u i+1,j ) -k y (u i,j+1 -2u i,j + u i,j+1 ) -k 2x [(u i-1,j -u i,j ) 2 -(u i,j -u i+1,j ) 2 ] -k 2y [(u i,j-1 -u i,j ) 2 -(u i,j -u i,j+1 ) 2 ] +k 3x [(u i-1,j -u i,j ) 3 -(u i,j -u i+1,j ) 3 ] + k 3y [(u i,j-1 -u i,j ) 3 -(u i,j -u i,j+1 ) 3 ] = -m eq ÿ = Y 0 cos(Ωτ ) (1)
For boundary conditions, we suppose that we have fixed identical systems in both directions, where we define 2(N + M + 1) extra variables and set them to zero as: u 0,j = u N +1,j = u i,0 = u 0,M +1 = 0.

Normalized equations

It proves convenient to define the following scaled dimensionless variables:

τ = ω 0 t x i,j = u i,j u D (2)
where u D = Y 0 2(cx+cy)ω 0 is the dynamic displacement of the associated linear system while neglecting the linear coupling with Q = meqω 0 2(cx+cy) . By replacing these variable into the equations of motion and dividing by Y 0 ω 0 meq 2(cx+cy) , we obtain the following nondimensional system: ẍi,j + x i + ξ x (-ẋi-1,j + 2 ẋi,j -ẋi+1,j ) + ξ y (-ẋi,j-1 + 2 ẋi,j -ẋi,j+1 )

-R x (x i-1,j -2x i,j + x i+1,j ) -R y (x i,j-1 -2x i,j + x i,j+1 ) -β 2x [(x i-1,j -x i,j ) 2 -(x i,j -x i+1,j ) 2 ] -β 2y [(x i,j-1 -x i,j ) 2 -(x i,j -x i,j+1 ) 2 ] +β 3x [(x i-1,j -x i,j ) 3 -(x i,j -x i+1,j ) 3 ] + β 3y [(x i,j-1 -x i,j ) 3 -(x i,j -x i,j+1 ) 3 ] = -2(ξ x + ξ y ) cos( Ω ω 0 t) (3)
The parameters appearing in Equation ( 3) are:

ξ x = cx meqω 0 , ξ y = cy meqω 0 , R x = kx keq+2kx+2ky , R x = kx keq+2kx+2ky , β 2x = k 2x Y 0 2(cx+cy)ω 3 0 meq , β 2y = k 2y Y 0 2(cx+cy)ω 3 0 meq , β 3x = k 3x Y 2 0 4(cx+cy)ω 4 0 meq and β 3y = k 3y Y 2 0 4(cx+cy)ω 4 0 meq (4)

Linear study

Equation (3) describes a system of Duffing oscillators linearly coupled by springs k x and k y in the x and y directions respectively and subjected to harmonic horizontal base excitation. To determine the natural frequencies and their eigenvectors, Equation (3) can be written in matrix form as

Ẍ + C Ẋ + K L X Linear part +F N L (X) Nonlinear part = F (t), (5)
where X = [x 1,1 , x 1,2 , . . . , x 1,M , . . . x i,1 , x i,2 , . . . , x i,M , . . . ] is the N×M-dimensional variable vector, C the damping matrix, F (t) = -2(ξ x + ξ y ) cos( Ω ω 0 t)[1, . . . , 1] T the excitation vector, K N L (X) the nonlinear stiffness vector and

K L =         K 1 -K 2 -K 2 K 1 -K 2 0 . . . . . . . . . 0 -K 2 K 1 -K 2 -K 2 K 1         (6)
where

K 1 =          1 + 2R x + 2R y -R y -R y 1 + 2R x + 2R y -R y 0 . . . . . . . . . 0 -R y 1 + 2R x + 2R y -R y -R y 1 + 2R x + 2R y          (7)
and K 2 = R x I M . The eigenvalues ω i,j are represented in the following form for few coupled oscillators in contact:

N = 1, M = 1 : ω 1,1 = 1 + 2R x + 2R y N = 2, M = 1 : ω 1,1 = 1 + R x + 2R y ω 2,1 = 1 + 3R x + 2R y (8)
We may express all normal frequencies relative to the same nondimensional reference frequency which is 1, so that

ω i,j = (1 + λ i R x + λ j R y ) (i = 1, . . . , N and j = 1, . . . , M ), (9)
We suppose that k x << k eq + 2k x + 2k y and k y << k eq + 2k x + 2k y . Consequently, R x and R y << 1 and:

ω i,j ≈ 1 + 1 2 λ i R x + 1 2 λ j R y (i = 1, . . . , N and j = 1, . . . , M ). (10)
This assumption leads to the creation of linear closed modes which permits to study the effect of the mode localization on the collective dynamics. In the following section, we will proceed to solve the normalized equivalent differential System (3) using a perturbation technique. Before doing so, we should note that the quadratic nonlinearities are much smaller than the Duffing coupling parameters, thereby their nonlinear effects are considered negligible.

Analytical-numerical approach

In this paper, we use the multiple scales method as an approximate analytical solution to solve nonlinear differential systems. The main advantages of the present approach are the capacity of handling weakly coupled nonlinear systems, which permits to visualize all physical responses branches and their properties in terms of modal interactions and bifurcation topology transfer.

We introduce the parameters listed in Equation ( 11) and set the external frequency an amount εω 0 Ω D away from the resonant frequency, whereby they contribute to the equations of amplitude.

ξ x = 1 2 εη x ξ x = 1 2 εη x R x = 1 2 εΓ x R y = 1 2 εΓ y and 2(ξ x + ξ y ) = ε 3/2 g (11)
We calculate the responses of the spheres lattice subjected to harmonic horizontal excitation using secular perturbation theory, by expanding x i,j (t) as a sum of standing-wave modes with slowly varying amplitude [START_REF] Lifshitz | Response of parametrically driven nonlinear coupled oscillators with application to micromechanical and nanomechanical resonator arrays[END_REF].

x i,j (t) = ε 1/2 N r=1 M p=1 (A r,p (T ) sin(nq r ) sin(dq p )e iωr,pt + c.c.)

+ε 3/2 x (1) i,j (t) + • • • , i = 1, ..., N and j = 1, ..., M (12)
where T = εt is a slow time variable, that authorizes the complex amplitude A r,p (T ) to vary slowly in time.

Since we proposed fixed boundary conditions, the possible wave components q r , d p can be given as

q r = rπ N +1 , r = 1, ..., N d p = pπ M +1 , p = 1, ..., M (13)
After replacing the proposed solution [START_REF] Yang | Interaction of highly nonlinear solitary waves with linear elastic media[END_REF] into the equation of motion, we can get at the order of ε 3 2 , N × M equations of the form:

ẍ(1) i,j + x (1)
i,j = r p (r, p th secular term)e iωr,pt + other terms [START_REF] Theocharis | Nonlinear Periodic Phononic Structures and Granular Crystals[END_REF] We must eliminate the N × M secular terms so that x i,j (t) remains finite which allows us to determine the equations for the slowly varying amplitudes A r,p (T ). To extract the equations for these amplitudes, we make use of the orthogonality of the modes, by multiplying the r, p th secular term by sin(nq r ) sin(dq p ) and summing over n and m. We also express all normal frequencies relative to 1, so that:

ω r,p = 1 + εΩ r,p (15)
We find that the equation of the r, p th amplitude A r,p (T ) is given by:

2iω r,p dAr,p dT + 2iω r,p (η x sin 2 ( qr 2 ) + η y sin 2 ( dp 2 ))A r,p + 2(Γ x sin 2 ( qr 2 ) + Γ y sin 2 ( dp 2 ))A r,p +3β 3x sin( qr 2 ) i,k,l j,s,o sin( q i 2 )sin( q k 2 )sin( q l 2 )A i,j A k,s A * l,o ∆ (1)
ikl;r ∆

jso;p e i(Ω i,j +Ω k,s -Ω l,o -Ωr,p)T +3β 3y sin(

dp 2 ) i,k,l j,s,o sin( d p 2 )sin( d r 2 )sin( d s 2 )A i,j A k,s A * l,o ∆ (2)
ikl;r ∆

jso;p e i(Ω i,j +Ω k,s -Ω l,o -Ωr,p)T

In solving Equations ( 16), we write the periodic steady state solution in the following form

A r,p = a r,p e i[Ω D -Ωr,p]T (17)
Substituting Equation ( 17) into ( 16), we obtain the required Equation ( 5) for the complex amplitudes a r,p .

-2ω r,p (Ω D -Ω r,p )a r,p + 2iω r,p (η x sin 2 ( qr 2 ) + η y sin 2 ( dp 2 ))a r,p + 2(Γ x sin 2 ( qr 2 ) + Γ y sin 2 ( dp 2 ))a r,p +3β 3x sin( qr 2 ) i,k,l j,s,o sin( q i 2 )sin( q k 2 )sin( q l 2 )a i,j a k,s a * l,o ∆ (1)
ikl;r ∆

(2) jso;p

+3β 3y sin( dp 2 ) i,k,l j,s,o sin( d p 2 )sin( d r 2 )sin( d s 2 )a i,j a k,s a * l,o ∆ (2) ikl;r ∆ (1) jso;p = 2 (N +1)(M +1) g n m sin (nq r ) sin (md p ) (18)
These complex algebraic equations which define the time independent mode amplitudes are the main result of the perturbation analysis applied to the normalized equivalent differential system. In order to study the collective dynamics of the weakly coupled 2D periodic structure, we start by solving the algebraic system for a set of design parameters listed in Table 1 which satisfies the modes localization assumption given by Equation [START_REF] Nesterenko | Propagation of nonlinear compression pulses in granular media[END_REF]. Table 1: Design parameters for the corresponding periodic structure depicted in Figure 1(c)

Results and discussions

Before starting our investigations we should note the following statements. First, the second member of the Equation ( 18) is proportional to the sum of standing waves, which is null for all i,j modes where at least one position is even. For a single direction (N = 1, or M = 1), the general periodic structure is reduced to an array of coupled oscillators, where a detailed study has been carried out [START_REF] Bitar | Collective dynamics of periodic nonlinear oscillators under simultaneous parametric and external excitations[END_REF].

Particularly, for two coupled complex algebraic equations, we used Mathematica to solve them numerically, where the stability has been performed based on the eigenvalues of the Jacobian matrix of the differential System ( 16) for each point. Figure 2 displays the intensity responses of two weakly coupled oscillators (N = 1, M = 2), as a function of the detuning parameter Ω D for the first design parameters listed in Table 1. Black, orange and gray curves indicate the Single Mode (SM), Double Mode (DM) and the unstable solution branches. SM solutions are generated by the null trivial solution a 1,2 = 0, where a 1,1 corresponds to a single forced Duffing oscillator. Remarkably, there are additional branches labeled DM, involving the excitation of both modes collectively. With two weakly coupled oscillators we can obtain up to four stable solutions for a given frequency. Hence, for 2D periodic structure we expect additional complex phenomena not obtained for one dimensional array.

Therefore, we study the case of 2 × 2 weakly coupled oscillators, for several design parameters. Starting by the case where we consider the same linear and nonlinear coupling parameters in both directions. Figure 4 shows the squares of the amplitudes of the different modes as a function of the detuning parameters for the second design parameters in Table 1. DM i,j represent the Double Mode solution branches generated by two different modes and QM is the Quadruple Mode involved by the excitation of all modes collectively.

Note that the modes localization have interesting effects on the response intensities with additional multimodal solutions that share the same frequency range with bifurcation topology transfer. In addition, as we choose identical coupling parameters, a 1,2 and a 2,1 become symmetric and interact in the same manner with a 1,1 . Moreover, the Triple Mode (TM) solution branches are missing in this case. In fact, as by replacing one of the non-excited modes a 1,2 , a 2,1 or a 2,2 by 0 into any of their equations, we obtain a null trivial solution, which corresponds to the case of a DM. For the fourth set of design parameters, the coupling in the y direction is highly decreased. The responses intensities plotted in Figure 5 show that the small solution branches join together in one DM 1,1;1,2 branch heading in the direction of the non-resonant branch of |A 1,1 | 2 , While the response intensity of |A 2,1 | 2 takes an elliptical form.

Conclusion

To conclude, a 2D periodic nonlinear structure composed of coupled two-beams-Mass structures under harmonic horizontal base excitation has been considered. The equations of motions of the equivalent Mass-Spring system has been normalized and solved using an analytical-numerical approach based on secular perturbation theory coupled with standing waves decomposition. Intensity responses for few coupled oscillators has been calculated for several design parameters that satisfies the modes localization condition. Complex nonlinear phenomena with high number of stable multimodal solution branches has been observed. Considering different linear and nonlinear coupling parameters in the different directions, the intensity responses topology can be modified with additional multimodal solution branches.

Figure 1 :

 1 Figure 1: (a) 2D periodic nonlinear lattice of spheres in contact under compression and subjected to horizontal base excitation. (b) A view of the system from the top. (c) An equivalent Mass-Spring system.

  r ) sin (md p )e i(Ω D -Ωr,p)T

  -3 8.10 -3 16.10 -3 16.10 -3 29.10 -3 29.10 -3 16.10 -3 2 8.10 -3 8.10 -3 16.10 -3 16.10 -3 15.10 -3 15.10 -3 16.10 -3 3 8.10 -3 4.10 -3 16.10 -3 10.10 -3 5.10 -3 19.10 -3 12.10 -3 4 8.10 -3 18.10 -4 16.10 -3 7.10 -3 2.10 -3 20.10 -3 98.10 -4

Figure 2 :

 2 Figure 2: Response intensities of two weakly coupled oscillators under harmonic base external excitation as a function of the detuning parameter Ω D , for the first set of design parameters listed in table 1.

Figure 3 :

 3 Figure 3: Response intensities of 2 × 2 weakly coupled oscillators under harmonic base external excitation as a function of the detuning parameter Ω D , for the second set of design parameters listed in table 1.

Figure 4 :

 4 Figure 4: Response intensities of 2 × 2 weakly coupled oscillators under harmonic base external excitation as a function of the detuning parameter Ω D , for the third set of design parameters listed in table 1.

Figure 5 :

 5 Figure 5: Response intensities of 2 × 2 weakly coupled oscillators under harmonic base external excitation as a function of the detuning parameter Ω D , for the fourth set of design parameters listed in table 1.

Acknowledgements

This project has been performed in cooperation with the Labex ACTION program (contract ANR-11-LABX-01-01).