Guillaume Cantin
email: guillaume.cantin@univ-nantes.fr

On the robustness of discontinuous patterns in degenerate reaction-diffusion systems with perturbed hysteresis

Keywords:

In this paper, we study the dynamics of a morphogenesis model which reproduces the formation of discontinuous patterns observed in numerous biological or ecological systems. This model is determined by a degenerate reaction-diffusion system with hysteresis in the non-diffusive equation. We focus on the robustness of discontinuous patterns under the action of a perturbation of the hysteresis process. We analyze the bifurcations of homogeneous stationary solutions to this nonlinear model and prove that the trivial solution is the only one to resist to a perturbation of strong intensity. We then prove an original result on the structural transformation of discontinuous patterns, which are seen to react by acquiring a supplementary discontinuity jump under the effect of a perturbation of small intensity. We interpret the supplementary jump for ecological systems as the possible emergence of an invasive ecosystem. Numerical simulations and animations are finally provided to guide intuition on the complex morphogenesis process of this dynamical system.

Introduction

The formation of patterns observed in numerous real-world biological and ecological systems, also known as morphogenesis, has become a very important research topic in mathematical modeling and in the study of nonlinear dynamical systems. Turing is now recognized as the first mathematician to propose a mechanistic model, determined by a system of reaction-diffusion equations, for reproducing morphogenesis in an activator-inhibitor chemical system [START_REF] Turing | The chemical basis of morphogenesis[END_REF]; the forms produced by his model are now well-known as Turing patterns. Nowadays, morphogenesis is a vibrant topic that continues to produce an active literature, in which reaction-diffusion systems play a central role. In these reactiondiffusion systems, Turing patterns have been intensely studied. The literature on this subject is so wide that it is impossible to cite every relevant paper. One can however refer to [START_REF] Li | Bifurcation analysis of an extended Klausmeier-Gray-Scott model with infiltration delay[END_REF][START_REF] Li | Interactions of time delay and spatial diffusion induce the periodic oscillation of the vegetation system[END_REF][START_REF] Liang | Nonlocal interactions between vegetation induce spatial patterning[END_REF][START_REF] Sun | Dynamic analysis of a plant-water model with spatial diffusion[END_REF], where Turing patterns are proved to emerge in vegetation models, or to [START_REF] Mondal | Emergence of Turing patterns and dynamic visualization in excitable neuron model[END_REF], where they are observed in a neuron model. The case of reaction-diffusion systems with a discrete diffusion operator also reveals the formation of remarkable patterns, as proved in [START_REF] Kim | Fully nonlinear Cucker-Smale model for pattern formation and damped oscillation control[END_REF]. Other related works on this subject can be found in [START_REF] Lacitignola | Turing pattern formation on the sphere for a morphochemical reaction-diffusion model for electrodeposition[END_REF][START_REF] Li | Hopf bifurcation and Turing instability in the reaction-diffusion Holling-Tanner predator-prey model[END_REF][START_REF] Mukherjee | Bifurcation scenario of Turing patterns in prey-predator model with nonlocal consumption in the prey dynamics[END_REF][START_REF] Murray | Parameter space for Turing instability in reaction diffusion mechanisms: a comparison of models[END_REF][START_REF] Ritchie | Turing instability and pattern formation on directed networks[END_REF][START_REF] Toole | Pattern formation in Turing systems on domains with exponentially growing structures[END_REF][START_REF] Wang | Vegetation pattern formation of a water-biomass model[END_REF] and the references cited therein. In particular, it is well explained how Turing patterns are continuous in space and it is often observed that they can take various forms such as spots, spirals or labyrinths. Nonetheless, another type of patterns seems to have been less investigated, namely, discontinuous patterns. Indeed, discontinuous patterns have been identified in several real-world systems, notably in microbiology [START_REF] Cygan | Stable discontinuous stationary solutions to reaction-diffusion-ode systems[END_REF][START_REF] Marciniak-Czochra | Reaction-difusion model of early carcinogenesis: The effects of influx of mutated cells[END_REF][START_REF] Marciniak-Czochra | Pattern formation in a diffusion-ODE model with hysteresis[END_REF], in forest ecosystems models [START_REF] Cantin | Non-existence of the global attractor for a partly dissipative reaction-diffusion system with hysteresis[END_REF][START_REF] Cantin | Mathematical modeling of forest ecosystems by a reactiondiffusion-advection system: impacts of climate change and deforestation[END_REF][START_REF] Kuznetsov | A cross-diffusion model of forest boundary dynamics[END_REF][START_REF] Huy | Stationary solutions to forest kinematic model[END_REF] or in competitive species models [START_REF] Aronson | A density-dependent diffusion system with stable discontinuous stationary solutions[END_REF][START_REF] Mimura | Stationary pattern of some density-dependent diffusion system with competitive dynamics[END_REF]. It turns out that the reaction-diffusion systems in which discontinuous patterns are proven to appear admit common properties. First, they are degenerate reaction-diffusion systems, that is, at least one equation of the model is non-diffusive. In other words, one species of the system is sedentary, such as trees in the forest model presented in [START_REF] Kuznetsov | A cross-diffusion model of forest boundary dynamics[END_REF] for instance, which are typically motionless individuals (it is interesting to note that the population of individuals of the sedentary species can paradoxically be capable of moving [START_REF] Einav | Indirect diffusion effect in degenerate reaction-diffusion systems[END_REF]). Furthermore, the non-diffusive equation admits a strong nonlinearity which is not invertible. This nonlinearity is often related to a hysteresis process that acts in the system (see for instance [START_REF] Aiki | A prey-predator model with hysteresis effect[END_REF][START_REF] Kopfovà | Hysteresis in biological models[END_REF][START_REF] Köthe | Hysteresis-driven pattern formation in reaction-diffusion-ODE systems[END_REF][START_REF] Staal | Hysteresis of tropical forests in the 21st century[END_REF]). Overall, discontinuous patterns seem to appear when the smoothing effect of the diffusion operator [START_REF] Cantrell | Spatial ecology via reaction-diffusion equations[END_REF][START_REF] Henry | Geometric theory of semilinear parabolic equations[END_REF] is perturbed by the combination of a degeneracy in the diffusion process and a non-invertible term. Note however that the case of non degenerate reaction-diffusion systems with discontinuous coefficients could also reveal discontinuous patterns. Although the existence of discontinuous patterns has been recently studied in the aforementioned papers [START_REF] Aronson | A density-dependent diffusion system with stable discontinuous stationary solutions[END_REF][START_REF] Cantin | How hysteresis produces discontinuous patterns in degenerate reaction-diffusion systems[END_REF][START_REF] Cygan | Stable discontinuous stationary solutions to reaction-diffusion-ode systems[END_REF][START_REF] Köthe | Hysteresis-driven pattern formation in reaction-diffusion-ODE systems[END_REF][START_REF] Huy | Stationary solutions to forest kinematic model[END_REF][START_REF] Marciniak-Czochra | Reaction-difusion model of early carcinogenesis: The effects of influx of mutated cells[END_REF][START_REF] Marciniak-Czochra | Pattern formation in a diffusion-ODE model with hysteresis[END_REF][START_REF] Mimura | Stationary pattern of some density-dependent diffusion system with competitive dynamics[END_REF], their characteristic properties are far from being completely understood. Indeed, very recent results have been established on the asymptotic behavior of the trajectories determined by these degenerate reaction-diffusion systems. Notably, a result of non-existence of the global attractor for a degenerate reaction-diffusion system with hysteresis has been proved in [START_REF] Cantin | Non-existence of the global attractor for a partly dissipative reaction-diffusion system with hysteresis[END_REF]. Meanwhile, the weak convergence of the solutions to a forest model towards heterogeneous stationary solutions has been established in [START_REF] Iwasaki | Asymptotic convergence of solutions to the forest kinematic model[END_REF] using the Lojasiewicz-Simon gradient inequality, and the weak convergence of the solutions to a simplified model towards discontinuous patterns has been proved in [START_REF] Cantin | How hysteresis produces discontinuous patterns in degenerate reaction-diffusion systems[END_REF] using a macroscopic mass effect under symmetry assumptions. In parallel, the local stability of the discontinuous patterns, that is, their behavior under a perturbation of the initial condition, has been studied in numerous papers (for instance in [START_REF] Cygan | Stable discontinuous stationary solutions to reaction-diffusion-ode systems[END_REF][START_REF] Härting | Stable patterns with jump discontinuity in systems with Turing instability and hysteresis[END_REF][START_REF] Köthe | Hysteresis-driven pattern formation in reaction-diffusion-ODE systems[END_REF][START_REF] Marciniak-Czochra | Instability of Turing patterns in reaction-diffusion-ode systems[END_REF]). But their structural stability, that is, their behavior under the effect of a perturbation of the nonlinear function involved in the system, has not been analyzed.

Contribution.

Therefore, to the best of our knowledge, the structural stability of the discontinuous patterns emerging in degenerate reaction-diffusion systems with respect to a perturbation has never been analyzed. In this paper, it is precisely our aim to bring a new contribution to a better understanding of the robustness of discontinuous patterns in this type of systems. To that aim, we consider a generic model for the formation of discontinuous patterns and we focus on the structural response of its solutions under the action of a smooth perturbation. The model we consider is written

