Luc Nguyen-Thê

Simone De Liberato

Motoaki Bamba

Cristiano Ciuti

Many-body physics of intersubband polaritons

Intersubband polaritons are light-matter excitations originating from the strong coupling between an intersubband quantum well electronic transition and a microcavity photon mode. In this paper we study how the Coulomb electron-electron interaction and the Pauli saturation of the electronic transitions affect the physics of intersubband polaritons. We develop a microscopic many-body theory for the physics of such composite bosonic excitations in a microcavity-embedded two-dimensional electron gas. As a first application, we calculate the modification of the depolarization shifts and the efficiency of intersubband polariton-polariton scattering processes.

I. INTRODUCTION

Intersubband polaritons are elementary excitations of microcavities embedding doped semiconductor quantum wells. In particular, they are the normal modes due to the strong coupling between a cavity photon mode and the transition between two conduction subbands in quantum wells containing a two-dimensional electron gas. Intersubband polaritons were experimentally demonstrated for the first time in 2003 1 and their physics is the object of many interesting investigations, both fundamental [2][3][4][5][6][7][8][9][10][11][12] and applied [13][14][15][16][17][18][19][20][21] .

So far, theoretical investigations have been based mostly on a bosonized Hamiltonian approach 2 . In the dilute regime, i.e., when the excitation density is much smaller than the electronic density in the doped quantum wells, intersubband polaritons behave as (composite) bosons. At higher excitation densities anyway, Pauli saturation effects due to electron fermionic statistics and deviations from the bosonic behavior are expected 22 .

An important aspect of the physics of such excitations is the role of the Coulomb interaction. In the case of bare intersubband transitions, Coulomb interaction is known not to modify dramatically the polaritonic lineshape 23 , being responsible mainly for a simple shift of the excitation frequency, the so-called depolarization shift. This is why, in many relevant situations, the presence of the Coulomb interaction can be in first approximation neglected, by using a renormalized (measured) value for the frequency of the intersubband transition instead of the bare one, calculated from the single-electron subband structure.

Very recently, the role of Coulomb interaction in intersubband polariton systems has started to attract interest. In Ref. [24] the part of the Coulomb interaction responsible for the depolarization shift has been studied in the bosonic excitation formalism, allowing to calculate the renormalization from first principles. Other works have studied the same problem using different gauges 25,26 and a Green functions approach 27 . As a matter of fact, up to now, only depolarization or static effects have been taken into account. Hence, the general problem of the many-body physics of intersubband polaritons is largely unexplored, in particular the physics of polariton-polariton interactions and polariton-polariton scattering.

In this paper we derive a microscopic many-body theory describing the intersubband polariton-polariton scattering due both to the Coulomb interaction and to the non-perfect bosonicity of the intersubband excitations. As a first application of the present theoretical formalism, we calculate how the depolarization shifts are modified and, most importantly, we determine the efficiency of the intersubband polariton-polariton scattering processes.

This article is organized as follows. In Sec. II, we present the physical system, introduce the fermionic electron-photon Hamiltonian including the different channels of the (screened) Coulomb interaction. Moreover, we write the intersubband excitations in second quantized formalism. In Sec. III, the composite boson commutator approach 28 to calculate many-body matrix elements is presented for the case of intersubband excitations. In Sec. IV, an effective bosonic Hamiltonian (up to the fourth order in the bosonic operators) is derived in such a way that the two-excitation properties are the same as those predicted by the microscopic composite boson approach. In Sec. V numerical applications of the present theory to the calculation of the interaction energy between intersubband excitations are presented. Sec. VI is devoted to the calculation of polariton-polariton scattering processes. Finally, conclusions are drawn in Sec. VII. Additional technical details are included in Appendix A, B and C.

II. THE SYSTEM AND ITS HAMILTONIAN DESCRIPTION A. Physical system

We will consider a photonic, planar microcavity embedding n QW symmetric quantum wells, each doped with an areal electron density of n el (see Fig. 1 for a schematic representation of such a system). As the quantum wells confine the electrons along the growth direction, the elec- tronic bands will be split into multiple, almost parallel subbands. The optically active transitions we are interested in will thus be the ones from the last occupied conduction subband to the first unoccupied one. The photonic microcavity width will be chosen to have one of the photonic modes almost resonant with the electronic transition between these two electronic conduction subbands (the first two for simplicity, even if more general situations are possible 29 ). In order to simplify the notation, we will consider that the quantum wells are all identically coupled to the relevant microcavity field.

The Hamiltonian describing the coupled electronphoton system can be written in the Coulomb gauge as

H = H C + H F , (1)
where H C includes only photonic degrees of freedom, while H F regroups all the terms presenting also electronic ones. Introducing the annihilation operators c µ,k for an electron in the subband µ with an in-plane wavevector k, and a q for a microcavity photon with an in-plane wavevector q, we can write the two Hamiltonians explicitly as

H C = q ω c,q a † q a q + ∆ q (a q + a † -q )(a -q + a † q ) H F = H Kin + H Rabi + H Coul , (2)
where

H Kin = µ={1,2},k ω µ,k c † µ,k c µ,k (3)
H Rabi = k,q χ q (c † 2,k+q c 1,k + c † 1,k+q c 2,k )(a q + a † -q ) H Coul = 1 2 k,k ′ ,q µ,µ ′ ,ν,ν ′ ∈{1,2} V µνν ′ µ ′ q c † µ,k+q c † ν,k ′ -q c ν ′ ,k ′ c µ ′ ,k ,
are respectively the electronic kinetic energy, the lightmatter coupling Hamiltonian and the electron-electron Coulomb interaction. Coefficients ω c,q and ω µ,k are the frequencies of a photon with in-plane wavevector q and of an electron with in-plane wavevector k in the subband µ, χ q is the light-matter coupling frequency, ∆ q the diamagnetic interaction coefficient and the q = 0 term has to be excluded from the sum in H Coul in order to compensate for the homogeneous positive background of the jellium model [START_REF] Giuliani | Quantum theory of the electron liquid[END_REF] .

In order to simplify the notation both spin and quantum well indexes have been omitted but one should keep in mind that there is an implicit sum over them. The explicit notation is given in Appendix A. We also neglect to explicitly write the photon polarization because, due to the intersubband selection rule, only the TM modes are coupled with matter.

The Coulomb matrix elements V µνν ′ µ ′ q are given by 31

V µνν ′ µ ′ q = e 2 2Sǫ 0 ǫ r q I µνν ′ µ ′ q (4)
I µνν ′ µ ′ q = dzdz ′ ψ µ (z)ψ ν (z)ψ ν ′ (z ′ )ψ µ ′ (z ′ )e -q|z-z ′ | ,
where S is the surface of the sample and ψ µ (z) is the real z-component wavefunction of an electron in the subband µ.

Due to the symmetry of the quantum wells, and thus of the wavefunctions, a certain number of matrix elements in Eq. ( 4) can be seen to be zero, in particular all matrix elements with an odd number of 1 and 2 indices

V 1112 q = V 1121 q = V 1211 q = V 2111 q = 0 (5) V 2111 q = V 2212 q = V 2122 q = V 1222 q = 0.
The other elements obey the relations

V 1122 q = V 1212 q = V 2121 q = V 2211 q (6) V 1221 q = V 2112 q ,
so that there are only four distinct matrix elements, V 1212 q ,V 1221 q ,V 1111 q and V 2222 q . In Fig. 2 we show a graphical representation of these four qualitatively different processes and the wavevector dependency of their respective coefficients for an infinite potential well. From such a figure we can clearly see a qualitative difference in the long wavelength behavior between V 1212 q and the other three coefficients. As explained in greater detail in Ref. [24], they indeed refer to qualitatively different physical processes. While V 1212 q describes an intersubband, coherent Coulomb scattering, that is responsible for the so-called depolarization shift, the others describe intrasubband scattering processes, in which electrons are scattered inside the same subband.

B. Screened Coulomb interaction

The dense two dimensional electron gas in the first electronic subband screens the Coulomb interaction. In order to take it into account we will thus replace the bare Coulomb interaction V µνν ′ µ ′ q defined by Eq. ( 4), with its static RPA-screened version Ṽ µνν ′ µ ′ q . From Refs. [31,32], these coefficients obey the Dyson equation

Ṽ µνν ′ µ ′ q = V µνν ′ µ ′ q (7) + α,β V µβαµ ′ q Π αβ (q, ω = 0) Ṽ ανν ′ β q ,
where Π αβ (q, ω) is the RPA polarization function for the α → β transition. The screened interactions thus take the form

Ṽ 1νν1 q = V 1νν1 q 1 -V 1111 q Π 11 (q, 0) (8)
Ṽ 2222

q = V 2222 q + [(V 1221 q ) 2 -V 1111 q V 2222 q ]Π 11 (q, 0) 1 -V 1111 q Π 11 (q, 0) Ṽ 1212 q = V 1212 q 1 -V 1212 q [Π 12 (q, 0) + Π 21 (q, 0)] ,
where ν = {1, 2}. Screening of intrasubband and intersubband processes can thus be encoded respectively in the dielectric functions ǫ(q) and ǫ 12 (q), where

ǫ(q) = 1 -V 1111 q Π 11 (q, 0), (9)
and

ǫ 12 (q) = 1 -V 1212 q [Π 12 (q, 0) + Π 21 (q, 0)]. ( 10
)
As we will see in the rest of this Section, only the longwavelength limit is relevant in the study of intersubband polaritons. We will thus always consider the q → 0 limit of the above experessions. In Fig. 2 the reader can find the screened and unscreened profiles of the different potentials as a function of the exchanged momentum normalized over the Fermi wavevevector q/k F .

C. Intersubband excitations

The theory in this paper will be developed using a zero temperature (T = 0) formalism. In order for this approximation to be valid, we need to have a vanishing thermal electronic population in the excited subband, that is ω 2,0 -E F ≫ k B T , where E F is the Fermi energy and k B the Boltzmann constant. Under such a condition the ground state of the uncoupled light-matter system is given by

|F = |k|≤kF c † 1,k |0 , (11)
where k F is the Fermi wavevector and |0 is the vacuum state, with no electrons in both subbands and no photons in the microcavity. In the ground state all electronic states up to the Fermi energy are occupied and all the others are empty.

If we neglect the incoherent intrasubband Coulomb processes (due to the coefficients V 2222 q , V 1221 q and V 1111 q ) in Eq. ( 3), the full Hamiltonian H can be solved, in the dilute regime, in terms of bosonic coherent excitations named intersubband polaritons 2 .

In order to do this, we first notice that the microcavity photon mode is coupled through H Rabi to only one particular linear superposition of electronic transitions

b † 0,q = 1 n QW N el k ν * 0,k c † 2,k+q c 1,k , (12)
where ν 0,k = Θ(k F -k) and Θ is the Heaviside function. The choice of the notation with the index 0 will be explained in the following. In Eq. ( 12), as well as in Eq. ( 3), the sum is implicit over the spin and quantum well indexes to simplify notations. This defines an intersubband collective excitation, which we will name bright intersubband excitation. Hamiltonian H Rabi can thus be written exactly

H Rabi = q Ω q (b 0,-q + b † 0,q )(a q + a † -q ), (13)
where Ω q = n QW N el χ q is the Rabi frequency. Eq. ( 13) makes it clear that bright excitations concentrate all the oscillator strength of the electron gas. Other similar collective excitations can be constructed by an orthonormalization procedure

b † i,q = 1 n QW N el k ν * i,k c † 2,k+q c 1,k , (14)
where index i runs from 1 to n QW N el -1 and the ν coefficients have support over the Fermi sea and satisfy the orthonormality relation

1 n QW N el k ν * i,k ν j,k = δ i,j . ( 15
)
0 1 2 3 0 2 4 0 1 2 3 0 0.5 1 1.5 q/kF V 1111 q (a.u.) V 1111 q (a.u.) 0 1 2 3 0 2 4 0 1 2 3 0 0.5 1 1.5 q/kF V 2222 q (a.u.) V 2222 q (a.u.) 0 1 2 3 0 2 4 0 1 2 3 0 0.5 1 1.5 q/kF V 1221 q (a.u.) V 1221 q (a.u.) 0 1 2 3 0 0.5 1 0 1 2 3 0 0.5 1 q/kF V 1212 q (a.u.) V 1212 q (a.u.)
FIG. 2: Panels in the first two lines show the wavevector dependency of the the four qualitatively different Coulomb processes. In each panel we plot both the bare potentials (dashed line) and the static-RPA screened ones (solid line). The first row is for a mid-infrared transition, with ω12 = 140meV and an electronic density in each quantum well n el = 5 × 10 11 cm -2 . The second row is instead for a THz transition, with ω12 = 15meV and n el = 2 × 10 11 cm -2 . In the third row we present schematic representation of the four relevant processes.

However, none of these new collective excitations is coupled to the microcavity photon field. This is why we call them dark intersubband excitations. We will consider them in Sec. III.

We then use the fact that the two electronic subbands are parallels, and we neglect the photonic exchanged momentum when compared to an electronic one, that is

ω 2,k = ω 1,k + ω 12 ω µ,|k+q| ≃ ω µ,k . (16)
We can thus rewrite the Hamiltonian as 2,24

