Emeric Bouin
email: bouin@ceremade.dauphine.fr

Guillaume Legendre
email: guillaume.legendre@ceremade.dauphine.fr

Yuan Lou
email: lou@math.ohio-state.edu

Nichole Slover
email: slover.9@osu.edu

EVOLUTION OF ANISOTROPIC DIFFUSION IN TWO-DIMENSIONAL HETEROGENEOUS ENVIRONMENTS

Keywords: Mathematics Subject Classification. 35K57, 92D15, 92D25. 1

We consider a system of two competing populations in two-dimensional heterogeneous environments. The populations are assumed to move horizontally and vertically with different probabilities, but are otherwise identical. We regard these probabilities as dispersal strategies. We show that the evolutionarily stable strategies are to move in one direction only. Our results predict that it is more beneficial for the species to choose the direction with smaller variation in the resource distribution. This finding seems to be in agreement with the classical results of Hasting [15] and Dockery et al. [11] for the evolution of slow dispersal, i.e. random diffusion is selected against in spatially heterogeneous environments. These conclusions also suggest that broader dispersal strategies should be considered regarding the movement in heterogeneous habitats.

1. Introduction 1.1. Background and motivation. In this paper, we consider populations of individuals that disperse in a bounded two-dimensional habitat, where the resources are distributed heterogeneously across the landscape. A natural question is how organisms should distribute themselves in the space to better match the available resources and, accordingly, what kind of dispersal strategies organisms should adopt to reach such distributions; see, e.g. [START_REF] Clobert | Dispersal[END_REF][START_REF] Fretwell | On territorial behavior and other factors influencing habitat selection in birds. I. Theoretical development[END_REF]. Most of previous studies on the evolution of dispersal assume that individuals move in two orthogonal directions with the same probability, which we refer as random dispersal; see, e.g. [START_REF] Hastings | Can spatial variation alone lead to selection for dispersal?[END_REF]. For spatially varying but temporally constant environments, Hastings considered a scenario in which the resident is randomly dispersing and reaches the equilibrium; some rare mutant, which is also randomly dispersing but differs from the resident only in the diffusion rate, is introduced to the habitat. He found that slower rates of dispersal will be selected, as the mixing of populations tends to reduce the growth [1,2]. Dockery et al. [START_REF] Dockery | The evolution of slow dispersal rates: a reaction diffusion model[END_REF] considered a system of two randomly diffusing competing populations in spatially varying but temporally constant environments, and two species are identical except their dispersal rates. They showed that the population with the smaller dispersal rate always drive the population with the larger dispersal rate to extinction, irrelevant of the initial data. This phenomenon is often termed as the evolution of slow dispersal, as any population with a positive dispersal rate will be replaced by a mutant with a smaller dispersal rate.

If we assume that individuals move, e.g. horizontally and vertically with two different probabilities (with the sum of probabilities equal to one) and regard these probabilities as dispersal strategies, what kind of strategies will be evolutionarily stable? Here, we are referring to the situation where the population moves east or west with probability p/2, and north or south with probability (1 -p)/2, respectively, with p being chosen in [0, 1]. Intuitively, random dispersal strategies, i.e. p = 1/2, might not be evolutionarily stable as the distribution of resources is generally not the same in the horizontal and vertical directions so that it could be more advantageous for the population to have a higher probability moving in one direction than the other. This is indeed the case, and one might attempt to conjecture that some mixed strategy p * in (0, 1), allowing the population to move in horizontal and vertical directions with different probabilities, would emerge as an evolutionarily stable dispersal strategy in this particular setting. A bit surprisingly, our results suggest that the only evolutionarily stable dispersal strategies are p * = 0 and/or p * = 1, i.e. it is more advantageous for the population to move in only one direction.

More specifically, we consider a system of two competing populations in two-dimensional heterogeneous environments. The populations are assumed to move horizontally and vertically with different probabilities, but are otherwise identical. We introduce a function F of the dispersal probability, which measures the difference between the spatial variations of the population distributions at equilibrium in horizontal and vertical directions: when it is positive, the species has more variations in the horizontal direction; when it is negative, it has more variations in the vertical direction. We show that F is monotone deceasing and that the evolutionarily stable dispersal strategies are to maximize the function F when it is positive and to minimize it when it is negative, i.e. the evolutionarily stable strategies are to move in one direction. As the population distribution is often positively correlated with the resource distribution, thus function F indirectly measures the difference between the resource variations in horizontal and vertical directions. Therefore, our results seem to predict that it is more favorable for the species to choose the direction with smaller variation in the resource distribution.

To explain these findings intuitively, consider a peculiar scenario in which the underlying habitat is a rectangular region and the resources are distributed inhomogeneously in the horizontal direction but homogeneously in the vertical one. For such case, as there is only spatial variations in the horizontal direction, the results of Hastings [START_REF] Hastings | Can spatial variation alone lead to selection for dispersal?[END_REF] and Dockery et al. [START_REF] Dockery | The evolution of slow dispersal rates: a reaction diffusion model[END_REF] for the evolution of slow dispersal suggest that it might be better for the population not to move horizontally, which is in agreement with our findings in this paper. These considerations also suggest that we should probably consider a broader set of dispersal strategies, e.g. those strategies which allow for condition-dependent movement [START_REF] Cantrell | Evolution of dispersal and the ideal free distribution[END_REF][START_REF] Cosner | Reaction-diffusion-advection models for the effects and evolution of dispersal[END_REF][START_REF] Gyllenberg | On the evolution of patch-type dependent immigration[END_REF][START_REF] Kisdi | Evolution of condition-dependent dispersal[END_REF][START_REF] Lam | Evolution of dispersal: evolutionarily stable strategies in spatial models[END_REF][START_REF] Lam | Evolutionarily stable and convergent stable strategies in reaction-diffusion models for conditional dispersal[END_REF][START_REF] Potapov | Evolutionarily stable diffusive dispersal[END_REF]. 1.2. Organization of the paper. Section 2 contains the formal derivations of the mathematical models at stage and the main results. We present numerical simulations in Section 3 to complement the analytical results and to provide some intuition and insights. In Section 4, we discuss the stability of semi-trivial equilibria and investigate properties of the invasion fitness. In Section 5, we further study the stability of the semi-trivial equilibria and identify all evolutionarily stable strategies. Section 6 is devoted to the classification of the global dynamics of the two-species competition model introduced in Section 2. In Section 7, we summarize our conclusions, and discuss possible extensions of the present work. Finally, some technical materials are given in the Appendix.

The models and analytical results

2.1. Formal derivation from random walks. In this section, we will closely follow the approach in [START_REF] Slover | The evolution of dispersal in two-dimensional habitats[END_REF]. Let the habitat be the discrete lattices of steps ∆x and ∆y in the full two-dimensional space R 2 . Assume that each individual moves horizontally to the left and to the right with probability θ 2 and vertically up and down with probability 1-θ 2 , with θ in (0, 1). Let N(t, x, y) denote the number of individuals of the population at time t and location (x, y) in R 2 . Then

N(t + ∆t, x, y) = θ 2 [N(t, x + ∆x, y) + N(t, x -∆x, y)] + 1 -θ 2 [N(t, x + ∆x, y) + N(t, x -∆x, y)] .
Using Taylor series expansions, we have

N(t + ∆t, x, y) -N(t, x, y) ∆t = (∆x) 2 2∆t θN xx + (∆y) 2 2∆t (1 -θ)N yy + (∆x) 2 ∆t • O(∆x) + (∆y) 2 ∆t • O(∆y),
where O(∆x) and O(∆y) denote terms which are bounded with respect to the ∆x and ∆y, respectively. Assuming that both (∆x) 2 2∆t → D and (∆y) 2 2∆t → D for some positive constant D as the lengths of the time step ∆t and of the space steps ∆x and ∆y tend to 0, we obtain, passing to the limit in the relation above, (2.1)

N t = D [θN xx + (1 -θ)N yy ] ,
which is the type of anisotropic diffusion operator to be considered in the present article.

The parameter θ can be regarded as a dispersal strategy. Namely, when θ = 0, the whole population will either move north or south with probability 1/2. Similarly, when θ = 1, the population will only move east or west. Most of previous studies assume that individuals are randomly diffusing, i.e. they move in two orthogonal directions with the same probability (θ = 1/2). Given arbitrarily distributed resources, is there some particular strategy θ in [0, 1] which can convey a competitive advantage? The main goal of the present article is to address this question.

2.2.

The single-species model. Incorporating the population dynamics into equation (2.1), we arrive at the following reaction-diffusion equation:

(2.2)      N t = D(θ)N xx + D(1 -θ)N yy + (a -N)N in Ω, t > 0, D(θ)N x ν x + D(1 -θ)N y ν y = 0 on ∂Ω, t > 0, N(0, •, •) = N 0 0 in Ω.
Here, the domain Ω is a smooth open bounded subset of R 2 , and we denote its boundary by ∂Ω. We assume without any further notice that Ω is strictly convex and ∂Ω is of class C 1 . The unit outward normal vector to ∂Ω is denoted by ν = (ν x , ν y ). Thus, the map ν is one-to-one and continuous from ∂Ω to S 1 .

For any θ in [0, 1], the function D is defined by

D(θ) := D + (D -D)θ, θ ∈ [0, 1].
Note that if D = 0, D(θ) = Dθ is reduced to the form of diffusion in equation (2.1), which is degenerate when θ = 0. To avoid such degeneracy, we assume in the remainder of the paper that D and D are positive constants satisfying

0 < D < D.
The equation has been completed with zero flux boundary conditions so that no individuals may escape the domain. The model is not mass conservative since individuals may reproduce according to monostable non-linearities. The free growth rate is given by the heterogeneous function a, that is assumed to satisfy the following assumption: (A1) Function (x, y) → a(x, y) is positive, Hölder continuous and non-constant in Ω.

By standard regularity theory for parabolic equations and a comparison argument, see, e.g. [START_REF] Cantrell | Spatial ecology via reaction-diffusion equations[END_REF], it can be shown that N is positive in Ω for all times and that N(t, •, •) tends to N θ uniformly in Ω as t tends to infinity, where N θ , the equilibrium distribution of the population, is the unique positive steady state of system (2.2), i.e. N θ satisfies

(2.3) D(θ)(N θ ) xx + D(1 -θ)(N θ ) yy + (a -N θ )N θ = 0 in Ω, (D(θ)(N θ ) x , D(1 -θ)(N θ ) y ) • ν = 0 on ∂Ω.
Note that when the function a is non-constant, so is N θ .

