K Gittins
email: katie.gittins@unine.ch

B Helffer

Courant-sharp Robin eigenvalues for the square and other planar domains

Keywords: MSC classification (2010): 35P99, 58J50, 58J37 Courant-sharp, Robin eigenvalues, square, planar domains

This paper is devoted to the determination of the cases where there is equality in Courant's nodal domain theorem in the case of a Robin boundary condition. For the square, we partially extend the results that were obtained by Pleijel, Bérard-Helffer, Helffer-Persson-Sundqvist for the Dirichlet and Neumann problems.

After proving some general results that hold for any value of the Robin parameter h, we focus on the case when h is large. We hope to come back to the analysis when h is small in a second paper.

We also obtain some semi-stability results for the number of nodal domains of a Robin eigenfunction of a domain with C 2,α boundary (α > 0) as h large varies.

1 Introduction.

Let Ω ⊂ R m , m ≥ 2, be a bounded, connected, open set with Lipschitz boundary and let h ∈ R, h ≥ 0. The case when h < 0 is mathematically interesting but less motivated by Physics. The Robin eigenvalues of the Laplacian on Ω with parameter h are λ k,h (Ω) ∈ R, k ∈ N, k ≥ 1, such that there exists a function u k ∈ H 1 (Ω) which satisfies

-∆u k (x) = λ k,h (Ω)u k (x) , x ∈ Ω , ∂ ∂ν u k (x) + h u k (x) = 0 , x ∈ ∂Ω ,
where ν is the outward-pointing unit normal to ∂Ω.

We recall that by the minimax principle, the Robin problem is associated with the quadratic form:

H 1 (Ω) u → Ω |∇u| 2 + h ∂Ω |u ∂Ω | 2 dσ ,
where u ∂Ω is the trace of u. So the spectrum is monotonically increasing with respect to h for h ∈ [0, +∞). That is, the Robin eigenvalues with h > 0 interpolate between the Neumann eigenvalues (h = 0) and the Dirichlet eigenvalues (h = +∞).

The Robin eigenvalues satisfy the celebrated Courant nodal domain theorem [12] stating that any eigenfunction corresponding to λ k,h (Ω) has at most k nodal domains. We consider the Courant-sharp Robin eigenvalues of Ω. We call a Robin eigenvalue λ k,h (Ω) Courant-sharp if it has a corresponding eigenfunction that has exactly k nodal domains. As for the Dirichlet and Neumann eigenvalues, λ 1,h (Ω) and λ 2,h (Ω) are Courant-sharp for all h ≥ 0.

An interesting question is whether it is possible to follow the Courant-sharp (Neumann) eigenvalues with h = 0 to Courant-sharp (Dirichlet) eigenvalues as h → +∞, or whether there are some critical values h * (k, Ω) after which the Robin eigenvalues λ k,h (Ω), h ≥ h * (k, Ω) become Courant-sharp or are no longer Courant-sharp. We note that throughout this paper, we denote the Dirichlet eigenvalues by λ D k and the Neumann eigenvalues by λ N k .

We consider the particular example where Ω is a square S in R 2 of sidelength π and the main question is: Is it possible to determine the Courant-sharp Robin eigenvalues of this square?

As λ 2,h (S) = λ 3,h (S) by a symmetry argument, it follows immediately that λ 3,h (S) is not Courant-sharp for any h ≥ 0. In addition, λ 4,h (S) is Courantsharp for all h ≥ 0, see Subsection 2.2.

It was asserted by Pleijel in [START_REF] Pleijel | Remarks on Courant's nodal line theorem[END_REF] that the only Courant-sharp Dirichlet eigenvalues of the square are for k = 1, 2, 4. This was shown rigorously in [START_REF] Bérard | Dirichlet eigenfunctions of the square membrane: Courant's property, and A. Stern's and A. Pleijel's analyses[END_REF]. The only Courant-sharp Neumann eigenvalues of the square are for k = 1, 2, 4, 5, 9, as shown in [START_REF] Helffer | Nodal domains in the square-the Neumann case[END_REF].

The first step to obtain the results of [START_REF] Bérard | Dirichlet eigenfunctions of the square membrane: Courant's property, and A. Stern's and A. Pleijel's analyses[END_REF][START_REF] Helffer | Nodal domains in the square-the Neumann case[END_REF] is to reduce the number of potential Courant-sharp eigenvalues by invoking an argument which was inspired by the founding paper of Pleijel [START_REF] Pleijel | Remarks on Courant's nodal line theorem[END_REF]. We employ a similar argument in Section 3 to reduce the possible cases that may give rise to Courant-sharp Robin eigenvalues. We have the following theorem.

Theorem 1.1. Let h ≥ 0. If λ k,h (S) is an eigenvalue of S with k ≥ 520, then it is not Courant-sharp.

We note that in the case of a Dirichlet boundary condition, the equivalent statement in [START_REF] Pleijel | Remarks on Courant's nodal line theorem[END_REF] gives k ≥ 34 and in the case of a Neumann boundary condition, [START_REF] Helffer | Nodal domains in the square-the Neumann case[END_REF], k ≥ 209. The strategies of [START_REF] Bérard | Dirichlet eigenfunctions of the square membrane: Courant's property, and A. Stern's and A. Pleijel's analyses[END_REF][START_REF] Helffer | Nodal domains in the square-the Neumann case[END_REF] are then either to re-implement the Faber-Krahn inequality, or to use symmetry properties of the corresponding eigenfunctions to further eliminate potential Courant-sharp eigenvalues. One is then reduced to the analysis of the nodal structure of very few families of eigenfunctions that belong to two-dimensional eigenspaces.

We will show that the Robin eigenfunctions satisfy analogous symmetry properties. We were not able to eliminate potential Courant-sharp cases via symmetry as it is possible that a Robin eigenvalue has multiplicity larger than 2 and the corresponding eigenfunctions have no common symmetries (see Subsection 7.2).

In addition, for a Robin eigenvalue λ k,h (S), we do not know how to take the relationship between k and h into account in an efficient way. Indeed, to prove Theorem 1.1 our arguments are independent of h as they rely on the monotonicity of the Robin eigenvalues and comparison to the corresponding Dirichlet and Neumann eigenvalues.

We note that the recent articles [START_REF] Freitas | Extremal domains and Polya type inequalities for the Robin Laplacian and union of rectangles[END_REF][START_REF] Gittins | Upper bounds for Courant-sharp Neumann and Robin eigenvalues[END_REF] also consider the Robin eigenvalues of Euclidean domains and make use of this monotonicity property. In [START_REF] Gittins | Upper bounds for Courant-sharp Neumann and Robin eigenvalues[END_REF], upper bounds are obtained for the Courant-sharp Neumann and Robin eigenvalues with h > 0 of a bounded, connected, open set Ω ⊂ R n with C 2 boundary. In [START_REF] Freitas | Extremal domains and Polya type inequalities for the Robin Laplacian and union of rectangles[END_REF], it is shown that the Robin eigenvalues with h > 0 on rectangles and unions of rectangles with prescribed area satisfy Pólya-type inequalities.

In addition, we treat the problem asymptotically as h → +∞. Hence we show that for h large enough the only Courant-sharp Robin eigenvalues are for k = 1, 2, 4.

Theorem 1.2. There exists h 1 > 0 such that for h ≥ h 1 , the Courant-sharp cases for the Robin problem are the same as those for h = +∞ (i.e. the Dirichlet case).

In order to prove this theorem, we follow the strategy due to Pleijel, [START_REF] Pleijel | Remarks on Courant's nodal line theorem[END_REF]. It is therefore necessary to estimate the number of nodal domains whose boundaries intersect the boundary of the square on at least a non-trivial interval. For such nodal domains, we cannot use the Faber-Krahn inequality for the Dirichlet problem. Nevertheless, there is a Faber-Krahn inequality for the Robin problem when h > 0 (see [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bucur | Faber-Krahn inequalities for the Robin-Laplacian: a free discontinuity approach[END_REF][START_REF] Daners | A Faber-Krahn inequality for Robin problems in any space dimension[END_REF]). We will see how this can be used for h sufficiently large in Subsection 3.3 and Section 4.

In Section 5, we analyse the number of nodal domains of Robin eigenfunctions in the general context of a planar domain with piecewise C 2,α boundary (α > 0). We obtain some semi-stability results for the number of nodal domains as the Robin parameter (h large) varies.

For the square, the results of Section 5 allow us to deal with the remaining case k = 5 which is not covered by Pleijel's strategy. In Section 6, we describe explicitly the situations where the eigenfunction corresponding to the fifth Robin eigenvalue has 2, 3, 4 nodal domains respectively (for h > 0 sufficiently large).

In a second paper, [START_REF] Gittins | Courant-sharp Robin eigenvalues for the square-Part II[END_REF], we hope to look at the situation where the Robin parameter h tends to 0 and to discuss the following conjecture.

Conjecture 1.3. There exists h 0 > 0 such that for 0 < h ≤ h 0 , the Courantsharp cases for the Robin problem are the same, except the fifth one, as those for h = 0 (i.e. the Neumann case) .

In light of the results of [START_REF] Pleijel | Remarks on Courant's nodal line theorem[END_REF][START_REF] Bérard | Dirichlet eigenfunctions of the square membrane: Courant's property, and A. Stern's and A. Pleijel's analyses[END_REF][START_REF] Helffer | Nodal domains in the square-the Neumann case[END_REF] and of the previous asymptotic results, a key question is to what extent is it possible to follow the Courant-sharp (Neumann) eigenvalues with h = 0 to Courant-sharp Robin eigenvalues as h → +∞? We prove a first general result concerning the possible crossings between curves corresponding to Robin eigenvalues.

We then focus on the cases k = 9, 25. For the case k = 9, we investigate if there exist critical values h * 9 (S), respectively h * 9 (S), after which the Robin eigenvalue λ 9,h (S) is not Courant-sharp, respectively before which it is Courantsharp, the next question being whether we have h * 9 = h * 9 . This question will be addressed in Subsection 7.2. In Subsection 7.3, for k = 25, we show that there exists h * 25 such that λ 25,h is not Courant-sharp for h < h * 25 , and we also investigate the structure of the nodal partitions of S for h sufficiently large.

Acknowledgements:

We are very grateful to Pierre Bérard, Dorin Bucur and Thomas Hoffmann-Ostenhof for useful discussions, and to Alexander Weisse for introducing us to the mathematics software system "SageMath" and helping us to produce some graphs of the Robin eigenvalues of the square. KG acknowledges support from the Max Planck Institute for Mathematics, Bonn, from October 2017 to July 2018.

All figures presented in this paper were created using SageMath.

2 Formulas for the eigenvalues and eigenfunctions of the Robin Laplacian for a rectangle.

Main formulas.

Here we follow the description given in [START_REF] Grebenkov | Geometrical structure of Laplacian eigenfunctions[END_REF] and we specialise to 2 dimensions. For rectangles Ω = (0, 1 ) × (0, 2 ) ⊂ R 2 and (x, y) ∈ Ω, an orthonormal basis for the Robin problem is given by

u p,q,h (x, y) = u p,h (x)u q,h (y), (2.1)
where, for p, q ∈ N (where N is the set of the non-negative integers), u p is the (p + 1)-st eigenfunction of the Robin problem in (0, 1 ):

u p,h (x) = sin(α p (h)x/ 1 ) + α p (h) h 1 cos(α p (h)x/ 1 ) ,
and similarly for u q,h (y) with y ∈ (0, 2 ). One should assume α p (h) = 0 (resp. α q (h) = 0), which holds for h = 0. For h = 0, the solution is trivial, hence not the right one! Here α p = α p (h) is the solution in [pπ, (p + 1)π) of

2α p h 1 cos α p + 1 - (α p ) 2 h 2 2 1 sin α p = 0 . (2.2)
The Robin eigenvalues are then given by

α p 1 2 + α q 2 2
.

So in 2 dimensions, the Robin eigenvalues correspond to pairs of non-negative integers (p, q). We analyse the 1D-situation in more detail and delete the reference to p, q, h. We note that the condition (2.2) reads (for h = 0 and α = 0),

α h = ± sin α + α h cos α .
In this way one understands the symmetry properties of the eigenfunctions better (see Lemma 2.1).

One also obtains the localisation of the eigenvalues in the following way.

If we consider the symmetric case,

α h = sin α + α h cos α , we get 2α h sin 2 α 2 = sin α , which leads to α tan α 2 = h . (2.3)
Similarly, if we consider the antisymmetric case, α h =sin α + α h cos α , we get 2α h cos 2 α 2 = -sin α ,

which leads to α h = -tan α 2 . (2.4)
With these formulas in mind, we get simpler expressions for the eigenfunctions. In the first case, we observe that

u(x) = sin(αx/ ) + α h cos(αx/ ) = sin(αx/ ) + cotan( α 2 ) cos(αx/ ) = 1 sin α 2 cos( αx -α 2 ) .
In the second case, we observe that

u(x) = sin(αx/ ) + α h cos(αx/ ) = sin(αx/ ) -tan( α 2 ) cos(αx/ ) = 1 cos α 2 sin( αx -α 2 ) .
In this way, we clearly see the symmetry properties of the eigenfunctions and we are closer to the Neumann case by considering x → cos( αxα 2 ) or x → sin( αxα 2 ) as eigenfunctions.

α 0 (h) α 1 (h) α 2 (h) Figure 1: Solutions α 0 (h), α 1 (h), α 2 (h) for h ≤ 100 .
The first case corresponds to p even. When h = 0, we have α = pπ and cos( αxα 2 ) = (-1) p 2 cos( pπx ). The second case corresponds to p = 2n + 1 odd (n ∈ N). When h = 0, we have α = pπ and sin( αxα

2 ) = ± cos( pπx ). By setting = π and then translating x → x + π 2 , we have that the Robin eigenfunctions (assuming h > 0) of the square S := (-π 2 , π 2 ) 2 are given by (2.1) with u p (x) = 1 sin

αp 2 cos α p x π , (2.5)
when p is even, and

u p (x) = 1 cos αp 2 sin α p x π , (2.6)
when p is odd. In Figure 1, we plot α 0 (h), α 1 (h), α 2 (h) for h ≤ 100 .

Particular cases

k = 1, 2, 3, 4.
We recall from the introduction that λ 1,h and λ 2,h (which for eigenfunctions of the form u p,q (x, y) correspond to (p, q) = (0, 0), (1, 0) respectively) are Courantsharp via Courant's nodal domain theorem and orthogonality of eigenfunctions. We note that λ 3,h (S) is not Courant-sharp since it corresponds to the case where (p, q) = (0, 1) so λ 3,h (S) = λ 2,h (S). Consider λ 4,h (S) with h > 0. Then p = q = 1 and the corresponding eigenfunction is

