Philippe Charron

Bernard Helffer

Pleijel's theorem for Schrödinger operators with radial potentials

In 1956, Å. Pleijel gave his celebrated theorem showing that the inequality in Courant's theorem on the number of nodal domains is strict for large eigenvalues of the Laplacian. This was a consequence of a stronger result giving an asymptotic upper bound for the number of nodal domains of the eigenfunction as the eigenvalue tends to +∞. A similar question occurs naturally for the case of the Schrödinger operator. The first significant result has been obtained recently by the first author for the case of the harmonic oscilllator. The purpose of this paper is to consider more general potentials which are radial. We will analyze either the case when the potential tends to +∞ or the 1 case when the potential tends to zero, the considered eigenfunctions being associated with the eigenvalues below the essential spectrum.

Introduction

The goal of this paper is to extend Pleijel's theorem for the Dirichlet Laplacian H(Ω) = -∆ in a bounded domain Ω to the case of the Schrödinger operator H V = -∆ + V in R d . We are interested in counting the number of nodal domains of an eigenfunction and to relate this number with the labelling of the corresponding eigenvalue. Throughout this paper, for any function f defined over a domain D ⊂ R d , µ(f ) will denote the number of nodal domains of f , namely the number of connected components of D \ f -1 (0). The starting point of the analysis is Courant's Theorem (1923) [START_REF] Courant | Ein allgemeiner Satz zur Theorie der Eigenfunktionen selbstadjungierter Differentialausdrücke[END_REF].

Theorem 1.1 (Courant) If φ n is an eigenfunction associated with the nth eigenvalue λ n of H(Ω) (ordered in non decreasing order and labelled with multiplicity), then µ(φ n ) ≤ n .