∂u ∂t = αw -q µ (u), ∂w ∂t = δ∆w -βw + αu. ( 1
)
As will be detailed below, the latter system is defined on a bounded domain Ω ⊂ R 1 ; it is supplemented with the Neumann boundary condition on ∂Ω and with a given initial condition (u 0 , w 0 ) defined in Ω. The first equation in [START_REF] Aiki | A prey-predator model with hysteresis effect[END_REF], describing the evolution of a sedentary species, is non-diffusive and contains a nonlinearity q µ (u) which is given by a non-invertible cubic function, perturbed by a smooth function parametrized by µ. The second equation describes the evolution of a mobile species which interacts with the first species. The unperturbed model has been studied, as a simplification of the forest model presented in [START_REF] Kuznetsov | A cross-diffusion model of forest boundary dynamics[END_REF], in [START_REF] Cantin | Non-existence of the global attractor for a partly dissipative reaction-diffusion system with hysteresis[END_REF] where it has been proved that it cannot admit the global attractor, and in [START_REF] Cantin | How hysteresis produces discontinuous patterns in degenerate reaction-diffusion systems[END_REF] where the weak convergence of a non trivial set of solutions towards discontinuous patterns has been established. Here, we focus on the effect of the perturbation on the stationary solutions to [START_REF] Aiki | A prey-predator model with hysteresis effect[END_REF]. We investigate the bifurcations of homogeneous stationary solutions and prove that the trivial solution is the only one to resist to a perturbation of strong intensity (Proposition 2). We then prove our main result, which states that the perturbation transforms discontinuous patterns, by adding a second discontinuity jump to the first discontinuity jump which is observed in the patterns of the unperturbed model (Theorems 4 and 5). Our analysis of the stationary solutions relies on geometric methods, with a careful analysis of a Hamiltonian system that admits a discontinuous right-hand side and is equivalent to the problem of finding heterogeneous equilibria of (1). Our new theoretical results are illustrated by numerical simulations and we provide additional animations along with this paper, so as to guide intuition on the behavior of the dynamical system under study.

Outline. Our paper is organized as follows. In Section 2, we present in detail the degenerate reaction-diffusion system we are studying and the basic assumptions on the perturbation introduced in the model. Local in time solutions and global solutions are then proven to exist. In Section 3, we analyze, by geometric methods, the distribution of homogeneous and heterogeneous stationary solutions to system (1) and their structural stability under the action of the perturbation. The emergence of two-jumps discontinuous patterns is interpreted for the case of a forest ecosystem. Finally, in Section 4, we show numerical experiments that support our theoretical statements.

Preliminaries and setting of the problem

In this section, we present the degenerate reaction-diffusion model under study, the assumptions on the perturbation introduced in the model, and basic results on the existence of local and global solutions.

Notations. Throughout this paper, the symbols k i , with i ∈ N, will denote positive constants. The space of continuous (respectively continuously differentiable) functions defined on an interval I ⊂ R with values in a Banach space X will be denoted C I, X) (respectively C 1 (I, X)). Let Ω denote an open, connected and bounded subset of R n , where n is a positive integer; we will assume that its boundary is sufficiently regular. Lebesgue spaces on Ω will be denoted L p (Ω) and Sobolev spaces on Ω will be denoted W k,p (Ω), with p ∈ [1, ∞] and k ∈ N. Those functional spaces are Banach spaces whose norms will be denoted • L p (Ω) and • W k,p (Ω) respectively. For p = 2, we simply note

H k (Ω) = W k,2 (Ω).

A generic model for the formation of discontinuous patterns

Various models reproducing the formation of discontinuous patterns in biological or ecological systems take the form of a degenerate reaction-diffusion system with hysteresis in the non-diffusive equation [START_REF] Kuznetsov | A cross-diffusion model of forest boundary dynamics[END_REF][START_REF] Marciniak-Czochra | Reaction-difusion model of early carcinogenesis: The effects of influx of mutated cells[END_REF][START_REF] Mimura | Stationary pattern of some density-dependent diffusion system with competitive dynamics[END_REF]. The following degenerate reaction-diffusion system unifies this variety of models:

                   ∂u ∂t = αw -q µ (u), (t, x) ∈ (0, +∞) × Ω, ∂w ∂t = δ∆w -βw + αu, (t, x) ∈ (0, +∞) × Ω, ∂w ∂ν = 0, (t, x) ∈ (0, +∞) × ∂Ω, u(0, x) = u 0 (x), w(0, x) = w 0 (x), x ∈ Ω. (2)
Here, Ω is a one-dimensional bounded domain. The unknown functions u = u(t, x), w = w(t, x) are defined in (0, +∞) × Ω, and the initial conditions u 0 , w 0 are defined in Ω. ∆ denotes the Laplace diffusion operator defined by ∆w = ∂ 2 w ∂x 2 . The outward normal vector to ∂Ω at point x is denoted by ν(x). The parameters α, β, δ are positive coefficients, and q µ (u) is a smooth function defined by a perturbation of a cubic:

q µ (u) = q(u) + µ p(u), q(u) = u(u 2 -1), u ∈ R, µ ≥ 0. ( 3
)
The unknowns u(t, x) and w(t, x) correspond to densities of biological individuals; u(t, x) is a sedentary species (e.g. trees) and w(t, x) is a mobile species (e.g. seeds).

Next, we assume that the perturbation p(u) is continuously differentiable in R and satisfies the following properties:

∀s ∈ R, p(s) + p(-s) = 0, ( 4
)
∀s ≥ 0, p(s) ≥ 0, (5)
∃M 1 > 0, ∀s ∈ R, |p(s)| + |p (s)| ≤ M 1 . ( 6
)
We also assume that q µ (0) increases when µ increases in [0, +∞), is greater than α β for µ sufficiently large, and that the function q µ (u) admits a unique minimum m µ on ûµ ∈ R + , such that m µ increases when µ increases in [0, +∞), ( 7)

m µ > α 2 β ûµ for µ sufficiently large. (8)
Finally, we assume that there exists µ 1 > 0 such that q µ (u) presents three monotone branches for µ < µ 1 and five monotone branches for µ > µ 1 .

With these assumptions, we easily prove the following lemma.

Lemma 1.

For each µ ≥ 0, the function q µ (u) satisfies

q µ (u) + q µ (-u) = 0, ( 9
)
for all u ∈ R. Furthermore, there exist µ 2 , µ 3 , µ 4 , such that

0 < µ 1 < µ 2 < µ 3 < µ 4
and for which the following properties hold.

(i) If 0 ≤ µ < µ 1 , then q µ (u) presents three monotone branches, admits three zeros and three intersection points with the line α 2 β u. (ii) If µ 1 < µ < µ 2 , then q µ (u) presents five monotone branches, admits five zeros but three intersection points with the line α 2 β u. (iii) If µ 2 < µ < µ 3 , then q µ (u) presents five monotone branches, admits five zeros and five intersection points with the line α 2 β u. (iv) If µ 3 < µ < µ 4 , then q µ (u) presents five monotone branches, admits a unique zero and five intersection points with the line α 2 β u. (v) If µ 4 < µ, then q µ (u) presents five monotone branches, admits a unique zero and a unique intersection point with the line α 2 β u.

A basic example of a perturbation p(u) satisfying the above assumptions is given by

p(u) = 20u (1 + 10u 2 ) 2 , u ∈ R. (10)
Note that all the numerical computations presented in this paper where produced with this example. We emphasize that the function q µ (u) derives from a potential Q µ (u) which can be written

Q µ (u) = u 4 4 - u 2 2 + µ u 0 p(s)ds, u ∈ R. (11)
The number of zeros of the function q µ (u) obviously determines the number of local minima of the potential Q µ (u). Hence, for 0 ≤ µ < µ 1 , Q µ (u) admits two local minima; for µ 1 < µ < µ 3 , it admits three local minima; for µ 3 < µ, it admits a unique global minimum. The effect of the perturbation parameter µ on the shapes of the function q µ (u) and of the potential Q µ (u) is illustrated in Figures 1 and2. We observe that the critical values µ 1 and µ 3 determine a structural modification of the potential Q µ (u). We will see below that the critical values µ 2 and µ 4 , which alternate with µ 1 and µ 3 , determine a structural modification of the homogeneous stationary solutions of the degenerate reaction-diffusion system (2). Now, our aim is to prove that the degenerate reaction-diffusion system (2) admits local solutions. To this end, we write it as a semi-linear equation in a proper Banach space.

µ = 0 (a) 0 0 u 0 < µ < µ1 (b) 0 0 u µ1 < µ < µ3 (c) 0 0 u µ3 < µ (d)

Abstract formulation and local solutions

Following [START_REF] Yagi | Abstract parabolic evolution equations and their applications[END_REF], we handle the degenerate reaction-diffusion system (2) in the Banach space X defined by

X = L ∞ (Ω) × L 2 (Ω), (12)
equipped with the product norm

