Youwa Mbetkwe

Laetitia

Jean Pierre Nguenang

Paul André Paglan

Thierry Dauxois

Andrea Trombettoni

Stefano Ruffo

Fractional Dynamics and Modulational Instability in Long-Range Heisenberg Chains

Keywords: Heisenberg spin chains, long-range interactions, fractional equations, modulational instability

We study the effective dynamics of ferromagnetic spin chains in presence of long-range interactions. We consider the Heisenberg Hamiltonian in one dimension for which the spins are coupled through power-law long-range exchange interactions with exponent α. We add to the Hamiltonian an anisotropy in the z-direction. In the framework of a semiclassical approach, we use the Holstein-Primakoff transformation to derive an effective long-range discrete nonlinear Schrödinger equation. We then perform the continuum limit and we obtain a fractional nonlinear Schrödinger-like equation. Finally, we study the modulational instability of plane-waves in the continuum limit and we prove that, at variance with the short-range case, plane waves are modulationally unstable for α < 3. We also study the dependence of the modulation instability growth rate and critical wave-number on the parameters of the Hamiltonian and on the exponent α.

I. INTRODUCTION

Simplified models of magnetic systems, like the Ising model, have allowed the understanding of complex magnetic phenomena and the theoretical interpretation of phase transitions [START_REF] Parisi | Statistical Field Theory[END_REF]. Specifically, the study of low dimensional spin systems has been at the center of a constant attention in the field of magnetic materials [START_REF] Mattis | The theory of magnetism. I, Statics and dynamics[END_REF][START_REF] Majlis | The Quantum Theory of Magnetism[END_REF][START_REF]Introduction to Frustrated Magnetism: Materials, Experiments, Theory[END_REF]. Applications of models of magnetism cover a very wide range of phenomena, from classical and quantum transport to nuclear magnetic resonance [START_REF] Feldner | Propriétés magnétiques des systèmes à deux dimensions[END_REF]. An important role in the study of magnetic systems is played by the description of strong correlations and of the nonlinearities induced by inter-spin interactions [START_REF] Kivshar | [END_REF][START_REF] Bar'yakhtar | Dynamics of Topological Magnetic Solitons[END_REF]. On the other hand, the effects of longrange interactions, where each unit is coupled to all others, motivated a remarkable activity in the last decades [START_REF] Campa | [END_REF][START_REF] Campa | Physics of Long-range Interacting Systems[END_REF][START_REF] Defenu | [END_REF]. First studies of models of classical spins with long-range interactions date back to Refs. [11,12]. A typical interaction which is relevant for a number of systems, ranging from gravitational ones to dipolar magnets and gases, is provided by the power-law decay 1/r α , where r is the distance among constituents. It is well known that a first criterion to determine the long-range nature of a system with power-law decaying interactions is based on the comparison of α with the space dimension d. If α < d and if the system is homogeneous, then it displays a diverging energy density. Thus, in order to obtain a well defined thermodynamic limit, the so-called Kac rescaling of the energy is required. In the following, we will refer to this region as the non-additive or strong long-range region. If α > d, the energy of the system is additive and it is useful to introduce a value of α, which we denote by α * , such that for α > α * the system behaves at criticality as a short-range system [13]. Since α * > d, there is a region of values of α, given by d < α < α * , in which the behavior of the system significantly differs from the one of the same system with short-range interactions, although the energy is additive: such a region is the weak long-range region. The actual value of α * depends on the specific model and the dimension d. As an example, for classical O(n) models in d = 1, it is known that α * = 2, see e.g. [14][15][16]. These studies have been also performed for quantum spin systems with long-range couplings [17,18].

A common feature for systems with long-range interactions is the behaviour of the dispersion relation as ∝ k σ where σ = α -d. This behaviour implies, in many cases, that the low-energy effective dynamics can be described by a fractional differential equation [START_REF] Tarasov | Fractional Dynamics: Applications of Fractional Calculus to Dynamics of Particles, Fields and Media[END_REF], which involves derivatives and/or integrals of non-integer order. Such equations are used to describe anomalous kinetics, anomalous transport phenomena, etc. [START_REF] Metzler | [END_REF]. Fractional dynamics is a field of study which investigates the behavior of microscopic and macroscopic systems that are either characterized by power-law non locality, or by power-law long term memory. Recently, models of coupled nonlinear oscillators with long-range couplings have been introduced with the aim of understanding the role played by nonlocal interactions in several phenomena: relaxation to equilibrium [21], anomalous transport [22], localized solutions [23]. Among various possible examples, we can mention the Fermi-Pasta-Ulam-Tsingou model [24] with long-range couplings [25,26], for which one can derive explicitly an effective fractional equation, the fractional Boussinesq equation, whose solutions depend on the long-range exponent α.

A natural extension of these studies is to consider other models, in particular those relevant for the study of magnetism. In this paper, we investigate the Heisenberg model with long-range couplings [27,28] and we derive for the first time in the large-S limit, S being the value of the spin, a fractional nonlinear Schrödinger (FNLS) equation for the effective dynamics. To this aim, we proceed by mapping the Heisenberg spin chain with long-range interactions into a continuum equation with the Riesz fractional derivatives using the Holstein-Primakoff transformation. The final equations of motion describe the dynamics of the Heisenberg model with long-range couplings in the continuum limit. In order to show an application of the equations that we obtain, we study modulational instabilities and we discuss how the modulationally unstable and stable regions depend on the range of the power-law exponent α. Modulational instabilities were introduced and studied in lattice models [29][30][31] to study the localization of energy and the stability of dynamical regimes [START_REF] Dauxois | Physics of Solitons[END_REF]. In particular, the modulational instabilities of the discrete nonlinear Schrödinger equations with long-range couplings were considered in Ref. [START_REF] Gori | [END_REF]. Here, the derivation of the effective fractional equation allows for a qualitative determination of modulationally stable regions in the space of parameters.

The paper is organized as follows. In Section II, we derive the effective Hamiltonian using bosonic operators. In Section III, we derive the FNLS equation and the corresponding integro-differential equation valid in the continuum limit. In Section IV, we study the modulational instability of plane-waves excitations. In Section V, we discuss the results and provide an outlook on future directions of investigation. In the Appendices, we present some details of the analytical computations leading to the results discussed in the main text.

II. THE MODEL AND THE HOLSTEIN-PRIMAKOFF TRANSFORMATION

In order to study the dynamics of the one-dimensional (1D) Heisenberg spin chain for spins S = (S x , S y , S z ) located on the sites of a chain, let us start by considering the Hamiltonian

H = -A N n=-N (S z n ) 2 - -N ≤n<m≤N J m,n S n • S m (1)
where the indices n and m denote the sites of the lattice and N is the total number of spins. S is the value of the spin.

In the first term of the Hamiltonian, A represents the anisotropy parameter and in the second term, the exchange interaction J m,n is supposed to decay algebraically with distance as a power-law

J m,n = J |m -n| α . (2)
We consider a ferromagnetic chain (J > 0) and choose α in the range 1 < α < 3, because for smaller values, 0 < α ≤ 1, the Hamiltonian diverges. This latter is the non additive strong long-range region, we therefore limit in this paper to analyze the weak long-range region. For larger values, α > 3, the system in the low-energy limit is expected to display a short-range behavior, which is exactly reached only when α → ∞.

The cross product of the spins at two different sites can be rewritten in terms of annihilation and creation operators as follows

S n • S m = 1 2 S + n S - m + S - n S + m + S z n S z m . (3)
We then use the Holstein-Primakoff transformation for spins operators in terms of bosonic operators in the framework of the low temperature approximation [START_REF] Ashcroft | Solid State Physics[END_REF][START_REF] Bishop | Solitons and Condensed Matter[END_REF][START_REF] Nguenang | [END_REF] as

S + n = √ 2ε 1 - ε 2 4 a + n a n + o ε 4 a n (4)
S - n = √ 2ε a + n 1 - ε 2 4 a + n a n + o ε 4 (5)
S z n = 1 -ε 2 a + n a n , (6)
where ε = 1/ √ S. Therefore, in the limit of large S, substituting Eqs. (4-6) into Eq. (3), we get

S n • S m = 1 + ε 2 a n a + m + a + n a m -a + n a n -a + m a m - ε 4 4 a + n a n a n a + m + a n a + m a + m a m + a + n a + n a n a m + a + n a + m a m a m -4a + n a n a + m a m + o ε 6 , (7)
while

(S z n ) 2 = 1 -2ε 2 a + n a n + ε 4 a + n a n a + n a n + o ε 6 . (8)
Hamiltonian (1) can thus be rewritten in terms of bosonic creation and annihilation operators as

H 1 = - -N ≤n<m≤N J m,n F a n a + m + a + n a m -a + n a n -a + m a m -G a + n a n a n a + m + a n a + m a + m a m + a + n a + n a n a m +a + n a + m a m a m -4a + n a n a + m a m - N n=-N -T a + n a n + Ia + n a n a + n a n + o ε 6 , (9)
where H 1 = H + N A + -N ≤n<m≤N J m,n and the coupling constants are F = ε 2 , G = F 2 /4, T = 2F A, and I = AF 2 . Equation ( 9) is a bosonic, Bose-Hubbard-like, Hamiltonian which will be used in the next Section.

III. EQUATIONS OF MOTION

The equations of motion are obtained from the Heisenberg evolution equation for a given bosonic operator a n as follows

i ȧn = [a n , H] , -N ≤ n ≤ N (10)
One has to compute the different commutators arising from the previous equations to get the following nonlinear discrete equation for the bosonic operators

i ȧn = - m =n J mn -F (a n -a m ) -G a n a n a + m + a + m a m a m + 4a + n a n a m -4a n a + m a m +T a n -I a n a + n a n + a + n a n a n . (11)
One can write mean-field, classical equations for the dynamics of the expectation values of the bosonic operators. One way to do it is to use the Glauber coherent states representation which are defined as a + n |f = f * n |f and a n |f = f n |f and apply it for the full equation of motion through

f | i . a n |f = f | [a n , H] |f .
We then get the following nonlinear equations of motion for the amplitude f n of the spin chain excitation's

i . f n = m =n J mn F (f n -f m ) + G f * m f 2 n + |f m | 2 f m + 4 |f n | 2 f m -4 |f m | 2 f n + T f n -2I |f n | 2 f n . (12)
Equation ( 12) is a discrete cubic nonlinear Schrödinger-like equation which is the subject of our analysis in the following.

A. The continuum limit

Following the different steps proposed by Tarasov [23,38], let us define the operation which transforms the above equation for f n (t) into a continuum medium equation for f (x, t). We assume that f n (t) are Fourier coefficients of some function f (k, t). Then, we define the field f (k, t)

on the interval [-K/2, K/2] as f (k, t) = +∞ n=-∞ f n (t)e -ikxn = F ∆ {f n (t)}, (13)
where x n = n∆x and ∆x = 2π/K is a distance between oscillators, and conversely

f n (t) = 1 K +K/2 -K/2 dk f (k, t)e ikxn = F -1 ∆ { f (k, t)}. (14)
These equations define a Fourier transform, which is obtained in the limit ∆x → 0 (K → ∞). In order to perform this limit, let us replace the discrete set of functions f n (t) = (2π/K)f (x n , t) with continuous function of two variables f (x, t), while letting

x n = n∆x = 2πn/K → x.
Then, change the sum into an integral, and Eqs. ( 13) and ( 14) become