(

Pleijel's theorem (1956) [START_REF] Pleijel | Remarks on Courant's nodal theorem[END_REF] says Theorem 1.2 (Pleijel's weak theorem) If the dimension is ≥ 2, there is only a finite number of eigenvalues of the Dirichlet Laplacian in Ω for which we have equality in (1.1).

Let us now give the strong form of Pleijel's theorem.

Theorem 1.3 (Pleijel's strong theorem) Let λ n the non decreasing sequence of eigenvalues associated to the Dirichlet realization of the Laplacian. For any d ≥ 2 for any orthonormal basis φ n of eigenfunctions φ n of H(Ω) (the Dirichlet Laplacian in Ω) associated with λ n (Ω),

lim sup n→+∞ µ(φ n ) n ≤ γ(d) = 2 d-2 d 2 Γ(d/2) 2 (jd 2 -1 ) d , (1.2)
where j ν denotes the first zero of the Bessel function J ν .

The theorem was proved by Pleijel [START_REF] Pleijel | Remarks on Courant's nodal theorem[END_REF] for d = 2 and then extended by Peetre [START_REF] Peetre | A generalization of Courant nodal theorem[END_REF] and Bérard-Meyer [START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF]. We recall from [START_REF] Bérard | Inégalités isopérimétriques et applications[END_REF]Lemma 9] that Pleijel's constant equals

γ(d) = w -1 d ω -1 d λ(B d ) -d/2 < 1 , (1.3)
where

• w d is the Weyl constant

w d := (2π) -d ω d , (1.4)
•

ω d := |B d | , (1.5)
where B d is the unit ball in R d and |D| denotes for an open set D ⊂ R d its volume;

• λ(B d ) is the Dirichlet ground state energy of the Laplacian in B d .

As γ(d) < 1, one recovers as a corollary that the inequality in Courant's theorem is strict for n large. The second point to notice is that the constant is independent of the open set. Complementary properties of γ(d) have been obtained by B. Helffer and M. Persson-Sundqvist [START_REF] Helffer | On nodal domains in Euclidean balls[END_REF]. In particular d → γ(d) is decreasing exponentially to 0. Finally note that this constant is not optimal (see [START_REF] Bourgain | On Pleijel's nodal domain theorem[END_REF], [START_REF] Steinerberger | A geometric uncertainty principle with an application to Pleijel's estimate[END_REF] and the discussion in [START_REF] Helffer | A review on large k minimal spectral k-partitions and Pleijel's Theorem[END_REF]).

The original proof of Pleijel's theorem is based on a combination between the Weyl formula [START_REF] Weyl | Über die asymptotische Verteilung der Eigenwerte[END_REF] and the Faber-Krahn inequality. Weyl's theorem reads, as λ → +∞,

N (λ) = |Ω| w d λ d 2 (1 + o(1)) , (1.6)
where N (λ) := #{λ j < λ} .

(1.7)

The Faber-Krahn inequality [START_REF] Faber | dass unter allen homogenen Membrane von gleicher Fläche und gleicher Spannung die kreisförmige die tiefsten Grundton gibt[END_REF][START_REF] Krahn | Über eine von Rayleigh formulierte minimal Eigenschaft des Kreises[END_REF] reads (

1.8)

There are a lot of Weyl's formulas available in the context of the Schrödinger operator H V := -∆ + V . The use of the Faber-Krahn inequality is more problematic, except of course for the case of bounded domains with bounded potential which can be treated like the membrane case. In 1989 Leydold [START_REF] Leydold | Knotenlinien und Knotengebiete von Eigenfunktionen[END_REF] obtained in his diploma thesis a weak Pleijel theorem for the isotropic harmonic oscillator (see also [START_REF] Bérard | On the nodal patterns of the 2D isotropic quantum harmonic oscillator[END_REF]). Two years ago Charron [START_REF] Charron | Théorème de Pleijel pour l'oscillateur harmonique quantique[END_REF] in his master thesis proved Pleijel's strong theorem also for the harmonic oscillator: Theorem 1.5 (Charron's theorem) Let (φ n ) n∈N be an orthogonal basis of eigenfunctions in L 2 (R d ) of the harmonic oscillator associated with (λ n ) n∈N . Then

lim sup n→+∞ µ(φ n ) n ≤ γ(d) .
(1.9)

The theorem is also proven in the case of the non-isotropic harmonic oscillator [START_REF] Charron | A Pleijel-type theorem for the quantum harmonic oscillator[END_REF]. The interesting fact is that the potential V (x) = i a i x 2 i (with a i > 0) does not appear on the right hand side of the upper bound. Note also that when there are no eigenvalue degeneracies a much stronger result is available in [START_REF] Charron | A Pleijel-type theorem for the quantum harmonic oscillator[END_REF].

A natural question to ask is whether the theorem can be extended to more general Schrödinger operators. We will answer positively this question under the additional assumption that the potential is radial. More precisely, we assume d ≥ 2

and we consider on R d a Schrödinger operator

H V = -∆ + V , where V is a radial potential: V (x) = v(r) , (1.10)
with |x| = r. We will assume that v ∈ C 1 (0, +∞) , (1.11) and that there exists R 0 > 0 such that v (r) > 0 , for r ≥ R 0 .

(1.12)

In order to allow some singularity at the origin, we assume either

v ∈ C 0 ([0, +∞)) , (1.13)
or that there exists s ∈ (0, 2) such that, as r → 0 ,

v(r) ≈ -r -s . (1.14)
Here, we say that a ≈ b is a/b and b/a are bounded, i.e. if there exists

C > 0 such that 1 C ≤ a b ≤ C .
We will study two cases according to the behavior of v at +∞ . Case A: v tends to +∞ as r → +∞ .

More precisely, we assume (1.11), (1.12) and either (1.13) or (1.14) and that there exists

1 m > 1 such that as r → +∞ v(r) ≈ r m , (1.15)
and v (r) ≈ r m-1 .

(1.16)

Case B: v tends to 0 as r → +∞ .

More precisely, we assume (1.11), (1.12) and either (1.13) or (1.14) and that there exists m ∈ (-2, 0) such that

v(r) ≈ -r m , (1.17)
and v (r) ≈ r m-1 .

(1.18)

In the two cases there is a natural selfadjoint extension starting from C ∞ 0 (R d ) (see Section 3). In Case A, the spectrum is discrete and consists of a non decreasing sequence of eigenvalues λ n tending to +∞ . In Case B the spectrum is divided in two parts, the essential spectrum: [0, +∞) and the discrete spectrum, which consists of an infinite sequence of negative eigenvalues (λ n ) n∈N tending to 0 (see for example Reed-Simon [START_REF] Reed | Method of Modern Mathematical Physics IV: Analysis of operators[END_REF], Vol. IV, Theorem XIII.6). Associated with this sequence (λ n ) n∈N , we can consider an orthonormal sequence of eigenfunctions φ n , where in Case A φ n is an Hilbertian basis of L 2 (R d ) and in Case B of the negative eigenspace.

Our analysis will contain two well-known potentials: the quantum harmonic oscillator (Case A) and the Coulomb potential (Case B). In both cases, we know the eigenvalues and an explicit basis of eigenfunctions but in the proof this property will not be used. Our aim is to prove the following result: Theorem 1.6 (Pleijel's theorem for Schrödinger) In Cases A or B, if (φ n ) n≥1 is an orthogonal sequence of eigenfunctions of H V associated with the above defined sequence λ n , then

lim sup n→+∞ µ(φ n ) n ≤ γ(d) . (1.19)
The paper is organized as follows. In Section 2, we discuss the general strategy and the methods used by Pleijel first and then by P. Charron. In Section 3, we review the general properties of the Schrödinger operator. In Section 4 we collect those Weyl-type results we need for the proof of Theorem 1.6.

In Section 5, we give the proof of our Pleijel's theorem in the two situations.

As recalled in the introduction, the original proof of Å. Pleijel was based on a tricky combination of Weyl's formula with the Faber-Krahn inequality. When considering the case of the Schrödinger operator in R d , Weyl's formula still exists but the use of Faber-Krahn is not easy: nodal domains could be unbounded and the variation of the potential inside a nodal domain could be very high. One has consequently to find an idea for proving that these two bad situations do not occur very often.

In the case of the harmonic oscillator Charron's proof relies on specific properties of the eigenfunctions and the potential. Namely, it used the fact that every eigenfunction is a linear combination of an exponential multiplied by polynomials whose degree can be controlled by a function of the labeling of the eigenvalue, that the hypersurfaces with constant potential are hyperspheres and the fact that the counting function N (λ) behaves nicely as λ → +∞ (Weyl's law).

In addition, it also used that, every nodal domain of an eigenfunction of a Schrödinger operator intersects the classically allowed region associated with the eigenvalue λ, i.e

V (-1) (-∞, λ) := {x ∈ R d | V (x) < λ} . (2.1)
This property is quite general and elementary.

The key was then to divide the classically allowed region in a finite number of annuli of the form V (-1) (a, b). Every nodal domain can either be contained in a single annulus or intersect more than one. To give an upper bound on the number of nodal domains not contained in one annulus, Charron uses properties of algebraic surfaces, as well as results arising in Morse theory adapted from Milnor [START_REF] Milnor | On the Betti Numbers of Real Varieties[END_REF].

Then, he used the Faber-Krahn inequality to give a lower bound on the volume of any nodal domain contained in a single annulus. Dividing the volume of each annulus by the volume of each nodal domain gave an upper bound on the number of nodal domains contained in that annulus.

The last step was to find an appropriate number of annuli to balance out both estimates.

To extend the methods of Charron's proof to more general potentials, we need to find Schrödinger operators such that:

(i) There are good lower bounds for the number of eigenvalues below any λ.

(ii) We can count the number of nodal domains that intersect a given energy hypersurface V (-1) (b).

(iii) We can give an upper bound on the number of nodal domains that are not contained in the classically allowed region.

In the case of the harmonic oscillator, the eigenvalues are known explicitly. However, for many potentials V , Weyl's law can be extended to the Schrödinger operator H V for estimating the number of eigenvalues. Hence, we need to find a suitable class of potentials where this law holds. So far, the only known method to give a suitable upper bound on the number of nodal domains that intersect an energy hypersurface are based on Milnor's theorem (see Subsection 3.7). Hence we need this hypersurface to be algebraic and the property that for any eigenvalue λ and any energy hypersurface (or at least a suitable family) the restriction of any associated eigenfunction u λ equals the restriction of a polynomial to this hypersurface. This is why we focus in this paper on the study of radial potentials. In this case, the energy hypersurfaces are hyperspheres {r = ρ} for some ρ > 0 and it can be shown that the restriction of an eigenfunction to a hypersphere is always a linear combination of hyperspherical harmonics, each one being the restriction to the hypersphere of a homogeneous harmonic polynomial. We will also have to control the degree of such polynomials by a function of λ or of its labelling. This last property will allow us to bound the number of nodal domains that are not contained in V (-1) (-∞, λ).

Another problem might arise when estimating the number of nodal domains contained in one annulus. In the case of the harmonic oscillator, summing over all annuli gives us an expression which can be compared directly with an integral. The error term that arises becomes negligible as λ → +∞ . It remains to show under which conditions on V the same method can be applied.

Finally, in the specific case of Coulomb-like potentials at the origin, we need to look at the behavior of the number of nodal domains near the origin.

3 On the spectral theory of the Schrödinger operators with radial potential

General theory

We first verify that our Schrödinger operator H V = -∆ + V is well defined by a Friedrichs procedure starting from its sesquilinear form defined on

C ∞ 0 (R d \ {0}) × C ∞ 0 (R d \ {0}) (u, v) → a(u, v) := R d ∇u(x) • ∇v(x) dx + R d V (x) u(x) v(x) dx .
Note that the left term has a meaning as soon as V ∈ L 1 loc (R d \ {0}) . In our case, this is a consequence of Assumption 1.11. Our operator, will be defined through a Friedrichs extension. This works as soon as q(u) := a(u, u) is bounded from below by -C ||u|| 2 L 2 . It is consequently enough to control the integral V <0 V (x) u(x) 2 dx from below. When d ≥ 3, we use the standard Hardy inequality (see [START_REF] Cowan | Optimal Hardy inequalities for general elliptic operators with improvements[END_REF] and references therein)

R d |∇u(x)| 2 dx ≥ d -2 2 2 R d 1 r 2 |u(x)| 2 dx , ∀u ∈ C ∞ 0 (R d ) , (3.1)
with r = |x| . This inequality extends to H 1 (R d ) by density. For d = 2, we can use the modified Hardy inequality (see also [START_REF] Cowan | Optimal Hardy inequalities for general elliptic operators with improvements[END_REF]) in a disk D(0, Ř) which reads

R 2 |∇u(x)| 2 dx ≥ 1 4 R 2 u(x) 2 |x| 2 log 2 (|x|/ Ř) dx , ∀u ∈ C ∞ 0 (D(0, Ř)) , (3.2)
which also extends to H 1 0 (D(0, Ř)) . Using these inequalities and a partition of unity, the semi-boundedness on

C ∞ 0 (R d \ {0}) follows immediately.
Let us now describe the form domain resulting from the Friedrichs extension procedure. We have:

Case A Q H = {u ∈ H 1 (R d ) | √ V u ∈ L 2 (D(0, R 1 ) c )} where R 1 is chosen such that v(r) ≥ 1 for r ∈ (R 1 , +∞) , Case B Q H = H 1 (R d ) .
We do not need to characterize the domain of the corresponding self-adjoint operator.

Nodal domains intersect the classically allowed region

We use a similar argument as in [START_REF] Leydold | Knotenlinien und Knotengebiete von Eigenfunktionen[END_REF] and [START_REF] Charron | A Pleijel-type theorem for the quantum harmonic oscillator[END_REF] . We assume that we are either in Case A or in Case B, but the result is much more general.

Proposition 3.1 Let λ be an eigenvalue below the essential spectrum, u λ be an eigenfunction of H V associated with eigenvalue λ and Ω be a nodal domain of u λ . Then

Ω ∩ V (-1) (-∞, λ) = ∅ . Proof. If for all x ∈ Ω, V (x) > λ , then λ = Ω |∇u λ (x)| 2 dx + Ω V (x)u λ (x) 2 dx Ω u λ (x) 2 dx ≥ Ω V (x)u λ (x) 2 dx Ω u λ (x) 2 dx > Ω λu λ (x) 2 dx Ω u λ (x) 2 dx = λ , (3.3)
hence a contradiction.

2

Therefore, any nodal domain is either contained in the classically allowed region {V < λ} or intersects the hypersurface V (-1) (λ) .

The radial Schrödinger operator

Although the exposition there is limited to the case d = 3, one can refer to Reed-Simon [START_REF] Reed | Method of Modern Mathematical Physics IV: Analysis of operators[END_REF], Vol. IV p. 90-91. The Laplace operator -∆ can be written as

-∆ = - ∂ 2 ∂r 2 - d -1 r ∂ ∂r + 1 r 2 (-∆ S d-1 ), (3.4)
where r = |x| is the radial variable and ∆ S d-1 is the Laplace-Beltrami operator, acting in L 2 (S d-1 ). The following proposition is standard (see for example [START_REF] Shubin | Asymptotic behaviour of the spectral function[END_REF], Theorem 22.1 and Corollary 22.1).

Proposition 3.2 Assume that d ≥ 2. The spectrum of -∆ S d-1 consists of eigenvalues ( + d -2) , ∈ N .
The multiplicity of the eigenvalue ( + d -2) is given by

Λ ,d := + d -1 d -1 - + d -3 d -1 ,
which coincides with the dimension of the space of homogeneous, harmonic polynomials of degree .

We denote by S d-1 ω → Y ,m (ω) an orthonormal basis of the Λ ,d -dimensional eigenspace associated with ( + d -2). We recall that each Y ,m is the restriction to S d-1 of a harmonic homogeneous polynomial of degree . We now consider the Schrödinger operator H V and assume

V (x) = v(r) . (3.5)
In this case, one can determine the spectrum by using polar coordinates. In the spherical coordinates, we can determine the spectrum by considering the (closure of the) union of the spectra of the family (indexed by ∈ N) of Sturm Liouville operators L defined by

L = - d 2 dr 2 - d -1 r d dr + ( + d -2) r 2 + v(r) , (3.6)
acting in L 2 ((0, +∞), r d-1 dr) , with a suitable Dirichlet like condition at 0 (see Reed-Simon [START_REF] Reed | Method of Modern Mathematical Physics IV: Analysis of operators[END_REF], p. 91, Proof of Lemma 1). Note that the "Dirichlet like" condition is expressed after the unitary transform u → r d-1

2 u sending L 2 ((0, +∞); r d-1 dr) onto L 2 ((0, +∞); dr) and becomes the standard Dirichlet condition for = 0. When > 0, no condition is given. The new operator is then:

L = - d 2 dr 2 + ( + d-1 2 )( + d-1 2 -1) r 2 + v(r) . Proposition 3.3 Let H V = -∆ + V , where V (x) = v(r) satisfies either Case A or Case B . Any eigenvalue of -∆ + V is of the form λ = λ n, , (3.7)
where λ n, is the n-th eigenvalue of L .

A corresponding basis of eigenfunctions has the form

u n, ,m (r, ω) = f n, (r)Y ,m (ω) , (3.8)
where Y ,m (ω) denotes an orthonormal basis of (hyper)spherical harmonics.

We recall that these functions form a basis of L 2 (R d ) in Case A (see [START_REF] Reed | Method of Modern Mathematical Physics IV: Analysis of operators[END_REF]) or a basis of the negative eigenspace in Case B.

Courant's nodal theorem and nodal behavior of eigenfunctions.

For the analysis of potentials with singularities it is worth to ask under which condition one can prove Courant's theorem or describe the local nodal structure of an eigenfunction. Under our assumptions the only point is the control at the origin. Outside the origin, (1.11) implies that the potential is C 1 and the structure of the nodal set is well known.

Looking at the proof of Courant's theorem, the only thing we need is the unique continuation theorem, i.e. we need to show that if an eigenfunction u λ is identically 0 in a non empty open set ω then it is zero in R d . The argument clearly works if there is only a singularity at 0, because ω \ {0} is an open set where u λ vanishes identically. See [START_REF] Hoffmann-Ostenhof | Local properties of solutions of Schrödinger equations[END_REF] for more properties. We will show in the next subsection that no nodal domain is contained in a sufficiently small ball around the origin. Hence the counting of nodal domains can start outside this ball.

No nodal domains in a small ball

In this subsection, we show that under our assumptions the nodal domain cannot be contained in a small neighborhood of the origin. We will start with Case B which is easier.

Case B

We have the following statement:

Proposition 3.4
If d ≥ 2 and in Case B there exists r d > 0 such that, if λ < 0 is an eigenvalue and u λ is a corresponding eigenfunction, there is no nodal domain ω of u λ contained in B(0, r d ) .

Proof

We first deduce from Assumption (1.14), because s < 2, that:

• For d ≥ 3 there exists r d such that v(r) + (d -2) 2 4r 2 > 0 , (3.9)
for r ∈ (0, r d ).

• For d = 2, there exists r 2 > 0 such that

v(r) + 1 4r 2 ln 2 (r/2r 2 ) > 0 , (3.10)
for r ∈ (0, r 2 ) .

We now use the identity

ω |∇u λ (x)| 2 dx + ω V (x)|u λ (x)| 2 dx = λ ω |u λ (x)| 2 dx ,
and get because λ < 0

ω |∇u λ (x)| 2 dx + ω V (x)|u λ (x)| 2 dx < 0 .
When d ≥ 3 , we use Hardy's inequality (3.1) (u λ is extended in R d by 0 outside ω and this extension is in H 1 (R d )) and get that

ω V (x) + (d -2) 2 4r 2 |u λ (x)| 2 dx < 0 .
This contradicts (3.9) if ω ⊂ B(0, r d ) .

When d = 2 , we use the modified Hardy inequality (3.2) with Ř = 2r 2 , (u λ is extended by 0 outside ω in D(0, Ř) and this extension is in H 1 0 (D(0, Ř))) and get a contradiction with (3.10) . 2

Case A

In Case A, with singularities, there is some difficulty because we consider λ large. When d ≥ 3 , the previous proposition will be true in a ball whose radius is r d (λ) ≈ λ -1 2 . For d = 2 , r 2 (λ) could be taken as r 2 (λ) ≈ λ -1 2for some > 0 . More precisely, we have Proposition 3.5 Under Assumption (1.14) and if d ≥ 3 , there exists a constant c d > 0 which depends on V and d only and λ 0 > 0 such that, for λ ≥ λ 0 and if u λ denotes an eigenfunction of H V , there are no nodal domains of u λ contained in B(0, c d λ -1/2 ) . If d = 2 , for any > 0 , there exists λ > 0 and c V that depends only on V such that, for λ ≥ λ and if u λ denotes an eigenfunction of H V , there are no nodal domains of u λ contained in B(0, c V λ -1/2-) .

Proof By (1.14), there exists C > 0 and r 0 > 0 , such that V > -Cr -s for 0 < r ≤ r 0 . As s ∈ (0, 2), there exists λ 0 such that for λ ≥ λ 0 , there exists

r d (λ) ∼ d-2 2 λ -1 2 such that (d -2) 2 4 r -2 -Cr -s > λ , ∀r ∈ (0, r d (λ)) .
This implies

(d -2) 2 4 r -2 + v(r) -λ > 0 , ∀r ∈ (0, r d (λ)) and λ ≥ λ 0 .
The proof is achieved by taking 0 < c < d-2 2 in the statement of the proposiiton and using the Hardy inequality as in the second part of the proof of Proposition 3.4. For the case d = 2 , a not optimal r 2 (λ) = λ -1 2 -(for some > 0) together with the modified Hardy inequality do the job for λ large enough.

Upper bound for the degree of the polynomials associated with the spherical harmonics

In this subsection, we prove the existence of a rather optimal upper bound on the degree of the polynomials associated with the spherical harmonics Y m appearing in the decomposition of an eigenfunction u λ .

Proposition 3.6

In Cases A or B, if λ is an eigenvalue of H V such that λ < lim inf r→+∞ v then there exists p λ such that, for any associated eigenfunction u λ and for any τ satisfying inf v < τ ≤ λ, we can find a polynomial P τ,λ of d variables of degree at most p λ such that on V (-1) (τ ) in R d the restriction of u λ is equal to the restriction of P τ,λ . Moreover, p λ satisfies

p λ ≤ max{ | ≥ 1 and inf r v(r) + ( + d -2) r 2 < λ } . (3.11)
Proof.

For given λ, u λ has the form (see (3.8))

u λ = λ=λ n, c n, ,m u n, ,m .
If we restrict u λ to the hypersphere of radius r τ = v (-1) (τ ), we get

u λ (r τ , ω) = λ=λ n, d λ,τ, ,m Y ,m (ω) .
Considering the property of Y ,m , we can choose for the proof of the proposition

P τ,λ = d λ,τ, ,m P ,m , (3.12)
where P ,m is the homogeneous harmonic polynomial of degree such that

(P ,m ) /{r=1} = Y ,m .
It remains, in order to prove (3.11), to determine the highest ≥ 1 such that λ n, = λ. By the minimax principle, we have

λ = λ n, ≥ m := inf v(r) + ( + d -2) r 2 .
We have indeed

λ n, = +∞ 0 |f n, (r)| 2 r d-1 dr + +∞ 0 v(r) + ( + d -2) r 2 |f n, (r)| 2 r d-1 dr . 2
The behavior of m should be analyzed but note that our assumptions imply that m > -∞ . Furthermore → m is strictly increasing, so we can set p λ := [p λ ] + 1 , where pλ is the solution of λ = m pλ and [x] means the integer part of x.

Application: Determination of an upper bound of p λ .

Case A We can assume that ≥ 1. This simply implies later a choice of p λ ≥ 1 If we consider v(r) = c r m as a model case for m > 1 and c > 0 , the infimum is obtained when

c m r m-1 -2 ( + d -2) r 3 = 0 , i.e. r = 2 ( + d -2) c m 1 m+2 . So we get inf r m + ( +d-2) r 2 = 2 ( +d-2) c m m m+2 + ( + d -2) 2 ( +d-2) c m -2 m+2 = 2 c m m m+2 c m+2 2 ( ( + d -2)) m m+2 .

This gives us pλ

∼ a m λ m+2 2m , (3.13)
with

a m = ( 2 c m ) -1 2 ( c m + 2 2 ) -2m m+2 .
For m = 2 , we recover what we got for the harmonic oscillator by direct computation.

To treat the general case, we use the lower bound:

v(r) ≥ c r m - C r s . (3.14)
We have to estimate

inf c r m - C r s + ( + d -2) r 2 .
We observe that:

inf c r m -C r s + ( +d-2) r 2 ≥ inf c r m + ( +d-2) 2r 2 + inf -C r s + ( +d-2) 2r 2 .
But, there exists (see below the computation in (3.17

) with m = -s) a constant C 0 > 0, such that, for all ≥ 1, inf - C r s + ( + d -2) 2r 2 ≥ -C 0 ,
and we can use the lower bound of the model case above to get:

inf r v(r) + ( + d -2) r 2 ≥ 1 2 2 c m m m+2 c m + 2 2 ( ( + d -2)) m m+2 -C 0 .
Hence we obtain like for the model case:

Corollary 3.7 In Case A, as λ → +∞ , p λ ≈ λ m+2 2m . (3.15)
Case B Let us now compute an example corresponding to Case B. If we take v(r) = -r m for m ∈ (-2, 0), the infimum is obtained when

mr m-1 + 2 ( + d -2) r 3 = 0 . i.e. r = 2 ( + d -2) -m 1 m+2
.

(3.16)

So we get inf -r m + ( +d-2)

r 2 = -2 ( +d-2) -m m m+2 + ( + d -2) 2 ( +d-2) -m -2 m+2 = -2 -m m m+2 m+2 2 ( ( + d -2)) m m+2 . (3.17) This gives us for λ → 0 , λ > 0 , pλ ∼ a m (-λ) m+2 2m , (3.18)
with

a m = ( 2 -m ) -1 2 ( m + 2 2 ) -2m m+2 . (3.19)
In the Coulomb case m = -1 and d = 3, we get

pλ ∼ a -1 (-λ) -1 2 , (3.20)
to compare with the direct computation which can be done for the Coulomb case.

In the general case, we can use

v(r) ≥ -C r -s , ∀r ∈ (0, R)
and v(r) ≥ -c r m , ∀r ∈ (R, +∞) .

We will use twice the analysis of the model, the first time with m replaced by -s. Let us first consider inf r∈(0,R)

v(r) + ( + d -2) r 2 ≥ inf r∈(0,R) -C r -s + ( + d -2) r 2 .
We observe (see (3.16) with m = -s) that for large enough the map r → -C r -s + ( +d-2)

r 2 is decreasing on (0, R). Hence inf r∈(0,R) -C r -s + ( + d -2) r 2 = -C R -s + ( + d -2) R 2 ≥ -C s .
For the second case, we can use

inf r∈(R,+∞) -C r m + ( + d -2) r 2 ≥ inf r∈(0,+∞) -C r m + ( + d -2) r 2 ,
and what we obtained for the homogeneous model.

Nodal domains on hyperspheres

Since the considered potentials V are radial, the energy hypersurfaces {V = αλ} are hyperspheres centered at 0. Also, the restriction of any eigenfunction u λ of H V to a hypersphere equals the restriction of some harmonic polynomial. We can use the following result proven in [START_REF] Charron | A Pleijel-type theorem for the quantum harmonic oscillator[END_REF], which is based on [START_REF] Milnor | On the Betti Numbers of Real Varieties[END_REF]: Proposition 3.9 Let P be a polynomial of degree k with d variables. Then its restriction to the hypersphere S d-1 admits at most 2 2d-1 k d-1 nodal domains.

We will combine this with the previous estimates obtained in Corollaries 3.7 and 3.8 to obtain an upper bound on the number of nodal domains on any hypersphere. [START_REF] Bourgain | On Pleijel's nodal domain theorem[END_REF] Weyl's formula

Preliminaries

For Schrödinger operators, Weyl's formula takes (under of course suitable assumptions to be discussed below) the form

N (λ) ∼ (2π) -d ξ 2 +V (x)≤λ dxdξ . (4.1)
After integration in the ξ variable, we get

N (λ) ∼ W (λ) , (4.2)
where

W (λ) := w d (λ -V ) d 2 + dx , (4.3)
with w d defined in (1.4).

This formula makes sense in case A (as λ → +∞) and in case B (as λ → 0 with λ < 0). Let us just compute the right hand side for the two toy models: the harmonic oscillator and the Schrödinger operator with Coulomb potential. For the harmonic oscillator, we simply get

W (λ) = w d (λ -r 2 ) d 2 + dx = h d w d λ d , (4.4)
with

h d := (1 -r 2 ) d 2 + dx > 0 .
More generally, if v(r) = r m for m > 0, we obtain, as λ → +∞,

W (λ) = w d h d,m λ d( 1 2 + 1 m ) . (4.5)
In the Coulomb case, we get, with λ < 0

W (λ) = w d λ + 1 r d 2 + dx = e d w d (-λ) -d 2 , (4.6)
with

e d := 1 |x| -1 d 2 + dx = |S d-1 | 1 0 (1 -r) d 2 r d 2 -1 dr < +∞ . (4.7)
More generally, if v(r) = -r m , for m ∈ (-2, 0), we obtain as λ → 0 (λ < 0),

W (λ) = w d h d,m |λ| d( 1 2 + 1 m ) . (4.8)
Observing that N (λ n ) = n -1 if λ n-1 < λ n , and assuming that the Weyl formula is proven (see below for the proof), we get conversly

λ n ∼ wd n 1 d with 1 = h d w d ( wd ) d , (4.9)
in the case of the harmonic oscillator and

λ n ∼ -v 3 n -2 3 , with 1 = e 3 w 3 (v 3 ) -3 2 , (4.10)
in the case of the Coulomb case.

More generally we have the proposition:

Remark 4.3
The theorem is still true if we consider the Dirichlet problem for H V in R d \B, where B is a ball centered at 0. In Case A, the assumptions of the theorem are satisfied in R d \ B. This follows of our assumption (1.16) on v .

For Case B, Theorem XIII.82 in [START_REF] Reed | Method of Modern Mathematical Physics IV: Analysis of operators[END_REF] reads:

Theorem 4.4 Let V be a measurable function on R d (d ≥ 2) obeying -c 1 (r + 1) -β ≤ V (x) ≤ -c 2 (r + 1) -β , (4.14)
and

|V (x) -V (y)| ≤ c 3 [1 + min{|x|, |y|}] -β-1 |x -y| , (4.15)
for some β < 2 and suitable constants c 1 , c 2 , c 3 > 0 .

Then

lim λ→+∞ N (λ)/W (λ) = 1 . Remark 4.5
The theorem is still true if we consider the Dirichlet problem for H V in R d \B, where B is a ball centered at 0. In Case B, the assumptions of the theorem are satisfied in R d \B. This is a consequence of our assumption on v in (1.18).

Treatment of the singularity

To cover the question of the treatment of the singularity at the origin we could think of using (Problem 132 in [START_REF] Reed | Method of Modern Mathematical Physics IV: Analysis of operators[END_REF]) to treat the singularity as a perturbation. Due to the use of the Cwickel-Lieb-Rozenblum inequality [10, ?, ?] in the argument, this approach works only under the condition d ≥ 3 . If we remember that we only need a lower bound for N (λ) one can proceed for d ≥ 2 in the following way. We can introduce a small ball B = B(0, ) around the singularity and look at the Dirichlet problem in R d \ B . We denote by N B (λ) the corresponding counting function. Because the eigenvalues are greater than the initial problem by monotonicity of the domain, the estimate of the N (λ) of the new problem will give the lower bound:

N B (λ) ≤ N (λ) .
The theorem in Reed-Simon [START_REF] Reed | Method of Modern Mathematical Physics IV: Analysis of operators[END_REF] can be applied in R d \ B (proof unchanged) and we get by Weyl's formula

N B (λ) ∼ W B (λ) , with W B (λ) = w d r≥ (λ -V ) d 2 + dx .
It remains to compare W B (λ) and W (λ) in our two cases.

Case B Here λ < 0 . It is enough to show that, for some > 0,

B(0, ) (λ -V ) d 2 + dx < +∞ .
We have

B(0, ) (λ -V ) d 2 + dx ≤ C 0 r d-1-s d 2 dr < +∞ ,
the finiteness resulting from the assumption s < 2 .

Case A. Here λ ≥ λ 0 > 0. We will show that B(0, ) (λ -V )

d 2
+ dx is relatively small in comparison with N (λ). In Case A, we have seen that

W (λ) ≈ (λ -r m ) d 2 + r d-1 dr ≈ λ d 2 + d m , (4.16)
We have

B(0, ) (λ -V ) d 2 + dx ≤ C 0 (λ + r -s ) d 2 r d-1 dr ≤ Ĉ ( Ĉ + λ d 2 ) .
This gives, as λ → +∞ ,

B(0, ) (λ -V ) d 2 + dx/W (λ) = O(λ -d m ) .
Hence in the two cases, we have shown that W A (λ) ∼ W (λ) . In conclusion, we have obtained the following theorem.

Theorem 4.6

In Cases A or B , if d ≥ 2 , we have

N (λ) ≥ W (λ)(1 + o(1)) , (4.17)
where the remainder o(1) is as λ → +∞ in Case A and as λ → 0 (λ < 0) in Case B .

5 Counting nodal domains

Preliminaries

We construct a radial partition of {V < λ} of cardinality ν(λ) with ν(λ) to be defined later. When v is increasing on (0, +∞), r λ := v -1 (λ) is well defined in (inf v, +∞) in Case A, and for any λ in (-∞, 0) in Case B. But we will only be interested λ → +∞ in Case A, and in λ → 0 (λ < 0) in Case B. Under the weaker Assumption (1.12), we can define in the two cases r λ by

r λ := sup v(r)=λ r ,
and obtain the asymptotics

r λ ≈ |λ| 1 m .
Note that in both cases r λ tends to +∞ . We recall that in this case, we assume λ ≥ λ 0 > 0 .

Analysis of

We introduce a partition of the classical region by introducing ν(λ) annuli, for i = 1, . . . , ν(λ),

D i (λ) := x ∈ R d | i -1 ν(λ) 1/d r λ < r < i ν(λ) 1/d r λ .
We note that each annulus has the same volume: The determination of λ 0 "large enough" will be given during the proof. If u λ denotes some eigenfunction, we denote by D(u λ ) the set of the nodal domains of u λ . We now introduce in D(u λ ) the following subsets. Definition 5.1

|D i (λ)| = ω d 1 ν(λ) r λ d ≈ λ d m ν(λ) . ( 5
A i (u λ ) = {Ω ∈ D(u λ ) | Ω ⊂ D i (λ)} .
Here, i can take the values 1, 2, . . . , ν(λ). Definition 5.2

B j (u λ ) = Ω ∈ D(u λ ) Ω ∩ r = j ν(λ) 1/d r λ = ∅ .
Here, j can take the values 1, 2, . . . , ν(λ) .

From Subsection 3.2, we know that every nodal domain is contained in at least one of these sets.

Remark 5.3 This partition will be refined by introducing, in the case of a singularity at the origin, a further partition of A 1 (u λ ).

Counting the nodal domains contained in one annulus of the partition

We first count in each of the annuli D i (λ) for i ≥ 2 . Except if there are no singularity, the treatment of A 1 (u λ ) will be done separately. Hence we first prove the Proposition 5.4

In Case A, if ν(λ) satisfies (5.2) and (5.3), we have the following inequality, as λ → +∞ ,

ν(λ) i=2 #A i (u λ ) ≤ γ(d) W (λ) 1 + O 1 ν(λ)
.

(5.4)

Proof

If Ω is a bounded nodal domain of u λ , the Faber-Krahn inequality (see Theorem 1.4) gives:

Ω |∇u λ (x)| 2 dx Ω u λ (x) 2 dx ≥ 1 |Ω| 2 d ω d 2 d λ(B d ) , (5.5)
where we recall that |Ω| denotes the volume of Ω .

We know that for a given bounded nodal domain Ω, we have

λ = Ω |∇u λ | 2 dx + Ω V u 2 λ dx Ω u 2 λ dx , which implies Ω |∇u λ | 2 dx Ω u 2 λ dx < λ -inf x∈Ω V (x) . (5.6)
For all Ω ∈ A i (u λ ), we can, under Condition (5.3), use (1.12) to obtain that for λ large enough,

Ω |∇u λ | 2 dx Ω u 2 λ dx < λ -v i -1 ν(λ) 1/d r λ , (5.7)
for i ≥ 2 . We can then combine (5.5) and (5.7) and obtain

|Ω| ≥ ω d λ(B d ) d 2 λ -v i-1 ν(λ) 1/d r λ d 2
.

(5.8)

Observing that

Ω∈A i (u λ ) |Ω| ≤ |D i (λ)| , we obtain that # A i (u λ ) ≤ 1 ω d λ(B d ) d 2   |D i (λ)| λ -v i -1 ν(λ) 1/d r λ d 2   .
(5.9)

Summing up for i = 2, . . . , ν(λ), we get

ν(λ) i=2 # A i (u λ ) ≤ 1 ω d λ(B d ) d 2   ν(λ) i=2 |D i (λ)| λ -v i -1 ν(λ) 1/d r λ d 2   .
(5.10)

We recognize on the right hand side a Riemann sum for the function

x → (λ -V (x)) d 2 in D(0, r λ ) \ (D ν(λ) (λ) ∪ D 1 (λ))
. More precisely, we can write, using the monotonicity of v ,

ν(λ) i=2 |D i (λ)| λ -v i-1 ν(λ) 1/d r λ d 2 ≤ |D 2 (λ)| λ -v 1 ν(λ) 1/d r λ d 2 + (λ -V (x)) d 2 + dx ≤ |D 2 (λ)| λ d 2 + (λ -V (x)) d 2 + dx .
Using the asymptotic behavior of v at +∞ given in (1.17), the computation of γ(d) in (1.3), the asymptotic behavior of W (λ) given in (4.11), and (5.1), we achieve the proof of Proposition 5.4.

Counting the nodal sets meeting the boundary of the annuli

Let us now turn to the study of the sets B i (u λ ). We have shown in Corollary 3.7 that p λ ≈ λ m+2 2m . Using Proposition 3.9, we obtain that the number of nodal domains in a given B i (u λ ) satisfies

#B i (u λ ) ≤ 2 2d-1 p d-1 λ ≤ C d λ (d-1)(m+2) 2m
.

(5.11)

Comparing with (4.11), we get that, for some C > 0 and λ ≥ λ 0 ,

i #B i (u λ ) ≤ C ν(λ) λ -m+2 2m W (λ) . If ν(λ) satisfies in addition, lim λ→+∞ ν(λ) λ -m+2 2m = 0 , (5.12)
we obtain

lim λ→+∞ i #B i (u λ ) W (λ) = 0 . (5.13) 5.2.4 Counting in D 1 (λ)
We still have to consider when there is a singularity around 0. We treat first the case d ≥ 3 . We know from Proposition 3.5 that we can replace D 1 (λ) by the annulus D 1 (λ, C) defined by To continue, we consider a partition of D 1 (λ) in two annuli: Hence, for λ ≥ λ 0 , we deduce from (5.6), that

D 1 (λ, C) := {x ∈ R d | 1 C λ -1 2 < r < r λ /ν(λ)
D 11 (λ) := {x ∈ R d | 1 C λ -1 2 < r < C} and D 12 (λ) := {x ∈ R d | C < r < r λ /ν(λ) 1 d } ,
Ω |∇u λ (x)| 2 dx Ω u λ (x) 2 dx ≤ 2λ .
As in the proof of (5.8) we obtain:

|Ω| ≥ ω d λ(B d ) d 2 2 -d 2 λ -d 2 .
( (5.20)

The case when d = 2 does not lead to new difficulties.

Conclusion for Case A

Summing all the upper bounds and having chosen ν(λ) satisfying (5.2), (5.3), and (5.12), we get, as λ → +∞ , µ(u λ ) ≤ γ(d) W (λ) (1 + o(1)) .

(5.21)

Using the asymptotic upper bound (4.17 Observing that N (λ n ) ≤ n -1 , we obtain Theorem 1.6.

Case B

The proof in case B follows the same lines. We define the sets D i (λ), A i (u λ ) and B i (u λ ) for i = 2, . . . , ν(λ) as in case A, with ν(λ) satisfying conditions equivalent to (5.2), (5.3) and (5.12), replacing λ → ∞ by λ → 0 (with λ < 0). Hence, we assume that ν(λ) satisfies the conditions lim )) .

(5.26)

The proof is the same as in case A. For W (λ), we can use the asymptotics (4.11).

For the cardinalities of the sets B i (u λ ), (5.13) holds in case B under condition (5.25) and we can use Corollary 3.8 together with (4.11).

The treatment of the singularity is slightly easier in this case. We use Proposition 3.4 to make a partition of D 1 (λ) in two annuli: Again, we choose C such that v is strictly increasing for |x| > C. Since λ < 0, there exists M > 0 such that #A 11 (u λ ) < M . To give an upper bound on #A 12 (u λ ), we follow the same steps as in case A.

The behavior of the number of eigenvalues below or equal to λ is given by (4.17), hence inequality (5.22) is verified. Since N (λ n ) ≤ n -1 , we obtain Theorem 1.6 in this case as well.

Theorem 1 . 4 2 d λ(D) ≥ ω 2 d

 1422 For any domain D ⊂ R d (d ≥ 2), we have |D| d λ(B d ) .

Corollary 3 . 8

 38 In Case B,p λ ≈ (-λ)

Case A 5 . 2 . 1 A

 521 first partition of the classical region.

. 1 )

 1 The cardinality ν(λ) satisfies a priori the conditionlim λ→+∞ ν(λ) = +∞ ,(5.2)but the condition limλ→+∞ ν(λ) -1/d r λ = +∞ ,(5.3)will appear along the proof.

  large C. The number of nodal domains µ 10 (λ) crossing the hypersphere {r = 1 C λ -

1 ω d λ(B d ) d 2 |D 12 . 18 ) 2 ,

 1212182 where we keep the liberty to choose C > 0 larger than the previous one. Again the number µ 11 (λ) of nodal domains crossing the hypersphere {r = C} is controlled by(5.11):lim λ→+∞ µ 11 (λ) W (λ) = 0 . (5.15)Control in D 12 (λ). The treatment of D 12 (λ) can be done for C large enough (in order to have the monotonicity of v) like the analysis of the D i (λ) for i ≥ 1. More precisely, we can replace (5.9) (for i = 1) by# A 12 (u λ ) ≤ (λ)| (λ -v(C)) u λ ) := {Ω ∈ D(u λ ) | Ω ⊂ D 12 (λ)} .Hence we get, for some constantC d > 0 , # A 12 (u λ ) ≤ C d ν(λ) 12 (u λ )/W (λ)) = 0 . (5Control in D 11 (λ).Note that in D 11 (λ), we have, for some constant C s > 0 , V (x) ≥ -C s λ s ∀x ∈ D 11 (λ) .

. 19 )

 19 Summing over the Ω's contained in D 11 (λ) and observing that the volume of D 11 (λ) is bounded, we get the existence of a constant Ĉd such that#A 11 (u λ ) ≤ Ĉd λ d 2 ,whereA 11 (u λ ) := {Ω ∈ D(u λ ) | Ω ⊂ D 11 (λ)} .In particular we get lim λ→+∞ #A 11 (u λ ) W (λ) = 0 .

  ), we get, as λ → +∞ ,µ(u λ ) ≤ γ(d) N (λ) (1 + o(1)) .(5.22)