Clearly, the parameter θ has a strong influence on N θ . As θ increases, the single species has more tendency to move horizontally than vertically, which may reduce the spatial variations of the population distributions in the horizontal direction and increase the variations in the vertical direction. In this connection, we have the following result: Theorem 2.1. Define, for θ in [0, 1], the function

(2.4) F (θ) := Ω ((N θ ) x ) 2 -((N θ ) y ) 2 dx dy.
Then F ′ (θ) < 0 for θ ∈ (0, 1). In particular, F is either strictly positive, strictly negative or sign-changing exactly once in (0, 1).

The function F , which plays a critical role in later analysis, can be regarded as a measurement of the difference between the variations of the population distributions in horizontal and vertical directions: when F is positive, we envision that the species at equilibrium has more spatial variations horizontally; when F is negative, it has more variations in the vertical direction. Theorem 2.1 implies that as the species increases the horizontal diffusion and reduces the vertical diffusion, then it tends to have more variations in the vertical direction than the horizontal direction.

As the population distribution is often positively correlated with the resource distribution, function F can also be viewed as an indirect measurement of the difference between the resource variations in the horizontal and vertical directions. Numerical results on the shape of the function F are presented in the next section (see Figure 1).

2.3.

The two-species competition model. Given arbitrarily distributed resources across the habitat, we may regard the parameter θ as a dispersal strategy and ask whether there is some value for θ which is evolutionarily stable. To address this question, we now move to the situation where two populations are competing for the same resources but adopt different dispersal strategies. We thus consider the following reaction-diffusion system for two competing species:

(2.5)

               U t = D(p)U xx + D(1 -p)U yy + (a -U -V )U in Ω, t > 0, V t = D(q)V xx + D(1 -q)V yy + (a -U -V )V in Ω, t > 0, (D(p)U x , D(1 -p)U y ) • ν = 0 on ∂Ω, t > 0, (D(q)V x , D(1 -q)V y ) • ν = 0 on ∂Ω, t > 0, U(0, •, •) = U 0 0, V (0, •, •) = V 0 0 in Ω,
in which the functions U and V represent the respective population densities of two competing species. By standard regularity theory and the maximum principle for parabolic equations, it can be shown that U and V are positive in Ω for all times. The competition for resources is neutral and independent of the dispersal strategy of the individuals, so that the death rate is given by U + V for both populations.

We are given two orthogonal space directions e 1 = (1, 0) and e 2 = (0, 1), so that (x, y) are the Cartesian coordinates in this basis. We may emphasise that after that choice, the problem is not rotationally invariant. As such, the two populations disperse with their own dispersal strategies, assimilated to the respective probabilities p and q to move in the direction e 1 , with p and q chosen in [0, 1]. As formally explained in Subsection 2.1, this way of dispersing results in a diffusion coefficient given by D(p) (resp. D(q)) in the direction e 1 and D(1 -p) (resp. D(1 -q)) in the direction e 2 for the first (resp. second) density.

We will adopt the viewpoint in the theory of adaptive dynamics. An important concept in adaptive dynamics is that of evolutionarily stable strategies (ESS). A strategy is said to be evolutionarily stable if a population using it cannot be invaded by any small population using a different strategy. In system (2.5), p and q are strategies for two populations. In terms of adaptive dynamics, we say that p in [0, 1] is an ESS if the semi-trivial steady state (N p , 0) is locally asymptotically stable for q = p, with q in [0, 1] and q close to p.

The following result characterizes the local stability of (N p 0) for p and q in [0, 1]. Theorem 2.2. There exists some continuous function q = q * (p), defined in [0, 1], satisfying 0 ≤ q * (p) ≤ 1 such that the following statements hold.

(i) If F is positive in [0, 1], then q * (p) > p and q * (p) ≡ 1 for p close to 1 such that (N p , 0) is stable for p < q < q * (p), unstable for q > q * (p) and q < p. In particular, if q * (p) ≡ 1, then (N p , 0) is stable for q > p, unstable for q < p.

(ii) If F has a unique root θ * in (0, 1), then q * (p) > p for p ∈ [0, θ * ) and q * (p) < p for p ∈ (θ * , 1], such that (N p , 0) is stable for min{q * (p), p} < q < max{q * (p), p}, unstable for q > max{q * (p), p} and q < min{q * (p), p}.

(iii) If F is negative in [0, 1], then q * (p) < p and q * (p) ≡ 0 for p close to 0 such that (N p , 0) is stable for q * (p) < q < p, and unstable for q < q * (p) and q > p. In particular, if q * (p) ≡ 0, then (N p , 0) is stable for q < p and unstable for q > p.

This result follows from Theorems 5.2, 5.4, and 5.6. In the next section, numerical results shed some insight into the stability of (N p , 0) and illustrate the conclusions of Theorem 2.2 (see Figure 5). Some biological intuition can also be gained from this Theorem as it provides a criterion for finding the ESS of system (2.5).

Corollary 2.3. The following conclusions hold.

(i) If F is positive in [0, 1], then p = 0 is the only ESS.

(ii) If F has exactly one root θ * in (0, 1) so that F is positive in [0, θ * ) and negative in (θ * , 1], then both p = 0 and p = 1 are ESS, and θ * is not evolutionarily stable.

(iii) If F is negative in [0, 1], and p = 1 is the only ESS.

Our remaining goals include understanding the global dynamics of system (2.5). This system possesses two semi-trivial steady states, given by (N p , 0) and (0, N q ), respectively. Theorem 2.2 addresses the local stability of (N p , 0) for arbitrary values of p and q, and the stability of (0, N q ) can be similarly determined. Furthermore, we shall show that there are only three alternatives for the global dynamics of system (2.5):

(i) the state (N p , 0) is globally stable;

(ii) the state (0, N q ) is globally stable;

(iii) the states (N p , 0) and (0, N q ) are both unstable, and there exists a unique positive steady state which is globally stable. We refer to the statements of Theorems 6.4, 6.5 and 6.6 for further details on the characterizations of the global dynamics of system (2.5). These analytical results on the dynamics of (2.5), complemented by numerical simulations in the next section for a free growth rate function of the form a(x, y) = λA(x) + (1 -λ)A(y), with λ in [0, 1], will help provide a more clear picture on the dynamics of system (2.5).

The numerical results

All the simulations presented here were achieved using the free and open-source software FreeFEM [START_REF] Hecht | New development in FreeFem++[END_REF]. The numerical approximation of the large-time solution to system (2.2) was based on a variational form of the problem and achieved using a spatial discretisation based on the finite element method, with P 1 Lagrange elements, combined with an implicit-explicit (IMEX) Euler scheme (see [START_REF] Asher | Implicit-explicit methods for time-dependent partial differential equations[END_REF] for instance) for the time-integration of the resulting ordinary differential equations. The linear terms in the reaction-diffusion equation are then treated implicitly in time, while the non-linear reaction term is dealt with explicitly, in order to enforce the stability of the scheme.

The mesh used to discretise the domain Ω realised as a disk of radius 2 was comprised of 3916 triangles and the length of the time step used was 0.01. The chosen initial state N 0 is the constant one, with value 0.5. Once the stationarity of the approximate solution was obtained in relative L 2 norm within a prescribed tolerance of 1e-15, the approximate steady state was used to compute an approximate value of F (p) and also for a finite element discretisation of the linear eigenvalue problem (4.1). The computation of an approximation of the smallest eigenvalue of (4.1), denoted by Λ(p, q), was done with the ARPACK package. Representations of an approximation to the nodal set of Λ(p, q) for different values of λ were then obtained by repeating the computation for numerous values of the parameters p and q taken in a discrete grid of the interval [0, 1].

Note that the state (N p , 0) is linearly stable when Λ(p, q) is positive, and unstable when Λ(p, q) is negative. Furthermore, Λ(p, q) vanishes whenever p = q, i.e. the nodal sets of Λ(p, q) always consist of the diagonal line p = q in the p -q plane. 3.1. The function F . Numerical approximations of the graph of the function F in the case of a free-growth function of the form a(x, y) = λA(x) + (1 -λ)A(y) are provided in Figure 1. For the simulations, we considered a disk of radius 2 centered at the origin for the domain Ω, anisotropic diffusion parameters D and D respectively equal to 0.1 and 10, and the function A(x) = 2 -sin(πx).

It is easily seen that, for all λ in [0, 1] and all θ in [0, 1], the value F (θ) for λ is equal to the value of -F (1 -θ) for 1 -λ. Due to this symmetry in the function F with respect to λ, we only plot the graph of F for different values of λ between 0 and 0.5, illustrating how the function goes from strictly negative, to sign-changing once, and to strictly positive as λ varies. As shown in Figure 1, the function F is strictly decreasing in θ, as predicted by Theorem 2.1.

Figures 2 and3 present the numerical approximations of the function a and of the steady state N θ for various values of θ in [0, 1] and λ respectively equal to 0.4 and 0.6, in the case where Ω is a disk of radius 2 centered at the origin, D = 0.1, D = 10, and A(x) = 2 -sin(πx). For such values of λ, the function does not appear more biased in the horizontal direction than in the vertical one. Yet, one can see clearly that the steady state shows very little variation in vertical direction for θ = 0, but changes as the value of the parameter θ grows to end with very little variation in the horizontal direction for θ = 1. This illustrates how the function F goes from being negative to positive as θ varies.

In Figure 4, we fix a value of θ and observe that F (θ), seen as a function of λ, is not necessarily monotone. In this case, we considered a disk of radius 2 centered at the origin for the domain Ω, anisotropic diffusion parameters D and D respectively equal to 0.1 and 10, and the function (c) The steady state for θ = 0.2. (e) The steady state for θ = 0.6. (e) The steady state for θ = 0.6. (f) The steady state for θ = 0.8. 3.2. Local stability of (N p , 0). Figure 5 presents numerical approximations of the nodal sets of Λ(p, q) for a free growth function of the form a(x, y) = λA(x)+(1-λ)A(y), the choice of the problem parameters being the same as for the graphs of F in Figure 1.

A(x) = 4 -1 4 x 2 .
Here, Λ(p, q) is positive if and only if (N p , 0) is linearly stable and Λ(p, q) < 0 if and only if (0, N q ) is linearly unstable. Again, for all λ in [0, 1] and all (p, q) in [0, 1] 2 , the value of Λ(p, q) for λ is equal to the value of Λ(1 -p, 1 -q) for 1 -λ. Due to this symmetry of Λ(p, q) with respect to λ, we only plot the nodal sets for values of λ between 0 and 0.5, providing a graphical illustration of how the nodal set of Λ(p, q) changes as λ varies, in connection with the results of Theorem 2.2, or more broadly, the conclusions of Theorems 5.2, 5.4, and 5.6.