U X = u L ∞ (Ω) + w L 2 (Ω) , U = (u, w) ∈ X.
We consider the differential operator Λ defined as the realization of -δ∆ + β in L 2 (Ω) with the Neumann boundary condition on ∂Ω. It is known that Λ is a positive definite self-adjoint and sectorial operator, of angle strictly less than π 2 , with domain

D(Λ) = H 2 N (Ω) = w ∈ H 2 (Ω) ; ∂w ∂ν = 0 on ∂Ω .
Hence, the diagonal operator A = diag {1, Λ} is also a sectorial operator in X, with angle strictly less than π 2 , and with domain

D(A) = L ∞ (Ω) × D(Λ).
Here, we emphasize that the domain of A is not compactly embedded in X, although the domain of Λ is compactly embedded in L 2 (Ω). This is due to the absence of diffusion in the first equation, and partly determines an original asymptotic behavior, as proved in [START_REF] Cantin | Non-existence of the global attractor for a partly dissipative reaction-diffusion system with hysteresis[END_REF].

Next, we consider an exponent η ∈ 3 4 , 1 . The sectorial operator Λ admits a fractional power Λ η whose domain is given by

D(Λ η ) = H 2η N (Ω) = w ∈ H 2η (Ω) ; ∂w ∂ν = 0 on ∂Ω ,
where H 2η (Ω) is the interpolation space W 2η,2 (Ω). We have the continuous embeddings

H 2η (Ω) ⊂ C Ω ⊂ L ∞ (Ω) ⊂ L 2 (Ω), (13)
and the norms u D(Λ η ) and Λ η u L 2 (Ω) are equivalent. The diagonal operator A also admits a fractional power A η and its domain is given by

D(A η ) = L ∞ (Ω) × H 2η N (Ω).
In this way, the degenerate reaction-diffusion system (2) can be written in an abstract form

dU dt + AU = F µ (U ), t > 0, ( 14
)
where F µ (U ) is the nonlinear operator defined by

F µ (U ) = αw -q µ (u) + u αu , U = (u, w) ∈ D(A η ). ( 15
)
Note that the domain of the nonlinear operator F µ is uniform with respect to the perturbation parameter µ. Next, elementary computations show that F µ satisfies the estimation

F µ (U ) -F µ ( Ũ ) X ≤ C 1 1 + U 2 X + Ũ 2 X A η (U -Ũ ) X + µC 2 A η (U -Ũ ) X , (16)
for all µ in R + and all U, Ũ in D(A η ), with positive constants C 1 , C 2 . Therefore, by virtue of Theorem 4.1 in [START_REF] Yagi | Abstract parabolic evolution equations and their applications[END_REF], we can state the existence of local in time solutions to the problem [START_REF] Iwasaki | Asymptotic convergence of solutions to the forest kinematic model[END_REF].

Theorem 1 (Local solutions). Let µ ≥ 0. For all U 0 ∈ X, the Cauchy problem defined by [START_REF] Iwasaki | Asymptotic convergence of solutions to the forest kinematic model[END_REF] and U (0) = U 0 admits a unique local solution U µ = (u µ , w µ ) defined on [0, T U0 ] with T U0 > 0 in the function space

u µ ∈ C [0, T U0 ], L ∞ (Ω) ∩ C 1 (0, T U0 ], L ∞ (Ω) , w µ ∈ C (0, T U0 ], H 2 N (Ω) ∩ C [0, T U0 ], L 2 (Ω) ∩ C 1 (0, T U0 ], L 2 (Ω) .

Global solutions and continuous dynamical system

We continue by showing that local in time solutions to the degenerate reaction-diffusion system (2) are global and determine a continuous dynamical system. The following proposition establishes an a priori estimate for any local solution.

Proposition 1 (Dissipative estimate). Let µ ≥ 0. There exists a positive exponent λ and a positive constant C 3 such that, for all U 0 ∈ X, the local in time solution U µ of the Cauchy problem defined by [START_REF] Iwasaki | Asymptotic convergence of solutions to the forest kinematic model[END_REF] and

U (0) = U 0 , defined on [0, T U0 ] with T U0 > 0, satisfies U µ (t) X ≤ C 3 e -λt U 0 X + 1 , ( 17
)
for all t ∈ [0, T U0 ].
Proof. Let U 0 ∈ X. We denote by U µ (t) the local in time solution to the degenerate reaction-diffusion system (2) stemming from U 0 , defined on [0, T U0 ] with T U0 > 0. We introduce the energy function L defined for t ≥ 0 by

L(t) = 1 2 Ω u 2 µ (t) + w 2 µ (t) dx.
The function L is continuously differentiable and we have

L(t) = Ω ∂u µ ∂t (t)u µ (t)dx + Ω ∂w µ ∂t (t)w µ (t)dx = α Ω u µ (t)w µ (t)dx - Ω u µ (t)q µ u µ (t) dx + δ Ω w µ (t)∆w µ (t)dx + α Ω u µ (t)w µ (t)dx -β Ω w 2 µ (t)dx.
We can apply the Green formula and the Neumann boundary condition to deduce that

Ω w µ (t)∆w µ (t)dx ≤ - Ω |∇w µ (t)| 2 dx ≤ 0.
Next, we easily prove that the following lower estimate holds for all s ∈ R:

sq µ (s) ≥ g 1 s 2 -g 2 ,
where g 1 is an arbitrarily large coefficient, provided g 2 is sufficiently large. Indeed, elementary computations lead to

sq µ (s) -(g 1 s 2 -g 2 ) = s 4 -s 2 + µsp(s) -g 1 s 2 + g 2 = s 4 -(1 + g 1 )s 2 + g 2 + µsp(s) ≥ s 4 -(1 + g 1 )s 2 + g 2 ≥ s 2 - 1 + g 1 2 2 + g 2 - (1 + g 1 ) 2 4 ≥ g 2 - (1 + g 1 ) 2 4 ,
where we have used assumptions (4) and ( 5) in order to minor µsp(s) by 0. Hence, if g 1 is chosen arbitrarily large, it suffices to choose g 2 such that g 2 -(1+g1) 2

4

≥ 0. With this lower estimate, we deduce

L(t) ≤ 2α Ω u µ (t)w µ (t)dx -g 1 Ω u 2 µ (t)dx -β Ω w 2 µ (t)dx + g 2 |Ω|.
Now we apply the generalized Young inequality ab

≤ 1 2ε a 2 + ε 2 b 2 , which is valid of all a, b ∈ R and all ε > 0, to write Ω u µ (t)w µ (t)dx ≤ 1 2ε Ω u 2 µ (t)dx + ε 2 Ω w 2 µ (t)dx, (18)
for any ε > 0. We choose ε = β 2α , which leads to

2α Ω u µ (t)w µ (t)dx ≤ 2α 2 β Ω u 2 µ (t)dx + β 2 Ω w 2 µ (t)dx. ( 19
)
We obtain

L(t) ≤ -g 1 - 2α 2 β Ω u 2 µ (t)dx - β 2 Ω w 2 µ (t)dx + g 2 |Ω|.
Finally, we choose g 1 so that g 1 > 2α 2 β and we denote

λ = min g 1 - 2α 2 β , β 2 .
We obtain

L(t) ≤ -2λL(t) + g 2 |Ω|.
Applying Gronwall lemma leads to

L(t) ≤ L(0)e -2λt + g 2 |Ω| 2 ,
which reduces to

Ω u 2 µ (t) + w 2 µ (t) dx 1/2 ≤ Ω u 2 0 + w 2 0 dx 1/2 e -λt + k 1 ,
which proves that

U µ (t) L 2 (Ω) 2 ≤ k 2 U 0 L 2 (Ω) 2 e -λt + 1 . ( 20
)
Afterwards, it remains to show how the dissipative estimate [START_REF] Huy | Stationary solutions to forest kinematic model[END_REF], which holds in L 2 (Ω) 2 , implies the stronger estimate (17) in X. Indeed, it suffices to apply similar arguments as in the proof of Proposition 11.1 in [START_REF] Yagi | Abstract parabolic evolution equations and their applications[END_REF]. Hence, we can show that w and u satisfy respectively

w µ (t) L ∞ (Ω) ≤ k 3 (1 + t -η )e -k4t U 0 L 2 (Ω) 2 + 1 , ( 21
)
u µ (t) L ∞ (Ω) ≤ k 5 e -k6t U 0 X + 1 , ( 22
)
for all t > 0. Finally, we combine the dissipative estimates ( 20) and ( 22) to obtain [START_REF] Köthe | Hysteresis-driven pattern formation in reaction-diffusion-ODE systems[END_REF], which completes the proof.

With the dissipative estimate [START_REF] Köthe | Hysteresis-driven pattern formation in reaction-diffusion-ODE systems[END_REF], we can directly state the existence of global solutions, that determine a continuous dynamical system.

Theorem 2 (Global solutions and continuous dynamical system). Let µ ≥ 0. For all U 0 ∈ X, the Cauchy problem defined by [START_REF] Iwasaki | Asymptotic convergence of solutions to the forest kinematic model[END_REF] and U (0) = U 0 admits a unique global solution U µ (t, U 0 ) = (u µ , w µ ) defined on [0, +∞) in the function space