5

 5 In Case B, if ν(λ) satisfies (5.23) and (5.24), then we have the following inequality, as λ → 0 (λ < 0 ),ν(λ) i=2 #A i (u λ ) ≤ γ(d) W (λ) (1 + O( 1 ν(λ)

D

  11 (λ, C) := {x ∈ R d | r d < r < C} , and D 12 (λ, C) := {x ∈ R d | C < r < r λ /ν(λ) 1 d } .

This condition appears when applying Weyl's formula given by Theorem 4.2 from[START_REF] Reed | Method of Modern Mathematical Physics IV: Analysis of operators[END_REF]. At least under stronger assumptions on the regularity of v for r → +∞, it should be possible to assume m > 0.

Acknowledgements.

Thanks to the ESI where the paper was initiated (B. H. and T. H.-O.). The authors would also thank I. Polterovich for helpful discussions at various stages of this work.

where the asymptotics is as λ → +∞ in Case A and as λ → 0 (λ < 0) in Case B.

Proof

Outside a ball we can use for estimating the integral defining W (λ) the comparison of v(r) with r m and then use the previous computations for the models. The control of the integral in a ball will be done in Subsection 4.3. 2

Weyl's formula under weak assumptions

There is vast literature on this subject: Reed-Simon [START_REF] Reed | Method of Modern Mathematical Physics IV: Analysis of operators[END_REF] and references therein (for the historics), D. Robert [START_REF] Robert | Propriétés spectrales d'opérateurs pseudo-différentiels[END_REF], H. Tamura [START_REF] Tamura | Asymptotic formulas with remainder estimates for eigenvalues of Schrödinger operators[END_REF], Tulovski-Shubin [START_REF] Shubin | On the asymptotic distribution of the eigenvalues of pseudodifferential operators in R n[END_REF], R. Beals [START_REF] Beals | A general calculus of pseudo-differential operators[END_REF], L. Hörmander [?, 15], A. Mohamed [START_REF] Mohamed | Comportement asymptotique, avec estimation du reste, des valeurs propres d'une classe d'opérateurs pseudo-différentiels sur R n[END_REF]. In the recent contributions the goal is to control the remainder but this is not important in the applications considered here. Here, we prefer to work under weaker asssumptions and can use Theorem XIII.81 in Reed-Simon (Vol. 4) [START_REF] Reed | Method of Modern Mathematical Physics IV: Analysis of operators[END_REF] for the case V → +∞ with a condition d ≥ 2 and m > 1, and for the case V → 0 as |x| → +∞, Theorem XIII.82. The treatment of the singularity is also explained (without detail) (see the discussion p. 277, lines -7 to -1, sending to Problem 132 therein). The idea there is to first prove a statement with V continuous and then to show that the addition of a potential W with compact support or in L Theorem XIII.81 in [START_REF] Reed | Method of Modern Mathematical Physics IV: Analysis of operators[END_REF] reads:


R Chosson

V Vandenberghe

A.-F Gourgues-Lorenzon

J Crépin

J.-C Brachet

A Cabrera

V Garat

CREEP BEHAVIOR AT HIGH TEMPERATURE OF THE OXYGEN STABILIZED ZIRCONIUM ALPHA PHASE OF FUEL CLADDING TUBES OXIDIZED IN LOCA CONDITIONS

During hypothetical Loss-Of-Coolant-Accident (LOCA) scenarii, zirconium alloy fuel cladding tubes are oxidized at high temperature (HT): zirconia and oxygenstabilized alpha, called alpha(O) are formed on the outer surface of the cladding. To model the behavior of the cladding during LOCA transients, the specific behavior of each layer must be known. This study focuses on the characterization of the viscoplastic behavior of the alpha(O) phase at HT. Model alpha(O) phase specimens enriched by 2 wt% in oxygen were produced from zirconium alloys cladding specimens. Creep tensile tests were performed in the tube axial direction under secondary vacuum between 900 and 1100°C for stresses between 2 and 31 MPa. The creep resistance of the alpha(O) phase is higher than that of the beta phase by three to four decades. The HT creep behavior of oxidized cladding was then modeled by finite element simulation.

I. INTRODUCTION

Understanding interactions between oxidation and mechanical behavior is necessary to predict the mechanical response of zirconium-based fuel claddings during a Loss-Of-Coolant-Accident (LOCA) hypothetical situation in Pressurized Water Reactors (PWR). In such a situation, fuel cladding tubes undergo a high temperature transient (potentially up to 1200°C) in an oxidizing environment and swell due to an internal pressure loading. Current safety requirements focus on cladding tube burst and swelling, which must be avoided or limited in order to guarantee fuel containment and core cooling [START_REF] Hache | The history of LOCA embrittlement criteria[END_REF] .