Subfigure 5(a) corresponds to statement (iii) of Theorem 2.2 with q * (p) ≡ 0, which shows that (N p , 0) is stable for p > q and unstable for p < q. In particular, p = 1 is the only ESS. In this case, the nodal set of Λ(p, q) consists precisely of the diagonal line q = p. These conclusions are consistent with Subfigure 1(a), in which F is strictly negative.

Subfigures 5(b) to 5(f) correspond to statement (ii) of Theorem 2.2, which shows that (N p , 0) is stable for min{q * (p), p} < q < max{q * (p), p}, unstable for max{q * (p), p} < q and q < min{q * (p), p}, and the nodal set of Λ(p, q) consists of the curve q = q * (p) and the diagonal line q = p. Hence, both p = 0 and p = 1 are the only ESS. These conclusions are in accordance with Subfigures 1(b) to 1(f), in which F changes sign exactly once in (0, 1), from positive to negative.

3.3. Global dynamics. Figure 6 presents numerical approximations of the nodal sets of Λ(p, q) and Λ(q, p) for a free growth function of the form a(x, y) = λA(x) + (1λ)A(y), which correspond to the stablity of semi-trival steady states (N p , 0) and (0, N q ), respectively. The domain Ω is again a disk of radius 2 centered at the origin, the anisotropic diffusion parameters D and D are respectively equal to 0.1 and 10, but the function A is now A(x) = 4 -1 4 x 2 . Due to a symmetry of Λ(p, q) and Λ(q, p) with respect to λ, we only plot the nodal sets for values of λ between 0 and 0.5.

For Figure 6, in the green colored region (N p , 0) is stable and (0, N q ) is unstable. By Theorems 6.4, 6.5 and 6.6, (N p , 0) is globally stable for (p, q) in the green region. Similarly, the one colored in red is the one in which (0, N q ) is globally stable. The white region is where both (N p , 0) and (0, N q ) are unstable, and there is a unique positive steady state which is globally stable. The white regions for λ = 0.4 and λ = 0.5 are substantially greater than those for smaller values of λ. Biologically, this suggests that if the spatial variations of the resource distribution in the vertical and horizontal directions become more comparable, the chances for the coexistence of the two competing populations could be greater.

The invasion fitness

In this section we consider the stability of (N p , 0), study some properties of the invasion fitness and establish Theorem 2.1. . Numerical approximations of the nodal sets of Λ(p, q) for λ taking the values 0, 0.1, 0.2, 0.3, 0.4, and 0.5, in the case where Ω is a disk of radius 2 centered at the origin, D = 0.1, D = 10, and A(x) = 2 -sin(πx). The subset colored in green is the one in which Λ(p, q) > 0, that is in which the steady state (N p , 0) is linearly stable. . Numerical approximations of the nodal sets of Λ(p, q) and Λ(q, p) for λ taking the values 0, 0.1, 0.2, 0.3, 0.4, and 0.5, in the case where Ω is a disk of radius 2 centered at the origin, D = 0.1, D = 10, and A(x) = 4 -1 4 x 2 . The subset colored in green is the one in which Λ(p, q) > 0 > Λ(q, p), i.e. in which (N p , 0) is stable, while that colored in red is the one in which Λ(q, p) > 0 > Λ(q, p), i.e. in which (0, N q ) is stable. The white region is where both (N p , 0) and (0, N q ) are unstable.

The linear stability of (N p , 0) is determined by the sign of the smallest eigenvalue, denoted by Λ := Λ(p, q), of the linear eigenvalue problem (4.1)

D(q)ϕ xx + D(1 -q)ϕ yy + (a -N p )ϕ + Λ(p, q)ϕ = 0 in Ω, (D(q)ϕ x , D(1 -q)ϕ y ) • ν = 0 on ∂Ω,
as seen in the following result.

Lemma 4.1. The semi-trivial equilibrium (N p , 0) is linearly stable if Λ(p, q) is positive and unstable if Λ(p, q) is negative. Similarly, the semi-trivial equilibrium (0, N q ) is linearly stable if Λ(q, p) is positive and unstable if Λ(q, p) is negative.

The proof of Lemma 4.1 is the same as that of Lemma 5.5 in [START_REF] Chen | Evolution of conditional dispersal: a reaction-diffusionadvection model[END_REF] and is thus omitted.

In the theory of adaptive dynamics [START_REF] Dieckmann | The dynamical theory of coevolution: a derivation from stochastic ecological processes[END_REF][START_REF] Diekmann | A beginner's guide to adaptive dynamics[END_REF][START_REF] Geritz | Evolutionarily singular strategies and the adaptive growth and branching of the evolutionary tree[END_REF], Λ(p, q) is termed as the invasion fitness or invasion exponent, which can be regarded as the payoff function for the mutant phenotype with trait q, when the resident phenotype with trait p is at the equilibrium. Namely, if Λ(p, q) is positive, the mutant with trait q can invade when rare; on the contrary, when Λ(p, q) is negative, the mutant with trait q cannot invade when rare. We shall now give some properties of the map (p, q) → Λ(p, q) that will be used later to describe more precisely the stability of both semi-trivial steady states.

4.1. The selection gradient. If p = q, that is when both phenotypes are identical, Λ(p, p) = 0 for any p in [0, 1], thus both semi-trivial steady states (N p , 0) and (0, N q ) are neutrally stable. In this section, we consider the stability of (N p , 0) for p and q sufficiently close to each other. The following result provides a criterion in determining the sign of Λ(p, q) in such case (see also [START_REF] Slover | The evolution of dispersal in two-dimensional habitats[END_REF]).

Lemma 4.2. There holds

∂Λ ∂q q=p = (D -D) Ω (N p ) 2 F (p)
, where F is the function defined by (2.4).

Proof. Consider the positive eigenfunction ϕ associated with Λ(p, q) and uniquely determined by Ω ϕ 2 = Ω (N p ) 2 . It is a smooth function of p and q, and, for simplicity of notation, we denote ∂ϕ ∂q by ϕ ′ . Differentiating system (4.1) with respect to q, we obtain (4.2)

         D(q)ϕ ′ xx + D(1 -q)ϕ ′ yy + (a -N p )ϕ ′ + Λ(p, q)ϕ ′ + ∂Λ ∂q (p, q)ϕ +(D -D)(ϕ xx -ϕ yy ) = 0 in Ω, (D(q)ϕ ′ x + (D -D)ϕ x , D(1 -q)ϕ ′ y -(D -D)ϕ y ) • ν = 0 on ∂Ω.
Multiplying the first equation in system (4.1) by ϕ ′ , integrating by parts the result over Ω and using the second equation in system (4.1), we have

- Ω (D(q)ϕ x ϕ ′ x + D(1 -q)ϕ y ϕ ′ y ) + Ω (a -N p )ϕϕ ′ + Λ(p, q) Ω ϕϕ ′ = 0.
Similarly, multiplying the first equation in system (4.2) by ϕ, integrating by parts the result over Ω and using the second equation in system (4.2), we obtain

- Ω (D(q)ϕ ′ x ϕ x + D(1 -q)ϕ ′ y ϕ y ) + Ω (a -N p )ϕ ′ ϕ -(D -D) Ω ((ϕ x ) 2 -(ϕ y ) 2 ) +Λ(p, q) Ω ϕ ′ ϕ + ∂Λ ∂q (p, q) Ω ϕ 2 = 0.
Subtracting the above two equalities then yields

(4.3) ∂Λ ∂q (p, q) Ω ϕ 2 = (D -D) Ω ((ϕ x ) 2 -(ϕ y ) 2 ).
Since Λ(p, p) = 0, it follows from the normalization of ϕ that ϕ |q=p = N p , which completes the proof.

In view of Lemma 4.2, it is critical to understand the sign of function F . The following result establishes Theorem 2.1. Lemma 4.3. For any θ in [0, 1], one has F ′ (θ) < 0.

Proof. For simplicity of notation, denote ∂N θ ∂θ by N ′ θ . Integrating by parts, one gets where the last equality follows from the boundary condition satisfied by N ′ θ . We therefore have

F ′ (θ) = - 2 D -D Ω [D(θ)((N ′ θ ) x ) 2 + D(1 -θ)((N ′ θ ) y ) 2 -(a -2N θ )(N ′ θ ) 2 ].
Let λ 1 denote the smallest eigenvalue of the linear problem

D(θ)ϕ xx + D(1 -θ)ϕ yy + (a -N θ )ϕ + λϕ = 0 in Ω, (D(θ)ϕ x , D(1 -θ)ϕ y ) • ν = 0 on ∂Ω.
It is well-known that λ 1 can be characterized by the variational formula (4.4)

λ 1 = inf ϕ =0, ϕ∈H 1 (Ω) Ω [D(θ)(ϕ x ) 2 + D(1 -θ)(ϕ y ) 2 -(a -N θ )ϕ 2 ] Ω ϕ 2 ,
and that λ 1 is the only eigenvalue such that its corresponding eigenfunction does not change sign in Ω. Using system (2.3) for N θ , we infer that λ 1 = 0 and that its corresponding eigenfunction is a scalar multiple of N θ . In particular, by choosing the test function ϕ = N ′ θ in (4.4), we have that

Ω [D(θ)((N ′ θ ) x ) 2 + D(1 -θ)((N ′ θ ) y ) 2 -(a -N θ )(N ′ θ ) 2 ] ≥ 0,
which yields

F ′ (θ) ≤ - 2 D -D Ω N θ (N ′ θ ) 2 ≤ 0.
This gives F ′ ≤ 0, with the equality if and only if

N ′ θ ≡ 0 in Ω. Finally, if N ′ θ ≡ 0 for some θ, then N θ satisfies (N θ ) xx -(N θ ) yy = 0 in Ω, ((N θ ) x , -(N θ ) y ) • ν = 0 on ∂Ω.
This, together with the boundary condition satisfied by N θ , implies that

(N θ ) x ν x = (N θ ) y ν y = 0 on ∂Ω.
According to Lemma 8.1 in the Appendix, N θ is then a positive constant function, and it follows from the first equation in (2.3) that the function a(x, y) is also a constant function, which contradicts assumption (A1). Consequently, one has