u µ ∈ C [0, +∞), L ∞ (Ω) ∩ C 1 (0, +∞), L ∞ (Ω) , w µ ∈ C (0, +∞), H 2 N (Ω) ∩ C [0, +∞), L 2 (Ω) ∩ C 1 (0, +∞), L 2 (Ω) .
Furthermore, the degenerate reaction-diffusion system (2) determines a continuous dynamical system S µ (t) defined in X by

S µ (t)U 0 = U µ (t, U 0 ), t ≥ 0. ( 23
)
Remark 1 (Lack of compactness). Using the dissipative estimate [START_REF] Huy | Stationary solutions to forest kinematic model[END_REF], it can be shown that the continuous dynamical system S µ (t) admits an absorbing set B µ ⊂ X which is bounded in D(A) (see [START_REF] Yagi | Abstract parabolic evolution equations and their applications[END_REF], Chapter 11, Section 4). However, we cannot show that the absorbing set B µ is compact, hence it turns out that the study of the asymptotic behavior of the dynamical system S µ (t) cannot be described by means of the global attractor.

Existence of a potential

To understand its dynamics, it is worth emphasizing that the behavior of the degenerate reactiondiffusion system (2) is partly governed by the potential H µ (u, w) defined by

H µ (u, w) = β w 2 2 -αuw + Q µ (u), u, w ∈ R,
so that system (2) can be written under the gradient form

∂u ∂t = - ∂H µ ∂u , ∂w ∂t = δ∆w - ∂H µ ∂w .

Analysis of homogeneous and heterogeneous stationary solutions

In this section, we focus on the effect of the perturbation parameter µ on the stationary solutions of the degenerate reaction-diffusion system (2). First, we analyze the bifurcations of homogeneous solutions. We exhibit a pitchfork bifurcation and two saddle-node bifurcations. Then, we study the heterogeneous stationary solutions. We show that the degenerate reaction-diffusion system (2) admits non trivial discontinuous patterns and that a bifurcation of heteroclinic orbits tests the robustness of these patterns.

Bifurcations of homogeneous stationary solutions

The homogeneous stationary solutions of the degenerate reaction-diffusion system (2) are the solutions of the following system αw -q µ (u) = 0, -βw

+ αu = 0, (24)
which is equivalent to

w = α β u, w = 1 α q µ (u).
The solutions of the latter system are obviously determined by the intersection points of the function q µ (u) with the line α 2 β u (or equivalently between the function 1 α q µ (u) and the line α β u). The effect of the perturbation parameter µ on the number of these intersection points is illustrated in Figure 4. β u (or equivalently between the function 1 α qµ(u) and the line α β u) determines the number of stationary homogeneous solutions for the degenerate reaction-diffusion system (2). (a) For 0 ≤ µ < µ2, the degenerate reaction-diffusion system (2) admits three homogeneous stationary solutions. (b) For µ2 ≤ µ < µ4, it admits five homogeneous stationary solutions. (c) For µ4 ≤ µ, it admits a unique homogeneous stationary solution.

0 0 u 0 ≤ µ < µ 2 1 α q µ (u) α β u (a) 0 0 u µ 2 < µ < µ 4 (b) 0 0 u µ 4 < µ (c)
The following proposition enumerates the stationary homogeneous solutions with respect to the perturbation parameter µ.

Proof. First, assertions (i), (ii) and (iii) are direct consequences of Lemma 1.

Next, let Ū denote a homogeneous stationary solution of the degenerate reaction-diffusion system (2). The stability of Ū can be analyzed by considering the spectrum σ( Ā) of the operator Ā defined by Ā = A -F ( Ū ), where F ( Ū ) denotes the Fréchet derivative of the function F . Indeed, Ū is linearly stable if and only if

σ( Ā) ⊂ {z ∈ C ; z > 0} . ( 25
)
Hence, we compute

Ā = A -F ( Ū ) = q µ (ū) -α -α -δ∆ + β .
We then consider the spectrum problem

(λ -Ā)U = P, λ ∈ C, U ∈ D( Ā), P ∈ X.
Denoting U = (u, w) and P = (p 1 , p 2 ) , we obtain

λ -q µ (ū) u + αw = p 1 , αu + λ + δ∆ -β w = p 2 , which leads to α 2 -λ -q µ (ū) (λ + δ∆ -β) w = αp 1 -λ -q µ (ū) p 2 . ( 26
)
If λ = q µ (ū), then equation ( 26) reduces to αw = p 1 .

Thus w does not belong to H 2 N (Ω) for all p 1 in L ∞ (Ω). We can deduce that q µ (ū) ∈ σ( Ā). As q µ (u - µ ) < 0 and q µ (-u - µ ) < 0, we obtain that U - µ and -U - µ are unstable for all µ ∈ (µ 2 , µ 4 ).

It remains to analyze the stability of O, U +

µ and -U + µ . Let us assume that λ = q µ (ū). We denote by (ω n ) n≥0 the sequence of eigenvalues of the Laplace operator -∆ in L 2 (Ω) equipped with Neumann boundary condition on ∂Ω; this sequence is strictly increasing and its first element is null. With these notations, λ ∈ σ( Ā) if and only if λ is a solution to one of the following equations:

α 2 -λ -q µ (ū) (λ -δω n -β) = 0, n ≥ 0, or equivalently λ 2 -q µ (ū) + δω n + β λ + q µ (ū)(δω n + β) -α 2 = 0, n ≥ 0. ( 27
)
The discriminant D n of the latter equation satisfies

D n = q µ (ū) -(δω n + β) 2 + 4α 2 > 0,
which guarantees that ( 27) admits two real roots λ n 1 and λ n 2 , for all n ≥ 0. These roots satisfy

λ n 1 + λ n 2 = q µ (ū) + δω n + β, λ n 1 × λ n 2 = q µ (ū)(δω n + β) -α 2 .
If ū is equal to u + µ or -u + µ , then q µ (ū) > 0. Hence, λ n 1 × λ n 2 > 0 if and only if

ω n < β δq µ (ū) × α 2 β -q µ (ū) . ( 28
)
Furthermore, if ū = ±u + µ , then we have q µ (ū) > α 2 β . Therefore, condition ( 28) is never satisfied. It follows that λ n 1 > 0 and λ n 2 > 0 for all n ≥ 0, which proves that U + µ and -U + µ are linearly stable. Similarly, if ū = 0, we have q µ (ū) > α 2 β for µ > µ 2 and q µ (ū) < α 2 β for µ < µ 2 . We can deduce that O is unstable for µ < µ 2 and linearly stable for µ > µ 2 . The proof is complete.

Without loss of generality, we can write Ω = (-, ) with > 0. It is convenient to write ∆w = ẅ and to introduce the variable v = ẇ, so that system (29) can be rewritten

               0 = αw(x) -q µ u(x) , x ∈ Ω, v(x) = β δ w(x) - α δ u(x), x ∈ Ω, ẇ(x) = v(x), x ∈ Ω, v(-) = v( ) = 0. (30)
It would be easy to solve the latter system if the function q µ (u) was invertible. However, as shown in Lemma 1, the function q µ (u) is not monotone, hence it is not invertible. At the contrary, there exist infinitely many functions ψ defined in R, which are inverses of q µ and thus satisfy

ψ • q µ (u) = q µ • ψ(u) = u,
for all u ∈ R. Among these inverses ψ, none is continuous everywhere on R; some are piecewise continuous and monotone, whereas others admit an infinite number of discontinuity points and are not monotone. Now, let ψ denote one possible inverse of q µ . Assume that ψ is piecewise continuous on R and consider the following system of ordinary differential equations:

v = β δ w - α δ ψ(αw), ẇ = v. ( 31
)
The latter system admits in its first equation a discontinuous right-hand side, thus should be solved by appropriate methods (see notably [START_REF] Filippov | Differential equations with discontinuous righthand sides: control systems[END_REF]). In particular, its orbits can admit singular points where they are continuous but not differentiable. Obviously, if (v, w) is a solution of system (31) defined on the domain Ω, which moreover satisfies the condition v(-

) = v( ) = 0, (32)
then it determines a solution (u, w) of system (30) by setting u = ψ(αw). As we shall see in the next section, the existence of such a solution depends on the perturbation parameter µ. Since there is an infinite number of choices for the inverse ψ, we can obtain an infinite number of heterogeneous solutions (u, w) of the degenerate reaction-diffusion system (2). We will now focus on particular solutions, which are monotone and discontinuous with respect to the space variable x ∈ Ω.

Existence and robustness of discontinuous patterns

Since the degenerate reaction-diffusion system (2) possibly admits infinitely many heterogeneous stationary solutions, we are interested in this section in finding those which are discontinuous. The following definition specifies what we call discontinuous pattern.

Definition 1. An heterogeneous stationary solution {u(x), w(x)} x∈Ω of the degenerate reactiondiffusion (2) such that u is discontinuous is called a discontinuous pattern.

To construct discontinuous patterns, we first observe that system (31) can be rewritten, for any inverse ψ µ of q µ , as a Hamiltonian system

v = ∂K µ ∂w , ẇ = - ∂K µ ∂v , (33)
where K µ (v, w) is the potential defined by