Oxidation at high temperature induces the formation of two layers in addition to the metallic beta phase: an outer layer of zirconia and a metallic sublayer of oxygen stabilized alpha phase, called alpha(O) (Figure 1), containing 2 -7 wt% of oxygen [START_REF] Corvalan-Moya | Numerical modeling of oxygen diffusion in the wall thickness of Low-Tin Zircaloy-4 fuel cladding tube during high temperature (1100-1200°C) steam oxidation[END_REF] .

The usual-phenomenological behavior laws used in the modeling of fuel claddings behavior in LOCA conditions are established for materials with low oxidation levels. Experimental works on highly oxidized materials will allow us to understand and quantify the role played by the oxidation. In the long term, a physics-based simulation tool modeling cladding's behavior in LOCA conditions should be built.

The creep behavior of zirconium and of its alloys have been extensively studied in an inert environment and are well understood, even though some interrogations remain in the alpha + beta domain [START_REF] Kaddour | Experimental determination of creep properties of zirconium alloys together with phase transformation[END_REF][START_REF] Kaddour | Microstructural influence on high temperature creep flow of Zr-1%NbO alloy in near-alpha, alpha+beta, and beta temperature ranges in a high vacuum environment[END_REF][START_REF] Massih | High temperature creep and superplasticity in zirconium alloys[END_REF] . Several studies were conducted on the effect of oxygen on the mechanical behavior of Zircaloy-2 at high temperature. Rizkalla et al. [START_REF] Rizkalla | Effect of oxygen on the deformation of zircaloy-2 at elevated temperatures[END_REF] and Choubey et al. [START_REF] Choubey | Flow stress of oxygenenriched zircaloy-2 between 1023 and 1873K[END_REF] performed compressive tests between 750°C and 1400°C and Choubey et al. [START_REF] Choubey | Flow stress of oxygenenriched zircaloy-2 between 1023 and 1873K[END_REF] showed that the compressive yield stress exponentially increases with increasing oxygen concentration (up to 1.2 wt%) in alpha and beta phases. Similar exponential influence of oxygen was found for the creep resistance of the alpha phase by Burton et al. [START_REF] Burton | The interaction of oxidation and creep in Zircaloy-2[END_REF] between 700 and 800°C but they observed a negligible effect on the creep behavior of the beta phase between 1050 and 1250°C. Comparable results were obtained with Zircaloy-4 by Chow et al. [START_REF] Chow | Creep behavior of oxidized Zircaloy-4 fuel sheathing[END_REF] between 650 and 1400°C. These studies were restricted to oxygen content up to 1.5 wt%. Only the dislocation creep regime was observed in both alpha(O) and beta phases and it was reported that this range of oxygen concentration does not significantly affect the stress dependence of the creep rate. One study 10 investigated oxygen concentrations up to 2 wt% by performing compressive tests. The creep behavior of the alpha phase enriched in oxygen for higher contents has not been reported yet.

This study aims at modeling with physic-based tools the mechanical response at high temperature of the oxidized cladding tube seen as a stratified material. The behavior of the alpha(O) phase at high temperature was studied first for the lower oxygen content, i.e. 2 wt%, in order to define a lower bound of the strengthening effect of the alpha(O) layer. Model alpha(O) specimens enriched by 2 wt% in oxygen were produced from zirconium alloys cladding specimens. To characterize their creep behavior at high temperature, tensile tests were performed along the axial direction under secondary vacuum. Viscoplastic flow rules have been identified. The creep behavior of the fuel cladding tube at high temperature and under internal pressure was then simulated by using finite element analysis (FEA) and the impact of the alpha(O) layer, for several thicknesses, has been estimated.

II. EXPERIMENTAL PROCEDURES

II.A. Creep tests in vacuum on model alpha(O) cladding specimens

Model alpha(O) phase specimens enriched by about 2 wt% of oxygen were produced from zirconium-based alloy cladding specimens thanks to a two-step procedure: an oxidation followed by a homogenization heat treatment. First, cladding specimens were oxidized in steam at 1100°C in order to get the required oxygen content as a zirconia layer. By using optical and scanning electron microscopy, the thicknesses of the resulting layers of zirconia and alpha(O) were found to be consistent with previous oxidation data [START_REF] Brachet | Hydrogen content, preoxidation, and cooling scenario effects on postquench microstructure and mechanical properties of Zircaloy-4 and M5 alloys in LOCA conditions[END_REF] .

Then, oxidized specimens were put in an electricalmechanical machine [START_REF] Kaddour | Experimental determination of creep properties of zirconium alloys together with phase transformation[END_REF] under load control and under secondary vacuum (10 -3 -10 -4 Pa). Two cement-pasted alumina rings delimited the 20 -25 mm gauge region over which axial elongation was continuously measured using a laser extensometer. Temperature was monitored using S-type thermocouples spot-welded onto the specimen surface at various locations. The temperature gradient along the gauge length was 10±5°C. Specimens underwent the homogenization heat treatment at 1200°C for 3 hours under a very low tensile load in order to prevent from buckling. During this treatment, zirconia dissolved and oxygen diffused through the specimen wall [START_REF] Stern | Metallurgical and mechanical behavoir of PWR cladding materials oxidized at high temperature[END_REF] .