F ′ < 0 in [0, 1].
Remark 1. If we allow Lipschitz domains with flat parts on the boundary, it is possible to construct domains such that F ≡ 0 (and thus F ′ ≡ 0) in [0, 1]. See the Appendix 8.2 for further discussions.

Concavity of Λ(p, q)

. The following result concerns the concavity of Λ(p, q) with respect to q.

Lemma 4.4. For any p in [0, 1], the function q → Λ(p, q) is concave on [0, 1]. Moreover, if Λ q (p * , q * ) = 0 for some couple (p * , q * ), then Λ qq (p * , q * ) < 0.

Proof. Let us fix p in [0, 1]. The concavity of q → Λ(p, q) follows from a standard argument based on the variational characterization of Λ(p, q), see [START_REF] Ni | The mathematics of diffusion[END_REF]. Nevertheless, we include here a proof of this result in order to facilitate the proof of the second statement of the Lemma. Consider the positive eigenfunction ϕ associated with Λ(p, q) such that Ω ϕ 2 = Ω (N p ) 2 . Differentiating (4.3) with respect to q yields Λ qq (p, q)

Ω ϕ 2 + 2Λ q (p, q) Ω ϕ ′ ϕ = 2(D -D) Ω (ϕ ′ x ϕ x -ϕ ′ y ϕ y ),
with the notations ϕ ′ = ∂ϕ ∂q , ϕ ′ x = ∂ϕx ∂q , and ϕ ′ y = ∂ϕy ∂q . Note that the second term in the left-hand side vanishes due to the normalization condition on ϕ, which implies that

Ω ϕ ′ ϕ = 0.
Multiplying the first equation in system (4.2) by ϕ ′ and integrating by parts the result over Ω then gives

(D -D) Ω (ϕ x ϕ ′ x -ϕ y ϕ ′ y ) = - Ω D(q)(ϕ ′ x ) 2 + D(1 -q)(ϕ ′ y ) 2 -(a -N p )(ϕ ′ ) 2 -Λ(p, q)(ϕ ′ ) 2 ≤ 0,
where we have used the second equation in system (4.2) and the fact that Ω ϕ ′ ϕ = 0, the inequality following from the variational characterization of Λ(p, q),

Λ(p, q) = inf ϕ =0, ϕ∈H 1 (Ω) Ω [D(q)(ϕ x ) 2 + D(1 -q)(ϕ y ) 2 -(a -N p )ϕ 2 ] Ω ϕ 2 .
It then holds that Λ qq (p, q) ≤ 0, where the equality holds if and only if the function ϕ ′ is a scalar multiple of ϕ. Since ϕ is positive in Ω, this implies that ϕ ′ ≡ 0 in Ω. Thus, if Λ q (p, q) = Λ qq (p, q) = 0, system (4.2) reduces to ϕ xx -ϕ yy = 0 in Ω, and (ϕ x , -ϕ y ) • ν = 0 on ∂Ω.

Due to the boundary condition for ϕ in system (4.1), we further have ϕ x ν x = ϕ y ν y = 0 on ∂Ω. As a consequence of Lemma 8.1, the function ϕ is constant, which implies, using the first equation in system (4.1), that the function a -N p is also constant. Integrating over Ω the main equation in system (2.3) with θ = p, we then obtain

Ω (a -N p )N p = 0,
which ensures that a -N p = 0 in Ω, so that system (2.3) for N p reduces to

D(p)(N p ) xx + D(1 -p)(N p ) yy = 0 in Ω, (D(p)(N p ) x , D(1 -p)(N p ) y ) • ν = 0 on ∂Ω.
It finally follows from the maximum principle that the function N p is constant. This contradicts the assumption of the function a being non-constant.

Lemma 4.5. If Λ(p, q) = 0 for some couple (p, q) such that p = q, then Λ q (p, q) = 0.

Proof. We argue by contradiction. Suppose that there exists a couple (p, q) such that p = q for which Λ(p, q) = Λ q (p, q) = 0. By Lemma 4.4, one has Λ qq (p, q) < 0, which implies that the function q → Λ(p, q) has a local maximum point at q = q. This contradicts the fact that Λ(p, p) = 0 and the concavity of Λ(p, q) in q.

A consequence of Lemma 4.5 is that the nodal set of Λ(p, q) within the region {(p, q) : 0 ≤ p, q ≤ 1} can be parameterized by a function q = q(p). However, we caution the reader that the domain of this function can be either [0, 1] or a strict subset of it.

5.

Local stability for (p, q) in [0, 1] 2 : nodal set of invasion fitness

We have previously considered the stability of the semi-trivial steady state (N p , 0) for p and q close to each other. In this section, we study the local stability of (N p , 0) for general p and q in [0, 1]. This is equivalent to giving a description of the nodal set of Λ(p, q) for p and q in [0, 1], and as well the sets where Λ(p, q) is positive or negative. By Lemma 4.3, the function F is strictly decreasing in [0, 1], so it suffices to consider three generic cases: F (θ) > 0 for all θ in [0, 1], F changes sign exactly once in (0, 1), and F (θ) < 0 for all θ in [0, 1].

5.1.

Stability when F (θ) > 0. First, we consider the case for which the function F is positive in [0, 1). Lemma 5.1. Suppose that F (θ) > 0 for all θ in [0, 1). If 0 ≤ p < q ≤ 1, then the steady state (0, N q ) is unstable. Conversely, if 0 ≤ q < p ≤ 1, then the steady state (N p , 0) is unstable.

Proof. We first consider the situation 0 ≤ p < q < 1. The stability of the steady state (0, N q ) is determined by the sign of the smallest eigenvalue, denoted by Λ(q, p), of the linear eigenproblem

D(p)ϕ xx + D(1 -p)ϕ yy + (a -N q )ϕ + λϕ = 0 in Ω, (D(p)ϕ x , D(1 -p)ϕ y ) • ν = 0 on ∂Ω.
By the variational characterization of Λ(q, p) we have

Λ(q, p) = inf ϕ =0, ϕ∈H 1 (Ω) Ω D(p)ϕ 2 x + D(1 -p)ϕ 2 y -(a -N q )ϕ 2 Ω ϕ 2 ≤ Ω [D(p)((N q ) x ) 2 + D(1 -p)((N q ) y ) 2 -(a -N q )(N q ) 2 ] Ω (N q ) 2
Multiplying the main equation in system (2.3) with θ = q by N q and integrating by parts the result over Ω yields

Ω D(q) ((N q ) x ) 2 + D(1 -q) ((N q ) y ) 2 -(a -N q )(N q ) 2 = 0.
Therefore, we have, using the assumptions that p < q and F > 0 in [0, 1],

Λ(q, p) ≤ [D(p) -D(q)] Ω ((N q ) x ) 2 + [D(1 -p) -D(1 -q)] Ω ((N q ) y ) 2 Ω (N q ) 2 = (D -D)(p -q) Ω ((N q ) x ) 2 -Ω ((N q ) y ) 2 Ω (N q ) 2 = (D -D)(p -q) F (q)
Ω (N q ) 2 < 0. If q = 1, we note that, as Λ(p, q) < 0 for 0 ≤ p < q < 1, Λ(p, 1) ≤ 0. Since Λ(1, 1) = 0, we see by Lemma 4.5 that Λ(p, 1) < 0 for all 0 ≤ p < 1. Hence, the steady state (0, N q ) is unstable for 0 ≤ q < p ≤ 1.

Similarly, we can show that if 0 ≤ q < p ≤ 1, then (N p , 0) is unstable.

The first main result of this section is the following.

Theorem 5.2. Suppose that the function F is positive in [0, 1]. Then, there exists some positive, continuous function q = q * (p), defined in [0, 1], satisfying p < q * (p) ≤ 1 for all p in [0, 1] and q * (p) ≡ 1 for p close to 1, such that (i) Λ(p, q) > 0 for 0 ≤ p < q < q * (p), (ii) Λ(p, q) < 0 for q * (p) < q ≤ 1 and 0 ≤ q < p ≤ 1. In particular, if q * (p) ≡ 1, then Λ(p, q) > 0 for p < q ≤ 1 and Λ(p, q) < 0 for 0 ≤ q < p.

Proof. Since Λ(p, p) = 0 for all p in [0, 1] and F is positive in [0, 1], there exists some positive real number δ such that Λ(p, q) > 0 for p and q in [0, 1] with 0 < q -p < δ.

Given any p in [0, 1], if Λ(p, 1) < 0, by Lemma 4.4 and the positivity of Λ(p, q) in the strip 0 < q -p < δ, there exists a unique q * = q * (p) in (p, 1) such that Λ(p, q * ) = 0, Λ(p, q) > 0 for p < q < q * and Λ(p, q) < 0 for q > q * . If Λ(p, 1) ≥ 0, we define q * (p) = 1. Again by Lemma 4.4, Λ(p, q) > 0 holds for p < q < q * = 1. This establishes statements (i) and (ii).

To show that q = q * (p) is a continuous curve, suppose that Λ(p, q) = 0 for some p < q < 1. By Lemma 4.5, there exists a smooth curve q = q * * (p) passing through (p, q) and such that Λ(p, q * * (p)) = 0, which can be extended to the left and to right until it reaches either p = 0 or q = 1, as by the choice of δ, q = q * * (p) and its extension can never enter the strip 0 < q -p < δ. For each p, there exists at most one q > p such that Λ(p, q) = 0.

Hence, q * * ≡ q * as long as these functions are strictly less than 1. Therefore, q * defines a continuous curve on [0, 1]. By the choice of δ, we see that q * (p) ≡ 1 for p close to 1.

Similarly to Theorem 5.2, the local stability of the semi-trivial steady state (0, N q ) can be determined as follows.

Theorem 5.3. Suppose that the function F is positive in [0, 1]. Then, there exists some positive, continuous function p = p * (q), defined in [0, 1], satisfying q < p * (q) ≤ 1 for q in [0, 1] and p * (q) ≡ 1 for q close to 1, such that Λ(q, p) > 0 for q < p < p * (q) and Λ(q, p) < 0 for p * (q) < p ≤ 1 and q < p.

The proof being the same as that for Theorem 5.2, we omit it. 5.2. Stability when F (θ) < 0. Next, we consider the case for which the function F is negative in [0, 1]. This case is similar to the previous one, so that we may state the following results without proofs.