f (k, t) = +∞ -∞ dx e -ikx f (x, t) = F{f (x, t)} (15) f (x, t) = 1 2π +∞ -∞ dk e ikx f (k, t) = F -1 { f (k, t)}. (16)
Note that f (k, t) is the Fourier transform of the field f (x, t), and that f (k, t) is the Fourier series of f n (t), defined by

f n (t) = (2π/K)f (n∆x, t).
The map of the discrete model into the continuum one can be defined by the chain of transformations T = F -1 L F ∆ , where F ∆ is the Fourier series transform F ∆ {f n (t)} = f (k, t), the passage to the limit ∆x → 0 is represented by the operator L{ f (k, t)} = f (k, t) and the inverse Fourier transform is

F -1 { f (k, t)} = f (x, t). One has therefore f n (t) F∆ --→ f (k, t) L -→ f (k, t) F -1 ---→ f (x, t) = T {f n (t)}. ( 17
)
After introducing the general formalism needed to perform the continuum limit, let us analyze, using these methods, the discrete description displayed in Eq. ( 12) and derive the corresponding continuum one. Let us begin by analytically evaluating each of the terms in this Equation. Since most of the analytic computations will be performed in the framework of the continuum limit, we will consider the limit N → ∞.

To this end, let consider the first term

A 1 = -N ≤n<m≤N J m,n . (f n (t) -f m (t)). ( 18
)
The computations are performed in Appendix B in the similar spirit of those in which fractional equations are derived in Refs. [START_REF] Tarasov | Fractional Dynamics: Applications of Fractional Calculus to Dynamics of Particles, Fields and Media[END_REF]23,[37][38][START_REF] Klafter | Fractional Dynamics[END_REF][START_REF] Herman | Fractional Dynamics[END_REF][41]. We obtain

A 1 = a α ∂ α-1 f (x, t) ∂ |x| α-1 . ( 19
)
where

a α = 2JΓ(1 -α) sin (πα/2) ( 20
)
where Γ is the Euler gamma function and the partial Riesz derivative is defined as

- ∂ α u (x, t) ∂ |x| α = 1 2π +∞ -∞ dp |p| α û (p, t)e ipx . ( 21
)
The following term reads

B = -N ≤n<m≤N J m,n f * m f 2 n (22) = 1 2π 3 +∞ -∞ -N ≤n<m≤N J |m -n| α f * (k, t)e -ikm dk +∞ -∞ f (k , t)e ik n dk +∞ -∞ f (k , t)e ik n dk . (23)
In the continuum limit N → ∞ and setting m -n = m , we get

B = 1 2π 3 +∞ -∞ +∞ m =-∞,m =0 Je -ikm |m | α f * (k, t)e -ikn dk +∞ -∞ f (k , t)e ik n dk +∞ -∞ f (k , t)e ik n dk (24) = 1 2π 3 +∞ -∞ Ĵα (-k) f * (k, t)e -ikn dk +∞ -∞ f (k , t)e ik n dk +∞ -∞ f (k , t)e ik n dk . (25)
It is useful to introduce the following function

Ĵα (x) = +∞ m =-∞,m =0 J e im x m α . (26)
Next, we use the so-called infrared limit approximation, which is helpful to derive the main relation that allows us to pass from the discrete medium to the continuum [23,37,38]. For 1 < α < 3, α = 2 and k -→ 0, the fractional power of |k| is a leading asymptotic term and

Ĵα (-k) = Ĵα (k) a α |k| α-1 + Ĵα (0). ( 27
)
This allows us to transform the nonlinear discrete equation into a fractional differential equation. In the range 1 < α < 3, provided α = 2, we get

B = 1 2π 3 a α +∞ -∞ |k| α-1 f * (k, t)e -ikn dk +∞ -∞ f (k , t)e ik n dk +∞ -∞ f (k , t)e ik n dk + Ĵα (0) |f (x, t)| 2 f (x, t). (28)
Therefore, the fractional derivative allows us to write

B = 1 2π 3 a α +∞ -∞ |k| α-1 f * (k, t)e -ikn dk +∞ -∞ f (k , t)e ik n dk +∞ -∞ f (k , t)e ik n dk + Ĵα (0) |f (x, t)| 2 f (x, t). (29)
Thereafter, we get

B = -a α f 2 (x, t) ∂ α-1 f * (x, t) ∂ |x| α-1 + Ĵα (0) |f (x, t)| 2 f (x, t). (30)
Following similar analytical steps, after the computations given in Appendix B, the three following terms are derived

C = -N ≤n<m≤N J mn f m |f m | 2 = -a α |f (x, t)| 2 ∂ α-1 f (x, t) ∂ |x| α-1 + Ĵα (0) |f (x, t)| 2 f (x, t), (31)
D = -N ≤n<m≤N J mn f m |f n | 2 = -a α |f (x, t)| 2 ∂ α-1 f (x, t) ∂ |x| α-1 + Ĵα (0) |f (x, t)| 2 f (x, t), (32)
E = -N ≤n<m≤N J mn |f m | 2 f n = -a α |f (x, t)| 2 ∂ α-1 f (x, t) ∂ |x| α-1 + Ĵα (0) |f (x, t)| 2 f (x, t). (33)
Next, introducing the following ansatz

f (x, t) = ψ(x, t) e -i2ε 2 At/ (34)
and the parameters U α = -F a α , B α = Ga α , D α = 2I -2G Ĵα (0), equation ( 12) now reads

i ∂ψ ∂t + U α ∂ α-1 ψ ∂ |x| α-1 + B α ψ 2 ∂ α-1 ψ * ∂ |x| α-1 + B α |ψ| 2 ∂ α-1 ψ ∂ |x| α-1 + D α |ψ| 2 ψ = 0. ( 35
)
This is a fractional cubic nonlinear Schrödinger-like equation that turns out to be the one governing the dynamics of a ferromagnetic Heisenberg spins chain involving long-range interactions in the framework of the approximation described in Appendix B, i.e. the fields f are slowly varying in space.

However, if the latter approximation is avoided, then new formulas of the terms C and E are retrieved by the computations done in Appendix C. It follows from this analysis that we get a fractional integro-differential cubic nonlinear Schrödinger equation given by

i ∂ψ ∂t + U α ∂ α-1 ψ(x, t) ∂ |x| α-1 + B α ψ 2 ∂ α-1 ψ * ∂ |x| α-1 + B α |ψ(x, t)| 2 ∂ α-1 ψ(x, t) ∂ |x| α-1 + V α |ψ| 2 ψ -Gb α +∞ b dy |ψ(x -y, t)| 2 ψ(x -y, t) + |ψ(x + y, t)| 2 ψ(x + y, t) |y| α -4b α +∞ b dy |ψ(x -y, t)| 2 ψ(x, t) + |ψ(x + y, t)| 2 ψ(x, t) |y| α = 0, (36)
with V α = -5G Ĵα (0) -2I. Equations ( 35) and ( 36) are the main results of the present paper. The approximation done to get the differential Eq. ( 35) corresponds to the excitations with small amplitudes, while Eq. ( 36) can also describe large amplitude excitations. It would be very interesting, but nontrivial, to get their analytical solutions. This task appears not to be straightforward and in the following we focus on the modulational instability of the extended linear solutions of Eq. [START_REF] Bishop | Solitons and Condensed Matter[END_REF].

IV. MODULATIONAL INSTABILITY FOR THE CONTINUUM MEDIUM

The modulational instability of Eq. ( 35) can be studied using the standard method described for example in Refs. [42,43]. Here we are interested in the stability of the homogeneous solution ψ(t) = ψ 0 e iDαψ 2 0 t/ for which the amplitude ψ 0 is a real quantity without loss of generality. We remind that we are in the region of parameters 1 < α < 3, with α = 2.

In the presence of a small perturbation a(x, t) in the system, one can write

ψ(x, t) = (ψ 0 + a)e iDαψ 2 0 t/ , (37)
in which a(x, t) ψ 0 . Substituting the ansatz (37) in Eq. ( 35), we obtain

i ∂a ∂t = -U α ∂ α-1 a ∂ |x| α-1 -B α ψ 0 (ψ 0 + 2a) ∂ α-1 a * ∂ |x| α-1 -B α ψ 0 (ψ 0 + a + a * ) ∂ α-1 a ∂ |x| α-1 -D α ψ 2 0 (a + a * ). ( 38
)
Next, splitting the perturbation term into real and imaginary parts, a(x, t) = u(x, t)+iv(x, t), and linearizing Eq. ( 38) with respect to u and v, it turns out that we get the following system of equations

∂u ∂t = -U α ∂ α-1 v ∂ |x| α-1 , (39)
∂v ∂t = U α ∂ α-1 u ∂ |x| α-1 + 2B α ψ 2 0 ∂ α-1 u ∂ |x| α-1 + 2D α ψ 2 0 u. ( 40
)
By introducing the following Fourier transforms

û(k, t) = +∞ -∞ u(x, t) exp(ixk)dx, (41) v
(k, t) = +∞ -∞ v(x, t) exp(ixk)dx, (42)
Equarion ( 39) is converted into a set of ordinary differential equations in the wavevector k domain,

∂ û ∂t = U α |k| α-1 v, ( 43
)
∂v ∂t = (-U α -2B α ψ 2 0 ) |k| α-1 û + 2D α ψ 2 0 û, (44)
that can be combined into

∂ 2 v ∂t 2 = 2D α ψ 2 0 -(U α + 2B α ψ 2 0 ) |k| α-1 U α |k| α-1 v, (45)
and the same equation for û. Consequently, perturbations can grow if and only if the prefactor of the right-hand-side is positive. In such a case, one can thus define the growth rate of the modulational instability as

G(k, α) = U α |k| α-1 2D α ψ 2 0 -(U α + 2B α ψ 2 0 ) |k| α-1 , (46)
which shows that perturbations with wavevectors in the range of 1) , are exponentially amplified. The growth rate, which depends on the parameters A and J, the wavevectors k and also the long-range interacting parameter α, is definitely the quantity allowing to characterize the regions of stability/instability of the system. In this specific case, the growth rate can be zero, an imaginary number or either a positive real number. If the growth rate is zero or an imaginary number, then, the system displays stability. However if the growth rate is a real and positive number, then, the system may display a modulational instability phenomenon. Therefore, as the modulational instability phenomenon is the main focus of this study, it does occur in specific regions of the values of the parameters. Henceforth, Fig. 1 presents in green the region in the (A, J) plane where the system displays the modulational instability phenomenon. Figure 2 shows the dependence of the modulational instability growth rate on the wavevector k and the fractional index α for particular values of the parameter J. Here, we realize that the larger is the exchange parameter J, the larger is the stable region. Our study is suitable only for non integer values of α between 1 and 3. While looking at different panels of Fig. 2 obtained for the same value of the anisotropy parameter A, it is realised that as the parameter J increases,the range of values of α for which the modulational instability occurs, reduces. Figure 3 presents the evolution of the critical wave-vector as a function of the exponent α for the three cases presented in the previous figure. Similarly, when the parameter J increases as seen in the panels of Fig. 3, the range of values of the critical wavevactor k c reduces. Figures 4 and5 show that when the parameter J is fixed and the parameter A varies, the phenomenon is reversed as compared to what is observed in Fig. 2 and Fig. 3. This means, as seen in the corresponding panels, that as the anisotropy parameter increases, the instability region increases. From this, one sees that if the parameter J tends to favor the stability of the system, the anisotropy parameter at variance induces the instability in such a system. Needless to mention is the fact that while looking at Fig. 2 to Fig. 5, it is realised that the power-law long-range exponent displays a minimal value below which the instability does not occur. This minimal value depends also on both the exchange parameter J and the anisotropy parameter A. For instance, if the anisotropy parameter A is fixed, the minimal value of α increases with increasing values of the exchange parameter J, as seen in panels a) and b) of Fig. 2 and Fig. 3. The dependence of the minimal value of α on the A and J is also observed when J is fixed and the parameter A varies as seen in panels a) and b) of Fig. 4 and Fig. 5 with a reversed effect.