K µ (v, w) = β 2δ w 2 - 1 2 v 2 - α δ Ψ µ (αw), (34) 0
-w + µ -a 0 a w + µ v w 0 ≤ µ < µ 2 , σ = 0 (a) - 0 -a 0 a Ω 0 ≤ µ < µ 2 , σ = 0 u(x) w(x) (b)
Figure 7: (a) Phase portrait of the Hamiltonian system (33) obtained with ψ 0 µ . In the region delimited by the heteroclinic orbits that join the saddles (0 -w + µ ) and (0, w + µ ), a continuous band of periodic orbits encircles the center (0, 0). The space Ω = (-, ) necessary to join two successive points of the axis {v = 0} along one of these orbits varies increasingly from 0 (for the stationary point (0, 0)) to +∞ (for the heteroclinic orbits). (b) By setting u(x) = ψ 0 µ αw(x) , each solution {v(x), w(x)} x∈Ω of the Hamiltonian system (33) determines a heterogeneous stationary solution {u(x), w(x)} x∈Ω of the degenerate reaction-diffusion system (2), such that u is piecewise continuous on Ω, with a unique discontinuity jump. [START_REF] Ritchie | Turing instability and pattern formation on directed networks[END_REF] obtained for ψ σ µ with σ > 0. The two saddles (0, -w + µ ) and (0, w + µ ) are not connected by heteroclinic orbits as in the case σ = 0. Now, only one saddle admits a homoclinic orbit that encircles a continuous band of periodic orbits. If σ is too large, then the existence of an orbit {v(x), w(x)} x∈Ω connecting -a < 0 to a > 0 is compromised.

0 -w + µ -a 0 a w + µ v w ≤ µ < µ 2 , σ > 0
However, Lemma 2 guarantees that for σ > 0 sufficiently small, an inverse ψ σ µ can be chosen arbitrarily close to ψ 0 µ . Hence, it is seen that two distinct discontinuous patterns {u 0 (x), w 0 (x)} x∈Ω and {u σ (x), w σ (x)} x∈Ω can be constructed and are arbitrarily close to each other. Therefore, we can conclude that the degenerate reaction-diffusion system (2) admits a continuum of discontinuous patterns. We obtain the following theorem.

Theorem 3 (Continuum of discontinuous patterns). For each µ ∈ [0, µ 2 ), the degenerate reactiondiffusion system (2) admits an infinite number of discontinuous patterns with a unique discontinuity Theorem 5. Assume that µ is sufficiently near µ 2 and that Ω is sufficiently large. Then for each ε > 0, one can find two distinct discontinuous patterns {u 0 (x), w 0 (x)} x∈Ω , {u σ (x), w σ (x)} x∈Ω such that u 0 -u σ L ∞ (R) + w 0 -w σ L ∞ (R) < ε, where u 0 (x) and u σ (x) admit two discontinuity jumps.

We skip the proof of the latter theorem, since it is very similar to that of Theorem 3. Now, if µ ∈ µ * (a), μ , then the existence of discontinuous patterns {u(x), w(x)} x∈Ω such that w(x) connects (0, -a) to (0, a) and u(x) admits two discontinuity jumps is not guaranteed, even if Ω is large (see Figure 11(a)). Finally, if µ ∈ (μ, µ 4 ), then these discontinuous patterns cannot exist, whatever the value of a and the size of the domain Ω are (see Figure 11(b)). Roughly speaking, if such an orbit would exist, then it should cross the homoclinic orbits that encircle the centers (0, -w + µ ) and (0, w + µ ), which is excluded, since orbits of an autonomous system cannot cross each other. Remark 3 (Interpretation in ecology). We recall that the degenerate reaction-diffusion system (2) has been obtained as a simplification of the forest model presented in [START_REF] Kuznetsov | A cross-diffusion model of forest boundary dynamics[END_REF]. The convergence of its orbits towards a discontinuous pattern admitting a unique discontinuity jump has been interpreted in [START_REF] Kuznetsov | A cross-diffusion model of forest boundary dynamics[END_REF] as the formation of an ecotone, that is, frontier between the forest and a neighbor ecosystem. In the simplified forest model, u(t, x) and w(t, x) correspond to the densities of trees and seeds respectively. Hence, the of a one-jump discontinuous pattern corresponds to the area by the forest, and lower level corresponds to a neighbor ecosystem without trees. The existence of discontinuous patterns admitting two discontinuity jumps suggests that the perturbation of the dynamics of a forest ecosystem can lead to the emergence of an intermediate ecosystem between the forest and the neighbor ecosystem. Furthermore, if the intensity of the perturbation increases, then the intermediate ecosystem can spread over the forest ecosystem, which reproduces the invasion of the forest ecosystem by a new ecosystem. The invading ecosystem can take the a as suggested in [START_REF] Cantin | On the degradation of forest ecosystems by extreme events: Statistical model checking of a hybrid model[END_REF]. Numerical simulations shown in Section 4 illustrate this ecosystem invasion process.

Case of a strong perturbation

We end our research of discontinuous patterns with the case of a strong perturbation, for which we assume that µ 4 < µ. In this case, the degenerate reaction-diffusion system (2) admits (0, 0) as a unique homogeneous stationary solution. By a reasoning similar to the cases of small perturbations, we can establish the phase portrait of the Hamiltonian system [START_REF] Ritchie | Turing instability and pattern formation on directed networks[END_REF], for each inverse ψ µ of q µ . If the inverse is chosen monotone and piecewise continuous, then it is easily seen that the Hamiltonian system (33) admits a saddle as unique critical point (see Figure 12). Therefore, no orbit connecting two points of the axis {v = 0} can be found, and no discontinuous pattern can be constructed.

Overall, the effect of a variation of the perturbation parameter µ on the dynamics of the Hamiltonian system (33) obtained with a monotone and piecewise continuous inverse ψ µ of q µ is schematized as a bifurcation diagram in Figure 13 (where we abstract from the singular points of the orbits of system (33)). A pitchfork bifurcation occurs at µ 2 ; particular patterns vanish at µ * (a); a bifurcation of heteroclinic orbits occurs at μ; finally, two simultaneous saddle-node bifurcations take place at µ 4 . Roughly speaking, discontinuous patterns admitting a unique discontinuity jump exist for a very small perturbation. If the intensity of the perturbation increases, then the discontinuous patterns survive and admit two discontinuity jumps. But if the intensity of the perturbation increases too much, these discontinuous patterns seem to vanish. Remark 4 (Restricted choice of inverses). We emphasize that in this section, the research of heterogeneous solutions has been realized by analyzing the dynamics of the Hamiltonian system [START_REF] Ritchie | Turing instability and pattern formation on directed networks[END_REF] with a restricted choice of an inverse ψ µ of q µ , since we have only considered monotone and piecewise continuous inverses. This restriction does not diminish the existence results of discontinuous patterns. However, an infinite number of other inverses ψ µ of q µ could be considered as well. For example, [START_REF] Ritchie | Turing instability and pattern formation on directed networks[END_REF] obtained with a monotone and piecewise continuous inverse ψµ of qµ. A pitchfork bifurcation occurs at µ2; particular discontinuous patterns connecting -a to a exist for µ < µ * (a) but vanish for µ > µ * (a) (hence they are depicted in green for µ < µ * (a) and in red for µ > µ * (a)); a bifurcation of heteroclinic orbits occurs at μ; finally, two simultaneous saddle-node bifurcations take place at µ4.

among these other inverses, completely discontinuous inverses could have been investigated. It is left as an open question for further research to determine if those inverses which have not been considered

here could also lead to the discovery of discontinuous patterns.

Numerical simulations

In this section, our aim is to illustrate our theoretical results by numerical simulations. Animations of the degenerate reaction-diffusion system (2), in a one-dimensional or two-dimensional domain, are also provided on the website https://pagesperso.ls2n.fr/ ˜cantin-g/robustness.html. We solve numerically the degenerate reaction-diffusion system (2) with four initial conditions and with various values of the perturbation parameter µ, so as to show that the orbits are attracted to discontinuous patterns. As proved in Theorems 3 and 5, these discontinuous patterns admit a unique discontinuity

Random initial condition

We end this section with a numerical simulation of the degenerate reaction-diffusion system (2) starting from a randomly generated initial condition, with µ = 1. The results are depicted in Figure 21. We observe that discontinuous patterns admitting two discontinuity jumps appear at places of the domain. The of these patterns seem to be very sensitive to the initial condition and very difficult to predict.

Conclusion

In this paper, we investigated of a degenerate reaction-diffusion system with hysteresis in its non-diffusive equation, perturbed by a smooth function. This degenerate reaction-diffusion system is a generic model for the formation of discontinuous patterns arising in biology or in ecology.

We have analyzed by geometric methods the effect of the perturbation on the stationary solutions. main results are the following:

• we studied the bifurcations of homogeneous stationary solutions and proved that the trivial equilibrium is the only homogeneous stationary solution capable of withstanding a strong perturbation;

• we studied the effect of the perturbation on the heterogeneous stationary solutions and showed how discontinuous patterns with a single discontinuity jump are transformed into patterns with two discontinuity jumps;