Theorem 5.4. Suppose that the function F is negative in [0, 1]. Then, there exists a positive continuous function q = q * (p), defined in [0, 1], satisfying 0 ≤ q * (p) < p for p in [0, 1] and q * (p) ≡ 0 for p close to 0 such that (i) Λ(p, q) > 0 for q * (p) < q < p ≤ 1, (ii) Λ(p, q) < 0 for 0 ≤ q < q * (p) and 0 ≤ p < q ≤ 1. In particular, if q * (p) ≡ 0, then Λ(p, q) > 0 for q < p and Λ(p, q) < 0 for q > p.

Theorem 5.5. Suppose that the function F is negative in [0, 1]. Then, there exists a positive continuous function p = p * (q), defined in [0, 1], satisfying p * (q) < q for q in [0, 1] and p * (q) ≡ 0 for q close to 0, such that Λ(q, p) > 0 for p * (q) < p < q and Λ(q, p) < 0 for 0 ≤ p < p * (q) and q < p.

5.3.

Stability when F changes sign once. We finally consider the case for which the function F possesses a unique root in (0, 1), denoted by θ * . This function being decreasing, this implies that it is positive in [0, θ * ) and negative in (θ * , 1]. Theorem 5.6. Suppose that the function F has a unique root θ * in (0, 1). Then, there exists a positive continuous function q = q * (p), defined in [0, 1], satisfying p < q * (p) ≤ 1 for p in [0, θ * ) and 0 ≤ q * (p) < p for p in (θ * , 1], such that (i) Λ(p, q) > 0 for p in [0, 1] and min{q * (p), p} < q < max{q * (p), p}, (ii) Λ(p, q) < 0 for p in [0, 1], max{q * (p), p} < q ≤ 1 and 0 ≤ q < min{q * (p), p}.

Proof. Since Λ(p, p) = 0 for all p in [0, 1] and F changes sign exactly once at some θ * in (0, 1), it follows from Lemma 4.4 and the implicit function theorem that there exist both a smooth curve q = q * (p), which passes through (θ * , θ * ), such that q * (p) > p for p less than and close to θ * , and q * (p) < p for p greater than and close to θ * , and a positive real number δ, such that Λ(p, q) = 0 in the stripe |q -p| < δ if and only if either q = p or q = q * (p). Using Lemma 4.5, we can extend the curve q = q * (p) to the left until it reaches either p = 0 or q = 1. Note that we can choose δ so small that this extension can never re-enter the strip 0 < q -p < δ once it leaves it. If it first reaches q = 1 at some p = p in (0, θ * ), one can proceed as in the proof of Theorem 5.2 to define q * (p) for p in [0, θ * ) and show that q = q * (p) is a continuous curve. If it never attains q = 1, note that it can only intersect the line q = p at q = p = θ * , thus it can be defined at p = 0 in such a way that p < q * (p) < 1 for p in (0, θ * ). From Lemma 4.5 and the implicit function theorem, the curve q = q * (p) is smooth in this scenario. Similarly, one can extend q * (p) to [θ * , 1] as a continuous curve.

These arguments also show that the nodal set of Λ(p, q) is contained in the line q = p and the curve q = q * (p), so that the conclusions in (i) and (ii) hold.

Similarly to Theorem 5.6, the local stability of (0, N q ) can be determined as follows.

Theorem 5.7. Suppose that the equation F = 0 has a unique root, denoted by θ * , in (0, 1). Then, there exists some positive, continuous function p = p * (q), defined in [0, 1], satisfying q < p * (q) ≤ 1 for q ∈ [0, θ * ) and 0 ≤ p * (q) < q for q ∈ (θ * , 1], such that (i) Λ(q, p) > 0 for q ∈ [0, 1] and min{p * (q), q} < p < max{p * (q), q}, (ii) Λ(q, p) < 0 for q ∈ [0, 1], max{p * (q), q} < p ≤ 1 and 0 ≤ p < min{p * (q), q}.

The proof being the same as that for Theorem 5.6, we omit it.

5.4. Proofs of Theorem 2.2 and Corollary 2.3. We are now in a position to prove some main results of the paper. Theorem 2.2 follows from Theorems 5.2, 5.4, and 5.6.

Proof of Corollary 2.3. If the function F is positive in [0, 1], statement (i) follows from Theorem 5.2. If F > 0 in [0, 1) with F (1) = 0, we can apply Lemma 5.1 to conclude that p = 1 is an evolutionarily singular strategy but not an evolutionarily stable one, and p = 0 is thus the only ESS. Statement (ii) can be proved similarly. Finally, statement (iii) is a direct consequence of Theorem 5.6.

Remark 2. While the function F plays a critical role in the analysis provided in the current section, it appears that it only captures some partial information on Λ(p, q) and cannot possibly determine entirely the nodal set of Λ(p, q). For instance, even the sign of Λ(0, 1) cannot be resolved using only the function F , as it depends on a, D and D in delicate manners. As an example, assume that a(x, y) = λA(x) + (1 -λ)A(y), where λ belongs to [0, 1) and A attains a strict global maximum. Then, for large D, choosing D sufficiently small, Λ(0, 1) < 0. However, for such choices of a, D and D, the function F changes from negative to sign-changing and to positive as λ varies from 0 to 1. We refer to the Appendix 8.3 for further discussions.

6. Full dynamics of the two-species model 6.1. Local stability of semi-trivial steady states. In this subsection, we investigate further the local stability of both semi-trivial steady states (N p , 0) and (0, N q ), for general p and q in [0, 1].

The next result shows that the states (N p , 0) and (0, N q ) cannot be simultaneously stable, i.e. bistability cannot occur. Lemma 6.1. The following assertions hold for any p and q in [0, 1].

(i) If Λ(p, q) > 0, then Λ(q, p) < 0.

(ii) If Λ(p, q) = 0, then either p = q or Λ(q, p) < 0.

Proof. It follows from the variational characterization of Λ(p, q) that

Λ(p, q) = inf ϕ =0, ϕ∈H 1 (Ω) Ω [D(q)(ϕ x ) 2 + D(1 -q)(ϕ y ) 2 -(a -N p )ϕ 2 ] Ω ϕ 2 ≤ Ω [D(q)((N q ) x ) 2 + D(1 -q)((N q ) y ) 2 -(a -N p )(N q ) 2 ] Ω (N p ) 2 = Ω (N p -N q )(N q ) 2 Ω (N p ) 2
, where the last equality follows from the equation of N q . If Λ(p, q) > 0, then

Ω (N q ) 3 < Ω N p (N q ) 2 ≤ Ω (N p ) 3 1/3 Ω (N q ) 3 2/3
, which implies that

Ω (N q ) 3 < Ω (N p ) 3 .
If we assume that Λ(q, p) ≥ 0, we have, by the same argument as above,

Ω (N p ) 3 ≤ Ω (N p ) 2 N q ≤ Ω (N p ) 3 2/3 Ω (N q ) 3 1/3
, from which we get

Ω (N p ) 3 ≤
Clearly, Γ 1 and Γ 2 correspond to the situation when (N p , 0) and (0, N q ) are neutrally stable. Therefore, we have, by Theorems 5.2 to 5.7, Γ 1 = {(p, q) ∈ [0, 1] 2 : q = p or q = q * (p)}, Γ 2 = {(p, q) ∈ [0, 1] 2 : p = q or p = p * (q)}.

Next, we define the sets

Σ 1 : = {(p, q) ∈ [0, 1] 2 : Λ(p, q) > 0 > Λ(q, p)}, Σ 2 : = {(p, q) ∈ [0, 1] 2 : Λ(p, q) < 0 < Λ(q, p)}, Σ 3 : = {(p, q) ∈ [0, 1] 2 : Λ(p, q) < 0, Λ(q, p) < 0}.
The sets Σ i (i = 1, 2, 3) are disjoint and

Σ 1 ∪ Σ 2 ∪ Σ 3 = [0, 1] × [0, 1]/(Γ 1 ∪ Γ 2 ∪ Γ 3 ).
Theorem 6.2. The following characterizations hold:

Σ 1 = {(p, q) ∈ [0, 1] × [0, 1] : (q -q * (p))(q -p) < 0}, Σ 2 = {(p, q) ∈ [0, 1] × [0, 1] : (p -p * (q))(p -q) < 0}, Σ 3 = {(p, q) ∈ [0, 1] × [0, 1] : (q -q * (p))(p -p * (q)) < 0}.
Proof. By Lemma 6.1, one has

Σ 1 = {(p, q) ∈ [0, 1] × [0, 1] : Λ(p, q) > 0}.
It then follows from Theorems 5.2, 5.3 and 5.4 that Σ 1 is determined by (q -q * (p))(qp) < 0. The proof for the characterization of Σ 2 is similar and thus skipped. From the new characterizations for Σ 1 and Σ 2 , it follows that (p, q) ∈ Σ 3 if and only if (q-q * )(q-p) > 0 and (p-p * )(p-q) > 0, which amounts to (q-q * (p))(p-p * (q)) < 0.

6.2. Stability of positive steady states of system (2.5). The following result shows that any positive steady state of system (2.5) is asymptotically stable. It is essentially due to He and Ni [START_REF] He | Global dynamics of the Lotka-Volterra competition-diffusion system: diffusion and spatial heterogeneity I[END_REF]. For the sake of completeness, we have included here a slightly different demonstration of this result.

Lemma 6.3. Suppose that the free growth rate function a is non-constant and that p is not equal to q. Then, any positive steady state of system (2.5) is linearly stable and thus locally asymptotically stable.

Proof. Let (U, V ) denote any positive steady state of system (2.5), i.e.

     D(p)U xx + D(1 -p)U yy + (a -U -V )U = 0 in Ω, D(q)V xx + D(1 -q)V yy + (a -U -V )V = 0 in Ω, (D(p)U x , D(1 -p)U y ) • ν = (D(q)V x , D(1 -q)V y ) • ν = 0 on ∂Ω.
The linear stability of this state is determined by the sign of the principal eigenvalue λ 1 of the linear problem

         D(p)ϕ xx + D(1 -p)ϕ yy + (a -2U -V )ϕ -ψV + λ 1 ϕ = 0 in Ω, D(q)ψ xx + D(1 -q)ψ yy -Uϕ + (a -U -2V )ψ + λ 1 ψ = 0 in Ω, (D(p)ϕ x , D(1 -p)ϕ y ) • ν = 0 on ∂Ω, (D(q)ψ x , D(1 -q)ψ y ) • ν = 0 on ∂Ω.