H = H C + j,k ω 12 b † j,q b j,q (17)
+ q Ω q (b 0,-q + b † 0,q )(a q + a † -q ) + N el 2 q Ṽ 1212 q b † 0,q + b 0,-q b † 0,-q + b 0,q ,
where H C is defined in Eq. ( 2) and the last line comes from the intersubband part of the Coulomb interaction only (we have omitted the constant ground state energy).

If the b 0,q operators are supposed to obey bosonic commutation rules

[b 0,q , b † 0,q ′ ] ≃ δ q,q ′ δ i,j , (18)
then the Hamiltonian in Eq. ( 17) is quadratic and can thus be solved by an Hopfield-Bogoliubov transformation 33 .

The different approximations that led to Eq. ( 17) (parabolicity, small exchanged momentum, perfect bosonicity, absence of intrasubband scattering) have been thorougly tested in multiple experiments 6,7,9,21,29,34 in the linear regime.

Still, these simplifying approximations oblige us to neglect a certain number of phenomena, like the interplay between nonparabolicity and Coulomb interaction 23 , or the scattering toward electron-hole pairs at large wavevectors, that is known to be an important factor in the thermalization and dynamics of exciton polaritons 35,36 . Moreover neglecting intrasubband Coulomb scattering does not allow us to consider polariton-polariton Coulomb interactions, that would ap-pear as an higher order term in the Hamiltonian of Eq. (17). Indeed, the bosonic approximation in Eq. ( 18) is satisfied only if its expectation value is taken between states corresponding to a weak excitation density. Otherwise, we will have a correction term in Eq. ( 18), roughly proportional to the ratio between the number of excitations, N exc , and of electrons, N el n QW , in such a way that this deviation can become important at higher excitation densities 22 .

To go beyond these limitations, we will derive an effective quartic bosonic Hamiltonian, taking care of both nonbosonicity and Coulomb scattering. In order to do so, we will have to calculate the matrix elements of the full Hamiltonian in Eq. ( 2) between states with multiple intersubband excitations. The algebra necessary to calculate the coefficients of such an effective Hamiltonian is rather heavy, and we will use the coboson formalism, originally developed for exciton polaritons 28,37 . The reader can refer to the next section for a brief summary of such a formalism, highlighting the aspects relevant for the algebra of the present paper.

III. COMPOSITE BOSON COMMUTATOR FORMALISM

Here we generalize the composite boson (coboson) approach, which was originally developed in Refs. [37] and [28] for the case of excitons. Such a formalism applied to intersubband excitations provides a systematic way to calculate many-body matrix elements. For simplicity, we consider the case where n QW = 1. The case with many quantum wells is more cumbersome, but relatively straightforward. Final results are obtained by replacing 1/N el by 1/(n QW N el ) in Eq. ( 27), ( 28) and (33).

We start by writing Hamiltonian H F in Eq. ( 2) in the appropriate form to use the composite boson commutator formalism. In particular it must annihilate the ground state. Our main interest will be to calculate matrix elements describing scattering events between intersubband excitations. We will thus limit ourselves to consider matrix elements of the form F | . . . b i,q Hb † i ′ ,q ′ . . . |F , where all intersubband excitation creation operators are on the right of the Hamiltonian and annihilation operators on the left. It is thus immediate to verify that terms in the Hamiltonian containing operators c † 1,k or c 1,k in normal order, with k outside the Fermi sea, do not contribute to such matrix elements, as they annihilate the state. We can thus remove them by truncating sums over k, relative to the first subband, to the Fermi sea (k ≤ k F ). We then rewrite H F in Eq. ( 2) in the following form (omitting the ground state energy and non-contributing terms)

H F = HKin + H Rabi + H Intra + H Depol , (19)
where

HKin = k ω2,k c † 2,k c 2,k - k ω1,k c 1,k c † 1,k , (20)
H Intra = 1 2 k,k ′ ,q Ṽ 2222 q c † 2,k+q c † 2,k ′ -q c 2,k ′ c 2,k + 1 2 k,k ′ ,q Ṽ 1111 q c 1,k+q c 1,k ′ -q c † 1,k ′ c † 1,k - k,k ′ ,q Ṽ 1221 q c 1,k+q c † 2,k ′ -q c 2,k ′ c † 1,k , (21)
H Depol = N el 2 q Ṽ 1212 q 2b † 0,q b 0,q +b † 0,q b † 0,-q + b 0,q b 0,-q , ( 22
)
all sums over the first subband are again truncated to the Fermi wavevector and

ω1,k = ω 1,k - |k ′ |<kF Ṽ 1111 |k-k ′ | (23) ω2,k = ω 2,k - |k ′ |<kF Ṽ 1212 |k-k ′ | .
In the previous equation we have thus explicitly separated the intersubband and intrasubband terms of the Coulomb interaction, putting them respectively into H Depol and H Intra . The first term will give the depolarization shift 24 , while the second one contains the electron-hole interaction which lowers the energy of bare intersubband excitations 27 and is also responsible for the scattering between intersubband excitations. Moreover, we commuted the electronic operators in such a way that now HKin and H Intra annihilate the ground state |F . Writing the Hamiltonian in such a form makes it also apparent the renormalization effect of the Coulomb interaction over the electronic dispersions, that is the same as the one that would be given by an Hartree-Fock approach 27,[START_REF] Giuliani | Quantum theory of the electron liquid[END_REF] .

The starting point of the composite boson formalism is the exact commutator between intersubband transitions defined in Eq. ( 14)

b m,q , b † i,q ′ = δ m,i δ q,q ′ -D mq,iq ′ , (24)
where D mq,iq ′ = D (1) + D (2) is the deviation from bosonicity operator defined by

D (1) = 1 N el k ν m,k ν * i,k+q-q ′ c 1,k+q-q ′ c † 1,k D (2) = 1 N el k ν m,k ν * i,k c † 2,k+q ′ c 2,k+q . (25)
To calculate many-body matrix elements involving intersubband excitations, we need to evaluate the following commutator

D mq,iq ′ , b † j,q ′′ = 1 N el n λ n,j m,i (q -q ′ ) (26) + λ n,i m,j (q -q ′′ ) b † n,q ′ +q ′′ -q , where λ n,j m,i (q) = 1 N el k ν m,k ν n,k+q ν * i,k+q ν * j,k . (27)
The term appearing in Eq. ( 26) is analogous to an effective scattering due to the nonbosonicity of the intersubband excitations (Pauli scattering): two intersubband excitations can exchange their electrons to create two excitations with different wavevectors, the total wavevector being conserved. With the identities above mentioned, it is possible to compute the scalar products of many-excitation states as

F | b n,q ′′′ b m,q b † i,q ′ b † j,q ′′ |F = δ mq,iq ′ δ nq ′′′ ,jq ′′ (28)
- 1 N el λ n,j m,i (q -q ′ ) + (m, q ↔ n, q ′′′ ).
From Eqs. ( 13) and (22) we see that matrix elements of H Rabi and H Depol can be written as scalar products of many-excitation states and can thus be calculated with the previous formula.

In order to evaluate the many-body matrix elements of the kinetic part of the Hamiltonian HKin , we need the commutator

HKin , b † i,q = ω12 b † i,q , (29)
where we have approximated, for simplicity, ω2,|k+q| with ω2,k . Many body matrix elements of HKin can thus be calculated by commuting it all the way to the right, leaving behind only terms in the form of scalar products between states with many intersubband excitations.

To complete the framework, we need to consider the matrix elements of H Intra . We start again by calculating its commutator with an intersubband excitation creation operator

H Intra , b † i,q = - j γ i,j b † j,q + V i,q , (30)
where γ i,j is a coupling between bright and dark excitations

γ i,j = q Ṽ 1221 q λ 0,0 j,i (q), (31)
and V i,q is called the creation potential operator,

V i,q = Q,m Ṽ 2222 Q δ m,i -Ṽ 1221 Q λ 0,0 m,i (Q) ( 32
)
×b † m,q+Q k c † 2,k-Q c 2,k + Q,m Ṽ 1111 Q λ 0,0 m,i (Q) -Ṽ 1221 Q δ m,i ×b † m,q+Q k c 1,k+Q c † 1,k .
The coefficient γ i,i in Eq. ( 30) is a renormalization of the intersubband transition energy due to the electron-hole interaction. Notice that the creation potential operator, V i,q , as well as HKin and H Intra , annihilate the ground state |F , a property which will be important when calculating matrix elements. The last commutator we need is the one between a creation potential operator and an intersubband excitation creation operator, that is

V i,q , b † j,q ′ = 1 N el e 2 n el 2ǫ 0 ǫ r (33) m,n,Q ξ n,j m,i (Q) Qǫ(Q) b † m,q+Q b † n,q ′ -Q , where ξ n,j m,i (Q) = δ m,i δ n,j I 2222 Q + λ 0,0 m,i (Q)λ n,j 0,0 (Q)I 1111 Q -δ m,i λ n,j 0,0 (Q)I 1221 Q -δ n,j λ 0,0 m,i (Q)I 1221 Q ,
is the so-called direct scattering.

In order to calculate an arbitrary matrix element concerning H Intra between many-excitations states, we can commute it all the way to the right, leaving behind terms that will have the form of a creation potential operator, V i,q , sandwiched between multiple intersubband excitation operators. We can then commute creation potential operators all the way to the right, leaving at the end only many-excitations states scalar products, that we can calculate using Eq. ( 28).

IV. EFFECTIVE BOSONIC HAMILTONIAN

Our aim is to derive a bosonic Hamiltonian capable of describing intersubband polariton-polariton interactions. We will thus introduce the bosonic ground state |G and B q , the annihilation operator for a bright bosonic intersubband excitation, with in-plane wavevector q satisfying the usual bosonic commutation rules

B q , B † q ′ = δ q,q ′ . ( 34
)
The index has been omitted as we are only interested in the dynamics of bright excitations. We will then calculate the coefficients of the bosonic Hamiltonian H B , imposing that it has the same expectation values that H F , for states with up to two excitations. While the most general Hamiltonian, up to n-body interactions, can be calculated with the same procedure, here we will limit ourselves to consider only resonant two-body interactions involving bright intersubband excitations. The actual calculation of the coefficients requires some rather complex operatorial algebra. For sake of simplicity, we have thus moved the bulk of the calculations in the Appendices, presenting here only the final form of the Hamiltonian and the parameter dependency of its coefficients. In Sec. III, we have provided the essential ingredients of the composite boson commutator approach 22,28,37 , that we then use extensively in Appendix B to calculate the actual coefficients of the bosonic Hamiltonian.

The bosonic Hamiltonian, describing the bright excitations, and using the quantities defined in Sec. III and in Appendix B, has thus the form

H B = q ( ω12 -γ 0,0 ) B † q B q +
q e 2 n el 4ǫ 0 ǫ r ǫ 12 (0)

I 1212 q q q=0 2B † q B q + (1 -ζ) B q B -q + B † q B † -q + q Ω q (B q + B † -q )(a -q + a † q ) - 1 n QW N el q,q ′ ,p Ω q 2 [1 + (1 -δ p,0 ) (1 -δ p,q ′ -q )] B † q+p B † q ′ -p B q ′ a q + h.c. - 1 n QW N el 1 2 q,q ′ ,p U B † q+p B † q ′ -p B q ′ B q , (35)
where the last nonlinear coefficient, describing the scattering of pairs of intersubband excitations, is U = e 2 n el ǫ 0 ǫ r ǫ 12 (0)

I 1212 q q q=0 (1 -ζ) - e 2 n el 2ǫ 0 ǫ r κ ξ 0,0 0,0 (0) + e 2 2ǫ 0 ǫ r n el 2π x(0). (36)
In the previous expressions n el denotes the electron density in each quantum well. A comparison with Eq. (17) shows that taking into account the nonbosonicity of intersubband excitations and the intrasubband terms of the Coulomb interaction resulted in the appearance of quartic terms in the Hamiltonian, that makes it possible to observe nonlinear polariton-polariton processes. In the first and fourth line in Eq. ( 35) and in Eq. (36) coefficients have no dependence over the wavevectors of the intersubband excitations. This is justified by the fact that these coefficients evolve significantly on a scale of the order of the Fermi wavevector k F , which is much larger than the photonic wavevectors of interest. There is also no spin dependence as one should expect from a naive comparison with excitons 38,39 . As explained in Sec. II C and Appendix A this is due to the fact that all intersubband of interest are made of two fermions with the same spin.

The first quartic term in Eq. ( 35) describes the saturation of the light-matter interaction. Its physical origin lies in the non-perfect bosonicity of the intersubband excitations and accounts for a renormalization of the Rabi frequency at high excitation densities. Its presence is easy to understand if one remembers 40 that the lightmatter interaction is proportional to the square root of the total electronic population in the first subband minus the one in the second subband. As the creation of one intersubband excitation corresponds, in average, to the promotion of one electron from the first to the second subband, we do expect that an increase in the intersubband excitations population will reduce the light-matter coupling, consistently with the minus sign in front of such a coefficient.

The second quartic term describes instead the scattering between pairs of intersubband excitations. It has three contributions as seen in Eq. ( 36). The first one comes from the intersubband Coulomb interaction. The second and third contributions instead are due to the intrasubband Coulomb interaction. It is important to notice that the two nonlinear coefficients behave differently with respect to the number of quantum wells in the cavity. While Ω q is proportional to the square root of the total number of electrons in the cavity, U depends only on the number of electrons in each well. This is due the fact that the photonic mode couples to all electrons in the cavity while the Coulomb interaction does not couple one well with another.