u 1,1 (x, y) = 1 cos 2 α1 2 sin α 1 x π sin α 1 y π ,
for (x, y) ∈ (-π 2 , π 2 ) 2 . We see that x = 0 and y = 0 are nodal lines of u 1,1 (x, y) which partition S into 4 nodal domains. There cannot be any further nodal lines of u 1,1 (x, y) as these would give rise to additional nodal domains so we would get a contradiction to Courant's nodal domain theorem. Therefore λ 4,h (S) with h ≥ 0 is Courantsharp. Hence, from this point onwards, we are only interested in the remaining eigenvalues, i.e. in the eigenvalues λ n,h (S) with n ≥ 4. Note that, due to the monotonicity of the Robin eigenvalues with respect to h, we have for n ≥ 4,

λ n,h (S) ≥ λ 4,h (S) ≥ λ 4,0 (S) = 2 .
(2.7)

Symmetry properties.

The use of symmetries was quite powerful in the context of the Neumann case, [START_REF] Helffer | Nodal domains in the square-the Neumann case[END_REF], via an argument due to Leydold, [START_REF] Leydold | Knotenlinien und Knotengebiete von Eigenfunktionen[END_REF]. That is, a Courant nodal theorem for eigenfunctions that satisfy certain symmetry properties. In addition, the number of nodal domains inherits some particular properties from these symmetries. The goal of this subsection is to show that this invariance by symmetry is common to all the Robin problems on the interval and the square.

Symmetry of Robin eigenfunctions in 1D.

We recall that h = 0 corresponds to the Neumann case and h = +∞ corresponds to the Dirichlet case. The Robin condition for [-2 , 2 ] reads

du dx (-/2) = hu(-/2) , du dx ( /2) = -hu( /2) .
We also observe the following invariance by symmetry.

Lemma 2.1. If u is an eigenfunction of the 1D-Robin problem, the function ũ(x) = u(-x) is also an eigenfunction of the same problem.

Hence, we necessarily have (using the conservation of the norm) u(-x) = ±u(x). Moreover, if u(0) = 0, we have u(-x) = u(x) and if u (0) = 0 we get u(-x) = -u(x). Therefore, the eigenfunctions u p (see (2.5) and (2.6)) are alternately symmetric and antisymmetric:

u p (-x) = (-1) p u p (x) , (2.8)
like in the Dirichlet or Neumann case.

Symmetry of Robin eigenfunctions in 2D.

In 2D, we now consider the possible symmetries of a general eigenfunction associated with the eigenvalues λ n,h of (-π 2 , π 2 ) 2 which reads,

u(x, y) = i,j:λ n,h (S)=π -2 (α 2 i +α 2 j ) a ij u i (x)u j (y) , (2.9)
where u p (or u p,h if we want to mention the reference to the Robin parameter) is the (p + 1) -st eigenfunction of the h-Robin problem in (-π 2 , π 2 ). By considering the transformation (x, y) → (-x, -y), we obtain

u(-x, -y) = i,j:λ n,h (S)=π -2 (α 2 i +α 2 j )
a ij (-1) i+j u i (x)u j (y) .

(2.10) Remark 2.2. We note that if (i+j) is odd for any pair (i, j) such that λ n,h (S) = π -2 (α 2 i +α 2 j ), then we get by (2.10), u(-x, -y) = -u(x, y) and as a consequence u has an even number of nodal domains. As we shall see later, other symmetries related to the finite group generated by the identity and the symmetries (x, y) → (-x, y) and (x, y) → (x, -y) can be considered (see [START_REF] Gittins | Courant-sharp Robin eigenvalues for the square-Part II[END_REF]).

In what follows, we obtain an upper bound for the number of Courant-sharp Robin eigenvalues of S via arguments that do not depend on the parameter h.

Upper bound for the number of Courant-sharp

Robin eigenvalues of a square.

In this section, we prove h-independent bounds for the number of Courant-sharp Robin eigenvalues. This was indeed the first step proposed by Pleijel [START_REF] Pleijel | Remarks on Courant's nodal line theorem[END_REF] in the Dirichlet case to reduce the analysis of the Courant-sharp cases to the analysis of finitely many eigenvalues. His proof was a combination of the Faber-Krahn inequality and the Weyl formula. In the Neumann case considered in [START_REF] Helffer | Nodal domains in the square-the Neumann case[END_REF], a new difficulty arises as it is not possible to apply the Faber-Krahn inequality to the elements of the nodal partition whose boundaries touch the boundary of the square at more than isolated points. In this section, we extend the analysis to the Robin case.

Lower bound for the Robin counting function.

Recall that for λ > 0, the Robin counting function for the corresponding eigenvalues of Ω is defined as

N R,h Ω (λ) := #{k ∈ N : k ≥ 1, λ k,h (Ω) < λ}. (3.1)
Similarly we have the Dirichlet counting function

N D Ω (λ) := #{k ∈ N : k ≥ 1, λ D k (Ω) < λ}, (3.2)
and the Neumann counting function

N N e Ω (λ) := #{k ∈ N : k ≥ 1, λ N k (Ω) < λ}. (3.3)
Due to the monotonicity of the Robin eigenvalues with respect to h ∈ [0, +∞), it is rather easy to have a lower bound for the N R,h Ω (λ). In particular, we have

N R,h Ω (λ) ≥ N R,+∞ Ω (λ) = N D Ω (λ)
. We also recall that for the Neumann counting function of S, we have

π 4 λ + 2 √ λ + 1 ≥ N N e S (λ) > π 4 λ, (3.4)
and for the Dirichlet counting function of S, if λ ≥ 2, we have by [START_REF] Pleijel | Remarks on Courant's nodal line theorem[END_REF],

N D S (λ) > π 4 λ -2 √ λ + 1 . (3.5)
Assume that λ ≥ 2 (this is true for n ≥ 4 by (2.7)). Then, by (3.5) and monotonicity of the Robin eigenvalues with respect to h,

N R,h S (λ) ≥ N D S (λ) > π 4 λ -2 √ λ + 1. (3.6)
With λ = λ n,h > λ n-1,h and Ψ an associated eigenfunction, (3.6) becomes

n > π 4 λ n,h -2 λ n,h + 2 . ( 3.7)
We now work analogously to the proof of Proposition 2.1 in [START_REF] Helffer | Nodal domains in the square-the Neumann case[END_REF]. Denote by Ω inn the union of nodal domains of Ψ whose boundaries do not touch the boundary of Ω (except at isolated points), and µ inn (Ψ) the number of nodal domains of Ψ in Ω inn . Similarly denote by Ω out the nodal domains in Ω \ Ω inn , and µ out (Ψ) the number of nodal domains of Ψ in Ω out . We have that

µ inn (Ψ) = µ(Ψ) -µ out (Ψ)
and we require an upper bound for µ out (Ψ).

Counting the number of nodal domains touching the boundary for the Robin problem.

We give a proof which holds for all the Robin problems in the square, except the Dirichlet case. We make use of the following theorem that is due to Sturm, [START_REF] Sturm | Mémoire sur une classe d'équations à différences partielles[END_REF][START_REF] Bérard | Sturm's theorem on zeros of linear combinations of eigenfunctions[END_REF].

Theorem 3.1 [START_REF] Sturm | Mémoire sur une classe d'équations à différences partielles[END_REF]. Let u = a m u m + • • • + a n u n be a non-trivial linear combination of eigenfunctions of the 1D-Robin problem in (-π 2 , π 2 ), with 1 ≤ m ≤ n, and {a j , m ≤ j ≤ n} real constants such that a 2 m + • • • + a 2 n = 0. Then, the function u has at least (m -1), and at most (n -1) zeros in (-π 2 , π 2 ). As observed originally by Å. Pleijel [START_REF] Pleijel | Remarks on Courant's nodal line theorem[END_REF], the analysis of the zeros of linear combinations of eigenfunctions appear in the following context. We observe that if an eigenfunction associated with λ n,h (see (2.9)) satisfies the Robin condition on the square, then its restriction to one side satisfies the Robin condition relative to the interval and is not zero (except of course in the Dirichlet case). In general, when the multiplicity is not one, this is no longer an eigenfunction but a linear combination of eigenfunctions on the segment (-π 2 , π 2 ). For example, the restriction to one side of the square, say x = π 2 , is a linear combination of eigenfunctions on the segment (-π 2 , π 2 ):

u(π/2, y) = i,j:λ n,h (S)=π -2 (α 2 i +α 2 j )
a ij u i (π/2)u j (y) .

We can then use Theorem 3.1 which gives a lower-bound on the number of zeros of u(π/2, y) in (-π 2 , π 2 ) by

i n (h) := min(i : λ n,h (S) = π -2 (α 2 i + α 2 j )) ,
and an upper-bound by

j n (h) := max(j : λ n,h (S) = π -2 (α 2 i + α 2 j )) . (3.8)
We have

λ n,h (S) = (α 2 in(h) + α 2 jn(h) )/π 2 ≥ i n (h) 2 + j n (h) 2 ≥ j n (h) 2 ,
which gives that

j n (h) ≤ λ n,h (S) .
We can argue in the same way for the other sides of the square. Therefore, the number of zeros of u(x, y) on the boundary of S is bounded from above by 4 λ n,h (S). Coming back to the number of "boundary" nodal domains, we have the following lemma.

Lemma 3.2. Let λ be a Robin eigenvalue of S with h < +∞. If Ψ is a Robin eigenfunction associated to λ, then

µ out (Ψ) ≤ 4 √ λ .
(3.9)

Remark 3.3. There are other proofs given in Pleijel [START_REF] Pleijel | Remarks on Courant's nodal line theorem[END_REF] and [START_REF] Helffer | Nodal domains in the square-the Neumann case[END_REF], but the one given above is much more general and not restricted to two-dimensional eigenspaces (and also not based on an explicit knowledge of the eigenfunctions).

On the other hand, the claim in [START_REF] Helffer | Nodal domains in the square-the Neumann case[END_REF] is much more involved. It says that taking the whole boundary into consideration, the number of points on the boundary in the nodal set of an eigenfunction cos θ u i (x)u j (y) + sin θ u j (x)u i (y) (i = j) is comparable with i + j (See Section 5 of [START_REF] Helffer | Nodal domains in the square-the Neumann case[END_REF]). The proof1 is restricted to eigenfunctions whose corresponding eigenvalues have multiplicity 2. It would be interesting to prove the same result for the Robin case for h < +∞ .

3.3

Upper bound for Courant-sharp Robin eigenvalues of a square.

By Lemma 3.2, we have