Right after the homogenization heat treatment, temperature was decreased down to that of tensile creep tests and a constant load was applied (Figure 2). As already done for unoxidized specimens [START_REF] Kaddour | Experimental determination of creep properties of zirconium alloys together with phase transformation[END_REF] , different load jumps were applied and thus several axial stresses and corresponding strain rates were identified from each test. Loading history effects were quantified in one test by using increasing, then decreasing load levels. The difference in true strain rate versus true stress curves was lower than 30%. Investigated stresses and temperatures respectively ranged from 2 to 31 MPa and from 900 to 1100°C. Fig. 2. Scheme of the thermal-mechanical history of the oxygen-enriched creep specimens.

II.B. Creep test in steam

One creep test in steam was performed on stressrelieved Zircaloy-4 cladding in the EDGAR test facility [START_REF] Forgeron | Experiment and modeling of advanced fuel rod cladding behavior under LOCA conditions : alpha-bêta phase transformation kinetics and EDGAR methodology[END_REF] . A sealed tubular specimen was heated by Joule effect (~20°C.s -1 ) up to 1100°C and temperature was monitored by using an optical pyrometer. Once the test temperature stabilized (~20 s), a constant internal pressure of 5 bar was imposed. The values of temperature and applied pressure were chosen to enhance oxidation and its possible effects on the viscoplastic response of the cladding. The increase in outer diameter was continuously measured by a laser extensometer at mid-length.

II.C. Metallurgical analysis

Using in-house experimental procedures [START_REF] Brachet | Hydrogen content, preoxidation, and cooling scenario effects on postquench microstructure and mechanical properties of Zircaloy-4 and M5 alloys in LOCA conditions[END_REF] , optical and SEM observations were performed on polished and etched samples. For oxygen quantification, electron probe microanalysis (EPMA) was systematically used on tested enriched specimens with a ±0.1 wt% typical accuracy [START_REF] Brachet | Hydrogen content, preoxidation, and cooling scenario effects on postquench microstructure and mechanical properties of Zircaloy-4 and M5 alloys in LOCA conditions[END_REF] . Texture was studied thanks to electron backscatter diffraction (EBSD) analysis, on a specimen sectioned perpendicular to the tube axis. In view of the microstructure of Figure 1, the size of the sampled region was about 2 mm², with a step size of 2 µm, covering about 500 grains.

III. RESULTS AND DISCUSSION

III.A. Creep tests in vacuum on model alpha(O) cladding specimens

III.A.1. The oxygen-enriched cladding microstructure

EPMA analysis confirmed that the oxygen concentration was uniform through the wall, ranging from 2 to 2.5 wt%. Figure 3 shows the outer part of the cladding, where there is no remaining zirconia. The material microstructure is composed of large alpha(O) grains, their size ranges from 30 to 200 µm. The alpha(O) phase predominates with a phase fraction of 85±5%. However, a 15±5% residual fraction of prior-beta phase remains close to the inner surface after the annealing treatment, whatever the subsequent creep test temperature. This suggests that this amount of beta phase was already present at the annealing temperature of 1200°C. This point is in agreement with previous studies on oxygen-enriched Zircaloy-4 9,14,15 (Figure 4). Based the curves of Figures 4, one expects a 10 -25% beta phase fraction at 1200°C for a global 2 -2.5 wt% oxygen concentration. On the other hand, diffusion might not have been sufficiently rapid to dissolve this beta phase during the creep tests performed at lower temperature (900 -1100°C). The corresponding fraction of beta phase was thus assumed to be present at high temperature during the entire creep test. It had to be taken into account in the interpretation of the mechanical results.

III.A.2. Creep ductility of the oxygen-enriched material in vacuum

From the creep tests, the 2 wt% oxygen-enriched cladding can sustain significant strain levels when tested in axial creep between 900 and 1100°C. Moreover, it was not possible to obtain the ductility at fracture for each specimen due to the limited elongation allowed by the extensometry device. A creep test at 900°C was stopped without failure after a 15% axial strain. Another test was performed at 1000°C up to failure. Laser profilometry of the broken specimen indicated a 22% diametral strain in the necking region while SEM observation revealed a local reduction in thickness of almost 100% at fracture, together with some dimples (Figure 5). On the one hand, such ductile behavior contrasts with the brittle behavior reported for zirconium-based alloys with similar and higher oxygen contents at lower temperatures, 25 and 135°C [START_REF] Stern | Metallurgical and mechanical behavoir of PWR cladding materials oxidized at high temperature[END_REF] . On the other hand, Tseng et al. [START_REF] Tseng | Deformation behavior of duplex zircaloy-4-oxygen alloys[END_REF] reported strain values higher than 8% during compression tests performed between 850 and 1100°C on 1.8 wt% oxygenenriched Zircaloy-4, so that high-temperature ductility of such specimens is in agreement with scarce available literature data.

III.A.3. Creep flow behavior of the oxygen-enriched material in vacuum

For every considered load level, only steady-state creep was observed. Creep results in the axial direction in vacuum are summarized in Figure 6. Two creep regimes are observed. For higher stresses (σ > 10 MPa), the stress exponent n is close to 5 and, thus, suggests a "dislocation creep" regime [START_REF] Frost | Deformationmechanism maps: the plasticity and creep of metals and ceramics[END_REF] . For lower stresses (σ < 5 MPa), data are consistent with a stress exponent n close to 1. These two macroscopic creep regimes are similar to the ones reported at lower temperatures in the alpha phase that is non-enriched in oxygen [START_REF] Kaddour | Experimental determination of creep properties of zirconium alloys together with phase transformation[END_REF] . However the grain size is rather coarse in the oxygen-enriched material and thus the linear regime observed here might result from another physical mechanism than a diffusional flow [START_REF] Frost | Deformationmechanism maps: the plasticity and creep of metals and ceramics[END_REF] : a Harper-Dorn mechanism for example In order to distinguish between these two possible mechanisms, materials with various grain sizes should be tested, which is out of the scope of the present study. Based on values of n respectively equal to 1 and 5, experimental results were found to be well described by the following power-law equation:

      - = RT Q T A n exp σ ε& (1)
The apparent values of activation energy Q are both in agreement with those previously reported in non-enriched alpha phase in Zircaloy-4 and Zr-1%Nb [START_REF] Kaddour | Experimental determination of creep properties of zirconium alloys together with phase transformation[END_REF] .

For comparison purposes, creep flow results for alpha and beta single-phase of typical zirconium alloys non-enriched in oxygen, also tested along the axis of the cladding tubes are also plotted in Figure 6. The alpha phase behavior was extrapolated to higher temperatures from reported data [START_REF] Kaddour | Experimental determination of creep properties of zirconium alloys together with phase transformation[END_REF] . Indeed, in this temperature range alpha single phase does not exist without oxygen addition. The creep resistance of the 2 -2.5 wt% oxygen enriched material is much higher than that of zirconium alloys nonenriched in oxygen (i.e. containing 0.12 -0.14 wt% of oxygen).