As we can see, Coulomb interaction and the nonperfect bosonicity also give a renormalization effect visible in the quadratic part of the Hamiltonian. The intersubband energy is renormalized both by Hartree-Fock terms and by the Coulomb electron-hole interaction. Moreover, the non-resonant terms of the depolarization shift are renormalized by a ζ < 1 coefficient as a consequence of the nonbosonicity. Again, the expressions of the coefficients can be found in Appendix B.

The previous Hamiltonian is adapted to describe optically pumped experiments, in which the populations of the dark states remain negligible. While we will not treat this case in the present work, we notice that, in the case of electrical pumping, the dark states population can become significant. Interaction between bright and dark modes is described by the Hamiltonian

H D = -2 i =0,q,q ′ ,p G i q-q ′ ,p n QW N el B † q-p B † i,q ′ -p B i,q ′ a q + h.c. - i =0,q,q ′ ,p U i n QW N el B † q+p B † i,q ′ -p B i,q ′ B q , ( 37
)
where the index i means that the mode is dark. In a mean field approach this Hamiltonian can be seen to describe energy shifts and saturation effects on the bright modes as a function of the dark state population. Coefficients in Eq. ( 37) can be calculated from a generalization of the formula given in Sec. III and in Appendix B.

Passing from the fermionic description H F to the bosonic one H B , we have kept only quartic terms in H B , discarding all higher order terms. This is effec-tively equivalent to limit ourselves to the lowest order in the excitation fraction N exc /(n QW N el ) 28 . We thus expect that our theory will correctly describe the polaritonpolariton interactions in the regime of moderate excitation densities 41 . In order to clarify this point we give here an explicit example not involving Coulomb interaction. The reader can refer to Appendix C for more complex and detailed examples.

We consider a process in which a state with one photon and N exc intersubband exitations is converted into a state with no photon and N exc + 1 intersubband excitations. The matrix element of the Coulomb part of the Hamiltonian vanishes in this case, as it conserves the number of photons, and we can thus perform the calculation exactly, both in the fermionic and in the bosonic formalism. We obtain respectively

F | a q b Nexc 0,q H F b † Nexc+1 0,q |F F | b Nexc 0,q b † Nexc 0,q |F F | b Nexc+1 0,q b † Nexc+1 0,q |F = Ω q N exc + 1 1 - N exc n QW N el = Ω q N exc + 1 1 - N exc 2n QW N el + O N 2 exc n 2 QW N 2 el ( 38
) G| a q B Nexc q H B B † Nexc+1 q |G (N exc + 1)!N exc ! = Ω q N exc + 1 1 - N exc 2n QW N el ,
which are indeed equal to first order in N exc /(n QW N el ).

V. NUMERICAL RESULTS

In this Section we provide the numerical values of the relevant coefficients that appear in the bosonic Hamiltonian in Eq. (35), highlighting the dependence of the main parameters.

A. Renormalization of the intersubband transition

We start by analyzing the renormalization effect due to the Coulomb interaction on the quadratic part of the effective Hamiltonian. Performing a Bogoliubov transformation on the first line of Eq. ( 35) we obtain the renormalized energy of the intersubband excitation (that we should call, more properly, intersubband plasmon [24][25][26][27] )

E ISBT = ω12 -γ 0,0 + e 2 n el 2ǫ 0 ǫ r ǫ 12 (0) I 1212 q q q=0 2 - e 2 n el 2ǫ 0 ǫ r ǫ 12 (0) I 1212 q q q=0 (1 -ζ) 2 . (39)
In Fig.

(3) we plot the dispersion of the intersubband transition energy considering a GaAs quantum well of length L QW = 11 nm (left panel) and L QW = 39 nm (right panel), corresponding to bare transitions of 140 meV 1,7 and 15 meV 25 respectively. The solid line depicts the result for the intersubband transition including all contributions in Eq. ( 35). The dashed line represents the result obtained while neglecting the Coulomb intrasubband interactions and the nonbosonicity. The dash-dotted line is the result obtained as the dashed line, but neglecting also the screening of the intersubband Coulomb interaction. Notice that the renormalized 36) corresponding to the intersubband Coulomb interaction. Thin red line: the absolute value of the second term in Eq. (36), namely the direct Coulomb interaction. Note that this term is negative, thus producing a red-shifted contribution. Dash-dotted line: the third term in Eq. (36), which is due to the exchange Coulomb interaction.

intersubband energy, E ISBT , in Eq. ( 39), converges to the bare transition energy ω 12 for vanishing doping. We see that the main effect is due to the screening of the intersubband Coulomb interaction. It appears clear that these renormalization effects are relevant in the THz regime (right panel), but rather moderate for mid-infrared transitions (left panel).

B. Interaction energy between intersubband excitations

Now, we consider the quartic part of the Hamiltonian H B , responsible for the Coulomb scattering of pairs of intersubband excitations. The interaction energy is given by N exc U/(n QW N el ). In Fig. 4 we plot the energy U (thick solid line) for a mid-infrared transition (left panel) and a THz transition (right panel). The other lines depict the individual contributions of the three terms in Eq. (36) (see caption for details). For the considered realistic parameters, the interaction energy grows with increasing electron doping density almost linearly for the mid-infrared case, while it saturates at high densities for the case of a THz transitions. Notice that in the case of excitons 28,38 , the direct scattering vanishes in the limit of small wavevectors, whereas it is not the case for intersubband excitations when the screened Coulomb interactions are considered. This is due the fact that electron-electron and hole-hole interactions are screened differently by the Fermi sea [START_REF] Lee | [END_REF] . Note that for a doping density in the range of a few 10 11 cm -2 , the energy U is of the order of a few meV both for the cases of THz and mid-infrared transitions. The interaction energy can thus reach values close to the meV when the density of excitations in the system becomes significant. This is rather promising, since as shown in the case of exciton-polaritons 42 , very interesting nonlinear polariton physics occurs when the interaction energy becomes comparable to the linewidth of the polariton modes. For THz polaritons, state-of-the-art samples 25 exhibits polariton width in the meV range.

VI. POLARITON HAMILTONIAN

In this Section we consider the interactions in the polariton basis. For simplicity we limit ourselves to the first mode of the cavity and we assume that the number of electrons is such that the system is in the strong coupling regime. The diamagnetic interaction in H C as well as all non-resonant terms in H B are thus neglected. A Bogoliubov transformation of the quadratic parts of H C and H B gives the expression of the polaritonic operators

p L q p U q = β q -α q α q β q B q a q (40)
where p U q and p L q are polaritonic operators of the upper (UP) and lower branch (LP) respectively, while α q and β q are the Hopfield coefficients. Note that results can be generalized to the ultrastrong coupling regime (vacuum Rabi frequency comparable to intersubband transitions) by considering the Bogoliubov transformation involving the antiresonant light-matter interaction terms 2 . Here, for simplicity, we consider the case of moderate values of the vacuum Rabi frequency. Because of the quartic terms in Hamiltonian H B polaritons interact with each other through their matter part and can scatter. In the following, we will focus on the lower branch. Replacing operators a q and B q by their expression, the two-body interaction between LP polaritons is found to be

H LP-LP = 1 2 q,q ′ ,p V q,q ′ ,p n QW N el p † L q+p p † L q ′ -p p L q ′ p L q ,( 41
)
where V q,q ′ ,p = β q+p β q ′ -p β q ′ (α q Ω q -β q U ) . (42)

As an example, we now consider the case where the lower branch is pumped at a wavevector q p so that the system is in the state p † Nexc L qp |G / √ N exc !. This Hamiltonian allows us to describe single-mode (Kerr) and multimode (parametric) coherent non-linearities 42 . For the parametric case, one has to consider interaction channel

p † Nexc L qp |G → p † L qp+p p † L qp-p p † Nexc-2 L qp |G , ( 43
)
where pairs of polaritons scatter from the pumped mode into signal-idler pairs. As for the case of excitonpolaritons we expect that the maximum efficiency of this parametric processes is achieved when the energy conservation condition is fulfilled 43 . A mean-field approach of the problem 42,44 shows that the matrix element between the initial and the final states is the relevant quantity to consider and has to be compared with the lifetime of the excitations. For high pump intensity i.e. N exc ≫ 1 this matrix element is

M qp,p = N exc n QW N el V qp,qp,p . (44)
As discussed in Sec. V B and shown in Fig. 4, for THz polaritons nonlinear interaction energies of the order of a few meV (thus comparable to THz polariton linewidths) can be achieved, thus paving the way to a very interesting coherent nonlinear physics for this kind of composite excitations.

VII. CONCLUSIONS

In this paper, we have presented comprehensively a microscopic theory for the manybody physics of intersub-band cavity polaritons. Using a composite boson commutator approach, we have derived the manybody interaction matrix elements involving (screened) Coulomb interaction processes and Pauli saturation effects. We have also derived an effective bosonic Hamiltonian reproducing the matrix elements in the two-excitation manyfold calculated by the composite boson approach. As a first application, we have calculated the renormalization of the intersubband transition frequency at the mean-field level and determined the strength of the interactions between intersubband excitations. We have also derived the renormalization of the light-matter coupling due to the Pauli saturation effects and considered the interaction physics in the polariton basis. Our results show that significant intersubband polariton-polariton interactions occur, especially for transitions in the THz regime. Our work paves the way to promising future studies of nonlinear quantum optics with intersubband cavity polaritons. and

b † 0,q = 1 n QW N el nQW j=1 k,σ ν * 0,j,k c (j) † 2,k+q,σ c (j) 1,k,σ ,(A2)
where ν 0,j,k = Θ(k F -k) for all j, k is a wavevector such that k < k F , σ, σ ′ ∈ {↓, ↑} and µ, µ ′ , ν, ν ′ ∈ {1, 2}. Bright intersubband excitations are thus linear superposition of pairs of fermions with the same spin. We could generalize Eq. (A2) by allowing the two fermions to have different spins but the resulting collective excitation would be dark and thus not relevant if we consider only polariton.

Coefficient d2

Only the light-matter part of H F contributes and we obtain d 2 = Ω q .

(B8)

Coefficient d3

This term comes from the intersubband Coulomb interaction and is given by

d 3 = e 2 n el 2ǫ 0 ǫ r ǫ 12 (0) I 1212 q q q=0 (1 -ζ) , (B9)
where

I 1212 q q q=0 × ζ = 1 N el Q I 1212 Q Q λ 0,0 0,0 (Q). (B10)
The coefficient Q is then d 3 (1 + 1/n QW N el ). As the number of electrons is large we can neglect the correction due to the normalization.

Coefficient d4

As for the coefficient d 2 , only the light-matter part of H F contributes, giving

d 4 = Ω q δ p,0 + δ p,q ′ -q - 2 n QW N el . ( B11
)
The last term, at the origin of the first quartic term in Eq. ( 35), is due to the nonbosonicity.

Coefficient d5

The expression of the coefficient is given by the following relations (to first order in 1/(n QW N el ))

d 5 = F | b 0,q+p b 0,q ′ -p -F | b 0,q+p b 0,q ′ -p b † 0,q ′ b † 0,q |F F | b 0,q b 0,q ′ H F b † 0,q ′ b † 0,q |F , (B12)
F | b 0,q+p b 0,q ′ -p H F b † 0,q ′ b † 0,q |F = 2 ( ω12 -γ 0,0 ) δ p,0 + δ p,q ′ -q - 2 n QW N el + e 2 n el ǫ 0 ǫ r ǫ 12 I 1212 q q q=0 δ p,0 + δ p,q ′ -q - 4 -2ζ n QW N el + 2 n QW N el e 2 n el 2ǫ 0 ǫ r κ ξ 0,0 0,0 (0) - e 2 2ǫ 0 ǫ r n el 2π x(0) , (B13)
where

x(p) = m,n,Q k F N el Qǫ(Q)
λ 0,n 0,m (Qp) ξ n,0 m,0 (-Q). (B14)

In Eq. (B13) the long wavelength limit has been taken because the variation of the coefficients is significant only on a scale of the order of the Fermi wavevector k F .

and final states to be proportional respectively to the kets

|ψ i ∝ b † Nexc 0,q |F (C1)
and

|ψ f ∝ b † 0,q+p b † 0,q-p b † Nexc-2 0,q |F . (C2)
The transition rate can then be written as where δ t converges to the Dirac delta function for long times. The orthogonal projection in Eq. (C4) can be calculated as

P ⊥ Hb † Nexc 0,q |F = -N exc ( ω12 -γ 0,0 ) i =0 γ 0,i b † i,q b † Nexc-1 0,q |F + N exc (N exc -1) 2n QW N el U b † Nexc 0,q |F + N exc (N exc -1) 2n QW N el e 2 n el 2ǫ 0 ǫ r m,n,Q ξ m,0 n,0 (Q) Qǫ(Q) b † m,q+Q b † n,q-Q b † Nexc-2 0,q |F - N exc (N exc -1) 2n QW N el e 2 n el ǫ 0 ǫ r ǫ 12 m,Q λ m,0 0,0 (Q -q)b † 0,Q b † m,2q-Q b † Nexc-2 0,q |F . (C6)
The previous expression is a first order in the pertubations. We can thus consider only the zero-th order in The second term can be calculated using the same method as for the transition rate and is found to contribute only to higher orders in the perturbation. After some rather cumbersome algebra, the first term is found to be of the form 1