Figure 1 :

 1 Figure 1: Effect of the perturbation parameter µ on the function qµ(u) defined by (3). (a)-(b) For 0 ≤ µ < µ1, qµ(u) presents three monotone branches and admits three zeros. (c) For µ1 < µ < µ3, qµ(u) presents five monotone branches and admits five zeros. (d) For µ3 < µ, qµ(u) presents five monotone branches and admits a unique zero.

Figure 2 :

 2 Figure 2: Effect of the perturbation parameter µ on the potential Qµ(u) defined by (11). (a)-(b) For 0 ≤ µ < µ1, Qµ(u) admits two local minima. (c) For µ1 < µ < µ3, it admits three local minima. (d) For µ3 < µ, it admits a unique global minimum.

Figure 4 :

 4 Figure 4:The number of intersection points between the function qµ(u) and the line α 2 β u (or equivalently between the function1 α qµ(u) and the line α β u) determines the number of stationary homogeneous solutions for the degenerate reaction-diffusion system (2). (a) For 0 ≤ µ < µ2, the degenerate reaction-diffusion system (2) admits three homogeneous stationary solutions. (b) For µ2 ≤ µ < µ4, it admits five homogeneous stationary solutions. (c) For µ4 ≤ µ, it admits a unique homogeneous stationary solution.

Figure 8 :

 8 Figure 8:Phase portrait of the system (33) obtained for ψ σ µ with σ > 0. The two saddles (0, -w + µ ) and (0, w + µ ) are not connected by heteroclinic orbits as in the case σ = 0. Now, only one saddle admits a homoclinic orbit that encircles a continuous band of periodic orbits. If σ is too large, then the existence of an orbit {v(x), w(x)} x∈Ω connecting -a < 0 to a > 0 is compromised.

Figure 12 :Figure 13 :

 1213 Figure 12:Phase portrait of the Hamiltonian system[START_REF] Ritchie | Turing instability and pattern formation on directed networks[END_REF] obtained with a monotone and piecewise continuous inverse ψµ of qµ. The only critical point is a saddle at (0, 0). No orbit connecting two points of the axis {v = 0} can be found, and no discontinuous pattern can be constructed.

Figure 20 :

 20 Figure 20: Numerical simulation of the degenerate reaction-diffusion system (2) stemming from the periodic initial condition (42), with µ = 0. The orbit converges to a form that admits several discontinuous patterns with two discontinuity jumps in which uµ and wµ are still synchronized.

Acknowledgments

The author is sincerely grateful to the anonymous reviewers for the useful comments which greatly improved the presentation of this paper.

If the diffusion coefficient δ is null, then the degenerate reaction-diffusion system (2) reduces to a system of ordinary equations, whose orbits are attracted to the local minima of the potential H µ (u, w). But when the diffusion coefficient δ is positive, the dynamics is completely modified. In particular, the orbits can weakly converge towards discontinuous patterns, as will be analyzed in the sequel.

Next, it is easy to prove that H µ (u, w) admits two local minima (u + µ , w + µ ), (-u + µ , -w + µ ) and a unique saddle (0, 0) if 0 ≤ µ < µ 2 , three local minima (u + µ , w + µ ), (0, 0), (-u + µ , -w + µ ) and two saddles (u - µ , w - µ ), (-u - µ , -w - µ ) if µ 2 < µ < µ 4 and a unique global minimum (0, 0) if µ 4 < µ. The effect of the perturbation parameter µ on the energy levels of the potential H µ (u, w) is illustrated in Figure 3.

Figure 3: Effect of the perturbation parameter µ on the potential Hµ(u, w). (a) For 0 ≤ µ < µ2, Hµ(u, w) admits two local minima (u + µ , w + µ ), (-u + µ , -w + µ ) and a unique saddle (0, 0). (b) For µ2 < µ < µ4, it admits three local minima (u + µ , w + µ ), (0, 0), (-u + µ , -w + µ ) and two saddles (u - µ , w - µ ), (-u - µ , -w - µ ). (c) For µ4 < µ, it admits a unique global minimum (0, 0).

We will show in Section 3 that the critical points of the potential H µ (u, w) correspond to the homogeneous stationary solutions of the degenerate reaction-diffusion system [START_REF] Aronson | A density-dependent diffusion system with stable discontinuous stationary solutions[END_REF]. Hence, when µ varies and crosses the critical values µ 2 and µ 4 , a structural modification of the potential H µ (u, w) is observed. Note that this modification occurs after a delay in µ compared with the critical values µ 1 and µ 3 that provoke the modifications of the potential Q µ (u), as shown in Figure 2. In this paper, we focus on the effect of a perturbation of small intensity on the dynamics of the degenerate reactiondiffusion system [START_REF] Aronson | A density-dependent diffusion system with stable discontinuous stationary solutions[END_REF], that is, a perturbation that generates a third local minimum in the potential H µ (u, w). The case of a strong perturbation will be briefly treated and we will show that it leads to trivial dynamics. Therefore, the limiting value between the case of a small perturbation and that of a strong perturbation is µ 4 .

Remark 2 (Lyapunov function

). Note that the potential H µ (u, w) can be used to construct a Lyapunov function L µ (u, w), defined as

The existence of a Lyapunov function highlights the dissipative nature of the degenerate reactiondiffusion system [START_REF] Aronson | A density-dependent diffusion system with stable discontinuous stationary solutions[END_REF]. In [START_REF] Iwasaki | Asymptotic convergence of solutions to the forest kinematic model[END_REF], such a Lyapunov function has been considered for proving the weak convergence of the orbits of a system similar to (2) towards stationary solutions.

Proposition 2.

The following assertions hold.

(i) If 0 ≤ µ < µ 2 , then the degenerate reaction-diffusion system (2) admits three homogeneous stationary solutions O(0, 0),

then it admits a unique homogeneous stationary solution O(0, 0). Furthermore, U + µ and -U + µ are linearly stable for all µ ∈ [0, µ 4 ), whereas U - µ and -U - µ are unstable for all µ ∈ (µ 2 , µ 4 ); the trivial solution O(0, 0) is unstable if µ < µ 2 and stable if µ > µ 2 .

We emphasize that the homogeneous stationary solutions of the reaction-diffusion system (2) exactly coincide with the critical points of the potential H µ (u, w) defined by (2.4). More precisely, the local minima of the potential H µ (u, w) correspond to the homogeneous stationary solutions of (2) which are linearly stable, whereas the saddles of H µ (u, w) correspond to the homogeneous stationary solutions of (2) which are unstable. A bifurcation diagram showing the distribution and the stability of these homogeneous stationary solutions is depicted in Figure 5.

Two saddle-node bifurcations occur simultaneously at µ = µ4, where the linearly stable solutions

We observe that a pitchfork bifurcation occurs at µ = µ 2 , where the trivial stationary solution O becomes linearly stable, while the two unstable stationary solutions U - µ , -U - µ appear. Afterwards, two saddle-node bifurcations occur simultaneously at µ = µ 4 , where the linearly stable solutions U + µ , -U + µ merge with the unstable solutions U - µ , -U - µ . Obviously, the pitchfork bifurcation and the saddlenode bifurcations occurring at µ = µ 2 and µ = µ 4 respectively are local bifurcations. It will turn out that two global bifurcations take place in background between µ 2 and µ 4 . Indeed, a bifurcation of heteroclinic orbits occurring at μ ∈ (µ 2 , µ 4 ) shall be discovered by analyzing the structure of the heterogeneous stationary solutions of the degenerate reaction-diffusion system [START_REF] Aronson | A density-dependent diffusion system with stable discontinuous stationary solutions[END_REF]. Furthermore, we shall see that the bifurcation of heteroclinic orbits precedes the vanishing of particular patterns joining two opposite values -a and a, over a critical threshold µ * (a) ∈ (µ 2 , μ), depending on the amplitude between these opposite values.

Existence of heterogeneous stationary solutions

Heterogeneous stationary solutions u(x), w(x) of the degenerate reaction-diffusion system (2) are the solutions of the following elliptic problem

with Ψ µ given by

For simplicity, we consider only those inverses ψ µ which are monotone and piecewise continuous and analyze three cases, distinguishing very small perturbations (obtained for 0 ≤ µ < µ 2 ), small perturbations (obtained for µ 2 < µ < µ 4 ) and strong perturbations (obtained for µ 4 < µ).

Case of a very small perturbation

We assume in this subsection that 0 ≤ µ < µ 2 . Let m µ = q µ (û µ ) denote the minimum of q µ on R + (note that m µ < 0). We have the following lemma.

Lemma 2. For each µ ∈ [0, µ 2 ), the function q µ possesses a one-parameter family

of piecewise continuous and monotone inverses, admitting a unique discontinuity jump. Furthermore, for each ε > 0, one can find σ > 0 such that 4(a), where the cubic function 1 α q µ (u) crosses the line α β u at three points). Now, for each