-k 0 < k < k 0 with k c = ±k 0 = ±(2D α ψ 2 0 /(U α + 2B α ψ 2 0 )) 1/(α-

V. CONCLUSION

We have investigated the effective dynamics of a Heisenberg ferromagnetic spin chain with algebraic long-range couplings using a semiclassical approximation. We considered a spin chain with long-range power-law interactions having a strength proportional to 1/(n-m) α in the regime 1 < α < 3 and α = 2. We have used the Holstein-Primakoff representation of the spins to derive a discrete nonlinear cubic Schrödinger equation. We have also shown how one can get, in the continuum limit, on the one hand, a fractional cubic nonlinear Schrödinger-like equation and, on the other hand, an integro-differential fractional cubic nonlinear Schrödinger-like equation, corresponding to excitations with small amplitudes and excitations with large amplitudes, respectively. This has been achieved after performing the analytical derivatives by firstly using the Riesz derivative of fractional calculus and using, secondly, a direct analysis of the Fourier spectrum in the k → 0 limit. A remarkable feature of this interaction is the existence of a transformation that replaces the set of coupled individual bosonic spin equations with a continuum medium equation with a fractional spatial derivative of order α. Such a transformation is an approximation and it appears in the infrared limit for k → 0. We also studied the appearance of modulational instability using the obtained fractional equation. The modulational stability has been studied for the fractional equation and the stable and unstable regions have been determined. We showed the shifting of the onset of the modulational instability region for α > 1 and we studied the dependence of the regions of instability on the parameters of the system. It clearly appears that the modulational instability is present only for α < 3 and that above α = 3 only stable regions are present, thus indicating that the system behaves as a short range system. It is also realized that the parameter J tends to favorise the stability of the system while the anisotropy parameter A tends to induce the instability in the system. As a future work, it would be interesting to compare the results obtained using the effective continuum equation with those of the original lattice model. In particular it would be rewarding to compare the modulational instability regions obtained with the lattice equation found for large-S with the findings presented here, since in general the lattice model is expected to have a larger region of instabilities. To what extent some unstable phases are missed by the study of the continuum fractional equation is a subject deserving future analytical and numerical studies. Finally, we observe that the approach presented here does not apply to α = 2 and it would be useful to study in detail the α → 2 limit.

Generalizing this one-dimensional study to 2D and 3D lattices would be very interesting and important. However, increasing the dimension of the lattice very seriously complicates the derivation of effective equations in the continuum and the search of the solutions. Problems solved in 1D have waited many years, if not decades, to be treated successfully in higher dimensions. However, this would be very helpful and important to check using numerical simulations whether modulational instability thresholds qualitatively depends on the dimension of the lattice.

Setting then m = m -n, we get the following expression

A 1 = 1 2π +∞ -∞ +∞ m =-∞,m =0 J |m | α f (k, t) e ikn -e ik(m +n) dk (B5) = 1 2π   +∞ -∞ +∞ m =-∞,m =0 J |m | α f (k, t)e ikn dk - +∞ -∞ +∞ m =-∞ Je ikm |m | α f (k, t)e ikn dk   (B6) = 1 2π   +∞ -∞ Ĵα (0) f (k, t)e ikn dk - +∞ -∞ Ĵα (k) f (k, t)e ikn dk)   (B7) = 1 2π +∞ -∞ ( Ĵα (0) -Ĵα (k)) f (k, t)e ikn dk, (B8)
where we set x = k∆x and Ĵα (x

) = +∞ m =-∞,m =0
Je im x / m α . It turns out that we get

A 1 = - a α 2π +∞ -∞ k α-1 f (k, t)e ikn dk. (B9)
Since the Fourier transform involving the absolute value of momentum |k| α is expressed by Riesz derivative in the real space as

1 2π +∞ -∞ |k| α f (k, t)e ikn dk - ∂ α f (x, t) ∂ |x| α (B10)
with x = n∆x, the term A 1 finally reads

A 1 = a α ∂ α-1 f (x, t) ∂ |x| α-1 . (B11)
The term C can be split as

C = n-1 m=-∞,m =n J |m -n| α |f m | 2 f m C1 + +∞ m=n+1,m =n J |m -n| α |f m | 2 f m C2 (B12)
that can be analytically computed separately. Then, we get

C 1 = 1 2π 3 +∞ -∞ n-1 m=-∞,m =n Je ikm |m -n| α f (k, t)dk +∞ -∞ f * (k , t)e -ik (m) dk +∞ -∞ f (k , t)e ik (m) dk (B13)
and next setting m -n = m , it turns out that

C 1 = 1 2π 3 +∞ -∞ -1 m =-∞,m =0 Je ikm |m | α f (k, t)e ikn dk +∞ -∞ f * (k , t)e -ik (m +n) dk +∞ -∞ f (k , t)e ik (m +n) dk (B14) = 1 2π 3 +∞ -∞ +∞ m =1,m =0 Je -i(k-k +k )m |m | α f (k, t)e ikn f * (k , t)e -ik n f (k , t)e ik n dkdk dk (B15) = 1 2π 3 +∞ -∞ +∞ m =1,m =0 Je -ikm ∆x |m | α f (k, t)e ikn dk +∞ -∞ f * (k , t)e -ik (-m +n) dk +∞ -∞ f (k , t)e ik (-m + dk (B16) C 1 = 1 2π +∞ -∞ +∞ m =1,m =0 Je -ikm |m | α f (k, t)e ikn dkf * n-m f n-m . (B17)
Before moving to the continuum limit, we should remind that there are terms with sums of type m f n±m . In this respect, we can proceed to the approximation based on the assumption that in these sums, the fields f n±m 's are slowly varying in space. Therefore, such terms as f n±m can be brought outside the sums over m [26]. In this case, we can now write

C 1 = f * n-m f n-m 1 2π +∞ -∞ +∞ m =1,m =0 Je -ikm |m | α f (k, t)e ikn dk. (B18)
Following the same later steps, we get

C 2 = 1 2π 3 +∞ -∞ +∞ m=n+1,m =n Je ikm |m -n| α f (k, t)dk +∞ -∞ f * (k , t)e -ik (m) dk +∞ -∞ f (k , t)e ik (m) dk (B19) = 1 2π 3 +∞ -∞ +∞ m =1,m =0 Je ikm |m | α f (k, t)e ikn dk +∞ -∞ f * (k , t)e -ik (m +n) dk +∞ -∞
f (k , t)e ik (m +n) dk , (B20) that can be rewritten as

C 2 = 1 2π +∞ -∞ +∞ m =1,m =0 Je ikm |m | α f (k, t)e ikn dkf * n+m f n+m . (B21)
Then, if we consider here also the same approximation used above to move from Eq. (B17) to Eq. (B18), we get

C 2 = f * n+m f n+m 1 2π +∞ -∞ +∞ m =1,m =0 Je -ikm |m | α f (k, t)e ikn dk. (B22)
It then follows that C is given by

C = f * n-m f n-m 1 2π +∞ -∞ +∞ m =1,m =0 Je -ikm |m | α f (k, t)e ikn dk + f * n+m f n+m 1 2π +∞ -∞ +∞ m =1,m =0 Je ikm |m | α f (k, t)e ikn dk (B23) As n ± m n, we get C = f * n f n 1 2π +∞ -∞ +∞ m =1,m =0 Je -ikm |m | α f (k, t)e ikn dk + f * n f n 1 2π +∞ -∞ +∞ m =1,m =0 Je ikm |m | α f (k, t)e ikn dk (B24)
which leads us to

C = f * n f n 1 2π +∞ -∞ +∞ m =1,m =0 Je -ikm + Je ikm ∆x |m | α f (k, t)e ikn dk (B25) = f * n f n 1 2π +∞ -∞ Ĵα (k) f (k, t)e ikn dk, (B26)
as

Ĵα (k) = +∞ m =1,m =0 Je -ikm + Je ikm |m | α . ( B27
)
Here we consider the continuum approximation done in Eq. ( 24) to get

C f * n f n 1 2π +∞ -∞ a α |k| α-1 + Ĵα (0) f (k, t)e ikn∆x dk (B28) = f * n f n 1 2π +∞ -∞ a α |k| α-1 f (k, t)e ikn dk + f * n f n 1 2π +∞ -∞ Ĵα (0) f (k, t)e ikn dk. (B29)
Considering that f n = f (x = n∆x, t), the continuum approximation using the Riesz derivative gives by Eq. ( 21) leads us to

C = -a α ∂ α-1 f (x, t) ∂ |x| α-1 |f (x, t)| 2 + Ĵα (0) |f (x, t)| 2 f (x, t). (B30)
Then splitting E as follows,

E = +∞ m=-∞,m =n J mn |f m | 2 f n = n-1 m=-∞,m =n J mn |f m | 2 f n E1 + +∞ m=n+1,m =n J mn |f m | 2 f n E2 , (B31)
one has to calculate each term.

E 1 = n-1 m=-∞,m =n J mn |f m | 2 f n (B32) = 1 2π 3 +∞ -∞ n-1 m=-∞,m =n J |m -n| α f (k, t)e ikm dk +∞ -∞ f * (k , t)e -ik m dk +∞ -∞ f (k , t)e ik n dk . (B33)
now, setting m -n = m , we get

E 1 = 1 2π 3 +∞ -∞ -1 m =-∞,m =0 Je ikm |m | α f (k, t)e ikn dk +∞ -∞ f (k , t)e -ik (m +n) dk +∞ -∞ f (k , t)e ik n dk (B34) = 1 2π 3 +∞ -∞ -1 m =-∞,m =0 Je i(k-k )m |m | α f (k, t)e ikn dk +∞ -∞ f (k , t)e -ik n∆x dk +∞ -∞ f (k , t)e ik n dk (B35) = 1 2π 3 +∞ -∞ +∞ m =1,m =0 Je -i(k-k )m |m | α f (k, t)e ikn dk +∞ -∞ f (k , t)e -ik n∆x dk +∞ -∞ f (k , t)e ik n dk (B36) = 1 2π 3 +∞ -∞ +∞ m =1,m =0 Je -ikm |m | α f (k, t)e ikn dk +∞ -∞ f (k , t)e -ik (n-m ) dk +∞ -∞ f (k , t)e ik n dk (B37) = 1 2π +∞ -∞ +∞ m =1,m =0 Je -ikm |m | α f (k, t)e ikn dkf * n-m f n . (B38)
Then, proceeding to the same approximation used to move from Eq. (B17) to Eq. (B18), we get

E 1 = 1 2π f * n-m f n +∞ -∞ +∞ m =1,m =0 Je -ikm |m | α f (k, t)e ikn dk. (B39)
Through the same approximation and very similar steps, one derives When we replace all these terms in the equation of motion (12), we get (C4)

E 2 = 1 2π f * n+m f n +∞ -∞ +∞ m =1,m =0 Je ikm |m | α f (k, t)e ikn dk (
i ḟ = F a α ∂ α-1 f (x, t) ∂ |x| α-1 -T f + 2I |f | 2 f +G -a α ∂ α-1 f * ∂ |x| α-1 + Ĵα (0) |f | 2 f -a α |f | 2 ∂ α-1 f ∂ |x| α-
Next, while introducing the following Ansatz:

f (x, t) = ψ(x, t) e -i2ε 2 At/ (C5)

We get the integro-fractional differential nonlinear Schrödinger equation given by Eq. ( 36).

FIG. 1 .FIG. 2 .

 12 FIG. 1. Here the green region corresponds the region in the (A, J) diagram where the modulational instability could occur.

FIG. 3 .

 3 FIG. 3. The critical frequency kc as a function of the exponent α with J = 5 (panel a), J = 15 (b) and ψ0 = 0.1. Both panels are obtained for A=160.

FIG. 4 .

 4 FIG. 4. The growth rate of the modulational instability as a function of α for differents values of A = 60 (panel a) and A = 160 (b) while the exchange parameter is fixed to J = 5 and ψ0 = 0.1. The blue region corresponds to the stable domain while the multicolor region corresponds to the unstable one.

FIG. 5 .

 5 FIG. 5. The critical frequency kc as a function of the exponent α with A = 60 (panel a), A = 160 (b) and ψ0 = 0.1. Both panels are obtained for J=5.

  k) α-1 f (k, t)e ikn dk + f * n f n f (k, t)e ikn dk.(B44)

1 - 1 - 1 + 5

 1115 1 + Ĵα (0) |f | 2 f + 4 Ĵα (0) |f | 2 f -4 Ĵα (0) |f | 2 f . (B45) Appendix C: The fractional integro-differential nonlinear equation |m | α f (k, t)e ikn dkf * n-m f n-m t)dk. (C1)In the continuum limit ∆x → 0, by replacing n and m respectively by x/∆x and y/∆x, we get |f(x -y, t)| 2 f (x -y, t) + |f (x + y, t)| 2 f (x + y, t) |y| α . (C2)where b = ∆x and p = k/b. In the same spirit of analytical computations, the term E will be given by -y, t)| 2 f (x, t) + |f (x + y, t)| 2 f (x, t) following fractional integro-differential cubic nonlinear equation given byi ḟ = F a α ∂ α-1 f (x, t) ∂ |x| α-Ga α f 2 ∂ α-1 f * ∂ |x| α-Ga α |f (x, t)| 2 ∂ α-1 f (x, t) ∂ |x| α-Ĵα (0)f |f | 2 -T f + 2I |f | 2 f +G   b α +∞ b dy |f (x -y, t)| 2 f (x -y, t) + |f (x + y, t)| 2 f (x + y, t) -y, t)| 2 f (x, t) + |f (x + y, t)| 2 f (x, t)