T = 1 2 p Γ p , ( C10
)
where it is implicitly assumed that the summation is restricted to small wavevectors. To obtain this result we made F t act at zero-th order and we used the fact that bright and dark intersubband excitations are not resonant. The counterintuitive 1/2 factor comes from the overcompleteness of the composite boson basis.

Bosonic case

In this paragraph we calculate the same quantities as in the previous one but we assume that intersubband excitations are bosons and that their dynamics is described by Hamiltonian H B . We can therefore use a traditional Fermi golden rule.

The initial and final states in this case are

|ψ i ∝ B † Nexc q |G (C11)
and

|ψ f ∝ B † q+p B † q-p B † Nexc-2 q |G .
The transition rate is thus given by the formula

Γ p = 2π N exc (N exc -1) n 2 QW N 2 el U 2 δ t (0), (C12)
which is the same as in Eq. (C8) up to third order in N exc /(n QW N el ). As all bright intersubband excitations are resonant the lifetime of the initial state is the sum of the transition rates

1 T = p Γ p , (C13)
where the summation is again restricted to small wavevectors. A comparison with Eq. (C10) shows that this method underestimates the true lifetime by a factor two. This is coherent with the results in Ref. [37], that show how an effective Hamiltonian approach that correctly calculates transition rates needs to take into account an ad hoc factor 1/2 when calculating lifetimes, due to the overcompleteness of the composite boson basis. This is of course simply implies a renormalization of the composite boson density of states.

FIG. 1 :

 1 FIG. 1: Schematic representation of the system under investigation. Panel (a): Multiple doped quantum wells (QWs) are embedded in a planar photonic microcavity. The image above shows also the fundamental TM photonic mode in the case of metallic mirrors. Due to intersubband selection rule, the intersubband transitions only couple to the z-component of the electric field. Panel (b): Energy profile of the quantum well along the z-direction and its two bound states (left). Electronic dispersion of the first two conduction subbands in each of the doped quantum wells, as a function of the inplane electronic wavevector (right). The Fermi wavevector, kF , the Fermi Energy, EF , and the intersubband gap, ω12, are highlighted.

nFIG. 3 :nFIG. 4 :

 34 FIG. 3: Renormalized energy EISBT of the intersubband transition for a bare energy ω12 = 140meV (left panel) and 15meV (right panel), as a function of the electron density. Solid line: result including all the contributions in Eq. (39), namely Coulomb interaction and nonbosonicity. Dashed line: result neglecting intrasubband Coulomb interactions and the non-bosonicity. Dash-dotted line: same as the dashed line, but with no screening of the intersubband Coulomb interaction.

F= ψt ψt -ψ i ψt 2 .

 2 t , i.e we do as if intersubband excitations were exact eigenstates of H F , obtainingF | b Nexc-2 0,q b 0,q-p b 0,q+p F t (H F -ψ i | H F |ψ i ) ≈ F t (0) F | b Nexc-2 0,q b 0,q-p b 0,q+p . (C7)Taking into account the normalization we thus obtain the transition rateΓ p = 2π N exc (N exc -1)With the same notations, the contribution of the many-body physics to the lifetime of the initial state is 37 t T

  37 tΓp = ψ f ψt F t (H F -ψ i | H F |ψ i )P ⊥ H F |ψ i . (C4)As in Appendix B, P ⊥ is the projector over the subspace orthogonal to |ψ i and F t verifies|F t (E)| 2 = 2πt δ t (E), (C5)

	2	,	(C3)
	where
	ψt =

Acknowledgments

We acknowledge partial support from the french ANR project THINQE-PINQE. MB is also supported by KAK-ENHI 24-632. We wish to thank C. Sirtori, Y. Todorov, A. Vasanelli for discussions.

Appendix A: Spin and quantum well number

In this appendix we give the explicit notations for Eq.

(3) and Eq. (12) with quantum well and spin indexes. Electronic wavefunctions are localized in quantum wells and we neglect electronic tunneling from one well to another. We thus neglect Coulomb interaction between electrons in different wells, which are sufficiently apart. Eq. ( 3) and (12) can then be written with all indexes

In order to determine the coefficients of the effective bosonic Hamiltonian in Eq. (35), we impose that it has the same matrix elements as the fermionic one in Eq. ( 2) over the Hilbert space spanned by one and two excitation states. That is we impose that

where T q ∈ [I, a q , B q ], t q ∈ [I, a q , b 0,q ], P ⊥ is the projector over the subspace orthogonal to the ket 41 and N q,q ′ and A q,q ′ are the normalizations of the orthogonalized states.

The constraints in Eq. (B1), together with the fact that we decided to consider, beyond quadratic terms, only the quartic ones, are enough to uniquely define H B . The calculation of the coefficients of H B thus reduces to the calculation of matrix elements of the fermionic Hamiltonian H F , that we can perform exploiting the commutator approach outlined in Sec. III.

We thus start from a general bosonic Hamiltonian

We also introduce the following five matrix elements (each one implicitly dependent on all the relevant wavevectors)

The relations between the coefficients of H B in Eq. (B2) and the matrix elements in Eq. (B3) can be found with some algebra to be

where, for consistency, terms under the square root have also to be developed to the first order in 1/(n QW N el ).

We have thus reduced the determination of the coefficients of H B to the calculation of the five matrix elements in Eq. (B3). Without loss of generality, in the following calculations, we will assume that the excitations wavevectors are small compared to the Fermi wavevector, which implies that the Pauli scattering involving two bright excitations, defined in Eq. ( 27), can be approximated by a Kroneker delta

The only contributing terms are the kinetic and Coulomb part of the Hamiltonian H F . The matrix element is

Appendix C: Transition rates and lifetimes

In this Appendix we calculate the transitions rates for scattering of pairs of intersubband excitations using both the effective bosonic Hamiltonian approach and the fermionic formalism developed by M. Combescot and coworkers 37 , showing that the two approaches give exactly the same results

Fermionic case

In this paragraph we use the composite boson commutator formalism to calculate the transition rate between an initial state |ψ i and a final state |ψ f , as well as the lifetime of the former, using a first-order, time-dependent perturbation theory. The calculation is equivalent to a Fermi golden rule as explained in Ref. [37].

The particular event we want to describe is the scattering of an initial pump beam of arbitrary intensity into a signal and an idler mode. We will thus consider initial


Ï C Lo

Le Treust
email: letreust@ceremade.dauphine.fr

Simona Rota Nodari
email: simona.rota-nodari@u-cergy.fr

SYMMETRIC EXCITED STATES FOR A MEAN-FIELD MODEL FOR A NUCLEON

In this paper, we consider a stationary model for a nucleon interacting with the ω and σ mesons in the atomic nucleus. The model is relativistic, and we study it in a nuclear physics nonrelativistic limit. By a shooting method, we prove the existence of infinitely many solutions with a given angular momentum. These solutions are ordered by the number of nodes of each component.

Introduction

This article is concerned with the existence of excited states for a stationary relativistic mean-field model for atomic nuclei in the nuclear physics nonrelativistic limit. To our knowledge, this model was first studied by Esteban and Rota Nodari; in two recent papers [START_REF] Esteban | Ground States for a Stationary Mean-Field Model for a Nucleon[END_REF][START_REF] Esteban | Symmetric ground states for a stationary relativistic mean-field model for nucleons in the non relativistic limit[END_REF], the authors showed the existence of so-called ground states (see [START_REF] Esteban | Ground States for a Stationary Mean-Field Model for a Nucleon[END_REF] for more details about the definition of ground states).

As the authors formally derived in [START_REF] Esteban | Symmetric ground states for a stationary relativistic mean-field model for nucleons in the non relativistic limit[END_REF], the equations of the model are given, in the case of a single nucleon, by (1.1) iσ • ∇χ + |χ| 2 ϕ -a|ϕ| 2 ϕ + bϕ = 0,

-iσ • ∇ϕ + 1 -|ϕ| 2 χ = 0,
with a and b two positive parameters linked to the coupling constants and the nucleon's and mesons' masses. This system is the nuclear physics nonrelativistic limit of the σ-ω relativistic mean-field model ( [START_REF] Walecka | A theory of highly condensed matter[END_REF][START_REF] Walecka | Theoretical Nuclear and Subnuclear Physics[END_REF]) in the case of a single nucleon. We remind that σ is the vector of Pauli matrices (σ 1 , σ 2 , σ 3 ), and ϕ, χ : R 3 → C 2 .

As in [START_REF] Esteban | Symmetric ground states for a stationary relativistic mean-field model for nucleons in the non relativistic limit[END_REF], we look for solutions of (1.1) in the particular form

(1.2) ϕ(x) χ(x) =       g(r) 1 0 if (r) cos ϑ sin ϑe iφ      
where f and g are real valued radial functions and (r, ϑ, φ) are the spherical coordinates of x. The system (1.1) then turns to a nonautonomous planar differential system which is

(1.3)    f + 2 r f = g(f 2 -ag 2 + b) , g = f (1 -g 2 ) .
In order to avoid solutions with singularities at the origin, we impose f (0) = 0, and, since we are interested in finite energy solutions of (1.1), we seek solutions of (1.3) that are localized i.e. solutions which fulfill (1.4) (f (r), g(r)) -→ (0, 0) as r -→ +∞ .

In [5, Proposition 2.1], Esteban and Rota Nodari showed that there is no nontrivial solution of (1.3) such that (1.4) is satisfied unless a -2b > 0. Hence, in what follows, we assume a -2b > 0.

For every given x, there exists a local solution (f x , g x ) of (1.5)

       f + 2 r f = g(f 2 -ag 2 + b) , g = f (1 -g 2 ) , f (0) = 0, g(0) = x.
The problem is to find x, such that the corresponding solution is global (i.e. defined for all r ≥ 0), and satisfies (1.4).

In [5, Proposition 2.1], Esteban and Rota Nodari proved that if (f x , g x ) is a solution of (1.5) satisfying (1.4) then g 2

x (r) < 1, for all r in [0, +∞). So, in particular, x = g x (0) must be chosen such that x 2 < 1. This creates additional difficulties to deal with.

Since the system of equations (1.3) is symmetric with respect to 0, we study the problem (1.5) with x ∈ [0, 1). Moreover, let us remark that if x = 0 then (f x , g x )(r) = (0, 0) for all r ≥ 0 is the unique solution of (1.5).

In [START_REF] Esteban | Symmetric ground states for a stationary relativistic mean-field model for nucleons in the non relativistic limit[END_REF], the authors proved the existence of a global localized solution (f x , g x ) of (1.5) such that f x (r) < 0 < g x (r) for all r ∈ (0, +∞). In this paper, we generalize this results by showing the existence of global localized solutions with any given number of nodes. Our main result is the following. Theorem 1.1. Assume a > 2b > 0. There exists an increasing sequence {x k } k≥0 in (0, 1) with the following properties. For every k ≥ 0,

(1) the solution (f x k , g x k ) of (1.5) is a global solution;

(2) both f x k and g x k have exactly k zeros on (0, +∞);

(3) (f x k , g x k ) converges exponentially to (0, 0) as r → +∞.

This theorem is the first result of existence of excited state solutions for the model studied in [START_REF] Esteban | Symmetric ground states for a stationary relativistic mean-field model for nucleons in the non relativistic limit[END_REF][START_REF] Esteban | Ground States for a Stationary Mean-Field Model for a Nucleon[END_REF] for which Esteban and Rota Nodari proved the existence of a ground state solution.

Our theorem is similar to the result obtained by Balabane, Cazenave, Douady and Merle ( [START_REF] Balabane | Existence of excited states for a nonlinear Dirac field[END_REF]) for a nonlinear Dirac equation. Our proof is based on a shooting method inspired by the one used by Balabane, Dolbeault and Ounaies ( [START_REF] Balabane | Nodal solutions for a sublinear elliptic equation[END_REF]).

In [START_REF] Balabane | Existence of excited states for a nonlinear Dirac field[END_REF], the authors proved the existence of infinitely many stationary states for a nonlinear Dirac equation. More precisely, they showed the existence of a bounded increasing sequence of positive initial data {x k } k such that the associated solutions are global and each component has k nodes.

In [START_REF] Balabane | Nodal solutions for a sublinear elliptic equation[END_REF], thanks to some estimations on the energy decay and the rotation speed, the authors proved the existence of infinitely many solutions for a sublinear elliptic equation. As in [START_REF] Balabane | Existence of excited states for a nonlinear Dirac field[END_REF], they showed the existence of an increasing sequence of initial data {x k } k such that the associated solutions are radial, compactly supported and have exactly k nodes.

As we remarked above, the first difficulty to deal with here is that, to obtain a localized solution, the initial condition x must be chosen in (0, 1). Moreover, we are looking for solutions such that each component has exactly k zeros on (0, +∞).