It is known (see [START_REF] Cantrell | Spatial ecology via reaction-diffusion equations[END_REF][START_REF] Smith | Monotone dynamical systems: an introduction to the theory of competitive and cooperative systems[END_REF]) that we may choose ϕ > 0 and ψ < 0 in Ω. Set W = ϕ/U and Z = -ψ/V so that W and Z are both positive in Ω and satisfy

     D(p)(U 2 W x ) x + D(1 -p)(U 2 W y ) y -U 3 W + U 2 V Z + λ 1 U 2 W = 0 in Ω, D(q)(V 2 Z x ) x + D(1 -q)(V 2 Z y ) y + UV 2 W -V 3 Z + λ 1 V 2 Z = 0 in Ω, (D(p)W x , D(1 -p)W y ) • ν = (D(q)Z x , D(1 -q)Z y ) • ν = 0 on ∂Ω.
Multiplying the first line of this system by W 2 and integrating the result over Ω yields

Ω D(p)U 2 W (W x ) 2 + D(1 -p)U 2 W (W y ) 2 + (UW ) 3 -(UW ) 2 (V Z) -λ 1 U 2 W 3 = 0.
Similarly, multiplying the second line of the system by Z 2 and integrating the result over Ω, we find that

Ω D(q)V 2 Z x (Z 2 ) x + D(1 -q)V 2 Z y (Z 2 ) y + (V Z) 3 -(UW )(V Z) 2 -λ 1 V 2 Z 3 = 0.
It suffices to show λ 1 > 0. We argue by contradiction by assuming that λ 1 ≤ 0. Then, one has (6.1)

Ω (UW ) 3 ≤ Ω (UW ) 2 (V Z),
and the equality in (6.1) holds if and only if λ 1 = 0, W is a positive constant, and UW is a positive scalar multiple of V Z. Similarly, if λ 1 ≤ 0, one has (6.2)

Ω (V Z) 3 ≤ Ω (UW )(V Z) 2 ,
and the equality in (6.2) holds if and only if λ 1 = 0, Z is a positive constant, and UW is a positive scalar multiple of V Z. Finally, it follows from the Hölder inequality that (6. 3)

Ω (UW ) 3 ≤ Ω (UW ) 3 2/3 Ω (V Z) 3 1/3 and (6.4) Ω (V Z) 3 ≤ Ω (V Z) 3 2/3 Ω (UW ) 3 1/3
, from which we see that inequalities (6.1), (6.2), (6.3) and (6.4) must all be equalities. As a consequence, λ 1 is zero, both W and Z are positive constants, and UW is a positive scalar multiple of V Z, i.e. U = cV for some positive constant c. Therefore, U satisfies

     D(p)U xx + D(1 -p)U yy + (a -(c + 1)U)U = 0 in Ω, D(q)U xx + D(1 -q)U yy + (a -(c + 1)U)U = 0 in Ω, (D(p)U x , D(1 -p)U y ) • ν = (D(q)U x , D(1 -q)U y ) • ν = 0 on ∂Ω.
Hence, we find that U ≡ N p /(c + 1) and U ≡ N q /(c + 1), which implies that N p ≡ N q . As the function a is non-constant, so is N p . Since p = q, by subtracting the equations of the systems respectively satisfied by N p and N q , we see that N p is a solution to system (8.1) and is therefore a non-constant function, which contradicts the assumption. 6.3. Global dynamics of system (2.5). As the two-species competition model (2.5) is strongly monotone, its global dynamics can be fully determined by the local stability of its equilibria in some cases, see [18, Chapter IV] for instance. Let us recall below some known facts. (a) If there is no positive steady state, then one of the semi-trivial equilibria is unstable and the other is globally asymptotically stable among non-negative and nonidentically zero initial data. (b) If there is a unique positive steady state and it is stable, then it is globally asymptotically stable. (c) If all positive steady states are asymptotically stable, then there is at most one of them. In particular, either (a) or (b) applies.

We are now ready to infer on the global stability of steady states.

Theorem 6.4. Suppose that the function F is positive in [0, 1] and let p * and q * be the functions introduced in Theorems 5.2 and 5.3, respectively. Then, one of the following statements holds.

(i) If p < q < q * , then the steady state (N p , 0) is globally asymptotically stable;

(ii) If q < p < p * , then the steady state (0, N q ) is globally asymptotically stable;

(iii) If either q * < q ≤ 1 or p * < p ≤ 1 holds, then system (2.5) has a unique positive steady state, which is also globally asymptotically stable among non-negative and not identically zero initial data.

Proof. We first establish statement (i). By Theorem 6.2, we see that under assumption p ∈ [0, 1] and p < q < q * , Λ(p, q) > 0 > Λ(q, p). Hence, (N p , 0) is stable and (0, N q ) is unstable. As system (2.5) is strongly monotone, by Lemma 6.3 and statements (a) and (c), (N p , 0) is globally stable. The proof of statement (ii) is similar to that of statement (i) and thus omitted.

For statement (iii), Λ(p, q) < 0 and Λ(q, p) < 0. Hence, both states (N p , 0) and (0, N q ) are unstable. As system (2.5) is strongly monotone, by Lemma 6.3 and statements (b) and (c) recalled above, there is a unique positive steady state which is globally asymptotically stable.

Note that p * ≡ 1 if and only if q * ≡ 1. For such scenario, alternative (iii) in Theorem 6.4 does not occur, the state (N p , 0) is globally stable when q > p and the state (0, N q ) is globally stable when q < p.

Similar to Theorem 6.4, if F < 0 in [0, 1], the global dynamics of system (2.5) can be characterized as follows. Theorem 6.5. Suppose that the function F is negative in [0, 1] and let p * and q * be the functions introduced in Theorems 5.4 and 5.5, respectively. Then, one of the following statements holds.

(i) If q * < q < p, then the steady state (N p , 0) is globally asymptotically stable.

(ii) If p * < p < q, then the steady state (0, N q ) is globally asymptotically stable.

(iii) If either 0 ≤ q < q * (p) or 0 ≤ p < p * (q) holds, then system (2.5) has a unique positive steady state, which is also globally asymptotically stable among non-negative and not identically zero initial data.

If the function F changes sign in (0, 1), the global dynamics of system (2.5) is given by the following result: Theorem 6.6. Suppose that the function F changes sign in (0, 1) and let p = p * (q) and q = q * (p) be the functions given in Theorems 5.6 and 5.7. Then, one of the following alternatives holds.

(i) If (q -q * )(p -p * ) < 0, then system (2.5) has a unique positive steady state, which is also globally asymptotically stable among non-negative and not identically zero initial data. (ii) If (q -q * )(q -p) < 0, then the steady state (N p , 0) is globally asymptotically stable. (iii) If (q -q * )(q -p) < 0, then the steady state (0, N q ) is globally asymptotically stable.

The proof of Theorem 6.6 is the same as that of Theorem 6.4 and follows from Theorem 6.2 and Lemma 6.3.

Discussions

In this paper, we considered a reaction-diffusion model for two competing populations, which disperse in a bounded two dimensional habitat by moving horizontally and vertically with different probabilities but are otherwise identical. We regard these probabilities as dispersal strategies and ask what strategies are evolutionarily stable.

Our main finding is that the only evolutionarily stable dispersal strategies are to move in one direction. In particular, when the resources are distributed inhomogeneously only in one direction, e.g. horizontally, our result implies that the evolutionarily stable strategy could simply be to move in the vertical direction, in which the resources are homogeneously distributed. More precisely, we introduced a function F of the dispersal probability, which measures the difference between the spatial variations of the population equilibrium distributions in horizontal and vertical directions: when it is positive, the species has more variations in the horizontal direction; when it is negative, it has more variations in the vertical direction. We show that function F is monotone decreasing and that the evolutionarily stable dispersal strategies are to maximize the function F when it is positive and to minimize it when it is negative. As the population distribution at equilibrium is often positively correlated with the resource distribution, function F also indirectly measures the difference between the resource variations in horizontal and vertical directions. Therefore, our results seem to predict that it is more favorable for the species to choose the direction with smaller variations in resource distributions. This finding seems to be in agreement with the classical results of Hasting [START_REF] Hastings | Can spatial variation alone lead to selection for dispersal?[END_REF] and Dockery et al. [START_REF] Dockery | The evolution of slow dispersal rates: a reaction diffusion model[END_REF] for the evolution of slow dispersal, i.e., random diffusion is selected against in spatially heterogeneous environments.

We further investigated the local and global dynamics of the two-species system and determined the dynamics of system (2.5) for three different cases of the selection gradient. We applied numerical simulations to illustrate how the shapes of function F , the local stability of the semi-trivial steady states and the global dynamics of the system sequentially change with respect to a certain parameter which measures the difference between the resource variations in the horizontal and vertical directions. Our numerical results suggest that if the spatial variations of resource distributions in vertical and horizontal directions become more comparable, the chances for the coexistence of two competing populations could be greater.

One of our future works is to extend the mathematical modelling and analysis to any dimensional habitats, and to continuous trait models. Another future work will be to include a temporal variation of the environment and ask how it affects the evolution of horizontal and vertical movement. For example, if we choose a(t, x, y) = λA(x) + (1λ)B(t, y), a natural question is when vertical movement will be selected as in [START_REF] Hutson | The evolution of dispersal rates in a heterogeneous time-periodic environment[END_REF].

8.2.

A remark about a possible degeneracy induced by the domain Ω. Throughout the paper, we have assumed that Ω is a strictly convex domain. We now comment on this point, showing with a very basic example that a domain with flat parts on its boundary may lead to a degeneracy of the function F .

Consider Ω = (0, 1) × (0, 1). Let f and D be given in Example 1. Set

a := -(D + D) W xx W + W.
It is easy to check that for each θ ∈ [0, 1], W also solves (2.3), i.e. N θ ≡ W for each θ ∈ [0, 1]. Since W is non-constant, a is also non-constant. Furthermore, for each p, q ∈ [0, 1], (2.5) has a continuum of positive steady states of the form (U, V ) = (sW, (1 -s)W ), s ∈ (0, 1). Moreover, in this case, F ≡ 0 in [0, 1] even though the function a is not constant. Indeed, since f is even, we have W (x, y) = W (y, x), which implies that F ≡ 0 in [0, 1].

8.3.

The sign of Λ(0, 1). In this subsection we construct an example to support the claim made in Remark 2, as shown by the following result:

Proposition 1. Assume that a(x, y) = λA(x) + (1 -λ)A(y), where λ belongs to [0, 1), and A is positive Hölder continuous function, which attains a strict global maximum at y = ŷ for some ŷ, i.e. A(y) < A(y) for every y = ŷ. Then, there exists some positive real number δ such that for D > 1/δ, then for sufficiently small D, Λ(0, 1) < 0.

We first establish some a priori estimate on N Set Ω := ∪ y<y<y Ω y for some y < y. In particular, there exists some δ > 0 such that for D > 1/δ, then for sufficiently small D, Λ(0, 1) < 0, provided that To see this, observe that for all y = ŷ and y ∈ (y, y),

1 |Ω y | Ωy Ωx A(z) dz |Ω x | dx < 1 |Ω y | Ωy Ωx A(ŷ) dz |Ω x | dx = A(ŷ).
The conclusion thus follows from Lemma 8.4.

5 Figure 1 .

 51 Figure 1. Numerical approximations of the graphs of the function F for λ taking the values 0, 0.1, 0.2, 0.3, 0.4, and 0.5, in the case where Ω is a disk of radius 2 centered at the origin, D = 0.1, D = 10, and A(x) = 2 -sin(πx).

  The function a.

  The steady state for θ = 0.

  The steady state for θ = 0.4.

  The steady state for θ = 0.8.

  The steady state for θ = 1.

Figure 2 .

 2 Figure 2. Isolines for numerical approximations of the function a and of the steady state N θ for θ taking the values 0, 0.2, 0.4, 0.6, 0.8, and 1, in the case where Ω is a disk of radius 2 centered at the origin, D = 0.1, D = 10, λ = 0.4, and A(x) = 2 -sin(πx).

  The function a.

  The steady state for θ = 0.

  The steady state for θ = 0.2.

  The steady state for θ = 0.4.

Figure 3 .

 3 Figure 3. Isolines for numerical approximations of the function a and of the steady state N θ for θ taking the values 0, 0.2, 0.4, 0.6, 0.8, and 1, in the case where Ω is a disk of radius 2 centered at the origin, D = 0.1, D = 10, λ = 0.6, and A(x) = 2 -sin(πx).

Figure 4 .

 4 Figure 4. Numerical approximations of the graphs of F (θ) as a function of λ for θ taking the values 0, 0.1, 0.2, 0.3, 0.4, and 0.5, in the case where Ω is a disk of radius 2 centered at the origin, D = 0.1, D = 10, and A(x) = 4 -1 4 x 2 .

Figure 5

 5 Figure 5. Numerical approximations of the nodal sets of Λ(p, q) for λ taking the values 0, 0.1, 0.2, 0.3, 0.4, and 0.5, in the case where Ω is a disk of radius 2 centered at the origin, D = 0.1, D = 10, and A(x) = 2 -sin(πx). The subset colored in green is the one in which Λ(p, q) > 0, that is in which the steady state (N p , 0) is linearly stable.

Figure 6

 6 Figure 6. Numerical approximations of the nodal sets of Λ(p, q) and Λ(q, p) for λ taking the values 0, 0.1, 0.2, 0.3, 0.4, and 0.5, in the case where Ω is a disk of radius 2 centered at the origin, D = 0.1, D = 10, and A(x) = 4 -1 4 x 2 . The subset colored in green is the one in which Λ(p, q) > 0 > Λ(q, p), i.e. in which (N p , 0) is stable, while that colored in red is the one in which Λ(q, p) > 0 > Λ(q, p), i.e. in which (0, N q ) is stable. The white region is where both (N p , 0) and (0, N q ) are unstable.

Lemma 8 . 3 .D

 83 0 . By definition, N 0 satisfies D(N 0 ) xx + D(N 0 ) yy + (a -N 0 )N 0 = 0 in Ω, (D(N 0 ) x , D(N 0 ) y ) • ν = 0 on ∂Ω.Define Ω x := {y : (x, y) ∈ Ω} and Ω y := {x : (x, y) ∈ Ω}, and let Ω := ∪ y * <y<y * Ω y .If Ω is strictly convex, we may express it asΩ = {(x, y) : x * < x < x * , y * (x) < y < y * (x)}for some x * < x * , y * (x) < y * (x) in (x * , x * ) and y * (x) = y * (x) at x = x * , x * . For any x ∈ (x * , x * ), clearly Ω x := {(x, y) : y * (x) < y < y * (x)}. Suppose that the function a satisfies (A1), and Ω is strictly convex, C 1 . For any ǫ > 0, there exists δ > 0 such that if D > 1/δ, then for sufficiently small D,1 |Ω x | Ωx a(x, y) dy -ǫ ≤ N 0 (x, y) ≤ 1 |Ω x | Ωx a(x, y) dy + ǫ holds for any (x, y) ∈ Ω. Proof. Given ǫ > 0 small, choose function ã ∈ C 2 (Ω) such that a + ǫ < ã ≤ a + 2ǫ in Ω and ã(x, y) is constant for x ∈ [x * , x * + ǫ] ∪ [x * -ǫ,x * ] and y * (x) < y < y * (x). For each x ∈ (x * , x * ), let Ñ(x, y) denote the unique positive solution of the equation  Ñyy + Ñ (ã -Ñ) = 0 in y * (x) < y < y * (x), Ñy (y * (x)) = -ǫ D , Ñy (y * (x)) = ǫ D .Then Ñ is independent of x for x ∈ [x * , x * +ǫ]∪[x * -ǫ, x * ], for which it satisfies Ñx = 0. Therefore, Ñ satisfiesD Ñx ν x + D Ñy ν y = ǫ|ν y | ≥ 0 for (x, y) ∈ ∂Ω and x ∈ [x * , x * + ǫ] ∪ [x * -ǫ, x * ]. For (x, y) ∈ ∂Ω and x ∈ [x * + ǫ, x * -ǫ], D Ñx ν x + D Ñy ν y ≥ -D Ñx ν x L ∞ + ǫ min x∈[x * +ǫ,x * -ǫ]|ν y | > 0 for sufficiently small D, since |ν y | is a strictly positive and continuous function for (x, y) ∈ ∂Ω and x ∈ [x * + ǫ, x * -ǫ]. This implies that, for sufficiently small D, D Ñx ν x + D Ñy ν y ≥ 0 holds on ∂Ω.In Ω, Ñ satisfiesD Ñxx + D Ñyy + Ñ(a -Ñ ) = D Ñxx + Ñ(a -ã) ≤ D Ñxx -Ñ ǫ < 0,provided that D is sufficiently small. This implies that Ñ is a super-solution for the equation of N 0 . Hence, if D is small, then N 0 ≤ Ñ in Ω. As D tends to infinity, Ñ converges uniformly to1 |Ωx| Ωx ã(x, y) dy in Ω. Hence, there exists some δ > 0 such that for D > 1/δ, if D is small, then Ñ ≤ 1 |Ω x | Ωx ã(x, y) dy + ǫ, which implies that N 0 ≤ Ñ ≤ 1 |Ω x | Ωx a(x, y) dy + 3ǫ holds in Ω. The lower bound of N 0 can be similarly established.

Lemma 8 . 4 .≤ min y≤y≤y 1

 841 The following result holds:|Ω y | Ωy 1 |Ω x | Ωx a(x, z) dz -a(x, y) dx.

Ωy 1

 1 |Ω x | Ωx a(x, z) dz -a(x, y) dx ≡ 0 for y ∈ [y, y].Proof. By Lemma 8.3, for any ǫ > 0, there exists some δ > 0 such that for D > 1/δ, then for sufficiently small D, N 0 (x, y) ≤ 1 |Ωx| Ωx a(x, z) dz + ǫ in Ω. Note thatΛ(0, 1) = inf ϕ =0, ϕ∈H 1 (Ω) Ω [D(ϕ x ) 2 + D(ϕ y ) 2 -(a -N 0 )ϕ 2 ] Ω ϕ 2.

Ω(N q ) 3 , which is a contradiction. Hence, statement (i) holds.If Λ(p, q) = 0, following the same argument as above, we see that Λ(q, p) ≤ 0. If Λ(q, p) = 0, the only possibility is that N p ≡ N q , that is p = q, which proves (ii).To describe the global dynamics of the two-species model, we first introduce the setsΓ 1 : = {(p, q) ∈ [0, 1] 2 : Λ(p, q) = 0},Γ 2 : = {(p, q) ∈ [0, 1] 2 : Λ(q, p) = 0}.

Acknowledgement. We thank Professor Yoshikazu Giga for the helpful discussions which motivated the study of anisotropic diffusion. We are also very grateful to Maxime Chupin for his help with the post-processing of our numerical simulations.

By choosing ϕ = ϕ(y) we have

holds for any ǫ > 0 and ϕ(y) = 0, where

Hence, given any x and taking the infimum of the right-hand side of the above inequality over all ϕ = ϕ(y) = 0, we have


Shyama Prasad Das

Emil J Hopfinger
email: emil.hopfinger@hmg.inpg.fr

Estelle Guyez

Parametrically

Parametrically forced, breaking gravity waves in a circular cylinder

Keywords: gravity waves, parametric instability, geyser formation

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Parametrically forced surface waves, known as Faraday waves, that are sub-harmonically excited, have been studied rather extensively [START_REF] Benjamin | The stability of plane free surface of a liquid in a vertical periodic motion[END_REF]Kumar & Tuckermann, 1994[START_REF] Edwards | Patterns and quasi-patterns in the Faraday experiment[END_REF]. However, far less is known about large amplitude waves and wave breaking conditions especially of gravity waves. Large amplitude sloshing in containers is encountered in rocket engine fuel tanks and ship tanks. In these applications it is important to be able to predict the forces exerted by the sloshing motion on the tank walls as well as the consequences of sloshing and wave breaking on interfacial heat and mass transfer. Experimental results on sloshing due to horizontal forcing of circular and square-base cylindrical containers have been reported by [START_REF] Royon-Lebeaud | Liquid sloshing and wave breaking in circular and square-base cylindrical containers[END_REF]. Here we present results on sloshing when the container is subjected to vertical (along the axis) motions. Vertical forcing allows exciting axi-symmetric as well as asymmetric wave modes depending on forcing frequency.

Parametrically forced gravity waves and the breaking scenarios in two-dimensions have been studied by [START_REF] Jiang | Period tripling and energy dissipation of breaking standing waves[END_REF]. Breaking occurs when the wave amplitude has reached a value such that the downward acceleration of the wave crest is equal or slightly larger than the gravity. In a circular cylinder, axi-symmetric wave breaking can be much more violent and geyser formation is possible [START_REF] Zeff | Singularity dynamics in curvature collapse and jet eruption on a fluid surface[END_REF]. [START_REF] Henderson | Single-mode Faraday waves in small cylinders[END_REF], from here on referred to as HM, determined the stability conditions and the steady-state wave amplitude of the axi-symmetric wave mode in a "small" circular cylinder. Small implies that dissipation is important.