µ inn (Ψ) ≥ µ(Ψ) -4 λ n,h . (3.10)
Now, Ω inn = i ω inn i is a finite union of nodal domains of Ψ. Assuming that Ω inn is not empty, we get, on each ω inn i , by Faber-Krahn (see [START_REF] Pleijel | Remarks on Courant's nodal line theorem[END_REF]), that

A(ω inn i ) πj 2 ≥ 1 λ n,h , (3.11)
where A(ω inn i ) denotes the area of ω inn i and j denotes the first positive zero of the Bessel function J 0 . Adding, and invoking (3.10), we find By invoking the upper bound of (3.4), we obtain an upper bound for n such that λ n,h (S) < 598. Indeed, suppose λ n,h (S) < 598, then We remark that the above arguments do not depend on the Robin parameter h. In the sections that follow, we consider the case where h is large and improve the result.

π j 2 = A(S) πj 2 > A(Ω inn ) πj 2 ≥ µ inn (Ψ) λ n,h ≥ µ(Ψ) -4 λ n,h λ n,h , from which we obtain π j 2 ≥ µ(Ψ) -4 λ n,h λ n,h . ( 3
n -1 = N R,h S (λ n,h (S)) = #{k ∈ N : k ≥ 1, λ k,h (S) < λ n,h (S)} ≤ π 4 λ n,h (
4 Analysis as h → +∞.

In this section we show that for h sufficiently large, the Courant-sharp Robin eigenvalues of the square are the same as those in the Dirichlet case, [START_REF] Pleijel | Remarks on Courant's nodal line theorem[END_REF][START_REF] Bérard | Dirichlet eigenfunctions of the square membrane: Courant's property, and A. Stern's and A. Pleijel's analyses[END_REF], that is the first, second and fourth, except possibly the fifth which we deal with in Section 5. We first briefly revisit the strategy that was used by Pleijel for the Dirichlet problem.

Pleijel's approach for Dirichlet.

Let us come back to Pleijel's argument. We recall from (3.5)

that if λ n ≥ 2 is Courant-sharp, then n > π 4 λ n -2 λ n + 2 . (4.1)
On the other hand, if λ n is Courant-sharp, the Faber-Krahn inequality gives the necessary condition

n λ n ≤ πj -2 < 0.54323 . (4.2)
Recall that j is the smallest positive zero of J 0 the Bessel function of order 0, and that πj 2 is the ground state energy of the disc of area 1. Combining (4.1) and (4.2), leads to the inequality

πj -2 > π 4 -2λ -1 2 n + 2λ -1 n , (4.3)
and to

λ n ≤ 50 . (4.4)
Then the proof is achieved in the following steps (see [START_REF] Bérard | Dirichlet eigenfunctions of the square membrane: Courant's property, and A. Stern's and A. Pleijel's analyses[END_REF] for the full details).

• By a direct computation of the quotient of n λn , it is possible to eliminate all the eigenvalues except for n = 1, 2, 4, 5, 7 and 9 .

• The eigenvalues for n = 7 and n = 9 are eliminated by symmetry arguments (analogously to Remark 2.2).

• The final step is to analyse the fifth eigenvalue for which a specific analysis of the nodal structure can be done (see [START_REF] Bérard | Dirichlet eigenfunctions of the square membrane: Courant's property, and A. Stern's and A. Pleijel's analyses[END_REF]).

In the subsections that follow, we work through these steps and investigate the extent to which they still work for h large.

4.2 Faber-Krahn for the Robin case.

We recall the result of Bossel-Daners [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Daners | A Faber-Krahn inequality for Robin problems in any space dimension[END_REF], which asserts that the Robin eigenvalues of the Laplacian satisfy the following Faber-Krahn inequality. For a Lipschitz domain ω ⊂ R 2 and h > 0,

λ 1,h (ω) ≥ λ 1,h (D ω ), (4.5)
where

D ω ⊂ R 2 is a disc such that A(D ω ) = A(ω).
We note that this analysis is only interesting for the nodal domains whose boundary meet the boundary of Ω along at least some arc. For the interior nodal domains, the approach via the standard Faber-Krahn inequality still applies. For the boundary domains, we have mixed boundary conditions with Robin on some arcs and Dirichlet on the remaining arcs. For a lower bound, by monotonicity, we can indeed use the Faber-Krahn inequality (with a Robin boundary condition on all the boundary).

Consider a scaling of the domain ω by t > 0, tω := {tx ∈ R 2 : x ∈ Ω}. As observed by Antunes, Freitas and Kennedy, [START_REF] Antunes | Asymptotic behaviour and numerical approximation of optimal eigenvalues of the Robin Laplacian[END_REF], the Robin eigenvalues satisfy the following scaling property.

λ n,h (ω) = t 2 λ n,h/t (tω), (4.6)
A serious issue here is that the scaling also affects the Robin parameter. So, in particular, replacing D by D 1 , the disc of area 1, we have

λ 1,h (D ω ) = λ 1,hA(ω) 1 2 (D 1 )/A(ω) . ( 4

.7)

When h = +∞, the reference is λ 1,+∞ (D 1 ). In the Robin case, if we start from h large, we will not necessarily have hA(ω)

12 large if we use this inequality with ω a "boundary" nodal domain. Hence we have to be careful in the application of the Faber-Krahn argument. This is actually the main difficulty. We recall the asymptotic behaviour of the first Robin eigenvalue as the Robin parameter tends to +∞ or to 0 (see, for example, [START_REF] Grebenkov | Geometrical structure of Laplacian eigenfunctions[END_REF]).

We recall that λ 1, h(D 1 ) → λ 1,+∞ (D 1 ) = πj 2 as h → +∞, and that there exists c > 0 such that, as h → +∞,

λ 1, h(D 1 ) = λ 1,+∞ (D 1 ) - c h + O 1 h2 . (4.8)
We also recall that there exists d > 0 such that as h → 0,

λ 1, h(D 1 ) = d h + O( h2 ) . (4.9)
We give the proof for completion. To determine the first eigenvalue for the disc of area 1 and radius π -1/2 , one looks for an eigenfunction of the form J 0 (απ 1/2 r) where the corresponding eigenvalue is πα 2 . For the asymptotic behaviour near h = 0 or h = +∞, we use the Taylor expansion of J 0 or J 0 at α = 0 and α = j. The Robin condition 2 reads

απ 1/2 J 0 (α) + hJ 0 (α) = 0 .
We recall that J 0 (0) = 0 and J 0 (0) < 0. We get for h ≥ 0, for the first solution α 2 π 1/2 J 0 (0) ∼ -hJ 0 (0) . Hence the corresponding eigenvalue satisfies as h → 0,

λ 1,h (D 1 ) = -(2π 1/2 J 0 (0))/(J 0 (0)) h + O(h 2 ) .
We also have J 0 (j) = 0 and J 0 (j) = 0. With τ = 1 h , we write τ απ 1/2 J 0 (α) + J 0 (α) = 0, and expanding at α = j, we obtain:

α = j -π 1 2 j τ + O(τ 2 ) ,
and

πα 2 = πj 2 -2π 3 2 j 2 τ + O(τ 2 ) .
The proof gives an explicit value for the constants c and d in (4.8) and (4.9).

We will apply the Faber-Krahn inequality to a nodal domain of a Robin eigenfunction u = u n,h associated with λ n,h . We observe that an eigenfunction u can be extended to all of R 2 as a solution ũ of -∆ũ = λũ (we have an explicit expression as a trigonometric polynomial). Hence the nodal sets of ũ have a nice local structure (see P. Bérard [START_REF] Bérard | Inégalités isopérimétriques et applications. Domaines nodaux des fonctions propres[END_REF] for a survey) and have the same properties as in the Dirichlet case. In particular, these nodal sets are locally Lipschitz domains (actually with piecewise analytic boundary). If we observe that a nodal set of u is the intersection of a nodal set of ũ with the square S, we immediately deduce that the ω inn i are Lipschitz domains.

The regularity of the "boundary domains" ω out j has to be analysed. By Lemma 3.2, the nodal set intersects the boundary finitely many times, so ∂ω out j consists of a finite number of arcs belonging either to S or to ∂S. So we can apply Theorem 4.1 of [START_REF] Bucur | Faber-Krahn inequalities for the Robin-Laplacian: a free discontinuity approach[END_REF]. Alternatively, we can use the strategy given in Section 3 of [START_REF] Kennedy | An isoperimetric inequality for the second eigenvalue of the Laplacian with Robin boundary conditions[END_REF] to obtain (4.5) for these domains (see also [30, p. 3620]). We will discuss the regularity of the nodal domains further in Section 5.

Note also that for a "boundary" domain ω out j , u| ω out j satisfies a mixed Robin-Dirichlet condition on its boundary but we can use the monotonicity with respect to the Robin parameter which leads to

λ n,h ≥ λ 1,h (ω out j ) , (4.10)
and then use the pure Robin Faber-Krahn inequality.

4.3 Pleijel's approach as h → +∞ .

In light of what was recalled in Subsection 4.1 for h = +∞, we now consider the different steps in the limit h → +∞. We first recall that the eigenvalues depend continuously on h until +∞, in particular ∀n ∈ N, lim

h→+∞ λ n,h = λ D n . (4.11)
We keep the notation of the previous section. If we are in the Courant-sharp situation, then µ(u) = n, where u is an eigenfunction associated with λ n,h . If there exists ω inn i such that A(ω inn i ) ≤ A(S)/n, we are done like in the Dirichlet case. We combine the latter inequality with inequality (3.11) to obtain (4.2). Together with (4.1), this gives λ n,h ≤ 50. In particular, for these eigenvalues n is finite and using (4.11) we get that for h sufficiently large, (4.2) is not satisfied for n = 1, 2, 4, 5, 7, 9.

If not, the situation is more delicate, but we can assume that there exists

ω out j such that A(ω out j ) ≤ A(S)/n , (4.12)
and we take one of smallest area with this property. Combining (3.7), (4.10), (4.5), (4.12) and (4.7), we find that

A(S) λ 1,hA(ω out j ) 1/2 (D 1 ) > π 4 - 2 λ n,h + 2 λ n,h . (4.13)
Here, comparing with (4.3), we need to have h := hA(ω out j ) 1/2 large enough if we want to arrive at the same conclusion as for the Dirichlet case. So we have to find a lower bound for A(ω out j ) 1/2 . This seems difficult, at least with explicit lower bounds. We will use our initial h-independent upper bound from the previous section. Hence, we can assume in this Courant-sharp situation, that n ≤ 520 . (4.14)

Below, we do not try to obtain explicit constants. The first claim is that, according to (4.9), there exist c 1 > 0 and h 1 > 0 such that

λ 1, h ≥ c 1 h if 0 ≤ h ≤ h 1 .
We now assume that n ≤ 520 and λ n,h is Courant-Sharp and get

λ 520,∞ ≥ λ n,h ≥ c 1 hA(ω out j ) -1 2 , if hA(ω out j ) 1 2 ≤ h 1 . This gives a contradiction if c 1 h 2 > h 1 λ 520,∞ . Hence, assuming h > h 1 2 1 c -1 2 1 λ 1 2
520,∞ , we can now assume that hA(ω out j )

1 2 > h 1 . Now, we have λ 520,∞ ≥ λ n,h ≥ A(ω out j ) -1 λ 1,h1 (D 1 ) , which implies A(ω out j ) ≥ λ 1,h1 (D 1 )/λ 520,∞ .
This gives the existence of c 2 > 0 such that A(ω out j ) ≥ c 2 (see also Lemma 5.3). Coming back to (4.13), we have

π 2 λ 1,c 1/2 2 h (D 1 ) > π 4 - 2 λ n,h + 2 λ n,h . (4.15)
Hence for h large enough, we also get in this case that λ n,h ≤ 50 (compare with inequality (4.3)).

We can now follow the proof of Pleijel for the Dirichlet case. The first step was to achieve (assuming h large enough) the restriction to the three cases left by Pleijel. This step now follows (using the continuity (4.12) of the eigenvalues with respect to h as h → +∞ as already observed in the previous case). The second step is to rule out the cases λ 7,h (S) and, for h sufficiently large, λ 9,h (S). Here the symmetry argument due to Leydold holds in the same way as for the Dirichlet case [START_REF] Bérard | Dirichlet eigenfunctions of the square membrane: Courant's property, and A. Stern's and A. Pleijel's analyses[END_REF] for the two cases corresponding to the seventh and the ninth Robin eigenvalues. We briefly recall the relevant particular case of the argument due to Leydold. Lemma 4.1. Let 0 ≤ h < +∞. Suppose that λ n,h (S) is a Robin eigenvalue with corresponding eigenfunction defined in (2.9). Suppose that n is odd and that the conditions of Remark 2.2 are satisfied. Then λ n,h (S) is not Courant-sharp.

We know indeed by the standard Courant nodal domain theorem that the number of nodal domains is not larger than n and by Remark 2.2 that it is even. Hence the number is less than n.

As an application, we observe that any eigenfunction corresponding to the seventh Robin eigenvalue is a linear combination of u 2,1 (x, y) and u 1,2 (x, y) (see Figure 3 and Appendix A) and that 1+2 is odd. So λ 7,h (S) is not Courant-sharp for any h ≥ 0.

Similarly, for h large, any eigenfunction corresponding to the ninth Robin eigenvalue is a linear combination of u 3,0 (x, y) and u 0,3 (x, y) (see Figure 3 and Appendix A) and 0 + 3 is odd.

Hence at this stage, we have proved the following proposition. Proposition 4.2. There exists h 1 > 0 such that for h ≥ h 1 , the Courant-sharp cases for the Robin problem are the same, except possibly for k = 5, as those for h = +∞ . So, having in mind what was done for the Dirichlet case [START_REF] Bérard | Dirichlet eigenfunctions of the square membrane: Courant's property, and A. Stern's and A. Pleijel's analyses[END_REF], in order to prove Theorem 1.2 for h large enough it remains to count the number of nodal domains of any eigenfunction corresponding to the fifth eigenvalue. This will be analysed in Section 6 as a direct consequence of Section 5.

5 A general perturbation argument.

Preliminary discussion.

We analyse a θ-dependent family Φ h,θ of eigenfunctions, more explicitly Φ h,θ,p,q (x, y) = cos θ u p,h (x)u q,h (y) + sin θ u p,h (y)u q,h (x) , for (x, y) ∈ (-π 2 , π 2 ) 2 . For most of the arguments in this section, we will not use the explicit expression of the eigenfunction, but only the property that Φ h,θ is a very smooth family of eigenfunctions (with respect to h and θ) where, for h ∈ (0, +∞], Φ h,θ is an eigenfunction of the h-Robin Laplacian associated with a smooth eigenvalue λ(h). The parameter θ, which above belongs to R/(2πZ), could also be thought of as belonging to some open neighbourhood of some point θ 0 in R. In addition, most of the arguments extend to more general domains. We consider the case of bounded, planar domains with piecewise C 2,α (α > 0) boundary.

For h = +∞ (or h = h 0 > 0) and θ = θ 0 , we assume that the number of nodal domains is known (for example, that the corresponding eigenvalue is not Courant-sharp). The aim of this section is to prove that by perturbation (i.e. for | 1 h -1 h0 | + |θ -θ 0 | small enough) the number of nodal domains cannot increase (see Proposition 5.7). The proof involves various general statements which are interesting in a more general context3 , hence not restricted to the case of the square.

Robin Faber-Krahn inequality revisited.

Proposition 5.1. Given h 1 > 0 and M > 0, we consider a smooth family Φ h,θ of h-Robin eigenfunctions on Ω, where Ω is a connected, bounded set with piecewise C 2,+ boundary4 , λ(h) ≤ M and h ∈ I ⊂ [h 1 , +∞) (I being a finite or infinite interval). Any nodal domain of Φ h satisfies the h-Faber-Krahn inequality.

Remark 5.2. We note that the square satisfies the assumptions of Proposition 5.1 but in this case there is a more direct proof. As in Subsection 4.2, we indeed observe that Φ h,θ admits an extension Φ h,θ to R 2 such that -∆ Φ h,θ = λ(h) Φ h,θ . This gives more information about the local nodal structure of Φ h,θ up to the boundary (actually in a neighbourhood of S ).

Proof. The proposition holds for an open set with C 2,+ boundary (hence without corners) as a direct application of Theorem B.1 in Appendix B. Hence Ω is a domain with rectifiable boundary of finite length and thus the Faber-Krahn inequality holds by [START_REF] Bucur | Faber-Krahn inequalities for the Robin-Laplacian: a free discontinuity approach[END_REF] (as mentioned in Subsection 4.2). The same is true for the nodal domains whose boundaries do not touch a corner.

It remains to treat the corners. The Dirichlet case was addressed by Helffer, Hoffmann-Ostenhof and Terracini in [START_REF] Helffer | Nodal domains and spectral minimal partitions[END_REF]. This argument involves a local conformal change of coordinates which leads to the analysis of an operator with higher singularities.

We do not know an appropriate reference for the Robin case. The guess is that the boundary of a nodal domain (whose closure touches the corner) consists of Lipschitz arcs of finite length, including the arcs for which one end touches the corner, which would allow us to use the Robin Faber-Krahn inequality for Lipschitz domains. Instead we use that according to [START_REF] Bucur | A variational approach to the isoperimetric inequality for the Robin eigenvalue problem[END_REF], the h-Faber-Krahn inequality holds for any open set with finite area. In this general case, the first eigenvalue is defined as in Definition 4.2 of [START_REF] Bucur | A variational approach to the isoperimetric inequality for the Robin eigenvalue problem[END_REF]. It is also proven in [START_REF] Bucur | A variational approach to the isoperimetric inequality for the Robin eigenvalue problem[END_REF] that with this choice of definition, this eigenvalue is not larger than any other definition given in a more regular situation. Proof. This follows directly from the h-Faber Krahn inequality. If ω is a nodal domain of Φ h satisfying the assumptions of the lemma, we have

M ≥ λ(h) ≥ λ(h 0 ) ≥ λ 1,h0 (D ω ) = λ 1,h0A(ω) 1 2 (D 1 )/A(ω) ∼ d h 0 /A(ω) 1 2 . (5.1)
This shows that as soon as we avoid the Neumann situation, the ground state energy in a domain ω tends to +∞ as the area of the domain tends to 0.

On the nodal set at the boundary.

Proposition 5.4. Under the assumptions of Proposition 5.1, there exists C > 0 such that, for any h ∈ I and any θ, the number of zeros of Φ h,θ at the boundary is less than C. Remark 5.5. In the case of the square the proposition follows from Sturm's theorem.

Proof. We will use the Euler formula with boundary. The conditions for its application are satisfied by using Theorem B.1 and it reads as follows (see, for example, [START_REF] Hoffmann-Ostenhof | Bounds on the multiplicity of eigenvalues for fixed membranes[END_REF]). Proposition 5.6. Let Ω be an open set in R 2 with C 2,+ boundary, u a Robin eigenfunction with k nodal domains, N (u) its zero-set. Let b 0 be the number of components of ∂Ω and b 1 be the number of components of N (u) ∪ ∂Ω. Denote by ν(x i ) and ρ(y i ) the numbers of curves ending at critical point x i ∈ N (u),

respectively y i ∈ N (u) ∩ ∂Ω. Then k = 1 + b 1 -b 0 + xi ν(x i ) 2 -1 + 1 2 yi ρ(y i ) . (5.2)
In our application, we immediately obtain that the number ρ(u) of boundary points (actually counted with multiplicity) in the nodal set of u satisfies

ρ(u) ≤ 2k -2 .
To achieve the proof, we observe that by Courant's nodal domain theorem, k is less than the minimal labelling of λ(h) and that this labelling is uniformly bounded if λ(h) is uniformly bounded. By monotonicity, this labelling is indeed bounded by the maximal labelling of an eigenvalue λ j (h 1 ) satisfying λ j (h 1 ) ≤ M .

It remains to treat what is going on in the neighbourhood of a corner x c . We first show that there cannot exist an infinite sequence of zeros of u in the boundary (outside the corner) tending to the corner x c . Indeed, by Proposition 5.1, similarly to the proof of Lemma 5.3, there exists some sufficiently small > 0 such that any line starting from one of these zeros (which necessarily belongs to the boundary of one nodal domain) should cross ∂D(x c , ) ∩ Ω transversally and only once. Hence the number of points is finite, and moreover not greater than the cardinality of N (u) ∩ D(x c , ) ∩ Ω. Observing that, by Lemma 5.3, the number of nodal domains of u in Ω is the same as the number of nodal domains of u in Ω \ D(x c ), we can apply the Euler Formula in Ω \ D(x c ) and get the same bound.

On the variation of the cardinality of the nodal domains by perturbation.

We assume that Ω is a bounded, planar domain with piecewise C 2,+ boundary. Our main result is the following proposition.

Proposition 5.7. Under the previous assumptions on Ω and the family Φ h,θ , let ρ(h, θ) denote the cardinality of the nodal domains of Φ h,θ . For any θ 0 , h 0 ∈ (0, +∞], there exists

η 0 > 0 such that if | 1 h -1 h0 | + |θ -θ 0 | < η 0 , then ρ(h, θ) ≤ ρ(h 0 , θ 0 ) .
We prove this proposition in the following subsections by analysing what is going on at the interior critical points and at the boundary points of the zero set.

Analysis in a neighbourhood of an interior point.

We treat what is going on at an interior point z 0 . We assume that z 0 is a critical point of Φ h0,θ0 associated with an eigenvalue λ(h 0 ). We choose 0 > 0 small enough such that • D(z 0 , 0 ) ⊂ Ω;

• Lemma 5.3 applies with M > λ(h 0 );

• the circle C(z 0 , 0 ) crosses the 2 half-lines emanating from z 0 transversally at 2 points z j (h 0 , θ 0 ) (j = 1, . . . , 2 ).

Here we have used the general results on the local structure of an eigenfunction of the Laplacian (see [START_REF] Bérard | Inégalités isopérimétriques et applications. Domaines nodaux des fonctions propres[END_REF] and Appendix B).

Lemma 5.8. With the previous notations and assumptions of Lemma 5.3, there exists

η 0 > 0 such that if | 1 h -1 h0 | + |θ -θ 0 | < η 0
, then the number of nodal domains of Φ h,θ intersecting the disc D(z 0 , 0 ) cannot increase.

Proof. If we look at the nodal structure inside D(z 0 , 0 ), we have 2 local nodal domains.

By local nodal domain of an eigenfunction Φ h,θ , we mean the nodal domains of the restriction of Φ h,θ to D(z 0 , 0 ). We note that any local nodal domain belongs to a global nodal domain but that two distinct local nodal domains can be included in the same global nodal domain.

In this case, there exists a path γ in Ω joining these two local domains on which Φ h,θ is positive (or negative), which necessarily will not be included in D(z 0 , θ 0 ).

Starting from (h 0 , θ 0 ) we now look at a small perturbation. By considering the restriction of Φ h,θ to the circle ∂D(z 0 , 0 ), we observe that the 2 zeros of Φ h,θ in ∂D(z 0 , 0 ) move very smoothly, we denote them by z j (h, θ).

We indeed observe that the tangential derivative of Φ h0,θ0 at each point z j (h 0 , θ 0 ) is not zero (again we use the general results for eigenfunctions, in particular the transversal property, see Appendix B). By perturbation, this condition is still true if we choose η 0 small enough. Hence the restriction of Φ h,θ changes sign at each point z j (h, θ). Moreover, there are 2 local domains ω j (h, θ) of Φ h,θ with the property that ∂ω j (h, θ) intersects ∂D(z 0 , 0 ) along the arc (z j (h, θ), z j+1 (h, θ)) (with the convention that j + 1 is 1 for j = 2 ).

We now observe that if ω j (h 0 , θ 0 ) and ω j (h 0 , θ 0 ) belong to the same nodal domain (j = j ), the property remains true for (h, θ) sufficiently close to (h 0 , θ 0 ) (i.e. for η 0 in the lemma sufficiently small). If, for (θ 0 , h 0 ), ω j (h 0 , θ 0 ) and ω j (h 0 , θ 0 ) do not belong to the same nodal domain, then there are two cases

• either the situation is unchanged by perturbation;

• or they belong after perturbation to the same nodal domain via a new path in D(z 0 , 0 ).

In the second case, the number of nodal domains touching ∂D(z 0 , 0 ) is decreasing.

On the other hand, by Lemma 5.3, any nodal domain that intersects D(z 0 , 0 ) crosses ∂D(z 0 , 0 ). This achieves the proof. Remark 5.9. If =2, Φ h0,θ0 is a Morse function whose Hessian has two nonzero eigenvalues of opposite sign. Then, for 0 small enough, Φ h,θ remains a Morse function for η 0 small enough and admits a unique critical point z h,θ in D(z 0 , η 0 ). Then there are four local nodal domains if Φ h,θ (z h,θ ) = 0 and three local nodal domains if Φ h,θ (z h,θ ) = 0 (see Subsection 6.3.1 for a detailed proof ).

Analysis at the boundary.

It remains to control what is going on at the boundary. We consider a point z 0 ∈ ∂Ω such that z 0 is a zero of Φ h0,θ0 which in addition is assumed to be critical when h 0 = +∞.

We first assume that we avoid the corners and successively consider three cases:

• h 0 = +∞, perturbation only in θ.

• 0 < h 0 < +∞, general perturbation.

• h 0 = +∞, general perturbation.

In the first case, the proof follows the same argument as that used in the proof of Lemma 5.8 and uses the local structure of a Dirichlet eigenfunction at the boundary (see [START_REF] Bérard | Inégalités isopérimétriques et applications. Domaines nodaux des fonctions propres[END_REF] and Appendix B).

For the second case, considering the proof of Lemma 5.8 once again, we choose 0 > 0 sufficiently small such that z 0 is the only boundary point in the nodal set. Then the proof goes in the same way.

In the third case, the situation is more delicate due to the complete vanishing of Φ +∞,θ0 on the boundary, which should not be the case for h < +∞. To deal with this, we need the following lemma.

Lemma 5.10. Let θ = θ 0 and Z bnd denote the intersection of the nodal set of Φ +∞,θ0 with the boundary. Then for any > 0 there exists h * such that the set {d(z, ∂S) < } ∩ {d(z, Z bnd ) > } does not meet the zero set of Φ h,θ for any h * ≤ h < +∞ and any θ such that |θ -

θ 0 | < 1 h * .
In other words we have some nodal stability up to the boundary as h → +∞.

Proof. We consider the following two cases. At a regular point of the boundary.

We consider a point z 0 of the boundary (or a closed interval I in the boundary) which is not a critical point for Φ +∞,θ0 . By perturbation, this is still true for |θ -θ 0 | small. In this case the normal derivative of Φ +∞,θ for z 0 ∈ I does not vanish, and to fix the ideas we can assume that

∂ ν Φ +∞,θ (z 0 , θ) > c > 0
(the other case would be treated similarly). By continuity, replacing c by c 2 , this is still true for Φ h,θ , z in a h-independent neighbourhood of I and 1 h small enough.

On the other hand, we know that Φ h,θ satisfies the Robin condition:

∂ ν Φ h,θ (z 0 , θ) + hΦ h,θ (z 0 , θ) = 0 . Hence Φ h,θ (z 0 , θ) = - 1 h ∂ ν Φ h,θ (z 0 , θ) < 0 .
This implies that there exists a neighbourhood of I and η > 0 such that, for

1 h + |θ -θ 0 | < η , Φ h,θ is negative (actually < -c 2h ).
At a corner.

After translation, we assume that the corner is at (0, 0). We also assume that (0, 0) does not belong to the nodal set of Φ +∞,θ0 and that Φ +∞,θ0 < 0 in Ω near the corner. We now use the previous argument outside of (0, 0). For 0 > 0 small enough we can take η > 0 small enough such that, for 1 h + |θ -θ 0 | < η, Φ h,θ (x, y) < 0 for {x 2 + y 2 = 2 0 } ∩ Ω. Suppose now that Φ h,θ (x, y) > 0 for some (x, y) ∈ D((0, 0), 0 ). Then there is a nodal domain inside D((0, 0), 0 ) and this is excluded by Lemma 5.3 provided that we have chosen 0 sufficiently small. Remark 5.11. We have not proven in full generality that Φ h,θ is negative at the boundary near the corner but this is not required. We do not know what occurs if the corner belongs to the zero set. If the corner is not in the zero-set of the Dirichlet eigenfunction, we can prove by the previous argument that this is still the case for h large enough. In the case of the square, we get immediately that ∂ 2

x,y Φ +∞,θ0 (0, 0) < 0 . We now estimate Φ h,θ (0, 0). Using the Robin condition, we obtain that

Φ h,θ (0, 0) ∼ h -2 ∂ 2
x,y Φ h,θ (0, 0) By perturbation, we also have ∂ 2

x,y Φ h,θ (0, 0) < 0 . This implies Φ h,θ (0, 0) < 0 .

This leads to the following result when z 0 ∈ ∂Ω. We assume that z 0 is a critical point of Φ +∞,θ0 associated with an eigenvalue λ(∞). We choose 0 small enough such that • Lemma 5.3 applies with M > λ(h 0 );

• C(z 0 , 0 ) ∩ Ω crosses the half-lines emanating from z 0 transversally at points z j (h 0 , θ 0 ) (j = 1, . . . , ).

Here we have used the general results for the local structure of an eigenfunction of the Dirichlet Laplacian (see [START_REF] Bérard | Inégalités isopérimétriques et applications. Domaines nodaux des fonctions propres[END_REF], see also [START_REF] Helffer | Nodal domains and spectral minimal partitions[END_REF] for the case with corners).

Lemma 5.12. With the previous notation and assumptions of Lemma 5.3, there exists

η 0 > 0 such that if | 1 h -1 h0 | + |θ -θ 0 | < η 0
, then the number of nodal domains of Φ h,θ intersecting the disc D(z 0 , 0 ) cannot increase. If = 1, the number of nodal domains equals two and remains fixed.

Application to the square.

We come back to the case of the square and prove Theorem 1.2. To this end, having in mind Proposition 4.2, it is sufficient to obtain the following. Proposition 5.13. There exists h 0 > 0 such that for any h > h 0 , any eigenfunction corresponding to 1 π 2 (α 0 (h) 2 + α 2 (h)) 2 has 2, 3, or 4 nodal domains (as in the Dirichlet case). Hence for h > h 0 , λ 5,h is not Courant-sharp.

Proof. The property is indeed true for h = +∞ and, by the results of the preceding sections, the number of nodal domains cannot increase and is necessarily > 1.

In the next section, we carry out a deeper analysis for the eigenfunction associated with the fifth eigenvalue, where we count the nodal domains case by case. For some cases, the proof will use the explicit properties of the eigenfunctions Φ h,θ (see below).

In relation to Proposition 5.13, we note that by choosing non-critical values of θ we can obtain that 2, 3 and 4 nodal domains are attained for h 0 large enough.

6 Particular case k = 5 .

Main statement.

Looking at the fifth eigenvalue corresponding to the pair (0, 2), which is Courantsharp for Neumann and not Courant-sharp for Dirichlet, we consider the family of eigenfunctions in (-π 2 , + π 2 ) 2 with θ ∈ (-π, π]:

Φ h,θ,0,2 (x, y) := cos θ cos(α 0 (h)x/π) cos(α 2 (h)y/π) + sin θ cos(α 2 (h)x/π) cos(α 0 (h)y/π) . (6.1)

Up to changing the sign of the eigenfunction, it is sufficient to consider θ ∈ [0, π). We prove the following proposition.

Proposition 6.1. There exists h 0 > 0 such that for any h > h 0 , any eigenfunction corresponding to 1 π 2 (α 0 (h) 2 + α 2 (h)) 2 has 2, 3, or 4 nodal domains (as in the Dirichlet case). More precisely, there are three critical values θ

* j (h) ∈ [0, π) (j = 1, 2, 3) such that θ * 1 (h) = arctan - 1 q 2 (h) , θ * 2 (h) = π 2 -θ * 1 (h) , θ * 3 = 3π 4 ,
where

q 2 (h) = cos α2 2 cos α0 2 ,
and such that Φ h,θ has:

• 3 nodal domains for θ ∈ [0, θ * 1 (h)]; • 2 nodal domains for θ ∈ (θ * 1 (h), θ * 2 (h)); • 3 nodal domains for θ ∈ [θ * 2 (h), θ * 3 );
• 4 nodal domains for θ = θ * 3 ;

• 3 nodal domains for θ ∈ (θ * 3 , π).

Note that for the whole family of eigenfunctions, we have symmetry with respect to the two axes. In addition, the corresponding eigenvalue 1 π 2 (α 0 (h) 2 + α 2 (h) 2 ) is the fifth eigenvalue for any h ∈ [0, +∞] (due to monotonicity of the Robin eigenvalues with respect to h and the table given in Appendix A). For h = 0, we have α 0 (0) = 0 and α 2 (0) = 2π.

6.2 The Dirichlet case.

= arctan(1/3), π 8 , π 4 , 3π 8 , θ * 2 = π 2 -arctan(1/3), π 2 , 5π 8 , θ * 3 = 3π 4 , 7π
8 correspond to the purple, magenta, blue, grey, green, black, orange, teal, red, navy curves respectively. For h = +∞, i.e. in the Dirichlet case, we have α 0 (+∞) = π and α 2 (+∞) = 3π. The figures of Pockel, [START_REF] Pockels | Über die partielle Differentialgleichung -∆u -k 2 u = 0 and deren Auftreten in mathematischen Physik[END_REF], give the various possibilities as a function of θ. We refer to [START_REF] Bérard | Dirichlet eigenfunctions of the square membrane: Courant's property, and A. Stern's and A. Pleijel's analyses[END_REF] for a more rigorous mathematical analysis but note that Pockel gives all the possible topologies. He also gives the pictures for the θ corresponding to transitions between these topologies. In Figure 2, we plot the fifth Dirichlet eigenfunction Φ +∞,θ,0,2 (x, y) = cos θ cos(x) cos(3y) + sin θ cos(3x) cos(y) for (x, y) ∈ (-π 2 , π 2 ) 2 and various values of θ. The critical values of θ corresponding to a change in the number of interior critical points or the number of boundary critical points in the nodal set are θ * 1 = arctan(1/3), θ * 2 = π 2 -arctan(1/3), and θ * 3 = 3π 4 .

As was proven in [START_REF] Bérard | Dirichlet eigenfunctions of the square membrane: Courant's property, and A. Stern's and A. Pleijel's analyses[END_REF] and can be seen in Figure 2, the fifth Dirichlet eigenfunction has either 2, 3 or 4 nodal domains. More precisely, we have for θ ∈ [0, π):

• 3 nodal domains for θ ∈ [0, θ * 1 ]; • 2 nodal domains for θ ∈ (θ * 1 , θ * 2 ); • 3 nodal domains for θ ∈ [θ * 2 , θ * 3 );
• 4 nodal domains for θ = θ * 3 ;

• 3 nodal domains for θ ∈ (θ * 3 , π). In what follows, we prove that this holds for h sufficiently large.

Application of Section 5.

For h large enough, we analyse ψ(θ, x, y) := Φ h,θ,0,2 (x, y) = cos θ cos α 0 x π cos α 2 y π

+ sin θ cos α 2 x π cos α 0 y π .
This solution has a double symmetry with respect to x → -x and y → -y.

Interior critical points.

We can look at the critical points of ψ as a function of θ. In the case of Dirichlet, the only possible critical point is for x = y = 0 and can only occur for cos θ + sin θ = 0 (we assume θ = Z π 2 ). For cos θ + sin θ = 0, x = ±y belong to the zero set of ψ. We show that the zero set is exactly given by x = ±y. We observe that the Hessian of ψ at (x, y) = (0, 0) is

H (x,y)=(0,0) = cos θ π 2 α 2 2 -α 2 0 0 0 α 2 0 -α 2 2 ,
which has negative determinant so (x, y) = (0, 0) is a non-degenerate critical point of ψ. We see that H (x,y)=(0,0) has one positive eigenvalue and one negative eigenvalue, so the Morse index of the critical point (0, 0) is 1. By the Morse Lemma, in a neighbourhood U of (0, 0), there is a diffeomorphism φ = (u, v) :

U → V ⊂ R 2 with φ(0, 0) = (0, 0) such that ψ := ψ • φ -1 has the form ψ(u, v) = ψ(0, 0) -u 2 + v 2 = cos θ + sin θ -u 2 + v 2 .
So we see immediately that the critical point (0, 0) is isolated. With the condition that cos θ + sin θ = 0, the zero set is given by u = ±v. Since φ is a bijection and x = ±y is contained in the zero set of ψ, the zero set of ψ is given by x = ±y. More generally, the same proof gives that the zero set of ψ(θ, •) -(cos θ + sin θ) is given near (0, 0) by x = ±y. We remark that in this case there are 4 nodal domains.

Boundary edge.

Considering the boundary edge x = π 2 , we have that either y = ± π 2 is in the nodal set, in which case there are 4 nodal domains by symmetry, or y = ± π 2 is not in the nodal set. In the latter case, Theorem 3.1 gives that there are at most two points on the boundary edge x = π 2 that are in the nodal set. If there are exactly two such points in the nodal set, then this corresponds to 3 nodal domains. If there are no boundary points in the nodal set, then this corresponds to 2 nodal domains. For example, see Figure 2.

Double point on the boundary.

We now analyse what is going on at the double point on the boundary. This occurs for Dirichlet when tan θ = 1 3 and for y = 0. Here the situation is simple (see [START_REF] Pockels | Über die partielle Differentialgleichung -∆u -k 2 u = 0 and deren Auftreten in mathematischen Physik[END_REF]). We observe that y = 0 is a double point for tan θ = -1 q2(h) . From Ψ(θ, π 2 , y) = 0, we have

cos α 2 y π + t cos α2 2 cos α0 2 cos α 0 y π = 0 .
The critical t = tan θ is defined by t = -1/q 2 (h) with

q 2 (h) = cos α2 2 cos α0 2 . Hence t = 1 3 + O( 1 h
), and we have near y = 0,

y 2 = c + O 1 h t + 1 q 2 (h)
.

Again, this is the perturbation of a Morse function depending on the parameters h and θ with the particularity that when ψ = 0 and y = 0, the critical point is always ( π 2 , 0). We remark that in this case there are 3 nodal domains.

Interior critical points for any h > 0.

In this subsection, we show that there are no other critical points than (0, 0) without any restriction on h > 0. It is immediate that (0, 0) is a critical point and we get the same condition as in the Dirichlet case. Writing ψ = 0 and ∇ψ = 0, we get as a necessary condition that

α 2 tan α 2 x π = α 0 tan α 0 x π , α 2 tan α 2 y π = α 0 tan α 0 y π . (6.2) Lemma 6.2. Let α 0 and α 2 satisfy (2.3). For x ∈ (-π 2 , + π 2 ), α 0 tan(α 0 x/π) = α 2 tan(α 2 x/π) if and only if x = 0. Proof. Let us look at the function [0, π 2 ] x → W (x) = α 0 sin(α 0 x) cos(α 2 x) -α 2 sin(α 2 x) cos(α 0 x) .
Up to some multiplicative renormalisation of the eigenfunctions, we recognise the Wronskian of the eigenfunctions u 0 and u 2 . But for the Wronskian, we have

W (x) = (λ 0 -λ 2 )u 0 (x)u 2 (x) .
Now we observe that W (0) = 0 and that by (2.3), W ( π 2 ) = 0. Moreover W has a unique critical point in (0, π 2 ) at the first zero of u 2 . Hence W (x) cannot vanish except at x = 0 and π 2 . It is clear that this implies that (0, 0) is the only possible critical point in (-π 2 , π 2 ) 2 . The condition that ψ(0, 0) = 0 implies cos θ + sin θ = 0.

7 Analysis of crossings.

In this section, we analyse the possible crossings of two curves h → λ p,q,h (S) and h → λ p ,q ,h (S) defined in an interval of [0, +∞). This is indeed quite important as we want to follow the labelling of these eigenvalues when h varies.

7.1 A general result.

Proposition 7.1. For distinct pairs (p, q) and (p , q ), with p ≤ q and p ≤ q , there is at most one value of h in [0, +∞) such that λ p,q,h (S) = λ p ,q ,h (S).

Proof Suppose that λ p,q,h (S) = λ p ,q ,h (S). Without loss of generality, suppose p < p ≤ q < q. Consider the variation of

(0, +∞) h → σ(h) := 1 π 2 α p (h) 2 + α q (h) 2 -α p (h) 2 -α q (h) 2 .
The zeros of σ correspond to the values of h for which the curves corresponding to (p, q), (p , q ) intersect. To analyse its variation, we note that

σ (h) = 2 π 2 α p (h)α p (h) + α q (h)α q (h) -α p (h)α p (h) -α q (h)α q (h) .
Now, we deduce from (2.3) and (2.4), that h → α k (h) satisfies the differential equation

α k α k hπ + α 2 k 2 + h 2 π 2 2 = π , (7.1)
which implies

α k α k hπ + α 2 k 2 + h 2 π 2 2 = πα 2 k . (7.2)
We introduce for h > 0 and k ∈ N,

a k (h) = hπ + α 2 k 2 + h 2 π 2 2 > 0 .
We deduce

σ (h) = 2 π α 2 p a p + α 2 q a q - α 2 p a p - α 2 q a q = - 4 π hπ + h 2 π 2 2 1 a p + 1 a q - 1 a p - 1 a q .
We now assume that σ(h) = 0, which implies

a p + a q = a p + a q .
This gives

σ (h) = - 4 π hπ + h 2 π 2 2
(a p + a q )(a p a q -a p a q ) (a p a q a p a q ) .

So the sign of σ (h) is the sign of a p a q -a p a q . For > 0, we can now write a p = a p -and a q = a q + , and compute a p a q -a p a q = (a p -)(a q + ) -a p a q = (a p -a q ) -2 < 0 .

Since the derivative of σ has constant sign, there can be at most one point of intersection.

Remark 7.2. The proof of Proposition 7.1 shows that if p < p ≤ q < q and λ p,q,h * = λ p ,q ,h * for some h * ≥ 0, then the map

h → π -2 (α p (h) 2 + α q (h) 2 -α p (h) 2 -α q (h) 2 )
is increasing for h > h * . Hence the curve π -2 (α p (h

) 2 + α q (h) 2 ) is below the curve π -2 (α p (h) 2 + α q (h) 2 ) for h > h * .

The eigenvalue λ 9,h (S).

The ninth eigenvalue of the Neumann Laplacian for the square is Courant-sharp, [START_REF] Helffer | Nodal domains in the square-the Neumann case[END_REF], and corresponds to the eigenvalue 2 2 + 2 2 = 8. This eigenvalue is simple and corresponds to the labelling (2, 2). The eigenfunction reads Φ 0,θ,2,2 (x, y) = cos 2x cos 2y , for (x, y) ∈ (-

π 2 , π 2 ) 2 .
It is easy to see that the Courant-sharp property is still true for h small enough. By deformation, the eigenfunction is

Φ h,θ,2,2 (x, y) = cos(α 2 (h)x/π) cos(α 2 (h)y/π),
with corresponding eigenvalue 2 π 2 (α 2 (h)) 2 . The nodal structure is given by

α 2 (h)x π = - π 2 , α 2 (h)x π = π 2 , α 2 (h)y π = - π 2 , α 2 (h)y π = π 2 ,
hence for this eigenfunction and for h ∈ [0, +∞), there are nine nodal domains as long as 2 π 2 (α 2 (h)) 2 is the ninth eigenvalue.

The issue is to follow its labelling and we observe that when h = +∞ the eigenvalue is 18 and, according to the ordered list of the Dirichlet eigenvalues, has minimal labelling 11 (see Appendix A). Because 9 < 11, this eigenfunction is NOT Courant-sharp for h sufficiently large.

On the other hand the eigenvalue 1 π 2 (α 0 (h) 2 + α 3 (h) 2 ) which has minimal labelling 10 for h = 0 arrives with labelling 9 at h = +∞. Hence some transition occurs for at least one h * 9 > 0 which satisfies

α 0 (h) 2 + α 3 (h) 2 = 2α 2 (h) 2 .
By Proposition 7.1, there is at most one point of intersection between the curves corresponding to (2, 2) and [START_REF] Bérard | Inégalités isopérimétriques et applications. Domaines nodaux des fonctions propres[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications. Domaines nodaux des fonctions propres[END_REF]. We recall that α 0 (0) = 0 , α 1 (0) = π , α 2 (0) = 2π , α 3 (0) = 3π and that α 0 (+∞) = π , α 1 (+∞) = 2π , α 2 (+∞) = 3π , α 3 (+∞) = 4π , so α 0 (h) 2 + α 3 (h) 2 is increasing from 9π 2 to 17π 2 when 2α 2 (h) 2 goes from 8π 2 to 18π 2 .

In order to show that the curves corresponding to the pairs (2, 2), (3, 0) do not intersect the curves corresponding to the other pairs, we consider the table in Appendix A.

From above, we see that the eigenvalues corresponding to the pairs (3, 3), (4, 2), [START_REF] Antunes | Asymptotic behaviour and numerical approximation of optimal eigenvalues of the Robin Laplacian[END_REF][START_REF] Bérard | Dirichlet eigenfunctions of the square membrane: Courant's property, and A. Stern's and A. Pleijel's analyses[END_REF] and so on are all larger than or equal to 18. So we need to consider the eigenvalues corresponding to the pairs (3, 1), (3, 2), (4, 0), (4, 1) and show that they do not correspond to the ninth, tenth or eleventh eigenvalues for any h > 0. Numerically we find that,

λ 3,1,h (S) = λ 1,3,h (S) ≥ 18 for h > 11.4225, λ 3,2,h (S) = λ 2,3,h (S) ≥ 18 for h > 2.6288, λ 4,0,h (S) = λ 0,4,h (S) ≥ 18 for h > 1.2668, λ 4,1,h (S) = λ 1,4,h (S) ≥ 18 for h > 0.4208.
So we are left to consider h ≤ 11.4225.

From below, we see that the eigenvalues corresponding to the pairs (0, 0), (1, 0), (0, 1), [START_REF] Adams | Sobolev Spaces[END_REF][START_REF] Adams | Sobolev Spaces[END_REF] are smaller than or equal to 8 for all 0 < h < +∞. So we need to consider the eigenvalues corresponding to the pairs (2, 0), (2, 1) and show that they do not correspond to the ninth, tenth or eleventh eigenvalues for any h > 0. Numerically we find that, λ 2,2,h (S) ≥ 13 for h > 2.9804 , λ 3,0,h (S) =λ 0,3,h (S) ≥ 13 for h > 3.5468 .

So we are left to consider h ≤ 3.5468 < 11.4225 .

With the table from Appendix A in mind, we now plot the Robin eigenvalues of the square (α m (h) 2 + α n (h) 2 )/π 2 for h ≤ 12 corresponding to the pairs (0, 0), (1, 0), (1, 1), (2, 0), (2, 1), (2, 2), (3, 0), (3, 1), (3, 2), (4, 0), (4, 1).

From Figure 3, we see that for h ≤ 12 the curves corresponding to the pairs (2, 2), (3, 0) do not intersect the curves corresponding to the other pairs. By Proposition 7.1, the curves corresponding to (2, 2) and (3, 0) intersect for a unique value of h = h * 9 > 0 . Since u 2,2 (x, y) is an eigenfunction corresponding to λ 9,h * 9 (S) that has 9 nodal domains, we have proved:

Proposition 7.3. There exists h * 9 > 0 such that λ 9,h is Courant-sharp for 0 ≤ h ≤ h *
9 and not Courant-sharp for h > h * 9 . By the bisection method, we compute h * 9 numerically and find that h * 9 ∼ 1.6967. By the above, λ 9,h is given by the pair (2, 2) for h ≤ h * 9 and the pair (3, 0) for h > h * 9 . Also, λ 10,h is given by the pair (0, 3) and λ 11,h is given by the pair (3, 0) for h ≤ h * 9 and the pair (2, 2) for h > h * 9 . This shows that whether the eigenfunction corresponding to a Robin eigenvalue of the square is an odd function or an even function depends on h (in the case where there are crossings).

. The intersection between the curves corresponding to (

For example, for λ 9,h with h ≤ h * 9 , we have that u 2,2 (-x, -y) = u 2,2 (x, y). On the other hand, for h > h * 9 , u 3,0 (-x, -y) = -u 3,0 (x, y) and u 0,3 (-x, -y) = -u 3,0 (x, y) .

So any linear combination of u 3,0 (x, y) and u 0,3 (x, y) is antisymmetric with respect to the transformation (x, y) → (-x, -y) . Hence λ 9,h is not Courantsharp for h > h * 9 (via Lemma 4.1). For h = h * 9 , any eigenfunction corresponding to λ 9,h * 9 (S) is a linear combination of u 2,2 (x, y), u 3,0 (x, y) and u 0,3 (x, y), so in general it is neither symmetric nor antisymmetric with respect to the transformation (x, y) → (-x, -y).

The eigenvalue λ 25,h (S).

Similarly there are crossings between λ 25,h and λ 27,h . As for the ninth eigenvalue, we first show that the curves corresponding to the pairs (4, 3), (5, 1) do not intersect the curves corresponding to the other pairs by considering the table in Appendix A.

From above, we see that the eigenvalues corresponding to the pairs (5, 4), (6, 3), (7, 0) and so on are all larger than or equal to 41. So we need to consider the eigenvalues corresponding to the pairs (5, 2), (4, 4), (5, 3), (6, 0), (6, 1), [START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF][START_REF] Antunes | Asymptotic behaviour and numerical approximation of optimal eigenvalues of the Robin Laplacian[END_REF] and show that they do not correspond to λ 25,h (S), λ 26,h (S), λ 27,h (S), λ 28,h (S) for any h > 0. Numerically we find that, λ 5,2,h (S) = λ 2,5,h (S) ≥ 41 for h > 12.6664, λ 4,4,h (S) ≥ 41 for h > 4.9398, λ 5,3,h (S) = λ 3,5,h (S) ≥ 41 for h > 3.4557, λ 6,0,h (S) = λ 0,6,h (S) ≥ 41 for h > 3.8230,

λ 6,1,h (S) = λ 1,6,h (S) ≥ 41 for h > 2.0624, λ 6,2,h (S) = λ 2,6,h (S) ≥ 41 for h > 0.4016.
So we are left to consider h ≤ 12.6664.

The eigenvalues corresponding to the pairs (3, 2), (3, 1) and below in the table in Appendix A are smaller than or equal to 25 for all 0 < h < ∞. So we need to consider the eigenvalues corresponding to the pairs (4, 0), (4, 1), (3, 3), (4, 2), (5, 0) and show that they do not correspond to λ 25,h (S), λ 26,h (S), λ 27,h (S), λ 28,h (S) for any h > 0. Numerically we find that, λ 4,3,h (S) = λ 3,4,h (S) ≥ 37 for h > 11.5497 , λ 5,1,h (S) = λ 1,5,h (S) ≥ 37 for h > 15.3826 .

So we are left to consider h ≤ 12.6664 < 15.3826 .

With the table from Appendix A in mind, we now plot the Robin eigenvalues of the square (α m (h) 2 + α n (h) 2 )/π 2 for h ≤ 16 corresponding to the pairs (4, 0), ( The Robin eigenvalues of the square (α m (h) 2 + α n (h) 2 )/π 2 for h ≤ 16 corresponding to the pairs (4, 0), (4, 1), (3, 3), (4, 2), (5, 0), (5, 1), (4, 3), (5, 2), (4, 4), (5, 3), (6, 0), (6, 1), (6, 2).

From Figure 4, we see that for h ≤ 16 the curves corresponding to the pairs (4, 3), (5, 1) do not intersect the curves corresponding to the other pairs. We also note that the curves corresponding to (4, 4) and (6, 0) give rise to the same Dirichlet eigenvalue as h → ∞ (see the table in Appendix A). There are only two crossings in Figure 4. By Proposition 7.1, there exists a unique value h = h * 25 at which the crossing occurs. So λ 25,h is given by the pair (4, 3) for h ≤ h * 25 and by the pair (5, 1) for h > h * 25 . Hence, we have obtained: Proposition 7.4. There exists h * 25 > 0 such that λ 25,h is given by the pair (4, 3) for h ≤ h * 25 and by the pair (5, 1) for h > h * 25 .

By the bisection method, we compute h * 25 numerically and find that h * 25 ∼ 3.1317. We note that u 4,3 (x, y) is antisymmetric with respect to the transformation (x, y) → (-x, -y), while u 5,1 (x, y) is symmetric with respect to this transformation. From the first observation, Lemma 4.1 gives that any eigenvalue with corresponding eigenfunction a linear combination of u 4,3 (x, y) and u 3,4 (x, y) has an even number of nodal domains. Hence λ 25,h is not Courant-sharp for h < h * 25 and λ 27,h is not Courant-sharp for h > h * 25 . For h * 25 < h < ∞ we investigate whether λ 25,h (S) is Courant-sharp or not by considering the corresponding eigenfunctions u 5,1 (x, y) and u 1,5 (x, y).

For (x, y) ∈ (-π 2 , π 2 ) 2 , we consider the function Φ h,θ,5,1 (x, y) = cos θ sin(α 5 (h)x/π) sin(α 1 (h)y/π) + sin θ sin(α 1 (h)x/π) sin(α 5 (h)y/π) . (7.3)
We also note that the lines x = 0 and y = 0 belong to the nodal set of Φ h,θ,5,1 for any θ.

It is known that λ 25,∞ (S) is not Courant-sharp, [START_REF] Bérard | Dirichlet eigenfunctions of the square membrane: Courant's property, and A. Stern's and A. Pleijel's analyses[END_REF]. In addition, by Theorem 1.2, we know that, for h sufficiently large, λ 25,h (S) is not Courant-sharp.

In Figure 5, we plot the twenty-fifth Dirichlet eigenfunction Φ +∞,θ,5,1 (x, y) = cos θ sin(6x) sin(2y) + sin θ sin(2x) sin(6y) for θ = 0, arctan(1/3), π 4 , arctan(3), π 2 , 5π 8 , 3π 4 , 13π 16 , 15π 16 .

We observe that for θ = arctan(1/3) there are double points at (0, 0), (± π 4 , π 2 ), and (± π 4 , -π 2 ). There are triple points at (± π 4 , 0). Similarly for θ = arctan(3), there are double points at (0, 0), (± π 2 , π 4 ), and (± π 2 , -π 4 ), and triple points at (0, ± π 4 ). The eigenfunction associated with the fifth Dirichlet eigenvalue on (-π 2 , π 2 ) 2 is Φ+∞,θ,2,0 (x, y) = cos θ sin(3x) sin(y) + sin θ sin(x) sin(3y). We see that Φ+∞,θ,5,1 (x, y) := Φ +∞,θ,5,1

x

2 - π 4 , y 2 - π 4 = Φ +∞,θ,2,0 (x, y).
As in the proof of Lemma 4.2 of [START_REF] Helffer | Nodal domains in the square-the Neumann case[END_REF], Φ +∞,θ,5,1 can be constructed by taking its values in the square (0, π 2 ) 2 and folding evenly over (-π 2 , π 2 ) 2 , that is with respect to the axes x = 0 and y = 0. Compare Figure 2 with Figure 5.

We now consider the case where h > 0. In order to make a numerical comparison to the Dirichlet case, we choose h = 20 in what follows. This value of h is small enough that we see some differences compared to the Dirichlet case and large enough that we keep the asymptotic structure.

To determine the critical points on the side y = For h = 20, we compute numerically that θ t ≈ 1.2492655.

Using (2.4), we obtain the following asymptotic expansions for α 1 (h) and α 5 (h) when h → ∞.

α 1 (h) = 2π - 4 h + O 1 h 2 , (7.11)
α 5 (h) = 6π - 12 h + O 1 h 2 . (7.12)
Substituting these expansions into (7.6) and solving for x c (h) gives that

x c (h) = π 4 + 1 2h + O 1 h 2 . (7.13)
Using the above asymptotic expansions for α 1 (h), α 5 (h) and x c (h), we obtain that as h → ∞,

tan θ m (h) = 1 3 + O 1 h , tan θ t (h) = 3 + O 1 h , and tan π 2 -θ t (h) = cot θ t (h) = 1 3 + O 1 h .
For h = +∞, we have

θ m (+∞) = π 2 -θ t (+∞) = arctan 1 3 ,
and we deduce from above that

θ m (h) = arctan 1 3 + O 1 h , and π 2 -θ t (h) = arctan 1 3 + O 1 h .
At this stage we do not get any information about the sign of θ m + θ t -π 2 . For this, we observe that by (7.7) and (7.10),

tan θ m tan( π 2 -θ t ) = tan θ m tan θ t = α 5 sin α1 2 α 1 sin α5 2 .
From (2.4), we have that for n odd,

α n sin αn 2 = - hπ cos αn 2 . Hence we obtain tan θ m tan( π 2 -θ t ) = cos α1 2 cos α5 2 .
From the asymptotic expansions (7.11) and (7.12), we obtain that as h → ∞,

cos α 1 2 ∼ -1 - 2 h 2 + O 1 h 3 , cos α 5 2 ∼ -1 - 18 h 2 + O 1 h 3 .
We deduce that as h → ∞,

tan θ m tan( π 2 -θ t ) = 1 + 16 h 2 + O 1 h 3 . (7.14) Let δθ = θ m + θ t -π 2 . Then, observing that θ m -arctan 1 3 = O( 1 h ) and π 2 -θ t - arctan 1 3 = O( 1 h ), we get tan θ m = tan π 2 -θ t + δθ 10 9 1 + O 1 h .

Dividing by tan(

π 2 -θ t ) leads to tan θ m tan( π 2 -θ t ) = 1 + δθ 10 3 1 + O 1 h .
From (7.14), we then get

δθ 10 3 1 + O 1 h = 16 h 2 .
From this we deduce that

δθ = θ m (h) + θ t (h) - π 2 ∼ 24 5h 2 .
This gives the strict positivity of θ m (h)+θ t (h)-π 2 for h large enough, a property which is numerically satisfied for h = 20. In Figure 6, we plot h → g(h) := α5 sin(

α 1 2 ) α1 sin( α 5
2 ) for 20 ≤ h ≤ 500 and note that it approaches 1 from above.

2 )

α1 sin( α 5
2 )

for 20 ≤ h ≤ 500.

In Figure 7, we plot Φ h,θ,5,1 for h = 20 and θ = 0, π 2 -

θ t , θ = 1 2 ( π 2 -θ t + θ m ), θ m , π 4 , π 2 -θ m , θ = 1 2 ( π 2 -θ m + θ t ), θ t , π 2 , 5π 8 , 3π 4 , 13π 16
. From this figure, we make the following observations. For 0 ≤ θ < π 2 -θ t , there are 12 boundary critical points, 5 interior critical points and 12 nodal domains.

For θ = π 2 -θ t , there are 12 boundary critical points, 3 interior critical points and 12 nodal domains.

For π 2 -θ t < θ < θ m , there are 12 boundary critical points, 1 interior critical point and 8 nodal domains (see, for example, Part (a) of Figure 7 (maroon)).

For θ = θ m , there are 8 boundary critical points, 1 interior critical point and 8 nodal domains. For θ m < θ < π 2 -θ m , there are 4 boundary critical points, 1 interior critical point and 8 nodal domains.

For θ = π 2 -θ m , there are 8 boundary critical points, 1 interior critical point and 8 nodal domains.

For π 2 -θ m < θ < θ t , there are 12 boundary critical points, 1 interior critical points and 8 nodal domains (see, for example, Part (c) of Figure 7 (deep sky blue)).

For θ = θ t , there are 12 boundary critical points, 3 interior critical points and 12 nodal domains.

For θ t < θ < 3π 4 , there are 12 boundary critical points, 5 interior critical points and 12 nodal domains.

For θ = 3π 4 , there are 16 boundary critical points, 5 interior critical points and 16 nodal domains.

For 3π 4 < θ < π, there are 12 boundary critical points, 5 interior critical points and 12 nodal domains.

By comparing Figure 7 with Figure 8 below, we see that for h = 20 and θ = 3π 4 , the nodal structure for the twenty-fifth Robin eigenfunction is not obtained from the nodal structure of the fifth Robin eigenfunction (by the aforementioned folding procedure on the square (-π 2 , π 2 ) 2 which holds for the Dirichlet case). For example, this can be seen by comparing the navy curve in Part (d) of Figure 7 with the navy curve in Figure 8.

= 0 (blue), θ = π 2 -θ t (orange), θ = 1 2 ( π 2 -θ t + θ m ) (maroon), θ = θ m (magenta), θ = π 4 (red), θ = π 2 -θ m (purple), θ = 1 2 ( π 2 - θ m + θ t ) (deep sky blue), θ = θ t (teal), θ = π 2 (grey), 5π 8 (lime), θ = 3π 4 (navy), θ = 13π
16 (gold).

The numerical experiment discussed above suggests that there are no new transitions and no new critical points appear as h increases from h = 20 to we see that the eigenvalues λ 84,h , . . . , λ 92,h are indeed given by the pairs (9, 4), (7, 7), (10, 0), [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF]

By the above, we have that there exist 0 < h a < h b < h c < h d < +∞ such that the following hold.

For 0 ≤ h ≤ h a , the curve corresponding to (9, 4) lies below that corresponding to (7, 7) which lies below (10, 0), which lies below [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF] which in turn lies below [START_REF] Bucur | Faber-Krahn inequalities for the Robin-Laplacian: a free discontinuity approach[END_REF][START_REF] Adams | Sobolev Spaces[END_REF]. So for 0 ≤ h ≤ h a , λ 84,h = λ 85,h is given by (9, 4), λ 86,h by (7, 7), λ 87,h = λ 88,h by (10, 0), λ 89,h = λ 90,h by [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF] and λ 91,h = λ 92,h by [START_REF] Bucur | Faber-Krahn inequalities for the Robin-Laplacian: a free discontinuity approach[END_REF][START_REF] Adams | Sobolev Spaces[END_REF].

Similarly for h a ≤ h ≤ h b , λ 84,h = λ 85,h is given by (9, 4), λ 86,h by (7, 7), λ 87,h = λ 88,h by (10, 0), λ 89,h = λ 90,h by (10, 1) and λ 91,h = λ 92,h by [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF].

For h b ≤ h ≤ h c , λ 84,h = λ 85,h is given by (9, 4), λ 86,h = λ 87,h by (10, 0), λ 88,h by (7, 7), λ 89,h = λ 90,h by (10, 1) and λ 91,h = λ 92,h by [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF].

For h c ≤ h ≤ h d , λ 84,h = λ 85,h is given by (10, 0), λ 86,h = λ 87,h by (9, 4), λ 88,h by (7, 7), λ 89,h = λ 90,h by (10, 1) and λ 91,h = λ 92,h by [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF].

For h ≥ h d , λ 84,h = λ 85,h is given by (10, 0), λ 86,h = λ 87,h by (9, 4), λ 88,h = λ 89,h by (10, 1), λ 90,h by (7, 7) and λ 91,h = λ 92,h by [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF].

We compute numerically that h a ∼ 2.1209, h b ∼ 2.1864, h c ∼ 3.7786, and h d ∼ 5.2167. In Figure 10, we plot the curves corresponding to the pairs [START_REF] Bers | Local behavior of solutions of general linear elliptic equations[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF], [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bérard | Sturm's theorem on zeros of linear combinations of eigenfunctions[END_REF], (9, 3), (9, 4), (7, 7), (10, 0), [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF], (10, 1), (10, 2), [START_REF] Bucur | A variational approach to the isoperimetric inequality for the Robin eigenvalue problem[END_REF][START_REF] Bérard | Sturm's theorem on zeros of linear combinations of eigenfunctions[END_REF], [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bers | Local behavior of solutions of general linear elliptic equations[END_REF], (11, 0) for h ≤ 10.

We note that the curves corresponding to [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF] and (10, 0) give rise to the same Neumann eigenvalue when h = 0. In addition, the curves corresponding to [START_REF] Bucur | A variational approach to the isoperimetric inequality for the Robin eigenvalue problem[END_REF][START_REF] Bérard | Dirichlet eigenfunctions of the square membrane: Courant's property, and A. Stern's and A. Pleijel's analyses[END_REF]

We see that it is possible that the labelling could switch more than once for a given eigenvalue (that is, the eigenvalue could be given by more than two pairs).

A Comparison Dirichlet-Neumann.

In this appendix, we recall from [START_REF] Bérard | Dirichlet eigenfunctions of the square membrane: Courant's property, and A. Stern's and A. Pleijel's analyses[END_REF][START_REF] Helffer | Nodal domains in the square-the Neumann case[END_REF] the Neumann eigenvalues λ N k of S (in the left-hand side below) and the Dirichlet eigenvalues λ D k of S (in the righthand side below) for k ≤ 129. We also recall the pairs (m, n) corresponding to these eigenvalues and the values m 2 + n 2 of the eigenvalues. The purpose is to illustrate the values of k ≤ 129 for which there are crossings between the curves corresponding to the Robin eigenvalues of S. We use colours to emphasise this. For example, the Robin eigenvalue λ 2,2,h (S) starts as the ninth eigenvalue when h = 0 but as h → +∞, it corresponds to the eleventh eigenvalue. In this appendix, we prove some well-known results for the nodal set of an eigenfunction of the Neumann problem and extend them to the Robin problem.

Neumann m n m 2 + n 2 k 0 0 0 1 1 0 1 2,3 0 1 1 2,3 1
Although used in various contributions, for example [START_REF] Helffer | Nodal sets for the groundstate of the Schrödinger operator with zero magnetic field in a non simply connected domain[END_REF], no detailed proofs seem to be published for the Neumann problem. For the Dirichlet problem, see [START_REF] Hoffmann-Ostenhof | Bounds on the multiplicity of eigenvalues for fixed membranes[END_REF] and [START_REF] Helffer | Nodal domains and spectral minimal partitions[END_REF] where the case with corners or cracks is also considered. In addition, we require these results under weaker regularity assumptions on the boundary.

B.1 Main statement. for some non-zero a ∈ R, where (r, ω) are polar coordinates of x around x 0 . The angle ω is chosen so that the tangent to the boundary at x 0 is given by the equation sin ω = 0.

3. The nodal set N (u) is the union of finitely many, C 2 -immersed circles in Ω, and C 1 -immersed lines which connect points of ∂Ω. Each of these immersions is called a nodal line. Note that self-intersections are allowed. The connected components of Ω \ N (u) are called nodal domains.

4. If u has a zero of order at a point x 0 ∈ Ω then exactly segments of nodal lines pass through x 0 . The tangents to the nodal lines at x 0 dissect the full circle of radius B(x 0 , α) (for α > 0 small enough) into 2 equal angles.

5. If u has a zero of order at a point x 0 ∈ ∂Ω then exactly segments of nodal lines meet the boundary at x 0 . The tangents to the nodal lines at x 0 are given by the equation cos ω = 0, where ω is chosen as in (B.2).

B.2 Proof of the theorem.

The C 2 -regularity of u up to the boundary is a consequence of standard Schauder estimates (see [START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Order[END_REF]). The proof now is in four steps.

It remains to verify that (ii) holds. We compare this condition with the property established in the previous theorem. By Theorem B.2, we get a control of ∇ψ in W 1,p in any ball B(0, r) hence by Sobolev's embedding theorem we have, as soon as p > 2, the control of ∇ψ in L ∞ (B(0, r)) (see, for example, Part II Case C' of Theorem 5.4 in [START_REF] Adams | Sobolev Spaces[END_REF]). It remains to control the constants appearing in the continuity of this injection. To do this, for r > 0, we introduce a cut-off χ(x/r) where χ = 1 on B(0, 1) and suppχ ⊂ B(0, 2), and apply the standard Sobolev embedding theorem to χ(x/r)∂ xi ψ and use the two estimates from Theorem B.2. We get This is sufficient to apply the lemma.

Remark B.5.

There is some controversy regarding Cheng's paper [START_REF] Cheng | Eigenfunctions and nodal sets[END_REF] when applied to a dimension larger than 2. The reason is that a harmonic homogeneous polynomial does not always satisfy Item (iv) when the dimension is larger than 2 (see Appendix E in [START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF]).

B.3 Remarks.

We note that all the proofs are local and the results can be obtained locally if we have the corresponding local regularity property. The proofs also work in the Dirichlet case (with a different condition on ω). Instead of the reflection argument, in order to construct ũ, we can introduce an extension via odd reflection: ũ(x 1 , x 2 ) = u(x 1 , x 2 ) for x 1 > 0 , -u(-x 1 , x 2 ) for x 1 < 0 . Analogously to the above, if u is an eigenfunction in C 2 (Ω) satisfying the Dirichlet condition, one can verify that ũ is in C 2 (R 2 ). Theorem B.6. Let Ω be an open set in R 2 with C 2,+ boundary and let u be a real-valued eigenfunction of the Laplacian with Dirichlet boundary conditions. Then u ∈ C 2 (Ω). Furthermore, u has the following properties: for some non-zero a ∈ R, where (r, ω) are polar coordinates of x around x 0 . The angle ω is chosen so that the tangent to the boundary at x 0 is given by the equation ω = 0.

3. The nodal set N (u) is the union of finitely many, C 2 -immersed circles in Ω, and C 1 -immersed lines which connect points of ∂Ω.

4. If u has a zero of order at a point x 0 ∈ Ω, then exactly segments of nodal lines pass through x 0 . The tangents to the nodal lines at x 0 dissect the full circle of radius B(x 0 , α) (for α > 0 small enough) into 2 equal angles.

5. If u has a zero of order at a point x 0 ∈ ∂Ω then exactly -1 segments of nodal lines meet the boundary at x 0 . The tangents to the nodal lines at x 0 are given by the equation sin ω = 0, ω = 0, π. We can, for example, refer to [START_REF] Helffer | Nodal domains and spectral minimal partitions[END_REF] for the Dirichlet case which gives the results (except C 2 regularity) under the weaker assumption that the boundary is piecewise C 1,+ .

3. 4

 4 Proof of Theorem 1.1.

  S) + 2 λ n,h (S) + 1 < 518.67. (3.14) Hence we have shown Theorem 1.1.

Lemma 5 . 3 .

 53 Let h 0 > 0 and M > 0. Then, under the same hypotheses as in Proposition 5.1, there exists 0 > 0 such that no nodal domain of an eigenfunction Φ h associated with λ(h) for the Robin problem with parameter h ≥ h 0 in some open set Ω and λ(h) ≤ M can have area less than 0 . (This includes the Dirichlet case).

Figure 2 :

 2 Figure 2: The fifth Dirichlet eigenfunction Φ +∞,θ,0,2 on the square (-π 2 , π 2 ) 2 , for various values of θ. The values θ = 0, θ * 1 = arctan(1/3), π 8 , π 4 , 3π 8 , θ * 2 =

Figure 3 :

 3 Figure3: The Robin eigenvalues of the square (α m (h) 2 + α n (h) 2 )/π 2 for h ≤ 12 corresponding to the pairs (0, 0), (1, 0), (1, 1), (2, 0), (2, 1), (2, 2), (3, 0), (3, 1),[START_REF] Bérard | Inégalités isopérimétriques et applications. Domaines nodaux des fonctions propres[END_REF][START_REF] Antunes | Asymptotic behaviour and numerical approximation of optimal eigenvalues of the Robin Laplacian[END_REF], (4, 0),[START_REF] Bérard | Dirichlet eigenfunctions of the square membrane: Courant's property, and A. Stern's and A. Pleijel's analyses[END_REF][START_REF] Adams | Sobolev Spaces[END_REF]. The intersection between the curves corresponding to (2, 2) and (3, 0) occurs at(1.6970, 11.4498).

Figure 4 :

 4 Figure4: The Robin eigenvalues of the square (α m (h) 2 + α n (h) 2 )/π 2 for h ≤ 16 corresponding to the pairs (4, 0), (4, 1),[START_REF] Bérard | Inégalités isopérimétriques et applications. Domaines nodaux des fonctions propres[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications. Domaines nodaux des fonctions propres[END_REF], (4, 2), (5, 0), (5, 1), (4, 3), (5, 2), (4, 4), (5, 3), (6, 0), (6, 1),[START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF][START_REF] Antunes | Asymptotic behaviour and numerical approximation of optimal eigenvalues of the Robin Laplacian[END_REF].

Figure 5 :

 5 Figure 5: The twenty-fifth Dirichlet eigenfunction Φ +∞,θ,5,1 for θ = 0, arctan(1/3), π 4 , arctan(3), π 2 , 5π 8 , 3π 4 , 13π 16 , 15π 16 (blue, magenta, red, purple, grey orange, navy, gold, olive respectively).

π 2 ,x π sin α 1 2 +sin θ sin α 1 x π sin α 5 2 . 6 ) 8 )

 2268 consider the function ψ(x, θ) := Φh,θ,5,1 (x, π/2) = cos θ sin α 5We have that ψ(x, θ) = 0 gives tan θ = -Let x c (h) denote a solution of(7.6). For h = 20, we compute numerically that x c (20) ≈ 0.8096522. Defineθ m (h) := arctan -sin( α5xc π ) sin( α1 2 ) sin( α1xc π ) sin( α5 2 ) . (7.7)For h = 20, we compute numerically that θ m (20) ≈ 0.3324691.To determine the critical points on x = 0, consider the function ϕ(x, y, θ) := Φh,θ,5,1 (x, y) = cos θ sin αWe note that ϕ(x, y, θ) = 0 if and only if tan θ = -

Figure 6 :

 6 Figure 6: The graph of h → g(h) := α5 sin( α 1

  θ = 0 (blue), θ = π 2 -θt (orange), θ = 1 2 ( π 2 -θt + θm) (maroon).

  θ = θm (magenta), θ = π 4 (red), θ = π 2 -θm (purple).

2 -

 2 θm + θt) (deep sky blue), θ = θt (teal), θ = π 2 (grey).

5 (d) 5π 8 (

 58 lime), θ = 3π 4 (navy), θ = 13π 16 (gold).

Figure 7 :

 7 Figure 7: The Robin eigenfunction Φ h,θ,5,1 for h = 20 and various values of θ.

Figure 8 :

 8 Figure 8: The fifth Robin eigenfunction Φ h,θ,0,2 for h = 20 and various values of θ. The values of θ areθ = 0 (blue), θ = π 2 -θ t (orange), θ = 1 2 ( π 2 -θ t + θ m ) (maroon), θ = θ m (magenta), θ = π 4 (red), θ = π 2 -θ m (purple), θ = 1 2 ( π 2θ m + θ t ) (deep sky blue), θ = θ t (teal), θ = π 2(grey), 5π 8 (lime), θ = 3π 4 (navy), θ = 13π 16 (gold).

Theorem B. 1 . 2 .

 12 Let Ω be an open set in R 2 with C 2,+ boundary. Let h ∈ [0, +∞) and let u be a real-valued eigenfunction of the Laplacian with h-Robin boundary conditions. Then u ∈ C 2 (Ω). Furthermore, u has the following properties:1. If u and ∇u vanish at a point x 0 ∈ Ω then there exists > 1, > 0 and a real-valued, non-zero, harmonic, homogeneous polynomial of degree such that:u(x) = p (x -x 0 ) + O(|x -x 0 | + ). (B.1) If u vanishes at x 0 ∈ ∂Ω, then (B.1) holds for some > 0 and u(x) = ar cos ω + O(r + ) (B.2)

2 .

 2 sup x∈B(0,r) |∇ψ(x)| ≤ C p r -2+d+ + 2 p , for p >For p > 2 sufficiently close to 2 (for example -1+ 2 p = 2 ), we get sup x∈B(0,r) |∇ψ(x)| ≤ C p r -1+d+ 2 , for p > 2 .

1 .. 3 ) 2 .

 132 If u and ∇u vanish at a point x 0 ∈ Ω then there exists > 1, > 0 and a real-valued, non-zero, harmonic, homogeneous polynomial of degree such that:u(x) = p (x -x 0 ) + O(|x -x 0 | + ). (BIf moreover x 0 ∈ ∂Ω, then u(x) = ar sin ω + O(r + ) (B.4)

  ,[START_REF] Bucur | Faber-Krahn inequalities for the Robin-Laplacian: a free discontinuity approach[END_REF][START_REF] Adams | Sobolev Spaces[END_REF].

	140 λ h
	130
	120
	90 100 110					(7, 6) (8, 5) (9, 3) (9, 4) (7, 7) (10, 0) (8, 6) (10, 1) (10, 2) (9, 5) (8, 7) (11, 0)
	0	5	10	15	20	25	h

Figure 9: The curves (α m (h) 2 + α n (h) 2 )/π 2 for h ≤ 29 corresponding to the pairs (7, 6), (

  and[START_REF] Bucur | Faber-Krahn inequalities for the Robin-Laplacian: a free discontinuity approach[END_REF][START_REF] Adams | Sobolev Spaces[END_REF] give rise to the same Dirichlet eigenvalue at h = +∞. (see Appendix A).Figure10:The Robin eigenvalues of the square (α m (h) 2 +α n (h) 2 )/π 2 for h ≤ 10 corresponding to the pairs[START_REF] Bers | Local behavior of solutions of general linear elliptic equations[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF], (

	λ h
	110 120 130	(7, 6) (8, 5) (9, 3) (9, 4) (7, 7) (10, 0) (8, 6) (10, 1) (10, 2) (9, 5) (8, 7) (11, 0)
	100
	90
	0	2	4	6	8	10	h

There is a small gap in the proof which can be repaired using Theorem 3.1 due to Sturm.

Note that there is a misprint in[START_REF] Grebenkov | Geometrical structure of Laplacian eigenfunctions[END_REF] after formula (3.9) for the Robin eigenvalue which is corrected here.

We thank T. Hoffmann-Ostenhof for the useful suggestion to establish and use Lemma 5.3. We also thank D. Bucur for his enlightening explanation of the results of[START_REF] Bucur | A variational approach to the isoperimetric inequality for the Robin eigenvalue problem[END_REF] and[START_REF] Bucur | Faber-Krahn inequalities for the Robin-Laplacian: a free discontinuity approach[END_REF].

This means C 2,α for some α > 0.

h = +∞ (see, for example, Figure 5).

We remark that in general, for h = h * 25 any eigenfunction corresponding to λ 25,h is a linear combination of u 4,3 (x, y), u 3,4 (x, y), u 5,1 (x, y) and u 1,5 (x, y). Such an eigenfunction might not have any common symmetries. We note that we have not shown that λ 25,h * 25 (S) is not Courant-sharp.

Multiple crossings: analysis of examples.

Although Proposition 7.1 asserts that the curves corresponding to two distinct pairs can cross at most once, it is possible that an eigenvalue λ n,h (S) is given by more than two distinct curves as h varies.

The situation for the eigenvalues λ 84,h , . . . , λ 92,h seems to be quite complicated. We claim that these eigenvalues are given by the curves corresponding to the pairs (9, 4), [START_REF] Bers | Local behavior of solutions of general linear elliptic equations[END_REF][START_REF] Bers | Local behavior of solutions of general linear elliptic equations[END_REF], (10, 0), [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF], [START_REF] Bucur | Faber-Krahn inequalities for the Robin-Laplacian: a free discontinuity approach[END_REF][START_REF] Adams | Sobolev Spaces[END_REF]. We first show that none of the curves corresponding to other pairs intersect these ones.

By considering Appendix A and monotonicity of the Robin eigenvalues with respect to h, we have that all curves corresponding to pairs (p, q) with π -2 (α p (h) 2 + α q (h) 2 ) ≥ 130 for all h ≥ 0 do not intersect the curves corresponding to the pairs (9, 4), [START_REF] Bers | Local behavior of solutions of general linear elliptic equations[END_REF][START_REF] Bers | Local behavior of solutions of general linear elliptic equations[END_REF], (10, 0), [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF], [START_REF] Bucur | Faber-Krahn inequalities for the Robin-Laplacian: a free discontinuity approach[END_REF][START_REF] Adams | Sobolev Spaces[END_REF]. That is [START_REF] Cheng | Eigenfunctions and nodal sets[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications. Domaines nodaux des fonctions propres[END_REF], [START_REF] Bucur | A variational approach to the isoperimetric inequality for the Robin eigenvalue problem[END_REF][START_REF] Bers | Local behavior of solutions of general linear elliptic equations[END_REF] and so on. From above, we must consider the curves corresponding to (10, 2), [START_REF] Bucur | A variational approach to the isoperimetric inequality for the Robin eigenvalue problem[END_REF][START_REF] Bérard | Sturm's theorem on zeros of linear combinations of eigenfunctions[END_REF], [START_REF] Bucur | Faber-Krahn inequalities for the Robin-Laplacian: a free discontinuity approach[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications. Domaines nodaux des fonctions propres[END_REF], [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bers | Local behavior of solutions of general linear elliptic equations[END_REF], [START_REF] Bucur | Faber-Krahn inequalities for the Robin-Laplacian: a free discontinuity approach[END_REF][START_REF] Bérard | Dirichlet eigenfunctions of the square membrane: Courant's property, and A. Stern's and A. Pleijel's analyses[END_REF], [START_REF] Bucur | A variational approach to the isoperimetric inequality for the Robin eigenvalue problem[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF], [START_REF] Cheng | Eigenfunctions and nodal sets[END_REF]0), [START_REF] Cheng | Eigenfunctions and nodal sets[END_REF][START_REF] Adams | Sobolev Spaces[END_REF], [START_REF] Cheng | Eigenfunctions and nodal sets[END_REF][START_REF] Antunes | Asymptotic behaviour and numerical approximation of optimal eigenvalues of the Robin Laplacian[END_REF], [START_REF] Bucur | Faber-Krahn inequalities for the Robin-Laplacian: a free discontinuity approach[END_REF][START_REF] Bérard | Sturm's theorem on zeros of linear combinations of eigenfunctions[END_REF], [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF].

Note that if q ≤ r, then α p (h) 2 + α q (h) 2 ≤ α p (h) 2 + α r (h) 2 for all h ≥ 0. So it suffices to show that the curves corresponding to [START_REF] Bucur | A variational approach to the isoperimetric inequality for the Robin eigenvalue problem[END_REF][START_REF] Bérard | Sturm's theorem on zeros of linear combinations of eigenfunctions[END_REF], [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bers | Local behavior of solutions of general linear elliptic equations[END_REF], [START_REF] Cheng | Eigenfunctions and nodal sets[END_REF]0) do not intersect those corresponding to [START_REF] Bucur | A variational approach to the isoperimetric inequality for the Robin eigenvalue problem[END_REF][START_REF] Bérard | Dirichlet eigenfunctions of the square membrane: Courant's property, and A. Stern's and A. Pleijel's analyses[END_REF], [START_REF] Bers | Local behavior of solutions of general linear elliptic equations[END_REF][START_REF] Bers | Local behavior of solutions of general linear elliptic equations[END_REF], (10, 0), [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF], [START_REF] Bucur | Faber-Krahn inequalities for the Robin-Laplacian: a free discontinuity approach[END_REF][START_REF] Adams | Sobolev Spaces[END_REF]. Numerically, we compute that λ 9,5,h (S) = λ 5,9,h (S) ≥ 130 for h > 26.9531, λ 8,7,h (S) = λ 7,8,h (S) ≥ 130 for h > 9.3456, λ 11,0,h (S) = λ 0,11,h (S) ≥ 130 for h > 7.3264.

So we need to consider h ≤ 27. We note that the curves corresponding to [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF] and [START_REF] Bucur | Faber-Krahn inequalities for the Robin-Laplacian: a free discontinuity approach[END_REF][START_REF] Antunes | Asymptotic behaviour and numerical approximation of optimal eigenvalues of the Robin Laplacian[END_REF] give rise to the same value at h = +∞.

Again by Appendix A and monotonicity of the Robin eigenvalues with respect to h, we have that all pairs (p, q) with π -2 (α p (h) 2 + α q (h) 2 ) ≤ 97 for all h ≥ 0 do not intersect the curves corresponding to the pairs (9, 4), (7, 7), (10, 0), [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF], [START_REF] Bucur | Faber-Krahn inequalities for the Robin-Laplacian: a free discontinuity approach[END_REF][START_REF] Adams | Sobolev Spaces[END_REF]. That is [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications. Domaines nodaux des fonctions propres[END_REF], [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Antunes | Asymptotic behaviour and numerical approximation of optimal eigenvalues of the Robin Laplacian[END_REF] and below. Hence, from below, we must consider the curves corresponding to (6, 6), [START_REF] Bers | Local behavior of solutions of general linear elliptic equations[END_REF][START_REF] Bérard | Sturm's theorem on zeros of linear combinations of eigenfunctions[END_REF], [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bérard | Dirichlet eigenfunctions of the square membrane: Courant's property, and A. Stern's and A. Pleijel's analyses[END_REF], (9, 0), (9, 1), (9, 2), [START_REF] Bers | Local behavior of solutions of general linear elliptic equations[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF], [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bérard | Sturm's theorem on zeros of linear combinations of eigenfunctions[END_REF], [START_REF] Bucur | A variational approach to the isoperimetric inequality for the Robin eigenvalue problem[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications. Domaines nodaux des fonctions propres[END_REF]. Similarly to the above, it suffices to consider [START_REF] Bers | Local behavior of solutions of general linear elliptic equations[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF], [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bérard | Sturm's theorem on zeros of linear combinations of eigenfunctions[END_REF], [START_REF] Bucur | A variational approach to the isoperimetric inequality for the Robin eigenvalue problem[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications. Domaines nodaux des fonctions propres[END_REF]. We compute numerically that λ 9,4,h (S) = λ 4,9,h (S) ≥ 117 for h > 17.5353, λ 7,7,h (S) ≥ 117 for h > 12.4168, λ 10,0,h (S) = λ 0,10,h (S) ≥ 117 for h > 28.8245, λ 8,6,h (S) = λ 6,8,h (S) ≥ 117 for h > 9.9784, λ 10,1,h (S) = λ 1,10,h (S) ≥ 117 for h > 16.9735. So we must consider h ≤ 29.

In Figure 9, we plot the curves corresponding to the pairs [START_REF] Bers | Local behavior of solutions of general linear elliptic equations[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF], [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bérard | Sturm's theorem on zeros of linear combinations of eigenfunctions[END_REF], (9, 3), (9, 4), (7, 7), (10, 0), [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF], (10, 1), (10, 2), [START_REF] Bucur | A variational approach to the isoperimetric inequality for the Robin eigenvalue problem[END_REF][START_REF] Bérard | Sturm's theorem on zeros of linear combinations of eigenfunctions[END_REF], [START_REF] Bossel | Membranes élastiquement liées: inhomogènes ou sur une surface: une nouvelle extension du théorème isopérimétrique de Rayleigh-Faber-Krahn[END_REF][START_REF] Bers | Local behavior of solutions of general linear elliptic equations[END_REF], [START_REF] Cheng | Eigenfunctions and nodal sets[END_REF]0) The second step is to use the double manifold as suggested in Donnelly-Feffermann, [START_REF] Donnelly | Nodal sets for eigenfunctions of the Laplacian on surfaces[END_REF][START_REF] Donnelly | Nodal sets of eigenfunctions: Riemannian manifolds with boundary, Moser volume, Analysis et Cetera[END_REF][START_REF] Donnelly | Nodal sets of eigenfunctions on Riemannian manifolds[END_REF]. As we only wish to prove a local result, by a diffeomorphism we can reduce to the case when the boundary is given by x 1 = 0. In these new coordinates, the operator reads

In addition, this diffeomorphism can be chosen as a conformal map (see [START_REF] Donnelly | Nodal sets for eigenfunctions of the Laplacian on surfaces[END_REF]), so more precisely, we have

Note that in the Neumann case, there are no linear terms. This would make the proof easier and would permit weaker assumptions. If u denotes the eigenfunction defined locally in x 1 > 0, we define ũ by ũ(x 1 , x 2 ) = u(x 1 , x 2 ) for x 1 > 0 , u(-x 1 , x 2 ) for x 1 < 0 .

We can then define the extension of the operator as H

where ρ, ã2 and c are the extensions of ρ, a 2 and c by reflection and ã1 is defined by odd reflection. So ρ, ã2 and c are Lipschitz and ã1 is only bounded. With this definition, we verify that ũ is an even function (with respect to x 1 ) that satisfies the Neumann condition, and a solution of

The other second derivatives match on x 1 = 0. Hence ũ is actually in C 2 (R 2 ). B.2.3 Nodal structure for solutions of a second-order elliptic operator with coefficients with less regularity.

The third step is to determine whether the local nodal structure that holds for the Laplacian still holds for this second-order elliptic operator which has coefficients with less regularity. This problem is analysed by Hardt-Simon in [START_REF] Hardt | Nodal sets for solutions of elliptic equations[END_REF] (at least in a weaker sense) and more precisely in [START_REF] Qi | Singular sets of solutions to elliptic equations[END_REF] (see Theorem 1.5 and Theorem 3.1). The following theorem is Theorem 3.1 of [START_REF] Qi | Singular sets of solutions to elliptic equations[END_REF] applied to ũ and

in the neighbourhood of a point in the zero set on the boundary, which is assumed to be (0, 0). From this point onwards, we omit the tildes.

Theorem B.2. Suppose that Lu = 0 and that u is not flat at (0, 0), that is, u has finite vanishing order at (0, 0). Then there exists a homogeneous harmonic polynomial P of degree d and, for any p > 1, an > 0 such that ψ := u -P satisfies:

Remark B.3. Note that to apply the theorem, we need to know that u is not flat. According to [START_REF] Qi | Singular sets of solutions to elliptic equations[END_REF] (p. 985, lines 7-9), this is the case under our assumptions and the reference is [START_REF] Garafolo | Monotonicity properties of variational integrals, A p weights, and unique continuation[END_REF].

This theorem gives a good indication of the nodal structure: it should be close to the zero set of the harmonic polynomial P whose structure is well known.

B.2.4 Cheng-Kuo's argument.

Hence the last step is to verify if Cheng's argument [START_REF] Cheng | Eigenfunctions and nodal sets[END_REF] applies (a former reference is [START_REF] Bers | Local behavior of solutions of general linear elliptic equations[END_REF]). We can apply the following lemma attributed by Cheng [START_REF] Cheng | Eigenfunctions and nodal sets[END_REF] to Kuo [START_REF] Kuo | On C ∞ -sufficiency of sets of potential functions[END_REF].

Lemma B.4. Suppose that u and p are smooth functions in R 2 such that, with ψ = u -p, we have for some d ≥ 1 and > 0, In [START_REF] Cheng | Eigenfunctions and nodal sets[END_REF], Cheng applies the lemma to C ∞ functions, but the regularity of u and p is not discussed there. The proof clearly holds for C 2 functions and this assumption is satisfied in our case.

To apply this lemma to the present situation, we observe that a homogeneous harmonic polynomial of degree d in dimension 2 satisfies (iii) and (iv) above. We note that (i) holds by Theorem B.2.