Usually in a shooting method, the localized solution with k nodes is obtained taking the solution whose initial data x is the supremum of a well-chosen open subset of {x : g x has k zeros}. Hence, the main difficulty of our shooting method is to prove that for any k ∈ N, there exists ε > 0 such that {x ∈ (0, 1) : g x has k zeros} ⊂ (0, 1 -ε).

To do this, we have to give some accurate estimations on the behavior of the solution when the initial condition x becomes close to 1. The presence of four rest points (± √ a -b, ±1) in the Hamiltonian system

(1.6) f = g(f 2 -ag 2 + b) g = f (1 -g 2 ) ,
associated with the system (1.3), makes this study difficult. Indeed, we would like to control the solutions (f x , g x ) thanks to the continuity of the flow comparing (f x , g x ) to (f 1 , g 1 ) whenever x is close enough to 1. The problem is that (f 1 , g 1 ) tends to the rest point (-√ a -b, 1) of the system (1.6). Thus, (f x , g x ) stay in a neighborhood of (-√ a -b, 1) a very long time if x is sufficiently close to 1. Since (f 1 , g 1 ) does not wind around (0, 0), it is hopeless to get estimations on the speed of rotations of (f x , g x ) around (0, 0) as in [START_REF] Balabane | Nodal solutions for a sublinear elliptic equation[END_REF]. Hence, we introduce another strategy to prove that (f x , g x ) winds around (0, 0).

First of all, we prove that (f x , g x ) exits the neighborhoods of (-√ a -b, 1) at finite time, possibly very large. Next, we want to control the position of (f x , g x ) when this occurs. To do this, we introduce the so-called Hamiltonian regularization. More precisely, we replace the system (1.3) by the Hamiltonian ones (1.6) in a neighborhood of the points (± √ a -b, ±1) (see Figure 1). Then, we can use the qualitative properties of the solutions of the Hamiltonian system (1.6) to know the position of the solution when it exits the neighborhood of (-√ a -b, 1). Finally, we iterate the reasoning to prove that if x is sufficiently close to 1, then g x has more than k zeros.

The idea of the Hamiltonian regularization is inspired by the proof of Le Treust in [START_REF] Le Treust | Existence of nodal solutions for Dirac equations with singular nonlinearities[END_REF]. In this paper, the author proved the existence of infinitely many compactly supported nodal solutions of a Dirac equation with singular nonlinearity. The main problem encountered is that the nonlinearity is singular and the main theorems of ODE fail to show local existence and uniqueness. To overcome this, Le Treust used a regularization by a Hamiltonian system whenever the problems occur. The advantage of such a regularization is that it gives a better control of the regularized solutions while keeping true some qualitative properties of the solutions of the nonautonomous system of equation.

In section 2, we introduce the regularized system and we prove the existence of nodal localized solutions of the regularized problem assuming some key lemmas. In the next section, we prove these lemmas. In section 4, we show that the localized nodal solutions of the original system (1.3) can be obtained as limits of nodal localized solutions of the regularized system. Finally, in the appendix, we give some useful properties of the Hamiltonian energy associated to the system.

The regularized problem and the shooting method

Hamiltonian System H ≤ 0 f g X 1 X 2 X 3 X 4 f = ± √ a -b -η 2 f = ± √ a -b -η ∂A √ a -b - √ a -b 1 -1 Figure 1. Regularized System 2.1. Construction of the regularized problem. Let ϕ : (η, f, g) ∈ 0, √ a -b - a 2 × [- √ a -b, √ a -b] × R → ϕ η (f, g) ∈ [0, 1] be a smooth function on 0, √ a -b - a 2 × (- √ a -b, √ a -b) × R such that for all (f, g) ∈ R 2 , all η ∈ 0, √ a -b -a 2 ϕ η (f, g) = ϕ η (|f |, |g|) = 0 if |f | ≥ √ a -b -η/2 1 if |f | ≤ √ a -b -η.
Consider the system of equations (2.1)

   f + 2ϕ η (f, g) r f = g(f 2 -ag 2 + b), g = f (1 -g 2 ),
and the Cauchy problem

(2.2)          f + 2ϕ η (f, g) r f = g(f 2 -ag 2 + b), g = f (1 -g 2 ), f (0) = 0, g(0) = x.
We denote by (f x,η , g x,η ) the solutions to problem (2.2).

Remark 2.1. When η > 0, in the neighborhood of the four points (± √ a -b, ±1), the system of equations (2.1) becomes the following autonomous one

(1.6) f = g(f 2 -ag 2 + b), g = f (1 -g 2 ).
This system is a Hamiltonian system associated with the energy

(2.3) H(f, g) = 1 2 f 2 (1 -g 2 ) + a 4 g 4 - b 2 g 2 .
Remark 2.2. The behavior of the solutions of (1.6) is easier to understand than the one of the solutions of (1.3). This is actually the reason why we introduce such a Hamiltonian regularization in the neighborhood of the saddle points (± √ a -b, ±1) of H.

2.2.

Properties of the regularized system. We fix η ∈ (0, √ a -b -a 2 ). We begin by studying the existence and the uniqueness of the solutions of (2.1). Lemma 2.3. Let x ∈ R. For any a, b > 0, there is τ η > 0 and (f x,η , g x,η ) ∈ C 1 [0, τ η ], R 2 unique solution of (2.1) satisfying f x,η (0) = 0, g x,η (0) = x. Moreover, (f x,η , g x,η ) can be extended on a maximal interval [0, R x,η ) with either R x,η = +∞ or R x,η < +∞ and lim r→Rx,η |f x,η | + |g x,η | = +∞. Furthermore, (f x,η , g x,η ) depends continuously on x and η, uniformly on [0, R] for any R < R x,η .

Proof. As in [START_REF] Cazenave | Existence of Localized Solutions for a Classical Nonlinear Dirac Field[END_REF], it is enough to write

f (r) = 1 r 2 r 0 s 2 g(s)(f 2 (s) -ag 2 (s) + b) + 2s(1 -ϕ η (f (s), g(s)))f (s) ds , g(r) = x + r 0 f (s)(1 -g 2 (s)) ds ,
and note that the right hand side of (2.1) is a Lipschitz continuous function of (f, g). The lemma follows from a classical contraction mapping argument.

Next, we define

A = {(f 0 , g 0 ) ∈ R 2 | 2f 2 0 -ag 2 0 -(a -2b) ≤ 0, g 2 0 ≤ 1}, the set of admissible points. Remark. If g 2 0 ≤ 1, (f 0 , g 0 ) ∈ A if and only if H(f 0 , g 0 ) ≤ H(0, 1) = a-2b 4
(see Figure 1).

A key property of the solutions of the system (2.1) is the behavior of the energy H along the trajectories of the solutions as stated in the following lemma.

d dr H(f x,η , g x,η )(r) = - 2ϕ η (f x,η , g x,η ) r f 2 x,η (1 -g 2 x,η ).
Remark 2.5. Let us remark that this property which is true for the non-regularized system (1.3) is also true for the regularized one thanks to our choice of regularization.

Proof. The proof is a straightforward calculation.

Next, we prove a result that ensures that for all x ∈ (0, 1) the solutions (f x,η , g x,η ) are global and live in A.

Lemma 2.6. Let a, b > 0 such that a -2b > 0 and let (f x,η , g x,η ) be the solution of (2.1) satisfying f x,η (0) = 0, g x,η (0

) = x. If x 2 ≤ 1, then g 2 x,η (r) ≤ 1 and f 2 x,η (r) ≤ a -b for all r ∈ [0, R x,η ) and R x,η = +∞. Moreover, (f x,η , g x,η )(r) ∈ A, for all r ∈ [0, +∞) and if x 2 < 1, then (f x,η , g x,η )(r) ∈ Å for all r ∈ [0, +∞).
Lemma 2.6 can be proved as in [START_REF] Esteban | Symmetric ground states for a stationary relativistic mean-field model for nucleons in the non relativistic limit[END_REF] using the monotonicity properties of the energy. For the reader's convenience, we rewrite the proof here.

Proof. First of all, we use the monotonicity of the function

F (x) = a 4 x 4 -b 2 x 2 , the facts that F (x) ≤ 0 in -2b a , 2b
a and 2b a < 1 to show that F (x) < F (1) for all x such that x 2 < 1.

Let g x,η (0) = x such that x 2 < 1 and suppose, by contradiction, that there exists r 0 such that g 2

x,η (r 0 ) = 1 and g 2 x,η (r) < 1 for all r ∈ [0, r 0 ). As a consequence of Lemma 2.4, the energy H(f x,η , g x,η )(r) is non-increasing on [0, r 0 ), that means

H(0, x) ≥ H(f x,η , g x,η )(r 0 ) , or equivalently, F (x) ≥ F (1).
The above inequality contradicts the properties of F . As a conclusion, g 2

x,η (r) < 1 for all r ∈ [0, R x,η ).

Then, applying Lemma 2.4, we obtain that the energy is non-increasing. Thus,

H(f x,η , g x,η )(r) ≤ H(0, x) < a -2b 4 , ∀r ∈ [0, R x,η ) ,
and by the remark following the definition of the set A, (f x,η , g x,η )(r) ∈ Å and

f 2 x,η (r) < a -b, for all r ∈ [0, R x,η ). In particular, R x,η = +∞. The case x = ±1 is straightforward.
Remark 2.7. As a consequence of Lemma 2.4 and Lemma 2.6, if x 2 ≤ 1, the energy

H(f x,η , g x,η )(r) is non-increasing on [0, +∞).
Then, we state the following perturbation result. Lemma 2.8. Let ( f , g) ∈ A. Let (f, g) be the solution of (1.6) with initial data ( f , g). Let ( f n , gn ) ∈ A and ρ n be such that

lim n→+∞ ρ n = +∞ and lim n→+∞ ( f n , gn ) = ( f , g). Let (f n η , g n η ) be a solution of    (f n ) + 2ϕ η (f n , g n ) ρ n + r f n = g n ((f n ) 2 -a(g n ) 2 + b) (g n ) = f n (1 -(g n ) 2 ) such that f n (0) = f n , g n (0) = gn . Then (1) (f n η , g n η ) is defined on [0, +∞), (2) (f n η , g n η )
converges to (f, g) uniformly on bounded intervals. The proof of this lemma is a straightforward modification of the proof of [START_REF] Esteban | Symmetric ground states for a stationary relativistic mean-field model for nucleons in the non relativistic limit[END_REF]Lemma 3.2].

Next, we introduce the following definition to count the number of times that the solutions cross the set {g = 0}. Definition 2.9. We say that a continuous function g defined on an interval I changes sign at r 0 ∈ I if g(r 0 ) = 0 and there exists ε > 0 such that [r 0 -ε, r 0 +ε] ⊂ I, g(r 0 -ε)g(r 0 + ε) < 0 and g(r 0 -r)g(r 0 + r) ≤ 0 for all r ∈ (0, ε).

Remark 2.10. Let (f x,η , g x,η ) be a nontrivial solution of (2.2) with x ∈ (0, 1). Then, g x,η changes sign at 0 ≤ r 0 < +∞ if and only if g x,η vanishes at 0 ≤ r 0 < +∞.

Indeed, since (f x,η , g x,η ) is a nontrivial solution and

r 0 < +∞, f x,η (r 0 ) = 0. Hence, g x,η (r 0 ) = f x,η (r 0 )(1 -g x,η (r 0 ) 2 ) = 0 and g x,η changes sign at r 0 .
Finally, we state the following lemma which gives us an important qualitative property of the solutions of the system (2.1).

Lemma 2.11. Let x ∈ (0, 1). If (f x,η , g x,η ) is a solution of (2.2) such that g x,η changes sign a finite number of times and

lim r→+∞ H(f x,η , g x,η )(r) ≥ 0, then, for all r ≥ 0, (2.4) |f x,η (r)| + |g x,η (r)| ≤ C exp(-K a,b r) with K a,b = min b 2 , 2a-b 2a
and C a positive constant. In particular, we get

(2.5) lim r→+∞ (f x,η , g x,η )(r) = (0, 0).
Proof. First of all, we remark that if lim

r→+∞ g x,η (r) = δ, then |δ| = 1. Moreover, if δ = 0 then (f x,η , g x,η
) tends either to (0, b a ) or to (0, -b a ) as r goes to +∞. Indeed, suppose by contradiction that δ = ±1, then lim r→+∞ H(f x,η , g x,η )(r) = H(0, 1), which contradicts the monotonicity of H (Lemma 2.4) since H(0, x) < H(0, 1) and x ∈ (0, 1). Hence, -1 < δ < 1. Next, suppose δ = 0 and let {r n } n be a sequence such that lim n→+∞ r n = +∞ and lim

n→+∞ f x,η (r n ) = k for some k ∈ R. Let (u, v)
be the solution of (1.6) with initial data (k, δ). It follows from Lemma 2.8 that (f x,η (r n + •), g x,η (r n + •)) converges uniformly to (u, v) on bounded intervals. Since, lim n→+∞ g x,η (r n + r) = δ for any r > 0, we have v(r) = δ ∈ (-1, 0) ∪ (0, 1) for any r ≥ 0. Hence, from the second equation of (1.6), we obtain u(r) = 0 for all r > 0. This means that (u, v) is an equilibrium point of (1.6), and, since δ ∈ (-1, 0)∪(0, 1), this implies k = 0 and δ = ± b a . As a conclusion, f x,η converges as r goes to +∞,