In this paper we present results on parametrically forced gravity waves in a circular cylinder of size similar to the one of HM but using a fluid with low kinematic viscosity and low surface tension so that dissipation remains small. We first compare the experimental stability space of forcing amplitude as a function of wave frequency with theoretical predictions taking into account surface tension and viscosity corrections according to HM. This stability space is then compared with wave breaking threshold in the forcing amplitude-wave frequency domain. The wave amplitude response curves, obtained for the same frequency domain, exhibit the existence of steady-state wave motion. In the unstable regime, of the axi-symmetric mode (0,1) a finite time singularity can occur with intense jet formation, a phenomenon demonstrated by [START_REF] Zeff | Singularity dynamics in curvature collapse and jet eruption on a fluid surface[END_REF] in fluids of high viscosity. Here this singularity is shown to exist in fluids of low viscosity and low surface tension; in water it is less pronounced because of the large surface tension causing capillary wave perturbations.

Experimental conditions

The experiments have been conducted in a circular, cylindrical test cell, 5 cm in diameter and 6 cm in height that is mounted on a vertically vibrating table that is used in the frequency range 10 to 14 Hz and vibration amplitudes A= 0 to 1 mm. In the container used, the sub-harmonically excited axi-symmetric gravity wave mode has a wave number k 01 R=3.832. The interesting aspect of parametrically forced waves is the close neighbourhood of the asymmetric wave modes (2,1), (3,1),

(1,2), (4,1); the wave motion will lock to the mode of frequency ω /2 with lowest dissipation. The container is filled to a depth h=3 cm with a low viscosity and low surface tension liquid, FC-72, of ν = 0.00406 cm 2 /s, σ =11.5 dyn/cm and ρ =1.69Kg/m 3 at 20°C. The liquid depth is sufficient to satisfy (asymptotically deepwater conditions). The container acceleration 01

( ) tanh k h ≅ 1 2 ( ) sin a t A t ω ω = used in the experiments is ≤ ) (t a 5 m/s 2
and the forcing frequency (Hz) covered is 10.8 ≤ ω/2π ≤13.6.

Because of the small size of the container measurements were made by visualizations and image analysis only. Generally, a digital camera with an acquisition speed of 60 frames/second (fps) was used. Near resonant conditions with high velocity geyser formation, images were also taken with a high-speed camera at 1000 fps.

Theoretical concepts

The theory of weakly non-linear axi-symmetric wave motion in a circular cylinder (taking into account viscous effects) has been developed by HM following [START_REF] Miles | Resonantly forced surface waves in a circular cylinder[END_REF]. The dispersion relation with surface tension added is of the form

) tanh( ) 1 ( 2 2 h k g k gk mn mn mn mn ρ σ ω + =
.

(1)

The boundary condition on the container wall is that the velocity potential

φ is 0 ) '( / = = ∂ ∂ = R k J r mn m R r
φ which gives for the axi-symmetric wave modes, n =1 and 2, k 01 R=3.8317 and k 02 R= 5.3314 and for the asymmetric modes k 21 R = 3.0542, k 31 R = 4.2012, k 41 R = 5.3176. The corresponding natural frequencies are given by (1). For FC-72 and the container of R=2.5cm filled to a depth, h=3cm, the natural frequency of mode (0,1) is 39.090 rad s -1 .

The natural frequency shift due to linear damping is where

) 1 ( ˆδ ω ω - = mn mn / nm δ κ ω =
is the damping ratio and the damping rate. Viscous dissipation in the Stokes boundary layer dominates over bulk dissipation. HM give the following analytical expression for deep water and a free surface without rigidity:

κ 2 1/ 2 2 1 ( / ) 1 ( ) 2 2 1 ( / ) mn mn mn m k R R m k R ν δ ω + = - (2)
Using the dispersion relation (1) and neglecting surface tension terms, (2) can be written in the form

, ( 3
) 2 3 1/ 1 ( / ) C R g δ ν = 4
where the coefficient , a function of the wave mode as indicated by (2). Its value predicted by ( 2) is about a factor of 2 less than the values determined from experiments (Section 4.1).

1

C

In the theory of HM the steady-state wave amplitude as a function of time (at a given spatial position) depends on a small parameter ε , that is for the axi-symmetric mode

, ( 4
γ δ ε ⎡ = - ⎣ ⎤ ⎦
Neutral stability is given by β γ = ± , with corresponding threshold forcing amplitude:

( )

1 2 2 2 2 01 2 2 2 01 ( 2) 1 3.832 tanh(3.832 / ) (2 ) A R hR ω ω δ ω ⎡ ⎤ - ⎢ ⎥ = + ⎢ ⎥ ⎣ ⎦ (7)
4 Results

Damping coefficient

For the determination of the instability threshold it is important to know the damping coefficient δ with good accuracy. When the wave amplitude was established, the forcing was stopped and the wave motion decayed freely in time. Since the damping rate is small, the decay is exponential and is of the form

0 exp( ) b b t δω = -
, where b 0 is the wave amplitude at the stop of the forcing. Taking images of the wave motion at 60fps, corresponds to a sampling frequency of about ten times the wave frequency. By interpolation the wave motion was reconstructed giving the successive wave amplitudes to an accuracy of about 3%. For the wave mode k 01 , the damping coefficient is δ ≅ 0.01 and corresponds to C 1 ≅1.75 in (3). This value was used in calculating the instability threshold from (7). Using the relation , the natural frequency of the axi-symmetric mode after viscosity

correction is rad s ) 1 ( ˆ01 01 δ ω ω - = 01 ˆ38.699 ω = -1
. From here on we drop the hat ( ).

01 01 ˆω ω ≡

Instability and wave breaking thresholds

Figure 1a shows the non-dimensional forcing amplitudes of the instability and wave breaking thresholds as a function of non-dimensional frequency . The solid line corresponds to (7), with δ = 0.01. There is good agreement between the experimental and calculated instability thresholds in the range 1.02< <0.96. For >1.02 the experimental values deviate substantially from the theoretical ones because the excited wave mode is k where the two overlap and breaking always occurs just above the instability threshold; in this range no stable waves exist. Figure 1a shows clearly that the limits of existence of stable wave motion depend on both, the forcing frequency and forcing amplitude related by the frequency-offset parameter

{ } ( ) 2 1 2 1 2 01 1 2 ) ( 1 2 + - - = ε ε δ ω ω (8)

Wave amplitude response

Jet formation

For container forcing amplitudes A > A d wave breaking occurs. Axi-symmetric wave mode breaking occurs via a high velocity geyser or jet formation. The geyser is formed by the collapse of the depression of the large amplitude axi-symmetric standing wave, a phenomenon referred to as finitetime singularity [START_REF] Zeff | Singularity dynamics in curvature collapse and jet eruption on a fluid surface[END_REF]. Images of the wave depression (lower part of image) and of geyser formation (upper part of image) are shown in Figures 2a and2b for two different forcing conditions; these conditions are indicated in figure 1a by the symbols and ◊. . This is shown by the solid lines in Figure 3.

Conclusions

The experimentally determined stability threshold of the parametrically forced axi-symmetric gravity wave mode in low viscosity and low surface tension fluid in a circular cylinder is shown to be in good agreement with the HM theory. Bifurcations to other wave modes take place when the forcing frequency approaches the natural frequency of these modes. Stable wave motions exist between the stability threshold and the wave-breaking threshold established in the present experiments. The amplitude response curve shows that in a certain frequency range these stable waves exhibit amplitude modulations. In the unstable regime, a finite time singularity occurs with intense geyser or jet formation, a phenomenon demonstrated by [START_REF] Zeff | Singularity dynamics in curvature collapse and jet eruption on a fluid surface[END_REF] in fluids of high viscosity. Here this singularity is demonstrated for a low viscosity and low surface tension liquid. A dimensionless expression for the jet velocity is proposed that takes into account the container size in addition to the fluid properties.

  general, there is a large difference between the instability threshold forcing amplitude and wave breaking threshold, except when

.

  The frequency domain of existence of wave breaking increases with forcing amplitude. The onset of wave breaking is at ε = 0.0322 giving β = -γ = 0.951. The corresponding dimensionless forcing frequency obtained from (8) is 0.969:

Figure

  Figure 1b shows the wave amplitude response b/R as a function of dimensionless frequency ω/2ω 01 for forcing amplitude A/R=0.006. All experimental points correspond to stable wave motions. The closed symbols indicate the existence of the axi-symmetric mode (0,1) and the open symbols refer to mode (3,1) when ω/2ω 01 >1.02 and to mode (2,1) when ω/2ω 01 <0.97. Stable axi-symmetric wave motions are indeed observed for γ β γ -< < but near γ there is an amplitude modulation.

  A composite of two photographs showing the axi-symmetric standing wave depression below and geyser formation above for (a), A/R=0.0167, ω/2ω 01 =0.96; (b), A/R=0.0235, ω/2ω 01 =1.017. The horizontal and vertical bars represent 1cm. The horizontal scale changes with location because of optical distortions.

Figure0

  Figure 3. Normalised jet velocity,

  b c is the stable wave amplitude above which bubble pinch-off (Figure2a) takes place. The value of the critical Weber number is not given but their experimental results suggest a value of We c ≈ 4000. The container radius should enter in the scaling because the critical wave amplitude b c must depend on R. Viscous cut-off prevents V g from becoming infinite as b→ b c . A non-smooth surface has a similar cut-off effect as viscosity. This is the case with water where perturbations in the form of parasitic capillary waves appear. Zeff et al. used a glycerin-water solution of ν=1.94 cm 2 s -1 to keep the surface smooth at large wave amplitudes. For the fluid used in the present experiments parasitic capillary waves are not a problem because of the low surface tension. Figure3shows the normalised velocity of the jet as a function of b/R. The ratio of the critical wave amplitude to radius is (b/R) c ≅1.02 and We=4000. When the jet velocity and wave amplitude are appropriately scaled, the Weber number relation is

Acknowledgements

The valuable help of Deepak Singh during preliminary experiments is gratefully acknowledged. The work was financially supported by contract CNS n° 60167 within the COMPERE program.