Creep properties established by Chow et al. [START_REF] Chow | Creep behavior of oxidized Zircaloy-4 fuel sheathing[END_REF] for alpha phase of Zircaloy-4 containing 0.14 -1.5 wt% of oxygen are also displayed in Figure 6. They are based on experimental creep results ranging from 650 to 1000°C and for strain rates between 5.10 -5 and 7.10 -2 s -1 . Therefore, the transition to a nearly linear creep regime could not be expected from these existing results. Nevertheless, the expected strengthening effect of the oxygen weight concentration c in the dislocation creep regime, modeled by Chow with an exponential term (Eq. 2 and 3) overestimates the strain rate by a factor of about 3 with respect to the present experimental values obtained on a biphasic material (85% alpha(O)) containing 2 -2.5 wt% of oxygen.

) exp( . ) ( 0 Bc c - = ε ε & & (2)       - × × = RT 320000 exp 10 03 . 2 43 . 5 5 0 σ ε& ; 8 . 2 = B
(3)

III.A.4. Estimate of the creep behavior of the alpha(O) phase layer alone

Our purpose is to describe the behavior of the alpha(O) phase layer during a LOCA transient -internal pressure at high temperature -in order to model the threelayered oxidized cladding. This approach will reduce the uncertainties of the modeling in high oxidation ranges. Performing tensile creep tests in vacuum allowed us to establish viscoplastic flow rules of the 2 wt% oxygenenriched material. These equations can be used to describe the alpha(O) layer of an oxidized cladding only after discussing the three following points:

• The mechanical influence of the small amount (15±5%) of beta phase present in the model material but not in the alpha(O) layer of claddings oxidized in steam. • The possible difference in texture (arising from difference in phase transformation history) and in mechanical loading (axial tension vs. internal pressure). • The investigated 2 wt% oxygen concentration lower than the 4 wt% mean oxygen content in the alpha(O) layer 2 only gives a lower bound of the strengthening effect of the alpha(O) layer on the three-layered structure.

Mechanical influence of the 15±5% beta phase A homogenization approach was used to evaluate the impact of the 15±5% remaining fraction of beta phase on the viscoplastic flow behavior of oxygen enriched specimens. The fact that the experimental creep rates are strongly lower than the ones for as-received beta and alpha phases (up to 4 orders of magnitude) indicates that the major alpha(O) phase is stronger than the beta phase (which cannot dissolve high amount of oxygen, Figure 4) and probably controls the global mechanical response of oxygen-enriched specimens.

With such a contrast in viscoplastic behavior between alpha(O) and beta phases, and by considering the high fraction of alpha(O) phase, an inverse approach based on a two-phase mean field model with a Taylor assumption [START_REF] Taylor | Plastic strain in metals[END_REF] -homogeneous strain rate -can give a relevant estimate of the viscoplastic flow behavior of the alpha(O) phase alone. The viscoplastic flow rule of the beta phase of Zircaloy-4 was taken from tensile creep results on cladding tubes reported by Kaddour et al. [START_REF] Kaddour | Experimental determination of creep properties of zirconium alloys together with phase transformation[END_REF] as in Figure 6. As a starting point, flow rule coefficients of the alpha(O) phase were taken from Table 1. Then, assuming that the strain rate was homogeneous, i.e. equal in the two phases (Eq. 4), stresses in the alpha and beta phases (respectively, σ α and σ β ) and the macroscopic stress, Σ, were numerically calculated (from above flow rules and by using Eq. 5) using a least squares iterative method.

β α ε ε = = Ε (4) β β α β σ σ f f + - = Σ ) 1 ( (5)
where β f is the fraction of beta phase. By these means, the influence of the beta phase on the viscoplastic flow rule was estimated to a factor two in the higher stress regime, and was negligible in the lower stress regime. This factor two was thus applied in the higher stress range only by dividing by two creep coefficient A of Eq. 1 -identified from the axial creep tests -leading to the curves of Figure 7. From this point, it appears that the model by Chow et al. overestimates the strain rate by a factor 6.

Influence of texture and mechanical loading mode

By considering the anisotropy in mechanical properties of the hexagonal alpha phase, transfer of the obtained viscoplastic flow rules to another mechanical loading mode and material texture has to be done carefully. EBSD analysis revealed that the textures of both our oxygen-enriched specimens and the alpha(O) layer of cladding pressurized in steam are comparable: in both, <c> axes were perpendicular to the cladding axis. However axial tension differs from biaxial loading resulting from prescribed internal pressure. Nevertheless, stress and strain anisotropy of oxygen-enriched alpha(O) at 900 -1100°C cannot be inferred from that of more conventional alpha phase, for which no data are available at such high temperature. Consequently, this effect could not be numerically taken into account in the present study and the alpha(O) layer was considered as isotropic.

Oxygen concentration effect

The oxygen concentration varies between 2 and 7 wt% within the alpha(O) layer of an oxidized cladding [START_REF] Corvalan-Moya | Numerical modeling of oxygen diffusion in the wall thickness of Low-Tin Zircaloy-4 fuel cladding tube during high temperature (1100-1200°C) steam oxidation[END_REF] . Therefore, the identified flow equations for 2 wt% give a lower bound of the strengthening effect of oxygen. Further studies are under way to assess the experimental creep behavior of model materials with higher oxygen content. An evaluation of the upper bound strengthening effect for 7 wt% oxygen content based on the extrapolated strengthening law of Chow et al. [START_REF] Chow | Creep behavior of oxidized Zircaloy-4 fuel sheathing[END_REF] is shown in part IV.

III.B. Creep test in steam

In a steam environment, at the beginning of the creep test, the observed strain rate is similar to the one observed under vacuum. After a short time (~40 s) the increase in diameter is strongly slowed down (Figure 8) -significant oxidation having already started during heating and temperature stabilizing stage. The predicted weight gain in oxygen according to the data from Brachet et al. [START_REF] Brachet | Hydrogen content, preoxidation, and cooling scenario effects on postquench microstructure and mechanical properties of Zircaloy-4 and M5 alloys in LOCA conditions[END_REF] is plotted in the same figure. From the tests carried out in vacuum [START_REF] Kaddour | Experimental determination of creep properties of zirconium alloys together with phase transformation[END_REF][START_REF] Kaddour | Microstructural influence on high temperature creep flow of Zr-1%NbO alloy in near-alpha, alpha+beta, and beta temperature ranges in a high vacuum environment[END_REF] , no "primary creep stage" is expected, nor observed, in this range of stress level and test temperature. This is a clear demonstration that the combined formation of zirconia and alpha(O) reinforces the oxidized cladding with respect to creep. Such strengthening was also estimated by modeling in the following part.

IV. MODELING OF THE CREEP BEHAVIOR OF THE OXIDIZED CLADDING SEEN AS A STRATIFIED MATERIAL

IV.A Constitutive equations

The Zircaloy-4 one-side oxidized cladding can be considered as a three-layer-stratified material at high temperature (1100°C): the inner beta layer, an intermediate alpha(O) layer and an outer zirconia layer. In the present study, the mechanical contribution of zirconia could not be estimated, its mechanical strength being unknown at the studied temperatures. This contribution can possibly be evaluated by inverse analysis once the creep strength of the alpha(O) layer with its gradient in oxygen content will be fully elucidated. In the present calculations that aim at providing a lower bound to environmentally-induced creep strengthening, only two layers were considered in the model, namely, an inner beta layer and an outer alpha(O) layer. Their creep flow rules are respectively given by Kaddour et al. [START_REF] Kaddour | Experimental determination of creep properties of zirconium alloys together with phase transformation[END_REF] for the beta phase of Zircaloy-4 and by corrected creep flow equations based on our experimental results for the 2 wt% alpha(O) or by the extrapolated Chow creep laws for the 2 -4 -7 wt% alpha(O) layer. Moreover the effect of oxygen atoms having diffused into the beta phase on its creep behavior was neglected, in agreement with the results of Burton et al. [START_REF] Burton | The interaction of oxidation and creep in Zircaloy-2[END_REF] and Chow et al. [START_REF] Chow | Creep behavior of oxidized Zircaloy-4 fuel sheathing[END_REF] . In absence of experimental data about creep anisotropy of this layer, an isotropic von Mises equivalent stress was used in the viscoplastic potential. No work-hardening and threshold stress were introduced, as experimentally observed in this range of temperature and stress levels. The test temperature of 1100°C was chosen for creep simulations according to experimental conditions.

IV.B Finite element model and numerical procedures

Using the finite element code Cast3M 2012 18 , the mechanical behavior of the oxidized cladding was modeled for creep tests under internal pressure in an inert atmosphere. Therefore, no growth of the alpha(O) layer was considered during the whole simulation. A small portion (1° sector) of the cladding was meshed with quadratic brick elements. The number of elements was constant for the beta layer regardless of it thickness. On the other hand, the number of elements for the alpha(O) layer was adapted to its thickness. Boundary conditions were chosen according to the cylindrical symmetry of the problem (Figure 9). The mesh geometry was updated at each time step. A 1-second time step was used. It was checked that a smaller mesh size and a smaller time step did not affect the macroscopic response of the structure.

IV.C Simulation results

The simulated diametral strain as a function of time for a creep test simulation under 5 bar internal pressure and at 1100°C in inert atmosphere is shown in Figure 10 for various alpha(O) phase layer thicknesses. At 1100°C, replacing an outer layer of 5, 10 and 50 µm of beta phase by alpha(O) phase divide the initial strain rate in the cladding by a factor 1.7, 2.9 and 11 respectively. The thin layer of strong alpha(O) phase thus strongly enhances the resistance to viscoplastic hoop deformation under internal pressure. For the range of alpha(O) layer thicknesses expected from oxidation in steam at 1100°C during the creep test, the experimental strain rate of the cladding is still significantly lower than the ones predicted by these simulations. As expected, the predicted strengthening effect due to a 2 wt% alpha(O) layer is a lower bound of the strengthening due to oxidation in steam. Extrapolating the law established by Chow et al. [START_REF] Chow | Creep behavior of oxidized Zircaloy-4 fuel sheathing[END_REF] for the highest oxygen content, i.e. 7 wt%, gives an another estimate closer to an upper bound of the effect of the alpha(O) layer (Figure 11). However extrapolating the same law 9 for 4 wt%, which is the mean oxygen content in the alpha(O) layer, leads to an insufficient prediction of reinforcement of the cladding. To improve this result, two assumptions can be made:

• On the one hand the increase in creep resistance of the alpha(O) phase for 4 wt% of oxygen can be underestimated by extrapolating the strengthening equation proposed by Chow et al. 9 based on experimental results until 1.5 wt%. The difference between our creep results on 2 wt% oxygen-enriched material and the prediction based on Eq. 2 supports this solution (Figure 6). • On the other hand, if the behavior of the alpha(O) layer is correctly described, the mechanical contribution of zirconia has also to be taken into account in order to predict the total reinforcement of the cladding.

Creep tests on 4 wt% oxygen-enriched materials will help answering this question in a near future. Despite its limitations, the model presented here is thus a first, promising step toward improvement of prediction of the high temperature mechanical behavior of cladding tubes.

V. CONCLUSIONS

This study aimed at modeling the high temperature behavior of zirconium-based cladding in LOCA conditions. In order to model the oxidized cladding as a layered material, the constitutive behavior of the beta layer was taken from literature, so that the experimental part of the study focused on characterizing the alpha(O) layer, in particular under low stress levels.

Zirconium-based specimens enriched in oxygen by up to about 2 wt% were obtained and microstructurally characterized. Then, tensile axial creep tests were performed in a vacuum atmosphere. It was observed that this material exhibits high ductility in the temperature range from 900 to 1100°C, a 15% axial strain being reached under tension without fracture. In addition, oxygen-enriched specimens showed a high creep resistance, in agreement with the known mechanical strengthening effect of the oxygen. For the first time to the authors' knowledge, two creep regimes (creep exponents of ~1 and ~5, respectively, for lower and higer stresses) were observed, similarly to the two creep regimes already reported in the alpha phase of as-received zirconium alloys.

Homogenization calculations showed the remaining 15±5% fraction of beta phase in the oxygen enriched sample has a limited, but non-negligible impact (a factor two) on the viscoplastic strain rate.

As a consequence, the 2 wt% oxygen enriched power laws were used to model the alpha(O) behavior in FE computations of a two-layered tube.

These calculations with a 2 wt% oxygen enriched alpha(O) phase highlight a significant strengthening effect, up to a decade in hoop strain rate under internal pressure, which is a lower bound. On the other hand calculations based on extrapolated strengthening effect for 7 wt% of oxygen overestimate strengthening as expected.

Additional experimental tests are planned on 4 wt% and 7 wt% oxygen enriched samples -if technically feasible-to estimate the impact of the alpha(O) layer, and thus, of the oxide layer on the behavior of the layered cladding in LOCA conditions more precisely.

Fig. 1 .

 1 Fig. 1. Optical micrograph of a Zy-4 cladding after creep test in steam at 1100°C for 10 min. Outer surface on the left, tube axis perpendicular to the image.

Fig. 3 .

 3 Fig. 3. Optical micrograph of an oxidized, then homogenized sample containing 2 wt% of oxygen. Same orientation as in Figure 1, outer surface on the left.

Fig. 4 .

 4 Fig.4. Zircaloy-4 -oxygen isopleth from the equilibrium phase diagram[START_REF] Chow | Creep behavior of oxidized Zircaloy-4 fuel sheathing[END_REF][START_REF] Chung | Pseudo-binary zircaloy-oxygen phase diagram[END_REF][START_REF] Woo | Transformation characteristics of rapidly heated and quenched Zircaloy-4 oxygen alloys[END_REF] . The oxygen content and homogenization temperature are indicated with a red ellipse.

Fig. 5 .

 5 Fig. 5. SEM micrograph of the fracture surface of the 2 -2.5 wt% oxygen enriched-material broken at 1000°C.

Fig. 6 .

 6 Fig. 6. Stress dependence of the strain-rate: present experimental data (symbols) for oxygen-enriched material, predictions at 1100°C for zirconium alloys which are non-enriched 3 (color lines) or enriched with oxygen 9 (grey lines).

Fig. 7 .

 7 Fig. 7. Creep test results (symbols) and viscoplastic flow rules (lines), all determined in uniaxial tension along the tube axis at 900°C.

Fig. 8 .

 8 Fig. 8. Predicted weight gain based on data from Brachet et al. 11 versus time, together with creep strain of Zy-4 at 1100°C and 5 bar of internal pressure in steam.

Fig. 9 .

 9 Fig. 9. Mesh, loading and boundary conditions used in FE simulations of internal pressure loading at 1100°C.

Fig. 10 .

 10 Fig. 10. Creep test in steam (red line) and creep simulations in inert atmosphere for various alpha(O) layer thicknesses (grey lines) -Zy-4 at 1100°C and 5 bar.

Fig. 11 .

 11 Fig. 11. Creep test in steam (red line) and creep simulations in an inert atmosphere for a 20-µm-thick alpha(O) layer with various oxygen contents (grey lines) -Zy-4 at 1100°C under 5 bar.

& & &

ACKNOWLEDGMENTS

The authors would like to thank J. Heurtel and A. Laurent from the Centre des Matériaux, and also D. Hamon, E. Rouesne, and S. Urvoy from CEA for their contribution to this work.


domain: math sdv

text:
Simon Garnotel¹³
email: simon.garnotel@u-picardie.fr

Stéphanie Salmon²

Olivier Balédent¹³

Numerical Cerebrospinal System Modeling in Fluid-Structure Interaction

Keywords: Cerebrospinal fluid, Intracranial pressure, Stroke volume, Fluidstructure interaction, Finite element method

Objective Cerebrospinal fluid (CSF) stroke volume in the aqueduct is widely used to evaluate CSF dynamics disorders. In healthy population, aqueduct stoke volume represents around 10% of the spinal stroke volume while intracranial subarachnoid space stroke volume represents 90%. Amplitude of the CSF oscillations through the different compartments of the cerebrospinal system is function of the geometry and the compliances of each compartments but we suspect that it could also be impacted be the cardiac cycle frequency. To study this CSF distribution, we have developed a numerical model of the cerebrospinal system taking into account cerebral ventricles, intracranial subarachnoid spaces, spinal canal and brain tissue in fluid-structure interaction. Materials and Methods A numerical fluid-structure interaction model is implemented using a finite element method library to model the cerebrospinal system and its interaction with the brain based on the fluid mechanics equations and the linear elasticity equations coupled in a monolithic formulation. The model geometry, simplified in a first approach, is designed in accordance to realistic volumes ratio of the different compartments: a thin tube is used to mimics the high flow resistance of the aqueduct. CSF velocity and pressure, and brain displacements are obtained as simulation results, whence CSF flow and stoke volume are calculated from this results.

Results Simulation results show significant variability of aqueduct stroke volume and intracranial subarachnoid spaces stroke volume in the physiological range of cardiac frequencies.

Conclusions Fluid-structure interactions are numerous in the cerebrospinal system and difficult to understand in the rigid skull. The presented model highlights significant variations of the stroke volumes under cardiac frequency variations only.

Introduction

Cerebrospinal fluid (CSF) stroke volume in the aqueduct (SVaq) is widely used to evaluate CSF dynamics disorders. In healthy population, SVaq represents around 10% of the spinal stroke volume (SVspi) while intracranial subarachnoid spaces (SAS) stroke volume (SVsas) represents around 90% of SVspi [START_REF] Balédent | Relationship between cerebrospinal fluid and blood dynamics in healthy volunteers and patients with communicating hydrocephalus[END_REF]. Amplitude of the CSF oscillations through the different compartments o the cerebrospinal system is function of the geometry and the compliances of each compartments but we suspect that it could also be impacted by the cardiac frequency. To study the CSF distribution, a numerical model of the cerebrospinal system has been developed taking into account cerebral ventricles, intracranial subarachnoid spaces, spinal canal and brain tissue in fluid-structure interaction.

Previous numerical models of the cerebrospinal system have already been studied, using fluid equations only [START_REF] Alperin | MRintracranial pressure (ICP): a method to measure intracranial elastance and pressure noninvasively by means of MR imaging: baboon and human study[END_REF] or fluid and porous media equations [START_REF] Linninger | Cerebrospinal fluid flow in the normal and hydrocephalic human brain[END_REF] but they are not relevant in this study, so a new numerical framework is introduced. Model geometry and parameters are detailed and numerical simulations are run. A first step of validation is done before results presenting and discussion.

Materials and Methods

Numerical framework

A numerical fluid-structure interaction model is implemented using the finite element method library FreeFem++ [START_REF] Hecht | New development in FreeFem++[END_REF] to model the cerebrospinal system and its interactions with the brain. This model is based on the fluid mechanics equations, the Navier-Stokes equations, and the solid mechanics equations, the linear elasticity equations, coupled in a monolithic formulation [START_REF] Sy | Algorithme for solving fluid-structure interaction problem on a global moving mesh[END_REF].

The Navier-Stokes equations, that describe the fluid behavior in a moving domain, read as follow:

𝜌 𝐹 𝜕𝒖 𝜕𝑡 + 𝜌 𝐹 ((𝒖 -𝑈). ∇)𝒖 -2𝜇 𝐹 ∇. 𝝐(𝒖) + ∇𝑝 = 𝒇 𝐹 ∇. 𝒖 = 0 Where 𝒖 is the velocity, 𝑝 the pressure, 𝑈 the domain velocity, 𝒇 𝐹 the volumetric external force (taken here to zero, gravity is neglected), 𝜌 𝐹 the density and 𝜇 𝐹 the dynamic viscosity for the fluid; and where:

𝝐(𝒖) = 1 2 (∇𝒖 + (∇𝒖) 𝑇 )
The linear elasticity equations, that describe the solid deformation, read as follow:

𝜌 𝑆 𝜕 2 𝒅
𝜕𝑡 2 -∇. 𝝈 𝑆 (𝒅) = 𝒇 𝑆 Where 𝒅 is the displacement, 𝒇 𝑆 the volumetric external force (taken here to zero, gravity is neglected), 𝜌 𝑆 the density for the solid; and where:

𝝈 𝑆 (𝒅) = 𝜆(∇. 𝒅)𝑰 + 2𝜇 𝑆 𝝐(𝒅) Where 𝜆 and 𝜇 𝑆 are the Lamé's coefficients and 𝑰 is the identity matrix. Lamé's coefficients are defined using the Young's modulus and the Poisson's ratio as follow:

𝜆 = 𝜈𝐸 (1 -2𝜈)(1 + 𝜈) 𝜇 𝑆 = 𝐸 2(1 + 𝜈)
Where 𝐸 is the Young's modulus and 𝜈 the Poisson's ratio. In order to couple this two equations, the continuity of velocity and stress is imposed at the fluid-structure interface as follow:

𝒖 = 𝜕𝒅 𝜕𝑡 𝜎 𝐹 (𝒖, 𝑝)𝒏 = -𝝈 𝑆 (𝒅)𝒏
Where 𝒏 is the normal vector and:

𝝈 𝐹 (𝒖, 𝑝) = -𝑝𝑰 + 2𝜇 𝐹 𝜖(𝒖) The monolithic formulation of this problem leads to the following matrix form, where fluid and structure problems are strongly coupled and solved at the same time:

( 𝐴 𝐵 𝑇 0 𝐵 𝜖 0 0 0 𝜖 ) ( 𝑈 𝑃 𝐹 𝑃 𝑆 ) = ( 𝐿 0 0 )
Where 𝜖 is a stabilization term, 𝐴 and 𝐵 are the representatives matrix of the problem, 𝑈 is the fluid velocity, 𝑃 𝐹 is the fluid pressure, 𝑃 𝑆 is the solid pressure (equal to zero) and 𝐿 represents external and coupling conditions.

Geometry

The model geometry, simplified in a first approach, is deigned in accordance to realistic volumes ratio of the different compartments of the cerebrospinal system. Brain, intracranial subarachnoid spaces and ventricles are designed using two dimensions' disks and a thin tube is used to mimics the high flow resistance of the aqueduct, Fig. 1.

Figure 1. Geometry of the model

The cranium has a radius of 10 cm, intracranial SAS have a thickness of 0.75 cm, ventricles have a radius of 3 cm, aqueduct has a radius of 0.2 cm and spinal canal has a radius of 1 cm.

Parameters

Fluid and structure parameters, provided from the literature, are taken in accordance to the CSF and brain mechanical physiological parameters, Tab. 1. The model input is located on the spinal canal end where a velocity profile, close to the physiological one, is imposed:

𝑣 = 𝐴 sin ( 2𝜋𝑡 𝑝 )
Where 𝐴 is the velocity amplitude and 𝑝 the period of the cardiac cycle. Remark 1. In physiological conditions, CSF oscillations come from brain expansion at each cardiac cycle. In our model, this is the inverse phenomena: CSF oscillations cause the brain deformations. This behavior has been chosen, for future use of physiological measures, because CSF flow is easily measurable by phase contrast MRI; contrariwise, brain deformation is hard to obtain.

Simulations are done with different values of the heart rate period (𝑝) keeping the same SVspi, Fig. 2. Twenty cardiac cycles are simulated and the results are extracted from the last cardiac cycle. CSF velocity and pressure, and brain displacements are obtained as simulation results, whence CSF flow and stroke volume are calculated from this results.

Results

Validation

Model validity is successfully performed on numerical benchmarks (specific cases) to ensure the correct agreement of the algorithm. In addition, the flow conservation is calculated in the used model geometry during a cardiac cycle: the net flow error is 1.4%±1.3SD.

Study case

As a second validation, geometry and flow are measured in a healthy subject by MRI to design the model and impose an input velocity. Simulation results, Fig. 3 and Fig. 4, show the good agreement with the PC-MRI measured data and with the physiological ICP curve shape and amplitude.

Cardiac cycle impact

Simulation results, Fig. 5 and Fig. 6, show significant variability of SVaq and SVsas in the physiological range of cardiac frequencies; an inversion of the distribution is even observed for slow cardiac frequencies. Pressure gradient (maximum pressure minus minimum pressure during a cardiac cycle) is equally modified by the cardiac frequency variations; it is exponentially increased as the cardiac frequency growth.

Discussion

The linear elasticity is justified in this application since small deformations are obtained in the simulation results. The proposed model is validated in numerical (benchmark) and physiological (flow conservation) points of view.

Numerical results in the study case are in accordance to PC-MRI measured data. In addition, pressure curve presents a physiological amplitude and shape (peaks and valleys).

Cardiac frequency has an impact on the CSF distribution between the two main intracranial compartments: ventricles and intracranial SAS; stroke volume distribution is in a physiological range for normal (around 60 -100 bpm) and high cardiac frequencies, contrariwise, the stroke volume distribution tends to be equal in the two main compartments as cardiac frequency decrease.

Cardiac frequency has equally an impact on the pressure gradient in the two compartments, pressure gradient grown as cardiac frequency increase. As the SVspi is conserved, high cardiac frequency causes a rapid inflow of CSF into the intracranial compartment that can explain this phenomenon.

As the absolute pressure term does not appear in the Navier-Stokes formulation (only his gradient), it is impossible to obtain an information about the absolute pressure of the CSF: only the gradient of the intracranial pressure is calculated.

Conclusions

A numerical model is presented taking into account the numerous fluid-structure interaction in the cerebrospinal system closed in the rigid skull. This model highlights significant variations of the stroke volume distribution under cardiac frequencies variations only.

As perspectives, spatial distribution of the intracranial pressure gradient will be studied. An improvement of the model including arterial and venous flow will be designed.

Another perspective could be the determination of the absolute intracranial pressure regarding Marmarou's studies and the pressure volume curve [START_REF] Marmarou | Compartmental analysis of compliance end outflow resistance of the cerebrospinal fluid system[END_REF][START_REF] Marmarou | A nonlinear analysis of the cerebrospinal fluid system and intracranial pressure dynamics[END_REF].

Figure 2 .

 2 Figure 2. Input of the model Remark 2. A lot of cardiac cycles are simulated to ensure the stabilization of the system.

Figure 3 .

 3 Figure 3. Simulation results and PC-MRI measurements of the flow.

Figure 4 .

 4 Figure 4. Simulation results of pressure.

Figure 5 .Figure 6 .

 56 Figure 5. Stroke volume function of cardiac frequency

Table 1 .

 1 Mechanical parameters

	Parameter	Symbol	Value	Unity
	Fluid density	𝜌 𝐹	1	g.cm -3
	Fluid viscosity	𝜇 𝐹	0.001	g.cm -1 .s -1
	Solid density	𝜌 𝑆	1.1	g.cm -3
	Solid Young's modulus	𝐸	4.10 3	g.cm -1 s -2
	Fluid Poisson's ratio	𝜈	0.35	-