lim r→+∞ (f x,η , g x,η )(r) = 0, ± b a and lim r→+∞ H(f x,η , g x,η )(r) = H 0, b a < 0.
Next, we claim that, if g x,η changes sign a finite number of times and

lim r→+∞ H(f x,η , g x,η )(r) ≥ 0,
then there exists R < +∞ such that

• or g x,η (r) > 0 and f x,η (r) < 0 for all r > R,

• or g x,η (r) < 0 and f x,η (r) > 0 for all r > R.

Indeed, by Remark 2.10, if g x,η changes sign a finite number of times, then g x,η vanishes a finite number of times and there exists R < +∞ such that g x,η (r) > 0 or g x,η (r) < 0 for all r > R. Thanks to the symmetries of the problem, we can suppose w.l.o.g. that g x,η (r) > 0 for all r > R.

Hence, it remains to prove that there exists R < R < +∞ such that f x,η (r) < 0 for all r > R. We proceed as follows : first we prove that we cannot have f x,η (r) > 0 for all r > R and second we show that f x,η vanishes at most once in [R, +∞).

Step 1. Suppose, by contradiction, that f x,η (r) > 0 for all r > R. This implies that g x,η (r) is increasing for all r > R and lim r→+∞ g x,η (r) = δ with 0 < δ ≤ 1.

Hence, as we proved above, we have lim

r→+∞ (f x,η , g x,η )(r) = 0, b a and lim r→+∞ H(f x,η , g x,η )(r) = H 0, b a < 0
which contradicts the fact that lim r→+∞ H(f x,η , g x,η )(r) ≥ 0 (see Remark 2.12 for an alternative proof). As a consequence, there exists R < R < +∞ such that f x,η ( R) = 0. Let us remark moreover that for such a R, we have f x,η ( R) < 0.

Step 2. Suppose next, by contradiction, that there exist a positive constant

R such that R < R < R < +∞, f x,η (R ) = 0 and f x,η (r) < 0 on ( R, R ). Since f x,η has to be nonnegative in a neighborhood of R , we can conclude that 0 < g x,η (R ) ≤ b a . Hence, lim r→+∞ H(f x,η , g x,η )(r) ≤ H (0, g x,η (R )) < 0 which contradicts the fact that lim r→+∞ H(f x,η , g x,η )(r) ≥ 0.
As a conclusion, there exists R < +∞ such that g x,η (r) > 0 and f x,η (r) < 0 for all r > R. This implies that g x,η (r) is decreasing for all r > R and lim r→+∞ g x,η (r) = δ with 0 ≤ δ < 1. We claim that δ = 0. Indeed, suppose by contradiction, δ = 0. As above, we obtain lim

r→+∞ (f x,η , g x,η )(r) = 0, b a and lim r→+∞ H(f x,η , g x,η )(r) = H 0, b a < 0
which contradicts the fact that lim r→+∞ H(f x,η , g x,η )(r) ≥ 0 (see Remark 2.12 for an alternative proof). Hence, lim r→+∞ g x,η (r) = 0 and g 2 x,η (r) ≤ b 2a for r large enough. Considering (2.1),

f x,η ≥ - 2ϕ η (f x,η , g x,η ) r f x,η + b 2 g x,η , g x,η ≤ 2a -b 2a f x,η .
Thus, for r large enough,

(g x,η -f x,η ) + K a,b (g x,η -f x,η ) ≤ 0, with K a,b = min b 2 , 2a-b 2a
. Integrating the above equation, we obtain

|f x,η (r)| + |g x,η (r)| ≤ C exp(-K a,b r)
for all r ≥ 0 with C > 0.

With exactly the same arguments, we treat the case g x,η (r) < 0 for all r > R.

Remark 2.12. The proof is very similar to the one of [START_REF] Esteban | Symmetric ground states for a stationary relativistic mean-field model for nucleons in the non relativistic limit[END_REF]Lemma 3.4]. With the same arguments of [5, Proof of Lemma 3.4], we can prove that if x ∈ (0, 1) and

(f x,η , g x,η ) is a solution of (2.2) such that lim r→+∞ (f x,η , g x,η )(r) = 0, ± b a ,
then f x,η has infinitely many zeros.

This property is equivalent to the fact that (f x,η , g x,η ) cannot tend to 0, ± b a , while being in one of the half-planes {f > 0} or {f < 0}.

This remark allows us to prove in an alternative way that if g x,η (r) > 0 for all r > R, then f x,η has to vanish at least once in (R, +∞) without using the fact that lim r→+∞ H(f x,η , g x,η )(r) ≥ 0 (Step 1 of the proof of Lemma 2.11).

Moreover, it proves also that if g x,η (r) > 0 and f x,η (r) < 0 for all r > R, then lim r→+∞ (f x,η , g x,η )(r) = (0, 0) without using lim r→+∞ H(f x,η , g x,η )(r) ≥ 0 (end of the proof of Lemma 2.11).

2.3. The shooting method. Following [START_REF] Balabane | Nodal solutions for a sublinear elliptic equation[END_REF], we define 2). Remark 2.13. By Lemma 2.6, we get that (f x,η , g x,η )(r) ∈ A for all r whenever x ∈ [0, 1]. Remark 2.7 ensures then that lim r→+∞ H(f x,η , g x,η )(r) exists for all x ∈ [0, 1].