In this way, ψ σ µ is monotone and piecewise continuous, with a unique discontinuity jump at s = σ. Finally, it is seen that ψ σ µ varies continuously with σ, which leads to [START_REF] Sun | Dynamic analysis of a plant-water model with spatial diffusion[END_REF]. Two possible inverses ψ 0 µ and ψ σ µ with σ > 0 are shown in Figure 6. Note that in the family {ψ σ µ }, ψ 0 µ is the only inverse of q µ which is symmetric with respect to 0, that is,

for all s ∈ R.

We can now begin our research of discontinuous patterns to the degenerate reaction-diffusion system (2), which are monotone and discontinuous. To that aim, we establish the phase portrait of the Hamitonian system [START_REF] Ritchie | Turing instability and pattern formation on directed networks[END_REF]. Next, we look in this phase portrait for orbits which join the axis {v = 0}, to ensure that condition [START_REF] Péron | Uncovering periodic patterns of space use in animal tracking data with periodograms, including a new algorithm for the Lomb-Scargle periodogram and improved randomization tests[END_REF] is fulfilled.

We first analyze the orbits of the Hamitonian system (33) obtained with the particular inverse ψ 0 µ of q µ . For convenience, we denote this system by (K 0 µ ), that is

Figure 6: Two possible piecewise continuous and monotone inverses ψ 0 µ and ψ σ µ of qµ in the case 0 ≤ µ < µ2. In the family {ψ σ µ }, ψ 0 µ is the only inverse of qµ which is symmetric with respect to 0.

The phase portrait of the (K 0 µ ) is shown in Figure 7(a). Recall that in this phase portrait, each orbit is parametrized by x ∈ Ω = (-, ). By evaluating the energy levels of the potential (34), we easily prove that the Hamiltonian system (K 0 µ ) admits three critical points: (0, 0) is a center encircled by periodic orbits, and (0, -w + µ ), (0, w + µ ) are saddles, connected by two heteroclinic orbits. To join the two saddles, an orbit {v(x), w(x)} x∈Ω needs an infinite space Ω, which is excluded by assumption. In the region delimited by the heteroclinic orbits, a continuous band of periodic orbits encircles the center (0, 0). The space Ω = (-, ) necessary to join two successive points of the axis {v = 0} along one of these orbits varies increasingly from 0 (for the stationary point (0, 0)) to +∞ (for the heteroclinic orbits). Hence, for each bounded domain Ω = (-, ), a continuous orbit {v(x), w(x)} x∈Ω can always be found that fulfills condition [START_REF] Péron | Uncovering periodic patterns of space use in animal tracking data with periodograms, including a new algorithm for the Lomb-Scargle periodogram and improved randomization tests[END_REF]. Moreover, for such an orbit, the component w(x) is always monotone with respect to x ∈ Ω = (-, ) and joins a negative value w(-) to a positive value w( ). Finally, by setting u(x) = ψ 0 µ αw(x) , the pair {v(x), w(x)} x∈Ω determines a heterogeneous stationary solution {u(x), w(x)} x∈Ω . Since ψ 0 µ admits a discontinuity jump at s = 0, the component u(x) also admits a discontinuity jump. We obtain the following proposition.

Then the degenerate reaction-diffusion system (2) admits a discontinuous pattern {u(x), w(x)} x∈Ω such that w is continuous and monotone on Ω, satisfies w(-) = b, w( ) = a and u is piecewise continuous on Ω, with a unique discontinuity jump.

We show in Figure 7(b) a discontinuous pattern {u(x), w(x)} x∈Ω obtained for b = -a with 0 < a < w + µ : its component w(x) joins -a to a and its component u(x) admits a unique discontinuity jump at x = 0.