ACKNOWLEDGMENTS

The authors thank CNRS for financial support through the project DyFraLonPo of the "Dispositif de soutien aux collaborations avec l'Afrique Subsaharienne".

Appendix A: Properties of the function J One thus has to compute the function

Introducing the poly-logarithm function defined as

in which ζ stands for the zeta function defined as

1/n α and z = ik∆x, it turns out that we get

Denoting Ĵα (0) = 2Jζ (α) and a α = 2JΓ(1 -α) sin (πα/2), we thus get

In the framework of the continuum approximation [START_REF] Herman | Fractional Dynamics[END_REF], we get the expression

For non-integer α > 3 [START_REF] Tarasov | Fractional Dynamics: Applications of Fractional Calculus to Dynamics of Particles, Fields and Media[END_REF], we have

Appendix B: Discussion of the different terms of the discrete nonlinear equation

We use the notation

where F{f (x, t)} is the Fourier transform of f (x, t) with respect to x.

In that case, the term A 1 defined in Eq. ( 18) reads


Thierry Daudé
email: thierry.daude@univ-fcomte.fr

Niky Kamran
email: niky.kamran@mcgill.ca

François Nicoleau
email: francois.nicoleau@univ-nantes.fr

Hölder Stability in the Inverse Steklov Problem for Radial Schrödinger operators and Quantified Resonances

Keywords: Inverse Steklov problem, Steklov spectrum, Weyl-Titchmarsh functions, moment problems, Hölder stability. 2010 Mathematics Subject Classification. Primaries 81U40, 35P25; Secondary 58J50. Contents

. It is noteworthy that with respect to the original Schrödinger operator, the type of perturbation being considered for the amplitude function amounts to the introduction of a finite number of negative eigenvalues and of a countable set of negative resonances which are quantified explicitly in terms of the eigenvalues of the Laplace-Beltrami operator on the boundary sphere.

Introduction 2 Introduction

In a recent paper, [START_REF] Daudé | Stability in the inverse Steklov problem on warped product Riemannian manifolds[END_REF], we have obtained a set of logarithmic stability estimates in the inverse Steklov problem for the Laplace-Beltrami operator on a class of warped product Riemannian manifolds defined on a d-dimensional closed ball. These manifolds can be thought of as deformations of the closed Euclidean dball in which the deformation is parametrized by the choice of radial warping function. The deformations considered in [START_REF] Daudé | Stability in the inverse Steklov problem on warped product Riemannian manifolds[END_REF] include both the regular and singular cases of warped product metrics [START_REF] Petersen | Riemannian Geometry, Third Edition[END_REF].

The approach taken in [START_REF] Daudé | Stability in the inverse Steklov problem on warped product Riemannian manifolds[END_REF] was based on expressing the warped product metric in a coordinate system in which the metric takes the form of a conformal rescaling of the flat Euclidean metric. This enabled us by using the transformation law of the Laplace-Beltrami operator under conformal changes of the metric to reformulate the inverse Steklov problem for the Laplace-Beltrami operator on the original deformed closed ball as the inverse Steklov problem for a Schrödinger operator on the Euclidean ball, with a potential determined by the warping function of the original deformed ball and its derivatives. The logarithmic stability estimates that we obtained were thus of the nature of the estimates obtained by Alessandrini [START_REF] Alessandrini | Stable determination of conductivity by boundary measurements[END_REF] and Novikov [START_REF] Novikov | New global stability estimates for the Gelfand-Calderón inverse problem[END_REF].

Our goal in the present paper is to improve the logarithmic stability results of [START_REF] Daudé | Stability in the inverse Steklov problem on warped product Riemannian manifolds[END_REF] in a significant way by obtaining instead a set of Hölder stability estimates and by highlighting the role of the resonances. In the approach we will take in the current paper, instead of starting from a deformed closed d-ball, we will reverse the initial step taken in [START_REF] Daudé | Stability in the inverse Steklov problem on warped product Riemannian manifolds[END_REF] and start with a Schrödinger operator on the Euclidean d-ball, endowed with a radial potential. One reason for doing so is that the set of warped product metrics for which Hölder stability estimates can be obtained for the inverse Steklov problem will be significantly more restricted than the ones for which logarithmic stability estimates were obtained in [START_REF] Daudé | Stability in the inverse Steklov problem on warped product Riemannian manifolds[END_REF]. Thus by considering Schrödinger operators, we will be broadening the range of inverse Steklov problems to which our Hölder stability results will apply. We emphasize nevertheless that the class of potentials for which we are able to establish Hölder stability for the inverse Steklov problem is still rather special and that any extension of our results to more general potentials may require stronger techniques than the ones we are using the present paper. We now proceed to describe the main results of our paper.

More precisely, on the d-dimensional closed Euclidean ball

M = (0, 1] × S d-1 , (1.1)
where d ≥ 3, we consider the Dirichlet problem for the Schrödinger operator with a potential q, given by

-△u + q u = 0, on M , u = ψ ∈ H 1 2 (∂M ), on ∂M , (1.2)
When q ∈ L ∞ (M ) and λ = 0 is not a Dirichlet eigenvalue of the above Schrödinger operator, the Dirichlet problem (1.2) has a unique solution u ∈ H 1 (M ). The Dirichlet-to-Neumann (DN) map Λ q is then (formally) defined as an operator from H 1/2 (∂M ) to H -1/2 (∂M ) by

Λ q ψ = (∂ ν u) |∂M , (1.3)
where u is the unique solution of (1.2) and (∂ ν u) |∂M is the normal derivative of u with respect to the outer unit normal vector ν on ∂M .

The DN map thus defined is a self-adjoint operator on L 2 (∂M, dS g ) where dS g denotes the metric induced by the Euclidean metric on the boundary sphere ∂M = S d-1 . Its spectrum (the so-called Steklov spectrum) is discrete and accumulates at infinity We shall thus denote the Steklov eigenvalues (counted with multiplicity)

0 = σ 0 < σ 1 ≤ σ 2 ≤ • • • ≤ σ k → ∞ . (1.4)
The Steklov spectrum will be the central object of study in this paper.

For the remainder of this paper, we shall assume that the potential is radial and we write q = q(r), where r denotes the Euclidean distance to the origin. It will be convenient to replace the radial coordinate r ∈ (0, 1] by a new radial coordinate x ∈ (0, ∞) defined by x = -log r, in which case the the boundary of M now corresponds to x = 0. The Euclidean metric then takes the form

g = f (x) 4 (dx 2 + dΩ 2 ) ,
where f (x) = exp(-x/2) and dΩ 2 denotes the round metric on the unit sphere S d-1 . The Dirichlet problem (1.2) gets transformed into

[-∂ 2 x -△ S + Q(x)]v = -(d-2) 2 4 v, on M , v = f d-2 ψ, on ∂M , (1.5)
where △ S denotes the Laplacian on the boundary sphere S d-1 and Q(x) := e -2x q(e -x ).

Thanks to the spherical symmetry of the potential Q, we can use separation of variables and reduce (1.5) to an infinite sequence of radial ordinary differential equations. We thus let {Y k , k ≥ 0} denote an orthonormal Hilbert basis of L 2 (S d-1 ) consisting of eigenfunctions of △ S ,

-△ S Y k = α k Y k , α k = k(k + d -2) . Writing v = k≥0 v k (x) Y k , (1.6)
we obtain an infinite sequence of ordinary differential equations on (0, ∞) given by

-v ′′ k + Qv k = -(α k + (d -2) 2 4 )v k = -κ 2 k v k , (1.7)
where

κ k := k + d -2 2 , k ≥ 0 .
As in [START_REF] Daudé | Stability in the inverse Steklov problem on warped product Riemannian manifolds[END_REF], we introduce the Weyl-Titchmarsh (WT) function M (z) associated to the Sturm-Liouville operator L on (0, +∞) given by

L = - d 2 dx 2 + Q . (1.8)
This function will play a central role in our subsequent analysis of the stability problem of the Steklov spectrum for our Schrödinger operator. We assume in this paper that

Q ∈ L 2 (0, ∞) . (1.9)
Under this assumption, it is well-known that L is of limit point-type at infinity, which means that for all z ∈ C \ [-β, ∞) with β >> 1, there exists, up to a zon-zero multiplicative constant, a unique solution

u(x, z) of -u ′′ + Qu = zu z ∈ C , (1.10)
which is L 2 at ∞. The Weyl-Titchmarsh function M (z) is then defined by

M (z) := u ′ (0, z) u(0, z) for all z ∈ C \ [-β, ∞) . (1.11)
Of course, the square-integrability hypothesis (1.9) we made on the potential Q does not necessarily imply that the initial potential q ∈ L ∞ (M ). Thus, the previous definition we gave for the DN map is not directly applicable in this L 2 setting. We overcome this difficulty by exploiting the separation of variables and follow the procedure used in Section 2 of [START_REF] Daudé | Stability in the inverse Steklov problem on warped product Riemannian manifolds[END_REF] to define the DN map; namely we expand the boundary data ψ in the Hilbert basis