I -1 = ∅ and, for k ∈ N and η ∈ (0, √ a -b -a 2 ) fixed, A k = {x ∈ (0, 1) : lim r→+∞ H(f x,η , g x,η )(r) < 0, g x,η changes sign k times on R + }, I k = {x ∈ (0, 1) : lim r→+∞ (f x,η , g x,η )(r) = (0, 0), g x,η changes sign k times on R + } (see Figure
f g x 1 f g x 2 f g x 3 H ≤ 0 Figure 2. x 1 ∈ A 1 , x 2 ∈ I 1 , x 3 ∈ A2
Remark 2.14. We want to find non trivial localized solutions of equations (2.1) with a given number of nodes that is to say, x ∈ (0, 1) such that lim r→+∞ (f x,η , g x,η )(r) = (0, 0) and g x,η changes sign k times on R + . To do this, we show by a shooting method that

I k = ∅ for all k ∈ N and all η ∈ (0, √ a -b -a 2
). The core of the shooting method is the following lemma which gives the main properties of the sets A k and I k . It is very similar to the properties stated in the proof of [2, Theorem 1] except that the sets A k and I k are always bounded since they are included in (0, 1). The good equivalent property which is adapted to our case is given by point (ii) of the next lemma.

Lemma 2.15. For all k in N and all η ∈ (0,

√ a -b -a 2 ) we have (i) A k is an open set, (ii) there is ε ∈ (0, 1) such that A k ∪ I k ⊂ (0, 1 -ε), (iii) if x ∈ I k , there exists ε > 0 such that (x -ε, x + ε) ⊂ A k ∪ I k ∪ A k+1 (iv) if A k is not empty, we have sup A k ∈ I k-1 ∪ I k , (v) if I k is not empty, we have sup I k ∈ I k ,
The proof of this lemma is given in Section 3. We are now able to prove the following proposition.

Proposition 2.16. There exists an increasing sequence {x k } k≥0 ⊂ (0, 1) such that

x k ∈ I k for all k ∈ N.
The proof is essentially the same as in [START_REF] Balabane | Nodal solutions for a sublinear elliptic equation[END_REF]. We write it down here for sake of completeness.

Proof. We prove by induction that for all k ∈ N,

A k = ∅, sup I k-1 < sup A k .
If this property is true for all k, then A k is not empty, sup A k ∈ I k by point (iv) of Lemma 2.15 and sup A k ≤ sup I k < sup A k+1 . Hence, if we choose x k = sup A k we get the proposition.

(1) Let k = 0. We have that for all x ∈ (0, (2) Let us assume now that for some k ∈ N, we have

A k = ∅, sup I k-1 < sup A k .
By point (iv) of Lemma 2.15, we get sup

A k ∈ I k which implies I k = ∅ and sup A k ≤ sup I k . Since I k = ∅, by point (v), we obtain that sup I k ∈ I k and, since sup A k ≤ sup I k , point (iii) ensures that there is ε > 0 such that (sup I k , sup I k + ε) ⊂ A k+1 .
As a conclusion, we have

A k+1 = ∅, sup I k < sup A k+1 .

Proof of Lemma 2.15

In this section, we fix η ∈ (0, √ a -b -a 2 ) and we prove Lemma 2.15. We show by induction that for all k, there is ε ∈ (0, 1) such that if x ∈ (1 -ε, 1) then g x,η has at least k + 1 changes of sign on R + . This implies

0≤i≤k A i ∪ I i ⊂ (0, 1 -ε)
and point (ii) follows.

Remark 3.2. The idea of the proof is that we can control the solutions (f x,η , g x,η ) thanks to the continuity of the flow on the parameter x (see Lemma 2.3) comparing (f x,η , g x,η ) to (f 1,η , g 1,η ) on an interval of the type [0, R] for R > 0. Moreover, (f 1,η , g 1,η ) tends to a stationary point (-√ a -b, 1) of the system (1.6). Thus, (f x,η , g x,η ) stay in a neighborhood of (-√ a -b, 1) a very long time if x is sufficiently close to 1. We also know thanks to Lemma 2.11 that (f x,η , g x,η ) exits this neighborhood at finite time, possibly very large. The problem is that we have to control the position of (f x,η , g x,η ) when this occurs. To do this, we replace the system (1.3) by the Hamiltonian ones (1.6) in this neighborhood. Then, we can use the conservation of the energy H along the trajectory of (f x,η , g x,η ) to know the position of (f x,η , g x,η ) when it exits the neighborhood of (-√ a -b, 1). After that, we can control the solutions (f x,η , g x,η ) thanks to the continuity of the flow comparing (f x,η , g x,η ) to a solution (f, g) of (1.6) that remains at all times on ∂A and tends to (-√ a -b, -1) at infinity. We get that if x is close enough to 1 then g x,η changes sign one time. We iterate this reasoning to obtain a solution for which g x,η changes sign more than k times on R + .

Step 1. Proof by induction

First of all, we take

f 0 = √ a -b -η/2 > a 2 > a-2b 2
and we define

X 1 := -f 0 , 2f 2 0 -(a -2b) a X 2 := (-f 0 , -1) X 3 := f 0 , - 2f 2 0 -(a -2b) a X 4 := (f 0 , 1) .
The points X i are on ∂A, for i = 1, . . . , 4. Furthermore, remind that ϕ η (f, g) = 0 whenever |f | ≥ f 0 (see Figure 1). Definition 3.3. Let k ∈ N and i ∈ {1, . . . , 4} be given. We denote by (H k i ) the following property: for all γ and R positive constants given, there exists ε > 0 such that for any x ∈ (1 -ε, 1), there exists a positive constant R > R which satisfies

(f x,η , g x,η )( R) ∈ B(X i , γ) ∩ A
and such that g x,η change k times of sign in [0, R].

In the second step, we show that the properties (H 0 1 ) is true. Next, in the third step, we prove that for k ∈ N given we have

       (H k 1 ) ⇒ (H k+1 2 ), (H k 2 ) ⇒ (H k 3 ), (H k 3 ) ⇒ (H k+1 4 ), (H k 4 ) ⇒ (H k 1 ) so that (H k 1 ) ⇒ (H k+2 1
).

As a consequence, we get by induction that the property (H 2k 1 ) is true for all k ∈ N. In particular, there is ε ∈ (0, 1) such that for all x ∈ (1 -ε, 1), g x,η changes at least 2k of sign on [0, +∞) so that

I i ∪ A i ⊂ (0, 1 -ε).
for all i ∈ {0, 1, . . . , 2k -1} and point (ii) of Lemma 2.15 is proved.

Step 2. Initialization: We prove that (H 0 1 ) is true. (1) Preliminary results. Let γ and R be positive constants given. First of all, remark that with the notation of Lemma A.1,

X 1 = (-f 0 , G 2 (H(0, 1))).
By continuity of G 2 , there exists δ > 0 such that H(0, 1) -δ > E c and

(-f 0 , G 2 (E)) -X 1 < γ
for all E ∈ (H(0, 1) -δ, H(0, 1)) where . is the Euclidean norm of R 2 . So, we have to prove that there exists ε > 0 such that for any x ∈ (1 -ε, 1), there exists a positive constant R 1 > R which satisfies

H(0, 1) -δ < H(f x,η , g x,η )(R 1 ) < H(0, 1), f x,η (R 1 ) = -f 0 , g x,η (R 1 ) = G 2 (H(f x,η , g x,η )(R 1 )) and g x,η does not change sign in [0, R 1 ]. (2) Control of the solutions of (2.2) in an interval [0, R] with R > 0.
We denote (f, g) the solution of Cauchy problem (2.2) with x = 1. It is easy to see that

H(f, g)(r) = H(0, 1), g(r) = 1, f (r) > - √ a -b for all r ∈ [0, +∞) and lim r→+∞ (f, g)(r) = (- √ a -b, 1).
As a consequence, there exists R > R such that for all r ≥ R

(f, g)(r) ∈ A ∩ {(f, g) : |f | > f 0 }.
Next, since H is continuous on R 2 , there exists 0 < δ < 1 such that for any (f , g) ∈ B((f, g)(R), δ ),

we have |H(f , g) -H(0, 1)| < δ
where B((f, g)(R), δ ) is the Euclidean ball of R 2 centered in (f, g)(R) of radius δ . Moreover, if we choose δ sufficiently small, we can assume that |u| > f 0 for all

(u, v) ∈ B((f, g)(R), δ ).
Finally, by Lemma 2.3, there exists ε ∈ (0, 1) such that for all x ∈ (1 -ε, 1),

(f x,η , g x,η ) -(f, g) ∞,[0,R] ≤ δ where . ∞,[0,R] is the uniform norm of C([0, R], R 2 ) so that g x,η is positive in [0, R]. (3) Control of the solutions of (2.2) in A ∩ {(f, g) ∈ R 2 : |f | > f 0 }. We define R 1 := inf{r > R : |f x,η (r)| ≤ f 0 }.
By Lemma A.2, we have that

H(f x,η , g x,η )(r) > 0, for all r ∈ [R, R 1 ). Moreover by Lemma A.1, since (f x,η , g x,η )(r) ∈ A ∩ {(f, g) ∈ R 2 : |f | > f 0 }, g x,η (r) ≥ 2 a (f 2 0 + b) -1 for all r ∈ [R, R 1 )
. Hence, by Lemma 2.11, we get that R 1 is well-defined, R < R 1 < +∞ and g x,η does not change sign in [0, R 1 ]. Furthermore, f x,η (R 1 ) = -f 0 and, since (f x,η , g x,η ) is solution of the Hamiltonian system of equation (1.6) on [R, R 1 ], we obtain

H(f x,η , g x,η )(R 1 ) = H(f x,η , g x,η )(R) ∈ (H(0, 1) -δ, H(0, 1)).
Hence, it remains to show that g

x,η (R 1 ) = G 2 (H(f x,η , g x,η )(R 1 )). Let R := inf{r > 0 : |f x,η (r)| > f 0 }, then f x,η ( R) = -f 0 = f x,η (R 1 ). Moreover, in [ R, R 1 ], (f x,η , g x,η
) is solution of the Hamiltonian system (1.6); this implies

H(f x,η , g x,η )( R) = H(f x,η , g x,η )(R 1 ). Finally, g x,η is decreasing on [ R, R 1 ]; in particular g x,η (R 1 ) < g x,η ( R).
Hence, by Lemma A.1, we deduce

g x,η ( R) = G 1 (H(f x,η , g x,η )(R 1 )), g x,η (R 1 ) = G 2 (H(f x,η , g x,η )(R 1 )).
Thanks to the remark we did in the preliminary results, we proved Step 2.

Step 3. Iteration: Let k ∈ N and suppose that property (H k 1 ) is true. We show that this implies property (H k+1

2

). The proof of this fact is similar to the one of Step 2 except that now (f, g) is a solution of autonomous system (1.6).

(1) Preliminary results. Let γ and R be positive constants given. First of all, remark that with the notation of Lemma A.1,

X 2 = (-f 0 , -G 1 (H(0, 1))).
By continuity of G 1 , there exists δ > 0 such that H(0, 1) -δ > E c and

(-f 0 , -G 1 (E)) -X 2 < γ
for all E ∈ (H(0, 1) -δ, H(0, 1)). So, we have to prove that there exists ε > 0 such that for any x ∈ (1 -ε, 1), there exists a positive constant R > R which satisfies

H(0, 1) -δ < H(f x,η , g x,η )(R ) < H(0, 1), f x,η (R ) = -f 0 , g x,η (R ) = -G 1 (H(f x,η , g x,η )(R )) and g x,η changes sign k + 1 times in [0, R ].
(2) Control of the solutions of (2.2) when the solutions exit a neighborhood of X 1 . We denote by (f, g) the solution of the following autonomous system

(1.6) (f, g)(0) = X 1 . It is clear that H(f, g)(r) = H(0, 1), -1 < g(r) < 1, f (r) > - √ a -b for all r ∈ [0, +∞),
and lim

r→+∞ (f, g)(r) = (- √ a -b, -1).
Hence, there is R > 0 such that for all r ≥ R

(f, g)(r) ∈ A ∩ {(f, g) : |f | > f 0 }.
Next, since H is continuous on R 2 , there exists δ > 0 such that for any

( f , g) ∈ B((f, g)(R), δ ),
we have |H( f , g) -H(0, 1)| < δ. Moreover, if we choose δ > 0 sufficiently small, we can assume that |u| > f 0 for all (u, v) ∈ B((f, g)(R), δ ).

By Lemma 2.8, there exist R > 0 and γ > 0 such that if ρ ≥ R and

( f , g) -X 1 < γ then (f f ,g,η , g f ,g,η )(• + ρ) -(f, g) ∞,[0,R] < δ .
Since by hypothesis, there are ε ∈ (0, 1) and for any

x ∈ (1-ε, 1) a constant R 1 > max(R, R) such that (f x,η , g x,η )(R 1 ) ∈ B(X 1 , γ) ∩ A
and g x,η changes sign exactly k times on [0, R 1 ], we get

(f x,η , g x,η )(• + R 1 ) -(f, g) ∞,[0,R] < δ .
In particular,

|(f x,η , g x,η )(R 1 + R) -(f, g)(R)| < δ and g x,η changes sign exactly k + 1 times on [0, R 1 + R]. (3) Control of the solutions of (2.2) in A ∩ {(f, g) ∈ R 2 : |f | > f 0 }. Let R 2 := inf{r > R + R 1 : |f x,η (r)| < f 0 }.
With the same arguments used in the proof of property (H 0 1 ), we prove that

(f x,η , g x,η )(R 2 ) ∈ B(X 2 , γ) ∩ A and g x,η changes sign exactly k + 1 times on [0, R 2 ]. We proved that (H k 1 ) ⇒ (H k+1 2
).

Thanks to the symmetry of the system, we also get

(H k 3 ) ⇒ (H k+1 4
).

The proof of the remaining implications

(H k 2 ) ⇒ (H k 3 ) (H k 4 ) ⇒ (H k 1 )
uses the same ideas.

3.2.

Proof of the remaining points of Lemma 2.15. In this part, we assume that η ∈ (0, √ a -b -a 2 ) is fixed. First of all, we remark that point (i) follows directly from Lemma 2.3. For the remaining points, we need the following preliminary lemma.

Lemma 3.4. There exists c

0 > 0 universal constant such that if (i) H(f x,η , g x,η )(R) < c0 R (ii) g x,η (R) ∈ 0, 2b a and f x,η (R) < 0 or g x,η (R) ∈ -2b a , 0 and f x,η (R) > 0, (iii) g x,η changes sign k times on [0, R]; for x ∈ (0, 1), R > 0, η ∈ (0, √ a -b -a 2 ) and k ∈ N, then x belongs to A k ∪ I k ∪ A k+1 .
Proof. We define

c 0 := 2b 9a(a -b) 80b 2 81a
.

We can assume thanks to the symmetries of the system that

H(f x,η , g x,η )(R) < c 0 R , g x,η (R) ∈ - 2b a , 0 and f x,η (R) > 0
for some x ∈ (0, 1) and R > 0. First of all, we remark that if there exists R such that H(f x,η , g x,η )( R) ≤ 0, then

ϕ η (f x,η , g x,η )( R) = 1
by Lemma A.2. Hence, by Lemma 2.4 and Remark 2.7, we deduce

H(f x,η , g x,η )(r) < 0 for all r > R. Since {(f, g) : g = 0} ∩ H -1 (-∞, 0) = ∅, g x,η
does not change sign anymore on ( R, +∞). Moreover, if g x,η has no changes of sign in (R, +∞) then either

lim r→+∞ H(f x,η , g x,η )(r) < 0, and x ∈ A k or lim r→+∞ H(f x,η , g x,η )(r) ≥ 0,
and x ∈ I k by Lemma 2.11. We assume, by contradiction, that x / ∈ A k ∪ I k ∪ A k+1 then g x,η changes sign at least once in (R, +∞). Next, we denote

R := inf{r > R : f x,η (r) ≤ 0} ∈ (R, +∞]. Since g x,η is increasing for all r ∈ [R, R], g x,η changes sign at most once before (f x,η , g x,η ) exits {(f, g) : f > 0}. Moreover, we claim that R < R < +∞. Indeed, if R = +∞, g x,η changes sign k or k + 1 times on R + . Then, we have either lim r→+∞ H(f x,η , g x,η )(r) < 0, and x ∈ A k ∪ A k+1 or lim r→+∞ H(f x,η , g x,η )(r) ≥ 0,
and x ∈ I k ∪ I k+1 by Lemma 2.11. Moreover, if x ∈ I k+1 , Lemma 2.11 ensures that g x,η decays exponentially to 0 this contradicts the fact that g x,η is positive and increasing between inf{r ≥ R :

g x,η (r) ≥ 0} ∈ (R, R) and R. Nevertheless, we assumed that x / ∈ A k ∪ I k ∪ A k+1 hence R < +∞. As a consequence, we get f x,η (R) = 0 and, since x / ∈ A k+1 , H(f x,η , g x,η )(r) > 0.
for all r ≤ R. Moreover, we have

g x,η (R) > 2b a since H(0, x) ≤ 0 for all x ∈ [-2b a , 2b a ]. Next, we denote R := sup{r ∈ (R, R), g x,η (r) ≤ 1 3 2b a } and R := inf{r > R : g x,η (r) ≥ 2 3 2b a },
and we remark that this quantities are well-defined. For all r ∈ (R, R ), we get

d dr H(f x,η , g x,η )(r) = - 2 r f 2 x,η (1 -g 2 x,η ) = - 4 r H(f x,η , g x,η )(r) - a 4 g x,η (r) 4 + b 2 g 2 x,η (r) , and
d dr r 4 H(f x,η , g x,η )(r) = r 3 g 2 x,η (r) ag 2 x,η (r) -2b ≤ 0 (3.1) since g 2 x,η (r) < 2b a for all r ∈ [R, R ]
. Moreover, we have

g x,η (r) ∈ 1 3 2b a , 2 3
2b a for all r ∈ [R , R ] and

1 3 2b a = g x,η (R ) -g x,η (R ) = R R f x,η (s)(1 -g 2 x,η (s))ds (3.2) ≤ √ a -b (R -R ).
Integrating inequality (3.1), we have, thanks to inequality (3.2),

(R ) 4 H(f x,η , g x,η )(R ) -(R ) 4 H(f x,η , g x,η )(R ) ≤ -c 1 (R ) 4 -(R ) 4 (3.3) ≤ -c 1 (R -R )(R 3 + R 2 R + R R 2 + R 3 ) ≤ -4c 1 2b 9a(a -b) R 3 = -c 0 R 3 for c 1 = 20b 2 81a , since c 0 = 4c 1 2b 9a(a -b) = 2b 9a(a -b) 80b 2 81a
.

Then, we obtain by inequalities (3.1) and

(3.3) (R ) 4 H(f x,η , g x,η )(R ) ≤ -c 0 R 3 + (R ) 4 H(f x,η , g x,η )(R ) ≤ R 4 - c 0 R + H(f x,η , g x,η )(R) < 0.
This is impossible since H(f x,η , g x,η )(R ) > 0.

3.2.1. Proof of point (iii) of Lemma 2.15.

Lemma 3.5. Let k ∈ N and η ∈ (0, √ a -b -a 2 ). If x ∈ I k then there is ε > 0 such that (x -ε, x + ε) ⊂ A k ∪ I k ∪ A k+1 .
Proof. By Lemma 2.11, there exists C, K > 0 such that

(|f x,η | + |g x,η |)(r) ≤ C exp(-Kr)
for all r and H(f x,η , g x,η ) converges exponentially to 0. We easily get that there is R such that the assumptions of Lemma 3.4 are fulfilled for x at R. Then, by Lemma 2.3, there is ε > 0 such that for all y ∈ (x -ε, x + ε), g y,η changes sign k times on [0, R],

H(f y,η , g y,η )(R) < c 0 R and f y,η (R) < 0, g y,η (R) ∈ (0, 2b a ) or f y,η (R) > 0, g y,η (R) ∈ (- 2b a , 0).
Thus, by Lemma 3.4, we have that

(x -ε, x + ε) ⊂ A k ∪ I k ∪ A k+1 .
3.2.2. Proof of point (iv) of Lemma 2.15.

Lemma 3.6. Let k ∈ N and η ∈ (0, √ a -b -a 2 ). If A k is non-empty, then sup A k ∈ I k ∪ I k-1 .
Proof. Thanks to points (i) and (ii) of Lemma 2.15,

x := sup A k ∈ (0, 1)\ ∪ n∈N A n .
Let {x i } ⊂ A k be such that lim

i→+∞ x i = x.
Suppose by contradiction that x / ∈ ∪ n∈N (I n ∪ A n ), then in particular H(f x,η , g x,η )(r) > 0 for all r > 0. Moreover, Lemma 2.11 ensures that g x,η changes sign an infinite number of times. Let R > 0 be such that g x,η changes sign more than k + 1 times in [0, R] at 0 < r 1 < • • • < r k+1 . Hence, there is ε > 0 such that for all j ∈ {1, . . . , k + 1}, all r ∈ (0, ε) we have g x,η (r j -r)g x,η (r j + r) ≤ 0 and g x,η (r j -ε)g x,η (r j + ε) < 0. Then, by Lemma 2.3, there is M > 0 such that if i ≥ M then g xi,η (r j -ε)g xi,η (r j + ε) < g x,η (r j -ε)g x,η (r j + ε)/2 < 0 for all j ∈ {1, . . . , k + 1}. Thus, for all i ≥ M and all j ∈ {1, . . . , k + 1}, there is a real number r i j ∈ (r j -ε, r j + ε) such that g xi,η (r i j ) = 0. Then, we get g xi,η (r i j ) = 0 so that g xj ,η changes sign more that k + 1 times at the points r i j . This is impossible because x i ∈ A k . Hence, we have that sup A k ∈ I m for some m ∈ N and by point (iii), we get the result. Proof. The proof follows the same ideas as the one of Lemma 3.6. We get that sup I k ∈ I j for some j ∈ N and by point (iii), we get the result.

Proof of theorem 1.1

We give now the proof of Theorem 1.1 by taking the limit when η tends to 0.

Proof. Let us fix k ∈ N. For all η ∈ (0, √ a -b -a 2 )
, there is x η ∈ (0, 1) such that lim r→+∞ (f xη,η , g xη,η )(r) = (0, 0), g xη,η has k changes of sign on (0, +∞) by Proposition 2.16. We also know that H(f xη,η , g xη,η )(r) > 0 for all r ∈ R + . Since for all x ∈ (0, 2b a ) H(0, x) < 0, we deduce that

{x η } η ⊂ 2b a , 1 .
Thus, there is a subsequence

{η n } n such that      lim n→+∞ η n = 0 lim n→+∞ x ηn = x 0 ∈ 2b a , 1 .
By Lemma 2.3, we get that for all R > 0, ε > 0, there exists

N > 0 such that if n ≥ N then (f x0,0 , g x0,0 ) -(f xη n ,ηn , g xη n ,ηn ) ∞,[0,R] ≤ ε
where . ∞,[0,R] is the uniform norm on the set C([0, R], R 2 ). Thus, (f x0,0 , g x0,0 ) is a solution of equations (1.3) such that (f x0,0 , g x0,0 )(0) = (0, x 0 ) H(f x0,0 , g x0,0 )(r) ≥ 0, for all r ∈ R + and Remark 2.7 ensures that x 0 ∈ 2b a , 1 .

To conclude, we have to show now that x 0 ∈ 2b a , 1 since (f 1,0 , g 1,0 ) is not a localized solution of (1.3).

Assume, by contradiction, that x 0 = 1; then

H(f x0,0 , g x0,0 )(r) = H(0, 1)
for all r ≥ 0. We denote H 0 := H(0, 1)/2 > 0 and we define

R n := inf{r > 0 : H(f xη n ,ηn , g xη n ,ηn )(r) ≤ H 0 } ∈ (0, +∞).
We have H(f xη n ,ηn , g xη n ,ηn )(R n ) = H 0 . We claim that {R n } n tends to +∞. Indeed, if there is a subsequence also denoted {R n } n and a real number R > 0 such that R n ∈ [0, R] then, we get the following contradiction

H(0, 1) = lim n→+∞ H(f xη n ,ηn , g xη n ,ηn )(R) ≤ H 0 . Next, since A is a compact set, there is (f 0 , g 0 ) ∈ A such that up to extraction, lim n→+∞ (f xη n ,ηn , g xη n ,ηn )(R n ) = (f 0 , g 0 )
and H(f 0 , g 0 ) = H 0 . We denote by T > 0 the period of the solution of the Hamiltonian system of equations (1.6) of energy equal to H 0 . Let us consider now the following Cauchy problem (4.1)

   f + 2ϕη(f,g) r+ρ f = g(f 2 -ag 2 + b), g = f (1 -g 2 ), (f, g)(0) = X.
Its solutions depend continuously on the parameters (X, ρ, η) on every interval

[0, R] just as in Lemma 2.3. So, (f xη n ,ηn , g xη n ,ηn )( . + R n ) tends uniformly on [0, (k+2)T 2 ] to a solution (f, g) of    f = g(f 2 -ag 2 + b), g = f (1 -g 2 ), (f, g)(0) = (f 0 , g 0 ).
Moreover, (f, g) is periodic of period T which implies that g has at least k + 1 changes of sign on 0, (k+2)T 2 i.e. there is ε > 0, and 0 < r 1 < • • • < r k+1 < (k+2)T 2 such that g(r i + r)g(r i -r) ≤ 0 for all r ∈ (0, ε) and g(r i + ε)g(r i -ε) < 0 for all i ∈ {1, . . . k + 1}.

Hence, there is N > 0 such that for all n > N and all i ∈ {1, . . . k + 1}, we get

g xη n ,ηn (R n + r i + ε)g xη n ,ηn (R n + r i -ε) < g(r i + ε)g(r i -ε)/2 < 0, which implies that there is r i n ∈ (r i -ε, r i + ε) such that g xη n ,ηn (R n + r i n ) = 0. As a conclusion, since (f xη n ,ηn , g xη n ,ηn
) is a solution of (2.1), g xη n ,ηn ((R n +r i n )) = 0 and g xη n ,ηn has at least k + 1 changes of sign at the points r i n . This is impossible because g xη n ,ηn has exactly k changes of sign. As a consequence x 0 < 1.

Moreover, with the same arguments used above, we prove that g x0,0 changes sign a finite number of times k 0 . Hence, Lemma 2.11 ensures that (f x0,0 , g x0,0 ) converge exponentially to (0, 0).

As a consequence, there is R > 0 such that g x0,0 changes sign k 0 times in [0, R],

H(f x0,0 , g x0,0 )(R) < c 0 R and f x0,0 (R) < 0, g x0,0 (R) ∈ (0, 2b a ) or f x0,0 (R) > 0, g x0,0 (R) ∈ (- 2b a , 0).
Hence, there is ε ∈ (0, 1) such that for all x ∈ (x 0 -ε, x 0 + ε) and all η ∈ (0, ε), g x,η changes sign k 0 of times in [0, R],

H(f x,η , g x,η )(R) < c 0 R and f x,η (R) < 0, g x,η (R) ∈ (0, 2b a ) or f x,η (R) > 0, g x,η (R) ∈ (- 2b a , 0).
By applying Lemma 3.4, we get, for η ∈ (0, ε) fixed,

(x 0 -ε, x 0 + ε) ⊂ A k0 ∪ I k0 ∪ A k0+1 .
Remark that the set A k and I k depends on η, hence, it is important to fix η before writing such a property. By definition of {x ηn } n , there is

N ∈ N such that for all n ≥ N x ηn ∈ (x 0 -ε, x 0 + ε), η n ∈ (0, ε).
As a conclusion, for n ≥ N fixed x ηn ∈ I k0 . This ensures that k = k 0 . Finally, by Remark 2.10, if g x0,0 changes sign k times on (0, +∞) then g x0,0 has k zeros on (0, +∞). To conclude, it remains to prove that f x0,0 has k zeros on (0, +∞) as well.

Let {r 1 , . . . , r k } be the zeros of g x0,0 on (0, +∞). First of all, we prove by induction that, for all i = 1 . . . , k, f x0,0 has i -1 zeros on (0, r i ). This property is true for i = 1. Indeed, suppose by contradiction that f x0,0 (r) = 0 for some r ∈ (0, r 1 ). Hence, using the first equation of (1.3), we get H(f x0,0 , g x0,0 )(r) < 0. That is impossible. Next, if the property holds true for i-1, then it holds true for i. Indeed, suppose that f x0,0 has i -2 zeros on (0, r i-1 ); we prove that f x0,0 has 1 zero on (r i-1 , r i ). By contradiction, if f x0,0 does not change sing on (r i-1 , r i ), the second equation of (1.3) implies that g x0,0 is monotone on (r i-1 , r i ). This contradicts the fact that g x0,0 changes sign at r i-1 and r i . Hence, f x0,0 has at least 1 zero at r ∈ (r i-1 , r i ). Now suppose that there exists r ∈ (r, r i ) such that f x0,0 (r) = 0. One of the following situations arise: g(r) < 0 and f (r) > 0 on (r, r) or g(r) > 0 and f (r) < 0 on (r, r). Using again the first equation of (1.3), in both cases we get g 2 x0,0 (r) ≤ b a which implies H(f x0,0 , g x0,0 )(r) < 0, a contradiction. Hence, for all i = 1 . . . , k, f x0,0 has i -1 zeros on (0, r i ). Finally, using the same arguments, we show that f x0,0 has 1 zero on (r k , +∞) and we conclude that f x0,0 has k zeros on (0, +∞).

Appendix A. Geometric properties of H.

We remind that A = {(f 0 , g 0 ) ∈ R 2 | 2f 2 0 -ag 2 0 -(a -2b) ≤ 0, g 2 0 ≤ 1} is the set of admissible points. Let us remark that (f, g) ∈ R 2 satisfies H(f, g) = H(0, 1) = i.e. if and only if g 2 = 1 or f 2 -a 2 (g 2 + 1) + b = 0. Å is the connected component of {(f, g) : H(f, g) = H(0, 1)} which contains (0, 0). Thanks to the symmetries of H, we can restrict our study of A to the set {(f, g) : f ≥ 0, g ≥ 0}.

G : 2f 2 0 -(a -2b) a , 1 → R g → H(f 0 , g).
Since G (g) = g(-b + ag 2 -f 2 0 ), we deduce that G is continuous, increasing in On the other hand, if g 2 ≤ 1 and H(f, g) ≤ 0, then g 2 ≤ 2b a . As a consequence, we get H(f , g) > 0.

Lemma 2 . 4 .

 24 Let x ∈ R and a, b > 0. Then for any r ∈ [0, R x,η ) we have

  sup I -1 < sup A 0 .

3. 1 .

 1 Proof of point (ii) of Lemma 2.15. Remark 3.1. This proof is the most technical point of the paper and contains the main novelties of our work. The introduction of the Hamiltonian regularization of subsection 2.1 allow us to control the behavior of the solution in the neighborhood of the stationary points (± √ a -b, ±1) of the autonomous system of equations (1.6).

3. 2 . 3 .Lemma 3 . 7 .

 2337 Proof of point (v) of Lemma 2.15. Let k ∈ N and η ∈ (0, √ a -b -a 2 ). If I k is non-empty, then sup I k ∈ I k .

1 4 2 f 2 ( 1 -g 2 ) + a 2 (g 4 - 1 )

 1221241 (a -2b) if and only if 0 = H(f, g) -H(0, 1) = 1 -b(g 2 -

Lemma A. 1 . 1 G 2 :∪ 1 Next, let f 0 ∈ a-2b 2 ,

 11212 Let f 0 ∈ a-2b 2 ,√a -b and E c = H f 0 , continuous functionsG 1 : [E c , H(0, 1)] → f 2 0 + b a , [E c , H(0, 1)] → 2 a (f 2 0 + b) -1, f 2 0 + b a such that, for i = 1, 2 and for all E ∈ [E c , H(0, 1)], H(f 0 , G i (E)) = E, {(f 0 , G 1 (E)), (f 0 , G 2 (E))} = {(f 0 , g) ∈ A : g ≥ 0} ∩ H -1 ({E}),andG 1 (E) ≥ G 2 (E).Proof. First of all, we observe thatA ∩ (R + ) 2 = {(f, g) ∈ [0, √ a -b] × [0, 1] : ag 2 ≥ 2f 2 -(a -(f, g) : f ∈ a -2b 2 , √ a -b , g ∈ 2f 2 -(a -2b) a ,√ a -b fixed, and define the function

. 1 E 2 aG 2 : 2 ,< 2f 2 -

 12222 < √ a -b. Moreover, we haveG 2f 2 0 -(a -2b) a = G(1) = H(0, 1)and we define E c := G As a consequence, G is a one-to-one function fromf 2 0 +b a , 1 onto [E c , H(0,1)] whose inverse is continuous and is denoted by G 1 . Therefore,G 1 : [E c , H(0, 1)] → G 1 (E) is a continuous function such that H(f 0 , G 1 (E)) = E. Similarly, we denote by G 2 the inverse of the restriction of the function G to the set [E c , H(0, 1)] → 2 a (f 2 0 + b) -1, f 2 0 + b a E → G 2 (E) is a continuous function such that H(f 0 , G 2 (E)) = E. Moreover, G 1 (E) ≥ G 2 (E) for all E ∈ [E c , H(0, 1)]. Lemma A.2. Let f 0 ∈ a √ a -b . If (f , g) ∈ A ∩ {(f, g) ∈ R 2 : |f | ≥ f 0 }, then H(f , g) > 0.Proof. First of all, since f 0 > a 2 > a-2b 2 , we have for(f , g) ∈ A ∩ {(f, g) ∈ R 2 : |f | ≥ f 0 } = {(f, g) ∈ [f 0 , √ a -b] × [-1, 1] : ag 2 ≥ 2f 2 -(a -2b)}, that 2b a (a -2b) a ≤ g 2 ≤ 1.

Acknowledgment. This work was partially supported by the Grant ANR-10-BLAN 0101 of the French Ministry of research. The authors would like to thank Maria J. Esteban for the helpful discussions.