We emphasize that the discontinuous pattern {u(x), w(x)} x∈Ω whose existence is established by Proposition 3 has been obtained by choosing ψ 0 µ as an inverse q µ . It is thus natural to ask how a choice of another inverse ψ σ µ with σ = 0 affects the construction of the discontinuous pattern. We show in Figure 8 the phase portrait of the Hamiltonian system (33) obtained for ψ σ µ with σ > 0. We observe that the two saddles (0, -w + µ ) and (0, w + µ ) are not connected by heteroclinic orbits as in the case σ = 0. Now, only one saddle admits a homoclinic orbit that encircles a continuous band of periodic orbits. In this continuous band of periodic orbits, we can as previously extract small orbits {v(x), w(x)} x∈Ω that join two successive points of the axis {v = 0}. But if σ is sufficiently far from 0, then the existence of an orbit {v(x), w(x)} x∈Ω connecting b < 0 to a > 0 is compromised as soon as b and a are separated by the homoclinic orbit. For example, we have depicted in red in Figure 8 an orbit that joins -a to a, which exists for σ = 0, but not any longer for σ sufficiently far from 0. jump. Furthermore, for each ε > 0, one can find two distinct discontinuous patterns {u 0 (x),

Case of a small perturbation

We continue with the case of a perturbation of small intensity, by assuming that µ 2 < µ < µ 4 . We recall that in this case, the degenerate reaction-diffusion system (2) admits five homogeneous stationary solutions O, -U - µ , U - µ , -U + µ , U + µ . It turns out that this case exhibits rich dynamics. As previously, we aim to solve the Hamiltonian system [START_REF] Ritchie | Turing instability and pattern formation on directed networks[END_REF] and extract from its phase portrait orbits which satisfy condition [START_REF] Péron | Uncovering periodic patterns of space use in animal tracking data with periodograms, including a new algorithm for the Lomb-Scargle periodogram and improved randomization tests[END_REF]. To this end, we identify again a one-parameter family {ψ σ µ } of inverses of q µ , which now admit two discontinuity jumps. Lemma 3. For each µ ∈ (µ 2 , µ 4 ), there exist σ min < σ max such that the function q µ possesses a one-parameter family ψ σ µ σmin≤σ≤σmax of piecewise continuous and monotone inverses, admitting two discontinuity jumps. Furthermore, for each ε > 0, one can find σ ∈ (σ min , σ max ) such that σ = 0 and

We skip the proof Lemma 3, since it is very similar to that of Lemma 2. We only indicate that the values min and σ max are determined by the components -u - µ , u - µ of the homogeneous stationary solutions -U - µ , U - µ respectively. We show in Figure 9 two inverses ψ 0 µ and ψ σ µ with σ > 0. Note that in the family {ψ σ µ }, ψ 0 µ is the only inverse of q µ whose both discontinuity jumps coincide with -u - µ , u - µ .

-

Figure 9: Two possible piecewise continuous and monotone inverses ψ 0 µ and ψ σ µ of qµ when µ2 < µ < µ4. In the family {ψ σ µ }, ψ 0 µ is the only inverse of qµ whose both discontinuity jumps coincide with the components - µ , u - µ of the homogeneous stationary solutions -U - µ , U - µ respectively.

We can now establish the phase portrait of the Hamiltonian system (33) obtained with the particular inverse ψ 0 µ of q µ . As in [START_REF] Turing | The chemical basis of morphogenesis[END_REF], we denote again this system by (K 0 µ ). By evaluating the energy levels of the potential (34), we easily obtain its phase portrait, which is given in Figures 10(a), 11(a) and 11(b) for different values of µ. We observe that it presents five critical points: (0, -w - µ ) and (0, w - µ ) are centers, encircled by continuous bands of periodic orbits; (0, -w + µ ), (0, 0) and (0, w + µ ) are saddles, which appear under the effect of the pitchfork bifurcation described in Section 3.1. Now, it is seen that a bifurcation of heteroclinic orbits occurs as µ increases from µ 2 to µ 4 . Indeed, if µ is sufficiently near µ 2 , then the saddle (0, 0) admits two homoclinic orbits that encircle the centers (0, -w - µ ), (0, w - µ ), as shown in Figures 10(a) and 11(a). When µ increases, these homoclinic orbits grow, until µ reaches a critical value μ, where they become heteroclinic orbits connecting the saddle (0, 0) to the saddles (0, -w + µ ) and (0, w + µ ). When µ increases over μ, the heteroclinic orbits become homoclinic again, but now they are attached to the saddles (0, -w + µ ) and (0, w + µ ), as shown in Figure This bifurcation of heteroclinic orbits is also schematized in Figure 13. We obtain the following proposition.

Proposition 4 (Bifurcation of heteroclinic orbits). There exists μ ∈ (µ 2 , µ 4 ) such that the Hamitonian system (K 0 µ ) defined by [START_REF] Turing | The chemical basis of morphogenesis[END_REF] admits a bifurcation of heteroclinic orbits at μ. If µ ∈ (µ 2 , μ), then the Hamitonian system (K 0 µ ) admits two homoclinic orbits attached the saddle (0, 0). If µ ∈ (μ, µ 4 ), then the Hamitonian system (K 0 µ ) admits two homoclinic orbits attached to the saddles (0, -w + µ ) and (0, w + µ ).

Figure 10: (a) Phase portrait of the Hamiltonian system (K 0 µ ) defined by [START_REF] Turing | The chemical basis of morphogenesis[END_REF], for µ2 < µ < µ * (a). The saddle (0, 0) admits two homoclinic orbits that encircle the centers (0, -w - µ ), (0, w - µ ). If the domain Ω is sufficiently large, then for each a ∈ (0, w + µ ), an orbit {v(x), w(x)} x∈Ω connecting (0, -a) and (0, a) exists. (b) Each orbit {v(x), w(x)} x∈Ω connecting (0, -a) and (0, a) determines a discontinuous pattern {u(x), w(x)} x∈Ω such that u(x) admits two discontinuity jumps. Now that a bifurcation of heteroclinic orbits is identified to occur at μ, we aim to highlight another structural modification of the dynamics of the Hamiltonian system (K 0 µ ) that occurs between µ 2 and μ. Indeed, let a ∈ (0, w + µ ) and assume that we are looking again for an orbit of the Hamiltonian system (K 0 µ ) that connects the points (0, -a) and (0, a). If µ is sufficiently near µ 2 and if the domain Ω is sufficiently large, then we can prove that such an orbit {v(x), w(x)} x∈Ω exists.

Theorem 4 (Existence of particular discontinuous patterns). For each a ∈ (0, w + µ ), there exists µ * (a) ∈ (µ 2 , μ) such that, for each µ ∈ µ 2 , µ * (a) , the Hamiltonian system (K 0 µ ) defined by (38) admits an orbit {v(x), w(x)} x∈Ω connecting the points (0, -a) and (0, a), provided Ω is sufficiently large.

Furthermore, each orbit {v(x), w(x)} x∈Ω of the Hamiltonian system (K 0 µ ) connecting the points (0, -a) and (0, a) determines a discontinuous pattern {u(x), w(x)} x∈Ω such that u(x) admits two discontinuity jumps.

Proof. Let µ ∈ (µ 2 , μ). The Hamiltonian system (K 0 µ ) admits two homoclinic orbits (o), (o ) of small amplitude attached to the saddle (0, 0) (the homoclinic orbits (o) and (o ) are shown in Figure 10(a)). These homoclinic orbits are encircled by two heteroclinic orbits (h) and (h ) of large amplitude that connect the saddles (0, -w + µ ) and (0, w + µ ) (the heteroclinic orbits (h) and (h ) are also shown in Figure 10(a)). Between the homoclinic orbits (o), (o ) and the heteroclinic orbits (h), (h ), there exists a continuous band of periodic orbits {v(x), w(x)}. Let us denote by S (-a, a) the space which is necessary for such an orbit to connect two points (0, -a) to (0, a). Obviously, S (-a, a) varies continuously with a. Moreover, we have lim

µ , the orbit {v(x), w(x)} coincides with one of the large heteroclinic orbits that connect the saddles (0, -w + µ ) and (0, w + µ ). Now, we denote by h + µ the maximum of the small homoclinic orbit that is attached to the saddle (0, 0) in the half-plane {v > 0}. It is seen that h + µ tends to 0 as µ tends to µ 2 . Hence, one can find µ sufficiently near µ 2 such that a > h + µ . As before, we have lim

Hence, S (-a, a) is bounded by below between h + µ and w + µ . We introduce S min = min

assume that Ω is sufficiently so that its length satisfies |Ω| ≥ S min . an orbit {v(x), w(x)} x∈Ω exists that connects the points (0, -a) and (0, a).

Finally, we set u(x) = ψ 0 µ αw(x) , x ∈ Ω. Since ψ 0 µ admits two discontinuity jumps, by virtue of Lemma 3, we can conclude that u(x) also admits two discontinuity jumps. The proof is complete. [START_REF] Turing | The chemical basis of morphogenesis[END_REF], for µ * (a) < µ < μ. The saddle (0, 0) admits two homoclinic orbits that encircle the centers (0, -w - µ ), (0, w - µ ), but no orbit connects the points (0, -a) and (0, a). (b) Phase portrait of the Hamiltonian system (K 0 µ ) defined by [START_REF] Turing | The chemical basis of morphogenesis[END_REF], for μ < µ < µ4. Two homoclinic orbits that encircle the centers (0, -w - µ ), (0, w - µ ) are now attached to the saddles (0, -w + µ ) and (0, w + µ ). No orbit connects the points (0, -a) and (0, a).

The discontinuous patterns whose existence is established by Theorem 4 have been constructed with a particular choice for an inverse ψ 0 µ of of q µ . By choosing other inverses ψ σ µ with σ = 0, arbitrarily close to ψ 0 µ , we can prove, as in Theorem 3, the existence of a continuum of discontinuous patterns.

jump for very small values of µ, and bifurcate to discontinuous patterns admitting two jumps when µ increases. The numerical integration of the degenerate reaction-diffusion system (2) has been performed using a splitting numerical scheme of Strang type [START_REF] Strang | On the construction and comparison of difference schemes[END_REF], with α = 1, β = 1, δ = 10, = 125 and with the perturbation function given by [START_REF] Einav | Indirect diffusion effect in degenerate reaction-diffusion systems[END_REF].

Linear initial condition

We begin with numerical simulations of the degenerate reaction-diffusion system (2) stemming from a simple linear initial condition (u 0 , w 0 ) given by

We experiment several values of the perturbation parameter: µ ∈ {0, 0.5, 5}. The results are depicted in Figures 141516.

Figure 14: Numerical simulation of the degenerate reaction-diffusion system (2) stemming from the linear initial condition [START_REF] Yagi | Abstract parabolic evolution equations and their applications[END_REF], with µ = 0. The orbit is attracted a admitting a unique discontinuity jump.

For µ = 0, we observe that the orbit is attracted to a discontinuous pattern admitting a unique discontinuity jump (Figure 14). For µ = 0.5, the orbit is attracted to a discontinuous pattern admitting two discontinuity jumps (Figure 15). Finally, for µ = 5, the discontinuous pattern vanishes and the orbit is attracted to the trivial equilibrium (Figure 16).

Remark 5 (Invading ecosystem).

As discussed in Remark 3, these numerical simulations illustrate the possible spreading process of an invading ecosystem over an originally existing ecosystem. If the perturbation if small, both the initial and the new ecosystems coexist. But if the perturbation increases and overcomes a certain threshold, the invasion can be complete.

Periodic-linear initial condition

Oscillations in time and space are of particular interest for better understanding the dynamics of biological and ecological systems (see for instance Hence, is to experiment the behavior of the degenerate reaction-diffusion system (2) with initial conditions admitting such oscillations. Here, we present numerical simulations of the degenerate reaction-diffusion system (2) stemming from an initial condition (u 0 , w 0 ) such that u 0 is periodic and w 0 is linear:

with θ > 0. We test two values of the perturbation parameter: µ ∈ {0, 1}. The results are shown in Figures 17,18.

-

17: Numerical simulation of the degenerate reaction-diffusion system (2) stemming from the periodiclinear initial condition (41), with µ = 0. The initial oscillations of u0 are modified and the orbit converges to a form that admits several patterns, each a discontinuity For µ = 0, we observe that the initial oscillations of u 0 are modified. Several oscillations completely vanish, whereas others resist and progressively take the form of a discontinuous pattern (on a total of 8 initial oscillations, 5 resist and become a discontinuous pattern). Therefore, orbit converges to a form that admits several discontinuous patterns, each with a unique discontinuity jump (Figure 17). Next, for µ = 1, the discontinuous patterns are perturbed. Several patterns now exhibit two discontinuity jumps (Figure 18).

Periodic initial condition

Next, we present numerical simulations of the degenerate reaction-diffusion system (2) stemming from an initial condition (u 0 , w 0 ) such that u 0 and w 0 are periodic:

with θ > 0, θ > 0 such that θ = θ . In this way, u 0 and w 0 do not have the same phase nor the same amplitude. We experiment two values of the perturbation parameter: µ ∈ {0, 1}. The results are shown in Figures 19,20.

For µ = 0, we observe that the orbit converges to a form that admits several discontinuous patterns with a unique discontinuity jump, in which u µ and w µ are still synchronized, although u 0 and w 0 do

Figure 18: Numerical simulation of the degenerate reaction-diffusion system (2) stemming from the periodiclinear initial condition (41), with µ = 1. The orbit converges to a form that admits several discontinuous patterns, some of them with one discontinuity jump, some others with two discontinuity jumps.

-

19: Numerical simulation of the degenerate system (2) stemming from the periodic initial condition (42), with µ = 0. Although u0 and w0 do not have the same phase, the orbit converges to a form that admits several discontinuous patterns with a unique discontinuity jump in which uµ and wµ are synchronized. A transient pattern is observed at the centre of the domain and vanishes after a long period, whereas other patterns are persistent.

not have the same phase. Next, for µ = 1, the discontinuous patterns are still synchronized, but they now admit two discontinuity jumps.

The final number and location of these discontinuous patterns seem to be very difficult to predict, thus represents an interesting perspective for further research. • we have proved that the discontinuous two-jumps patterns vanish when the perturbation exceeds a certain threshold.

Several numerical simulations and animations have been provided to support our theoretical results. However, these new results were obtained under a restrictive assumption. Indeed, the existence of discontinuous patterns with two discontinuity jumps, proved in this paper (Theorems 4 and 5), was established by considering a strict subset of the huge set of inverses of the hysteresis process. Hence, it remains many choices of inverses to explore. Therefore, we believe that the study of this morphogenesis model will produce in a near future a great number of original results and will further reveal rich dynamics.