{Y k } k≥0 of L 2 (S d-1 ) as ψ = k≥0 ψ k Y k ,
and define the DN map Λ q as a sum of operators Λ k q by

Λ q ψ = ∞ k=0 (Λ k q ψ k )Y k .
(1.12)

The "diagonalized" DN operators Λ k q are then computed using the separation of variables to give

Λ k q ψ k = - (d -2) 2 v k (0) -v ′ k (0) , (1.13)
where

-v ′′ k + Qv k = -κ 2 k v k v k (0) = ψ k , v k ∈ L 2 at + ∞ .
As was proved in [START_REF] Daudé | Stability in the inverse Steklov problem on warped product Riemannian manifolds[END_REF], the operators Λ k q , k ≥ 0 can be further simplified to take the form of multiplication operators by making use of the Weyl-Titchmarsh function M evaluated at the points -κ 2 k :

Λ k q ψ k = - (d -2) 2 -M (-κ 2 k ) ψ k , (1.14)
and the Steklov spectrum {σ k , k ≥ 0} is given in terms of the Weyl-Titchmarsh function M as

σ k = - (d -2) 4 -M (-κ 2 k ) . (1.15)
There is an important representation formula first obtained in [START_REF] Simon | A new approach to inverse spectral theory, I. Fundamental formalism[END_REF] for the Weyl-Titchmarsh function in terms of a the Laplace transform of a unique amplitude function A, under the hypothesis that Q ∈ L 1 (0, ∞):

M (-κ 2 ) = -κ - ∞ 0 A(α)e -2κα dα , ∀κ > 1 2 ||Q|| 1 . (1.16)
We shall use a slightly refined version of this formula which applies in our L 2 -setting and which will serve as the starting point of our formulation of the stability problem for the Steklov spectrum.

Let us now explain what we mean precisely by local stability estimates. The starting point is a fixed potential Q ∈ L 2 (0, ∞) which we perturb through the addition of a certain exponential series to its corresponding amplitude A. The set of admissible exponential series parametrizing the perturbations lies in an infinite-dimensional space. Using powerful results of Killip -Simon [START_REF] Killip | Sum rules and spectral measures of Schrödinger operators with L 2 potentials Ann[END_REF] we then show that Q ∈ L 2 (0, ∞).

Given now such a pair of potentials Q, Q, we assume that the difference between their corresponding Steklov spectra is uniformly bounded in absolute value by a small error ǫ > 0,

|σ k -σk | < ǫ , ∀k ≥ 0 . (1.17)
Our main goal is to estimate the difference Q -Q of these potentials. Again, our result is local: for any fixed parameter T > 0, we get stability estimates in the space L 2 (0, T ). Roughly speaking, this means that

|| Q -Q|| L 2 (0,T ) ≤ C T g(ǫ) , (1.18)
where g(ǫ) → 0 when ǫ → 0, and C T is a constant only depending on T . Now, we can state our main result in this paper :

Theorem 1.1. Let Q ∈ L 2 (0, ∞)
be a square-integrable potential with amplitude function A, and let δ ≥ 3 -d be any fixed parameter. Set µ k := λ k + δ where λ k = 2k + d -3 + δ and let {c k , k ≥ 0} be a sequence of real numbers such that

• i) c k ≤ 0 for all k ≥ 0.
• ii) The power series k≥0 c k t λ k has a radius of convergence R > 1.

Then the function à defined by

Ã(α) = A(α) + k≥0 c k e -µ k α , α > 0 , (1.19)
is the amplitude function of a potential Q ∈ L 2 (0, ∞).

Moreover, under the hypothesis (1.17), for any fixed T > 0, there exists a positive constant C T such that

|| Q -Q|| L 2 (0,T ) ≤ C T ǫ θ , (1.20)
where the Hölder exponent θ ∈ (0, 1 2 ) is independent of T and is given by

θ = 1 2 min(1, log R log( 9M0 2 ) ), M 0 = max{2, 4(d -3 + δ) + 1} . (1.21)
Note that when the radius of convergence R ≥ 9M0 2 , we can take as Hölder exponent θ = 1 2 . Note also that when the initial potential is the trivial potential Q = 0, the perturbed potentials Q can be seen as a generalization of the so-called Bargmann potentials (see section 5 for details).

Finally, it is noteworthy that with respect to the original Schrödinger operator, the type of perturbation being considered for the amplitude function A amounts to the introduction of a finite number of negative eigenvalues - 2 which are equally spaced on the negative real axis (for k greater than some k 0 ). These resonances are quantified explicitly in terms of the parameter δ and the eigenvalues of the Laplace Beltrami operator ∆ S on the boundary sphere.

Notation and set-up of the model

On the d-dimensional closed Euclidean ball

M = (0, 1] × S d-1 , (2.22)
where d ≥ 3, we consider the Dirichlet problem for the Schrödinger operator with potential q ∈ L ∞ (M ), given by -△u + q u = 0, on M , u = ψ ∈ H 1 2 (∂M ), on ∂M .

(2.23)

Although we shall shortly make several assumptions about the potential q (including the requirement that it be radial), we begin by recalling a few general facts concerning the Dirichlet problem (2.23) for a general potential q ∈ L ∞ (M ). These will lead us to the definition of the Dirichlet-to-Neumann (DN) map and the associated Steklov spectrum, which will be central objects of study in this paper.

We first recall (see for example Theorem 8.3 in [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF]) that if q ∈ L ∞ (M ) and zero is not a Dirichlet eigenvalue of the above Schrödinger operator (which is the case for example if q ≥ 0), then the Dirichlet problem (2.23) has a unique solution u ∈ H 1 (M ). The Dirichlet-to-Neumann (DN) map Λ q is then defined as an operator from H 1/2 (∂M ) to H -1/2 (∂M ) by

Λ q ψ = (∂ ν u) |∂M , (2.24)
where u is the unique solution of (2.23) and (∂ ν u) |∂M is the normal derivative of u with respect to the outer unit normal vector ν on ∂M . Here (∂ ν u) |∂M is generally defined in the weak sense as an element of H -1/2 (∂M ) by

Λ q ψ|φ = M du, dv dV ol ,
for any ψ ∈ H 1/2 (∂M ) and φ ∈ H 1/2 (∂M ) such that u is the unique solution of (2.23) and v is any element of H 1 (M ) such that v |∂M = φ. It is easily checked that if ψ is sufficiently smooth, we have

Λ q ψ = g(ν, ∇u) |∂M = du(ν) |∂M = ν(u) |∂M ,
so that the expression in local coordinates for the normal derivative is then given by

∂ ν u = ν i ∂ i u . (2.25)
It is well known that the DN map is a pseudo-differential operator of order 1 which is self-adjoint on L 2 (∂M, dS g ) where dS g denotes the metric induced by the Euclidean metric on the boundary ∂M = S d-1 . Therefore, the DN map Λ q has a real and discrete spectrum accumulating at infinity, known as the Steklov spectrum. We shall denote the Steklov eigenvalues (counted with multiplicity) by

0 = σ 0 < σ 1 ≤ σ 2 ≤ • • • ≤ σ k → ∞ . (2.26)
We refer the reader to [START_REF] Girouard | Spectral geometry of the Steklov spectrum[END_REF] and references therein for an excellent survey of the known results on the Steklov spectrum.

As of now and for the remainder of this paper, we shall assume that the potential is radial and write q = q(r), where r denotes the Euclidean distance to the origin. It will be convenient for our subsequent analysis to replace the radial coordinate r ∈ (0, 1] by a new radial coordinate x ∈ (0, ∞) defined by x = -log r, in which case the the boundary of M now corresponds to x = 0. The Euclidean metric

g = dr 2 + r 2 dΩ 2 ,
where dΩ 2 denotes the round metric on the unit sphere S d-1 , then takes the form

g = f (x) 4 (dx 2 + dΩ 2 ) ,
where f (x) = exp(-x/2), and the Dirichlet problem (2.23) gets transformed into

[-∂ 2 x -△ S + Q(x)]v = -(d-2) 2 4 v, on M , v = f d-2 ψ, on ∂M , (2.27)
where △ S denotes the Laplacian on the boundary sphere S d-1 and Q(x) := e -2x q(e -x ).

Separation of variables

Our first step in the analysis of (2.27) is to exploit the spherical symmetry of the potential Q in order to separate variables and reduce (2.27) to an infinite sequence of radial ordinary differential equations. We thus let {Y k , k ≥ 0} denote an orthonormal Hilbert basis of L 2 (S d-1 ) consisting of eigenfunctions of △ S ,

-△ S Y k = α k Y k , α k = k(k + d -2) ,
where the eigenvalue α k is of multiplicity k+d-1 d-1

-k+d-3 d-1 , and let

κ k := k + d -2 2 , k ≥ 0 .
We separate variables by letting

v = k≥0 v k Y k , (2.28)
where v k = v k (x) depends only on the radial variable x ∈ (0, ∞), which gives rise upon substitution of (2.28) into (2.27) to the infinite sequence of ordinary differential equations on (0, ∞) given by

-v ′′ k + Qv k = -(α k + (d -2) 2 4 )v k = -κ 2 k v k . (2.29)
For our subsequent analysis of the stability problem, it will be very useful to complexify the spectral parameter κ k and consider instead the ordinary differential equation

-v ′′ + Qv = zv z ∈ C , (2.30)
where z is now to be thought of as a complex spectral parameter.

The Weyl-Titchmarsh function

We now introduce the Weyl-Titchmarsh function M (z) associated to the Sturm-Liouville operator

L = - d 2 dx 2 + Q , (2.31)
defined by the left-hand side of (2.30). This function will play a central role in our subsequent analysis of the stability problem of the Steklov spectrum for our Schrödinger operator.

We first recall that in order for the Weyl-Titchmarsh function to be well-defined, we need to assume that L is of limit point-type at infinity, meaning that for all z ∈ C \ [-β, ∞) with β >> 1, there exists, up to a zon-zero multiplicative constant, a unique solution u(x, z) of (2.30) which is L 2 at ∞. The Weyl-Titchmarsh function M (z) is then defined by

M (z) := u ′ (0, z) u(0, z) for all z ∈ C \ [-β, ∞) . (2.32)
It is easy to show that this property will be guaranteed if

Q ∈ L 2 (0, ∞) , (2.33)
a property which we will require Q to satisfy from now onwards. Indeed, from (1.5) in [START_REF] Simon | A new approach to inverse spectral theory, I. Fundamental formalism[END_REF], we know that L will be of limit point-type at infinity if which shows that (2.33) ensures indeed the property that L is of limit point-type at infinity,

β 2 := sup y>0 y+1 y max{Q(x), 0}dx < ∞ . ( 2

The DN map and the Weyl-Titchmarsh function

We begin by remarking that the square-integrability hypothesis (2.33) made on the potential Q appearing in the radial Sturm-Liouville operator L does not necessarily imply that the potential q of the initial Schrödinger operator considered in the set-up of Section 2 will satisfy the condition q ∈ L ∞ (M ), which we assumed in order to define the DN map Λ q . In other words, the definition we gave for the DN map may not be applicable in the L 2 setting defined by the condition (2.33).

We may nevertheless circumvent this difficulty by exploiting the separation of variables worked out in Section 2.1 and following the procedure used in Section 2 of [START_REF] Daudé | Stability in the inverse Steklov problem on warped product Riemannian manifolds[END_REF] to define the DN map. Indeed, by expanding the boundary data ψ in the Hilbert basis {Y k } k≥0 of L 2 (S) as

ψ = k≥0 ψ k Y k , the DN map Λ q is represented as a sum of operators Λ k q by Λ q ψ = ∞ k=0 (Λ k q ψ k )Y k , (2.36)
where the "diagonalized" DN operators Λ k q are computed using the separation of variables to give (see (2.16) in [START_REF] Daudé | Stability in the inverse Steklov problem on warped product Riemannian manifolds[END_REF])

Λ k q ψ k = - (d -2) 2 v k (0) -v ′ k (0) , (2.37)
where

-v ′′ k + Qv k = -(µ k + (d -2) 2 4 )v k = -κ 2 k v k v k (0) = ψ k , v k ∈ L 2 at + ∞ .
As shown in the calculation leading from (2.14) to (2.15) in [START_REF] Daudé | Stability in the inverse Steklov problem on warped product Riemannian manifolds[END_REF], the form given in (2.37) for the operators Λ k q , k ≥ 0 can be further simplified to obtain an expression of these operators as multiplication operators:

Λ k q ψ k = - (d -2) 2 -M (-κ 2 k ) ψ k . (2.38)
Finally, we recall from Lemma 2.3 in [START_REF] Daudé | Stability in the inverse Steklov problem on warped product Riemannian manifolds[END_REF] the explicit formula giving the Steklov spectrum {σ k , k ≥ 0} in terms of the Weyl-Titchmarsh function M ,

σ k = - (d -2) 4 -M (-κ 2 k ) . (2.39)
The latter formula will serve as the starting point of our formulation of the stability problem for the Steklov spectrum.

The amplitude A of a radial potential and the Weyl-Titchmarsh function

As stated in the Introduction, central to our analysis of the stability problem lies a remarkable representation formula first obtained in [START_REF] Simon | A new approach to inverse spectral theory, I. Fundamental formalism[END_REF] (under the hypothesis that Q ∈ L 1 (0, ∞)) for the Weyl-Titchmarsh function the Sturm-Liouville operator L in terms of a the Laplace transform of an amplitude function A. We shall be using a slightly refined version of this formula which applies to the class of L 2 potentials considered in our paper.

In order to state this formula, we first recall from Theorem 2.1 in [START_REF] Simon | A new approach to inverse spectral theory, I. Fundamental formalism[END_REF] that if Q ∈ L 1 (0, ∞), then the Weyl-Titchmarsh function M may be expressed in the form of the Laplace transform of an amplitude function A by

M (-κ 2 ) = -κ - ∞ 0 A(α)e -2κα dα , ∀κ > 1 2 ||Q|| 1 .
(2.40)

It was proved in [START_REF] Gesztesy | A new approach of inverse spectral theory, II. General potentials and the connection to the spectral measure[END_REF] (in the remark following (1.17)) that the above equality also holds for all κ ∈ C such that Re κ >

M (-κ 2 ) = -κ - ∞ 0 A(α)e -2κα dα , ∀ Re κ > 2 max{ 2β 2 , eβ 2 } .
(2.41)

3 The problem of stability

Statement of the problem and strategy

The stability problem may be stated as follows in general terms: Given a pair of potentials Q, Q such that the difference between their corresponding Steklov spectra is uniformly bounded in absolute value by a small error ǫ > 0,

|σ k -σk | < ǫ , ∀k ≥ 0 , (3.1)
what can we say about the difference Q -Q of these potentials? As a first step, we can use the expression (2.39) of the Steklov spectrum in terms of the Weyl-Titchmarsh function M and Simon's representation formula (2.41) for M in terms of the Laplace transform of the amplitude function A to reformulate the condition (3.1) in terms of A. We have

|σ k -σk | = |M (-κ 2 k ) -M (-κ 2 k )| = | ∞ 0 (A(α) -Ã(α))e -2κ k α dα| ,
and making the change of variables α = -log t, our hypothesis (3.1) on the difference of the Steklov spectra takes the form

| 1 0 t -δ A(-log t) -Ã(-log t) t 2k+d-3+δ dt| < ǫ , (3.2)
where δ will be a fixed real parameter that will be properly chosen later. We can see that (3.2) is effectively a Hausdorff moment problem, and thus one should not expect better stability results than the logarithmic stability estimates of the type obtained in [START_REF] Daudé | Stability in the inverse Steklov problem on warped product Riemannian manifolds[END_REF] for the Steklov spectra of deformed balls, or in [START_REF] Alessandrini | Stable determination of conductivity by boundary measurements[END_REF] and [START_REF] Novikov | New global stability estimates for the Gelfand-Calderón inverse problem[END_REF] for the Steklov spectra of certain Schrödinger operators. Nevertheless, as we shall explain in Section 3.2 below, one can approach the stability problem from a different starting point by working directly with perturbations of Simon's amplitude function A by a certain families of exponential series obtained from power series of Müntz type. We shall see that his leads in turn to Hölder type stability results which are significantly stronger than the logarithmic stability results mentioned earlier, albeit at the cost of restricting the class of potentials to a somewhat small subset of the set of square integrable potentials. More precisely, given the amplitude function A associated to a square-integrable potential function Q by (2.41), our strategy will consist in defining a perturbed amplitude function à as in (3.5) and then use an important stability estimate due to Still (Theorem 2 in [START_REF] Still | On the approximation of Müntz series by Müntz polynomials[END_REF]) to obtain a Hölder estimate on à -A. This will be the substance of Section 3.2. The next step, worked out in Section 4, will consist in using the powerful results of [START_REF] Killip | Sum rules and spectral measures of Schrödinger operators with L 2 potentials Ann[END_REF] to construct an L 2 potential Q associated to the perturbed amplitude Ã. Finally, we shall use the methods of boundary control theory of [START_REF] Avdonin | The boundary control approach to the Titchmarsh-Weyl m-function. I. The response operator and the A-amplitude[END_REF] to estimate the difference Q -Q.

It is noteworthy that with respect to the original Schrödinger operator, the type of perturbation being considered for the amplitude function A amounts to the introduction of a finite number of negative eigenvalues (corresponding to the choice of a negative real parameter δ) and of a countable set of resonances on the negative real axis, which admit a precise quantitative expression through to the eigenvalues µ k of the Laplace Beltrami operator ∆ S on the boundary sphere.

Improved stability by Still's method -first main result

Let Q ∈ L 2 (0, ∞) be a given fixed potential and let A be the corresponding amplitude function, given by the representation formula (2.40).

Having in mind the inequality (3.2), we set for k ≥ 0,

λ k := 2k + d -3 + δ , (3.3)
where δ ≥ 3 -d is an arbitrary fixed real parameter (so that λ k ≥ 0) and

h(t) = t -δ Ã(-log t) -A(-log t) (3.4)
We define formally a new amplitude à by adding to A a power series

Ã(α) = A(α) + k≥0 c k e -(λ k +δ)α , (3.5)
or equivalently

h(t) = k≥0 c k t λ k . (3.6)
We assume that the series defining h(t) has a radius of convergence R > 1, so that h ∈ C 0 ([0, 1]). Furthermore we assume that h is such that the estimate (3.2) holds, that is

| 1 0 h(t) t λ k dt| ≤ ǫ , ∀k ≥ 0 . (3.7)
Our goal for this section is to obtain a good approximation of Hölder or Lipschitz type for ||h|| 2 2 , under the above assumptions. We shall do so by using Theorem 2 in the paper [START_REF] Still | On the approximation of Müntz series by Müntz polynomials[END_REF] by Still and the procedure used in Section 4.3 of [START_REF] Daudé | Stability in the inverse Steklov problem on warped product Riemannian manifolds[END_REF]. In order to do so we first recall some of the notation used in [START_REF] Daudé | Stability in the inverse Steklov problem on warped product Riemannian manifolds[END_REF].

Theorem 3.1. Given ǫ > 0 and R > 1 as above and letting M 0 = max{2, 4(d -3 + δ) + 1}, we have, for some universal constant B > 0, the estimate

||h|| 2 2 ≤ B 2 ǫ + R 1-d ǫ log R log( 9M 0 2
) .

(3.8)

We note that the estimate (3.8) is generally a Hölder type estimate for ||h|| 2 2 , but that if R > 9M0 2 , this estimate is Lipschitz.

Proof. Given a sequence Λ

∞ := (λ n ) n≥0 of integers such that 0 ≤ λ 0 < λ 1 < • • • and λ k → ∞ as k → ∞,
we define for fixed n ≥ 1 the finite sequence

Λ n := 0 ≤ λ 0 < λ 1 < ... < λ n , (3.9)
giving rise to the vector space M(Λ) of "Müntz polynomials of degree λ n ":

M(Λ n ) = {P : P (t) = n k=0 a k t λ k } . (3.10)
Recall that according to the Müntz-Szász's Theorem, if Λ ∞ is a sequence of positive real numbers as above, then span {t λ0 , t λ1 , ...} is dense in

L 2 ([0, 1]) if and only if ∞ k=1 1 λ k = ∞. (3.11)
We remark that if λ 0 = 0, the denseness of the Müntz polynomials in C 0 ([0, 1]) in the sup norm is also characterized by (3.11).

Given now a function f in C 0 ([0, 1]) or in L 2 ([0, 1]), the error of approximation of f with respect to

M(Λ n ) is defined by E p (f, Λ n ) := inf P ∈M(Λn) ||f -P || p , (3.12)
where p = 2 or p = ∞ depending on whether

f ∈ C 0 ([0, 1]) or f ∈ L 2 ([0, 1]
). For our application, we have λ k := 2k + d -3 + δ, giving λ k+1 -λ k = 2 > 0, so by Theorem 2 of [START_REF] Still | On the approximation of Müntz series by Müntz polynomials[END_REF], we know that

E ∞ (h, Λ n ) ≤ CR -λn+1 , (3.13)
for some positive constant C.

We have, denoting by π n the orthogonal projection onto the subspace M(Λ n ),

||h|| 2 2 = ||π n (h)|| 2 2 + ||h -π n (h)|| 2 2 . (3.14)
Our next step is to combine the estimate (4.57) from [START_REF] Daudé | Stability in the inverse Steklov problem on warped product Riemannian manifolds[END_REF] Section 4.3 of and the estimate (3.13) to obtain an estimate for the norm of π n (h). In order to do so, we use the Gram-Schmidt process to obtain polynomials (L m (t)) with L 0 (t) = 1, and for m ≥ 1,

L m (t) = m j=0 C mj t λj , (3.15)
where we have set

C mj = 2λ m + 1 m-1 r=0 (λ j + λ r + 1) m r=0,r =j (λ j -λ r ) . (3.16)
The family (L m (t)) defines an orthonormal Hilbert basis of L 2 ([0, 1]). We may now recall the estimate (4.57) from [START_REF] Daudé | Stability in the inverse Steklov problem on warped product Riemannian manifolds[END_REF],

||π n (h)|| 2 2 ≤ ǫ 2 n k=0 k p=0 |C kp | 2 , (3.17
||h -π n (h)|| 2 2 = E 2 (h, Λ n ) ≤ E ∞ (h, Λ n ) .
Now, according to the estimate (4.72) of [START_REF] Daudé | Stability in the inverse Steklov problem on warped product Riemannian manifolds[END_REF], we have

||π n h|| 2 2 ≤ B 2 ǫ 2 g(n) 2 , (3.19)
where B is a positive constant and g : [0, +∞[ is a monotone increasing function defined for t ∈ [0, +∞) by

g(t) = 3 2 1 9M0 2 2 -1 √ 2t + 1 9M 0 2 t+1 , M 0 = max{2, 4(d -3 + δ) + 1} . (3.20)
Now, repeating the steps that lead from the inequalities (4.72) to (4.73) in [START_REF] Daudé | Stability in the inverse Steklov problem on warped product Riemannian manifolds[END_REF], we choose n as a function of ǫ so as to control the norm of the projection ||π n h|| 2 2 of h and thus set n(ǫ

) := [ (g -1 ( 1 √ ǫ ))]
where square brackets denote the integral part function. Since g is a monotone increasing function, we have

g(n(ǫ)) ≤ 1 √ ǫ , (3.21
n(ǫ) = log( 1 √ ǫ ) log( 9M0 2 ) . (3.23) Plugging this into (3.18) gives ||h|| 2 2 ≤ B 2 ǫ + CR -λ n(ǫ)+1 . (3.24)
Now, using the expression λ k = 2k + d -3 + δ and (3.23), we have

R -λ n(ǫ)+1 = R 1-d-δ R -2n(ǫ) ∼ R 1-d-δ R - 2 log 1 √ ǫ log( 9M 0 2 ) ∼ R 1-d-δ e log ǫ log( 9M 0 2 ) log R ∼ R 1-d-δ ǫ log R log( 9M 0 2
) . (3.25) Substituting (3.25) into (3.24), we obtain

||h|| 2 2 ≤ B 2 ǫ + R 1-d-δ ǫ log R log( 9M 0 2
) .

(3.26)

In terms of the amplitude function A in the variable α ∈ (0, ∞), using the relation

||h|| 2 2 = 1 0 t -2δ (A(-log t) -Ã(-log t)) 2 dt , we obtain ∞ 0 e (2δ-1)α (A(α) -Ã(α)) 2 dα ≤ B 2 ǫ + R 1-d-δ ǫ log R log( 9M 0 2
) .

(3.27)

4 From the perturbed amplitude à to a potential Q ∈ L 2 (0, ∞)

Statement of the main result

Our objective for Section 4 is to establish a result on the existence of square-integrable potentials Q associated to perturbed amplitudes à as defined in (3.5). As we shall see, this will require a few additional hypotheses on the perturbation of the amplitude function A given by (3.5) and thus on the perturbation of the starting potential Q.

We set for k ≥ 0,

µ k := λ k + δ = 2k + d -3 + 2δ (4.1) so that Ã(α) = A(α) + k≥0 c k e -µ k α , α > 0 . (4.2)
Thus, in dimension d greater than 3, µ k may be a negative real; so we split the series in (4.2) as

Ã(α) = A(α) + N -1 k=0 c k e -µ k α + k≥N c k e -µ k α , α > 0 , (4.3)
such that for k = 0, ..., N -1, µ k < 0 and for k ≥ N + 1, µ k ≥ 0, with the convention that the first sum in (4.3) does not appear if all the µ k 's are positive (i.e if N = 0). Now, it is convenient to rewrite (4.3) as:

Ã(α) -A(α) = N -1 k=0 2c k sinh(|µ k |α) + k≥0 c k e -|µ k | α . (4.4)
In this form, we note that the perturbation Ã(α)-A(α) has exactly the same expression as the amplitude given in [START_REF] Gesztesy | A new approach of inverse spectral theory, II. General potentials and the connection to the spectral measure[END_REF], Eq. (11.9), modulo the fact we consider a convergent series instead a finite sum, and that the coefficients |µ k | for k = 0, ..., N -1 appear simultaneously in the first finite sum and also in the convergent series.

As we shall see in the next section, the first finite sum in the (RHS) of (4.4) corresponds to the introduction of negative eigenvalues - 2 for k ≥ 0. We state this result in the form of a theorem, that will be proved in Section 4.3 below by applying Theorem 1.2 from the paper [START_REF] Killip | Sum rules and spectral measures of Schrödinger operators with L 2 potentials Ann[END_REF] by Killip and Simon: Theorem 4.1. Let Q ∈ L 2 (0, ∞) be a square-integrable potential with amplitude function A, let {c k , k ≥ 0} be a sequence of real numbers such that

• i) For all k ≥ 0, c k ≤ 0.
• ii) the power series k≥0 c k t λ k has a radius of convergence R > 1.

Then the function à defined by

Ã(α) = A(α) + k≥0 c k e -µ k α , α > 0 , (4.5)
is the amplitude function of a potential Q ∈ L 2 (0, ∞).

The spectral measure

In order to proceed with the proof of Theorem 4.1, we need to first compute the difference of the spectral measures of Q and Q in terms of the data contained in the perturbed amplitude Ã. In the first instance we will work heuristically so as to set the stage for the mathematical objects at play.

Let us first recall from [START_REF] Gesztesy | On local Borg-Marchenko uniqueness results[END_REF] that the Weyl-Titchmarsh function M is a function of Herglotz type, meaning that for all z ∈ C such that Imz > 0 we have Im M (z) > 0. This implies that we have the representation formula

M (z) = c + R 1 λ -z - 1 1 + λ 2 dρ(λ) ,
where dρ(λ) is the (positive) spectral measure associated to (2.30). It can be constructed by taking the following weak limit (in the distributional sense)

dρ(E) = w-lim ǫ↓0 1 π Im (M (E + iǫ)) dE . (4.6)
We denote by M (-κ 2 ) the putative Weyl-Titchmarsh function associated with the amplitude Ã(α), so thanks to (2.41), we have (formally and for suitable κ),

M (-κ 2 ) -M (-κ 2 ) = - ∞ 0 Ã(α) -A(α) e -2κα dα. (4.7)
Now, using (4.4), we easily write the difference M (-κ 2 ) -M (-κ 2 ) as

M (-κ 2 ) -M (-κ 2 ) = -2 N -1 k=0 c k |µ k | 4κ 2 -µ 2 k - k≥0 c k 2κ + |µ k | . (4.8)
We note that the above series is indeed convergent since by hypothesis the power series k≥0 c k t λ k has radius of convergence R > 1 so that in particular k≥0 |c k | < ∞. Now, using (4.6), we are able to define the difference of the spectral measures dρ(E) -dρ(E) (for more details, we refer the reader to ([9], Eqs. (11.7)-(11.9), p.637).

For E ≥ 0,

dρ(E) = dρ(E) - 2 π k≥0 c k √ E 4E + µ 2 k dE, (4.9)
and for E < 0,

dρ(E) = dρ(E) - 1 2 N -1 k=0 c k |µ k | δ(• + µ 2 k 4 ) dE, (4.10)
where δ(• -a) stands for the usual delta distribution, centered at the point a.

As was explained in ( [START_REF] Gesztesy | A new approach of inverse spectral theory, II. General potentials and the connection to the spectral measure[END_REF], Section 11), this corresponds to the introduction of a finite number of negative eigenvalues - As earlier, we now apply Theorem 1.2 from the paper [START_REF] Killip | Sum rules and spectral measures of Schrödinger operators with L 2 potentials Ann[END_REF] by Killip and Simon in order to prove Theorem 4.1 establishing existence of a potential Q ∈ L 2 (0, ∞) associated to the perturbed amplitude Ã. Besides the positivity of the perturbed measure dρ, there are four conditions stated in the theorem of Killip and Simon that we need to verify on dρ in order for their theorem to apply. In what follows, we state these conditions and show they are satisfied under the hypotheses of Theorem 4.1.

• Positivity of the measure dρ(E):

This follows immediately on account of (4.9), (4.10), and the hypothesis c k ≤ 0 , ∀ k ≥ 0.

• Weyl condition: The Weyl condition on dρ states that the support of dρ should decompose as

supp dρ = [0, ∞) ∪ { Ẽj } N j=1 , with Ẽ1 < Ẽ2 < • • • < 0 , with Ẽj → 0 if N = ∞ .
Again this is immediate from the fact that the measure dρ associated to Q satisfies the Weyl condition, from the identities (4.9), (4.10) and from the fact that our perturbation is only adding real resonances and a finite number of negative eigenvalues eigenvaluesµ k 2 4 for k = 0, ..., N -1. • Normalization: We need to verify that dρ satisfies a certain estimate whose formulation requires the introduction of the Hardy-Littlewood maximal function of a measure. The argument is more elaborate and we present it therefore is greater detail. Following [START_REF] Killip | Sum rules and spectral measures of Schrödinger operators with L 2 potentials Ann[END_REF], we introduce a measure dν on (1, ∞) parametrized by k, with E = k 2 ,

dν dk = Im M (k 2 + i0) -k ,
which gives on account of (4.7) with κ = -ik,

dν dk = dν dk -2 n≥0 c n k 4k 2 + µ 2 n .
We then define the Hardy-Littlewood maximal function M s of the measure ν by

(M s ν)(x) = sup 0<L≤1 1 2L |ν|([x -L, x + L]) .
and compute |ν|([x -L, x + L]):

|ν|([x -L, x + L]) = |ν|([x -L.x + L]) -2 n≥0 c n x+L x-L k 4k 2 + µ 2 n dk . Since x+L x-L k 4k 2 + µ 2 n dk = 1 8 log (x + L) 2 + µ 2 n (x -L) 2 + µ 2 n , we obtain |ν|([x -L, x + L]) = |ν|([x -L.x + L]) -2 n≥0 c n 8 log (x + L) 2 + µ 2 n (x -L) 2 + µ 2 n .
Taking x = k >> 1, the normalization condition we need to verify is

+∞ 1 log 1 + (M s ν)(k) k 2 k 2 dk < ∞ . (4.11)
We have

log (k + L) 2 + µ 2 n (k -L) 2 + µ 2 n = log 1 + 4kL (k -L) 2 + µ 2 n = log(1 + O( L k )) ,
uniformly in n and L ∈ (0, 1]. But by condition ii) in Theorem 4.1, we know that the series

n≥0 |c n | is convergent, so it follows that |ν|([x -L, x + L]) = |ν|([x -L.x + L]) + O( L k ) , k → ∞ ,
uniformly in L ∈ (0, 1]. Therefore we have

(M s ν)(k) = (M s ν)(k) + O( 1 k
) . (4.12)

Now, using the inequality

log 1 + (x + y) 2 ≤ C x 2 + log(1 + y 2 ) ,
for some constant C > 0, we obtain using (4.12) the estimate

log 1 + (M s ν)(k) k 2 = log 1 + (M s ν)(k) + O( 1 k ) k 2 ≤ C O( 1 k 4 ) + log 1 + (M s ν)(k) k 2 ,
which implies the normalization condition (4.11) after integration over R, using the fact that ν is associated to a potential Q ∈ L 2 (0, ∞).

• Lieb-Thirring condition: The Lieb-Thirring condition j | Ẽj | 3/2 < ∞ is trivially satisfied here as was the case for the Weyl condition since all we are doing is to add a finite number of negative eigenvalues.

• Quasi-Szegö condition: This condition states that if dρ 0 is the free spectral measure, that is the spectral measure associated to the zero potential Q ≡ 0, then

∞ 0 log 1 4 dρ dρ 0 + 1 2 + 1 4
dρ 0 dρ √ E dE < ∞ . (4.13)
Again, the verification of this condition for the perturbed measure dρ is more elaborate and we therefore present it in greater detail. The spectral spectral measure dρ 0 has the expression

dρ 0 (E) = 1 π χ (0,∞) (E) √ E dE ,
so that using (4.9), we have

dρ(E) -dρ 0 (E) = dρ(E) -dρ 0 (E) - 2 π n≥0 c n √ E 4E + µ 2 n dE . (4.14)
We now use (4.6) to express the spectral measure dρ in terms of the Jost function ψ(x, κ) associated to the potential Q, where we let E = -κ 2 . Recall that the Weyl-Titchmarsh function M is given by

M (-κ 2 ) = ψ ′ (0, κ) ψ(0, κ) , (4.15)
so that using (4.15), we obtain

Im M (-κ 2 ) = 1 2i M (-κ 2 ) -M (-κ 2 ) = 1 2i ψ ′ (0, κ) ψ(0, κ) - ψ ′ (0, κ) ψ(0, κ) = 1 2i W (ψ, ψ) |ψ| 2 (0, κ) , (4.16)
where W denotes the Wronskian. But ψ and ψ are solutions of the same linear second-order ODE since the potential Q and the spectral parameter E = -κ 2 are both real. It follows that the Wronskian W (ψ, ψ) is independent of x. Now since ψ is the Jost function, we have, in terms of the parameter k = iκ introduced in the normalization condition the asymptotics

ψ(x, κ) ≃ e -κx = e ikx , for x → ∞ , so that W (ψ, ψ) = 2ik .
Substituting the latter into (4.16), we obtain

Im M (-κ 2 ) = iκ |ψ(0, κ)| 2 , which plugged into in (4.6) gives for E > 0 dρ(E) = 1 π √ E |ψ(0, √ E)| 2 dE . ( 4
(E) = 1 |ψ(0, √ E)| 2 -2 n≥0 c n 4E + µ 2 n . (4.18)
We now analyze the asymptotics of d ρ dρ0 (E) in the limit E → ∞. On the one hand we have

1 4E + µ 2 n = 1 4E 1 + O( µ n 2 E ) ,
and on the other hand we know that since the radius of convergence R of the series n≥0 c n t λ n satisfies R > 1 and since λ n = O(n), the series n≥0 |c n |µ 2 n is convergent. The identity (4.18) implies

dρ dρ 0 (E) = 1 |ψ(0, √ E)| 2 - 1 2 n≥0 c n 1 E + O( 1 E 2 ) .
Using the asymptotics on the modulus of the Jost function given by

|ψ(0, √ E)| = 1 + a E + O 1 E 2 ,
where a is a real constant, we obtain that

dρ dρ 0 (E) = 1 + b E + O( 1 E 2 ) , (4.19)
for some real constant b, which implies in turn that

dρ 0 dρ = 1 - b E + O( 1 E 2 ) . ( 4
dρ 0 dρ = 1 + O( 1 E 2 ) , (4.21)
which implies that the Quasi-Szegö condition (4.13) is satisfied since (4.21) implies that log 1 4

dρ dρ 0 + 1 2 + 1 4
dρ 0 dρ = O( 1 E 2 ) .

A few examples : Bargmann potentials

In this section, we consider perturbations of the potential Q(x) = 0, (with the associated amplitude function A(α) = 0), and we give examples of amplitudes Ã(α) for which we can calculate explicitly the associated potentials. These examples are borrowed from ([9], section 11).

More precisely, we assume here that the potential Q ∈ L 2 (0, ∞), (and thus its associated amplitude A), is fixed, and that Q ∈ L 2 (0, ∞) belongs to the infinite dimensional class, denoted C Q , defined above, that is we assume that the associated amplitude to Q has the form Ã(α) = A(α) + k≥0 c k e -µ k α , α > 0 , (

where µ k = 2k + d -3 + 2δ and δ ≥ 3 -d. Moreover, one assumes that c k ≤ 0 for all k ≥ 0 and the power series k≥0 c k t λ k has a radius of convergence R > 1.

To obtain these local stability estimates, we shall make intensive use the local version of the classical Gel'fand-Levitan equations, (see for instance ( [START_REF] Avdonin | Inverse problem on the semi-axis: local approach[END_REF], Eq. (2.24)) which we recall here. For 0 ≤ x ≤ t ≤ T , we consider the integral equation V (x, t) + T x K(t, s)V (x, s) ds = -K(x, t), (

where the integral kernel K(t, s) is given by K(t, s) = p(2T -t -s) -p(|t -s|), (6. Proof. This is an immediate application of (3.27) and the Cauchy-Schwartz inequality.

Now, let us introduce some notation to simplify the presentation below. In what follows, the parameters x and T are assumed to be fixed and t is a variable lying in the interval [x, T ]. We denote by K the integral operator on L 2 (x, T ) with kernel K(t, s), Kf (t) =

in the sense of the operator norm on L 2 (x, T ). As a consequence for ǫ > 0 sufficiently small, the operator I + (I + K) -1 ( K -K) is invertible, and using again the resolvent identity, one obtains easily

(I + K) -1 = I + (I + K) -1 ( K -K) -1 (I + K) -1 . (6.15)
It follows that, for ǫ << 1, the operator norm of (I + K) -1 is uniformly bounded:

||(I + K) -1 || ≤ 2 ||(I + K) -1 ||. (6.16) Thus, thanks to (6.12) , (6.14) and (6.16), one has:

|| Ṽ -V || L 2 (x,T ) ≤ C T f (ǫ) (6.17)

In the same way, differentiating the integral equation (6.2) with respect to x, one obtains:

∂V ∂x (x, t) + T x K(t, s) ∂V ∂x (x, s) = -p ′ (t -x) + p ′ (2T -x -t) + K(t, x)V (x, x), (6.18) and by the same argument, we get immediately

|| ∂ Ṽ ∂x - ∂V ∂x || L 2 (x,T ) ≤ C T f (ǫ). (6.19)
Finally, using (6.6), (6.17 Then taking x = 0 and using (6.5), we see that

|| Q -Q|| L 2 (x,T ) ≤ C T f (ǫ), (6.21)
and the proof of Theorem 1.1 is complete.

µ 2 k 4

 24 for k = 1, ..., N , (corresponding to the case where µ k is negative), and of a countable set of real resonances -|µ k |

  .34) But by the Cauchy-Schwarz inequality, we have sup

2 +

 2 CR -λn+1 , (3.18) using (3.13) and the inequality

µ 2 k 4

 24 for k = 0, ..., N -1, whereas the second sum corresponds to the introduction of real resonances -|µ k |

µ 2 k 4

 24 for k = 0, ..., N -1, and of real resonances -|µ k | 2 for k ≥ 0. Now, we can explain precisely our strategy : in the next section we show, using the Killip-Simon conditions, that under the hypotheses of Theorem 4.1, there exists a potential Q ∈ L 2 (0, ∞) associated to the above spectral measure dρ(E) allowing us to define rigourously the associated Weyl-Titchmarsh function M (z) for z ∈ C\[-β, +∞[ for β >> 1. The amplitude function associated to M (z) is automatically given by Ã(α) thanks to the uniqueness of the inverse Laplace transform and analytic continuation.

4. 3

 3 The Killip-Simon conditions and the proof of Theorem 4.1

4 )Lemma 6 . 1 .

 461 These integral equations are uniquely solvable for all x ∈ (0, T ) and we can recover the underlying potential using the relation:Q(T -x) x, x)) = p(2T -x)V (x, x) + 2p ′ (2T -2x) -T x (p(2T -x -s) -p(s -x)) ∂V ∂x (x, s) ds -T x (p ′ (2T -x -s) -p ′ (s -x)) V (x, s) ds.(6.6)Let us begin with an elementary result: Under the hypotheses of Theorem 1.1, there exists a constant C T depending only onT such that ||p -p|| (C 0 (0,2T ),||•||∞) ≤ C T f (ǫ),(6.7)wheref (ǫ) = B 2 ǫ + R 1-d ǫ

  ) := p(t -x) -p(2T -x -t). (6.10) Thus, the solution V (x, .) of the integral equation (6.2) can be written asV := V (x, .) = (I + K) -1 d. (6.11)Using (6.11) and the usual resolvent identity, one obtainsṼ -V = (I + K) -1 d -d + (K -K)(I + K) -1 d (6.12)By Lemma 6.1, one has the uniform estimate for t, s ∈ [0, T ],| d(t) -d(t)| ≤ C T f (ǫ) , | K(t, s) -K(t, s)| ≤ C T f (ǫ),(6.13)thus using Schur's lemma, one gets || K -K|| ≤ C T f (ǫ),

  ) and (6.19), mimicking the above arguments, one has for all 0 ≤ x ≤ T ,|| d dx Ṽ (x, x) -d dx (V (x, x)) || L 2 (x,T ) ≤ C T f (ǫ). (6.20)

  1 2 ||Q|| 1 . In[START_REF] Avdonin | The boundary control approach to the Titchmarsh-Weyl m-function. I. The response operator and the A-amplitude[END_REF] (see Section 5, Algorithm 1, point 2), it is proved that if β 2 < ∞, where β 2 is defined in (2.34), then the integral in (2.40) is absolutely convergent for Re κ > 2 max{ √ 2β 2 , eβ 2 }. But we saw in (2.35) that for square-integrable potentials Q, one has the estimate β 2 ≤ ||Q|| 2 < ∞. It therefore follows from (2.40) that the Weyl-Titchmarsh function M admits the representation

* Research supported by the French National Research Projects AARG, No. ANR-12-BS01-012-01, and Iproblems, No. ANR-13-JS01-0006 † Research supported by NSERC grant RGPIN 105490-2018

First example

We define for α ≥ 0, Ã(α) = 2(γ 2 -β 2 ) e -2γα , (

where β > 0 and γ ∈ [0, β[. Of course, it corresponds to a Müntz series of the type (4.2) (with a single term) with c 0 = 2(γ 2 -β 2 ) < 0 and µ 0 = 2γ ≥ 0.

Thus, we take δ = γ + 3-d 2 ≥ 3 -d since d ≥ 3 by hypothesis. It is known that

(1 + β-γ β+γ e -2βx ) 2

(5.

2)

The associated Jost function is given in the variable κ = -ik by

(see [START_REF] Gesztesy | A new approach of inverse spectral theory, II. General potentials and the connection to the spectral measure[END_REF], case 2, p. 636) and is holomorphic in Re κ > -β. The unique root of the Jost function is given by κ = -γ which is a real resonance.

Second example

For α > 0, we define the amplitude

where c 1 > 0 is a normalization constant and κ 1 > 0. It corresponds to a Müntz series with two terms and with two µ k of different sign. The associated potential is given by

ahe Jost function has the form in the κ variable

(see [START_REF] Gesztesy | A new approach of inverse spectral theory, II. General potentials and the connection to the spectral measure[END_REF], case 1, p. 635). We note that the Jost function is vanishing at κ = κ 1 wich correspond to the single negative eigenvalue -κ 2 1 .

Gel'fand-Levitan equations and local stability estimates

In this section, we deduce from the estimates for the difference of the amplitudes A -Ã obtained in Section 4 a set of new Hölder local stability estimates for the difference of the associated potentials Q -Q. By local stability, we mean that we are able to control the norm ||Q -Q|| L 2 (0,T ) with respect to ǫ, if the Steklov spectra of the underlying Schrödinger operators are close up to ǫ as in (1.17), (T being any fixed positive parameter).