Dimitrios Mitsotakis

Denys Dutykh
email: denys.dutykh@univ-savoie.fr

John D Carter
email: carterj1@seattleu.edu

On the nonlinear dynamics of the traveling-wave solutions of the Serre system

Keywords: Solitary waves, cnoidal waves, stability, finite element method MSC: [2010]76B15 (primary), 76B25, 76M10 (secondary)

We numerically study nonlinear phenomena related to the dynamics of traveling wave solutions of the Serre equations including the stability, the persistence, the interactions and the breaking of solitary waves. The numerical method utilizes a high-order finite-element method with smooth, periodic splines in space and explicit Runge-Kutta methods in time. Other forms of solutions such as cnoidal waves and dispersive shock waves are also considered. The differences between solutions of the Serre equations and the Euler equations are also studied.

Introduction

The Serre equations (also known as the Green-Naghdi or Su-Gardner equations) [19,[START_REF] Serre | Contribution à l'étude des écoulements permanents et variables dans les canaux[END_REF][START_REF] Su | Korteweg-de Vries equation and generalizations. III. Derivation of the Korteweg-de Vries equation and Burgers equation[END_REF] approximate the Euler equations of water wave theory and model the onedimensional, two-way propagation of long waves. If a denotes a typical amplitude of a wave, d the mean depth of the fluid, and λ a typical wavelength, then the Serre equations are characterized by the parameters ε ≐ a d = O(1) and σ ≐ d λ ≪ 1, contrary to the Boussinesq equations which model the propagation of small-amplitude, long waves, i.e. ε ≪ 1 and σ ≪ 1, when the Stokes number is S ≐ ε σ 2 = O [START_REF] Ablowitz | Interactions and asymptotics of dispersive shock waves -Korteweg-de Vries equation[END_REF]. The Boussinesq equations are often called weakly nonlinear, weakly dispersive equations while the Serre equations are often called fully-nonlinear shallow-water equations. In dimensionless and scaled variables, the Serre equations take the form:

η t + u x + ε(ηu) x = 0 , u t + η x + εuu x -σ 2 3h [h 3 (u xt + εuu xx -ε(u x ) 2 )] x = 0 , (1.1)
for x ∈ R, t > 0, along with the initial conditions η(x, 0) = η 0 (x) , u(x, 0) = u 0 (x) .

Here η = η(x, t) is the free surface displacement, while

h ≐ 1 + εη , (1.3)
is the total fluid depth, u = u(x, t) is the depth-averaged horizontal velocity, and η 0 , u 0 are given real functions, such that 1 + εη 0 = h 0 > 0 for all x ∈ R. In these variables, the location of the horizontal bottom is given by y = -1. For a review of the derivation and the basic properties of this system we refer to [START_REF] Barthélémy | Nonlinear shallow water theories for coastal waves[END_REF].

The Euler equations along with the model system (1.1) admit traveling wave solutions, i.e. waves that propagate without change in shape or speed [3,[START_REF] Bona | Decay and analyticity of solitary waves[END_REF][START_REF] Lavrentiev | On the theory of long waves[END_REF]. Solitary waves form a special class of traveling wave solutions of these systems. The other important class of traveling wave solutions is the class of cnoidal wave solutions which can be thought of as the periodic generalization of solitary waves. Many Boussinesq-type equations are known to possess solitary wave and periodic solutions, but do not admit nontrivial, closedform solutions. In contrast, the Serre equations admit closed-form solitary and cnoidal (periodic) wave solutions. The solitary wave solutions of the Serre system traveling with constant speed c s are given by h s (ξ) = (a 0 + a 1 sech 2 (K s ξ)) σ, u s (ξ) = c s 1 -a 0 σh s (ξ) ǫ , (1.4) where ξ = xc s t, K s = 3a 1 4σa 2 0 c 2 s , c s = (a 0 + a 1 ) σ, a 0 > 0, and a 1 > 0. By taking a 0 = σ and a 1 = εσA s the formulas for the classical solitary waves that are homoclinic to the origin are obtained.

The cnoidal waves of the Serre system traveling with constant speed c c are given by

h c (ξ) = (a 0 + a 1 dn 2 (K c ξ, k)) σ , u c (ξ) = c c 1 - h 0 h(ξ) ǫ , (1.5)
where h 0 = a 0 + a 1 E(m) K(m), K c = √ 3a 1 2 a 0 (a 0 + a 1 )(a 0 + (1k 2 )a 1 ), c c = a 0 (a 0 + a 1 )(a 0 + (1k 2 )a 1 ) σh 2 0 , k ∈ [0, 1], m = k 2 , a 0 > 0, and a 1 > 0. Here K and E are the complete elliptic integrals of the first and second kind respectively. Note that (1.4) are the k → 1 limit of (1.5).

Another fundamental property of the Serre system is the conservation of the energy which plays also the role of the Hamiltonian,

H(t) = 1 2 ∫ ∞ -∞ (εhu 2 + εσ 2 3 h 3 u 2
x + εη 2 )dx, in the sense that H(t) = H(0) for all t > 0 up to the maximal time T of the existence of the solution.

In this paper we study the problem of the nonlinear stability (orbital and asymptotic) of the traveling waves of the Serre system by using numerical techniques. We provide numerical evidence of stability with respect to certain classes of perturbations. Phenomena such as perturbations of the traveling waves, perturbations of the Serre system and interactions of traveling waves are studied analyzing the stability properties of the waves at hand. We also study the interactions of dispersive shock waves (DSWs) in the Serre system. The physical relevance of the Serre equations is addressed whenever possible.

The paper is organized as follows. The numerical method is presented briefly in Section 2. The compatibility of the solitary waves of the Serre and the Euler systems is examined in Section 3. The head-on collision of solitary waves is studied in Section 3.2. A number of issues related to the stability of the traveling waves are discussed in Sections 4. The interaction of DSWs is presented in Section 5.

The numerical method

The numerical method of preference is a high-order Galerkin / Finite element method (FEM) for the spatial discretization combined with the classical fourth-order explicit Runge-Kutta method in time. In some cases adaptive time-stepping methods, such as the Runge-Kutta-Fehlberg, the Cash-Karp and the Dormand-Prince methods [20], were employed to verify that there are no spurious solutions or blow-up phenomena. This numerical scheme has been shown to be highly accurate and stable since there is no need for a restrictive condition on the step-size but only mild conditions of the form ∆t ≤ C∆x cf. [START_REF] Mitsotakis | On the Galerkin/Finite-Element Method for the Serre Equations[END_REF]. The conservation of the Hamiltonian was monitored and was usually conserved to within 8 to 10 significant digits. In order to ensure the accuracy of the numerical results obtained with the FEM we compared most with the analogous results obtained with the pseudo-spectral method described and analyzed in [START_REF] Dutykh | Finite volume and pseudo-spectral schemes for the fully nonlinear 1D Serre equations[END_REF]. The experiments presented in this paper also serve as numerical benchmarks for the efficacy of the numerical scheme.

We consider (1.1) with periodic boundary conditions and, for simplicity, assume ε = σ = 1. We rewrite (1.1) in terms of (h, u) rather than (η, u). This is done by using (1.3) and yields the initial-boundary value problem

h t + (hu) x = 0 , u t + h x + uu x -1 3h [h 3 (u xt + uu xx -(u x ) 2 ] x = 0 , ∂ i x h(a, t) = ∂ i x h(b, t), i = 0, 1, 2, . . . , ∂ i x u(a, t) = ∂ i x u(b, t), i = 0, 1, 2, . . . , h(x, 0) = h 0 (x) , u(x, 0) = u 0 (x) , (2.1)
where x ∈ [a, b] ⊂ R and t ∈ [0, T ]. Considering a spatial grid x i = a + i ∆x, where i = 0, 1, ⋯, N, ∆x is the spatial mesh length, and N ∈ N, such that ∆x = (ba) N. We define the space of the periodic cubic splines

S = φ ∈ C 2 per [a, b] φ [x i ,x i+1 ] ∈ P 3 , 0 ≤ i ≤ N -1 ,
where

C 2 per = f ∈ C 2 [a, b] f (k) (a) = f (k) (b)
, 0 ≤ k ≤ r and P k is the space of polynomials of degree k. The semi-discrete scheme is reduced to finding h, ũ ∈ S such that

( ht , φ) + ( hũ) x , φ = 0 , B(ũ t , φ; h) + h( hx + ũũ x ), φ + 1 3 h3 (ũũ xx -(ũ x ) 2 ), φ x ) = 0 , (2.2)
where B is defined as the bilinear form that for fixed h is given by B(ψ, φ; h) ≐ ( hψ, φ) + 1 3 ( h3 ψ x , φ x ) for φ, ψ ∈ S .

(2.3)

The system of equations (2.2) is accompanied by the initial conditions

h(x, 0) = P{h 0 (x)} , ũ(x, 0) = P{u 0 (x)} , (2.4)
where P is the L 2 -projection onto S satisfying (Pv, φ) = (v, φ) for all φ ∈ S. Upon choosing appropriate basis functions for S, (2.2) is a system of ODEs. For the integration in time of this system, we employ the classical, four-stage, fourth-order explicit Runge-Kutta method.

Solitary waves

In this section we study how close are the solitary waves of the Euler equations to those of the Serre equations. In other words we verify the consistency of the Serre equations and the ability to approximate well the basic solitary wave dynamics of the Euler equations.

Consistency of solitary waves

The Serre system and the Euler equations both possess solitary waves that decay exponentially to zero at infinity. Although the justification of the Serre equations ensures that its solutions will remain close to Euler solutions, it is not known how close remain an Euler solitary wave to a Serre solitary wave when it is used as initial condition to the Serre system. While the Serre system admits solitary wave solutions of the form given in (1.4), there are no known closed-form solitary wave solutions of the Euler system but only Fenton's asymptotic solution [START_REF] Fenton | A ninth-order solution for the solitary wave[END_REF]. Although this solution is an accurate approximation, modern numerical techniques enable us to compute solitary waves of the Euler equation even more accurately. For this reason we compute Euler solitary waves numerically. The numerical method is a Petviashvili iteration applied to the Babenko equation, [9, [START_REF] Dutykh | Efficient computation of steady solitary gravity waves[END_REF][START_REF] Petviashvili | Equation of an extraordinary soliton[END_REF]. In order to integrate the full Euler equations in time, we employ the method of holomorphic variables. This formulation was first coined by L. [START_REF] Ovsyannikov | To the shallow water theory foundation[END_REF] [START_REF] Ovsyannikov | To the shallow water theory foundation[END_REF] and developed later by A. [START_REF] Dyachenko | Nonlinear dynamics of the free surface of an ideal fluid[END_REF] [START_REF] Dyachenko | Nonlinear dynamics of the free surface of an ideal fluid[END_REF] in deep waters. The extension to the finite depth case was given in [START_REF] Li | A Numerical Study of the Exact Evolution Equations for Surface Waves in Water of Finite Depth[END_REF]. The resulting formulation is discretized in the conformal domain using a Fourier-type pseudo-spectral method. For the time integration we employ an embedded Runge-Kutta scheme of 5(4)th order along with the integrating factor technique to treat the dispersive linear part.

In order to demonstrate the ability of the Serre equations to approximate the Euler equations, we first compare the characteristics of two solitary waves with speeds c s = 1.1 and c s = 1.2. The solitary waves are not identical but their differences are small and more pronounced at the higher speed. For example, a 10% increase in the speed leads to an increase in solitary wave amplitude of almost 50% while the normalized difference between the Euler and Serre solitary waves increased by more than a factor of two. The amplitude of several Euler and Serre solitary waves are presented in Table 1. Next, we examine how the solitary waves of the Euler system propagate when they are used as initial conditions to the Serre system. Specifically, we use the numerically generated solitary wave solutions of the Euler equations and the exact formula u = c s η (1+η) to define the initial conditions η 0 and u 0 for the Serre equations. Then, we numerically integrate the Serre system. Figures 1 and2 contain plots of the solutions at t = 150 obtained using the Euler solitary waves with c s = 1.1 and c s = 1.2. These figures demonstrate that the difference between the Euler solitary wave and numerical Serre solution is greater when c s = 1.2 than when c s = 1.1. We note that the value c s = 1.2 is a relatively large value since the largest value we can use to generate an Euler solitary wave is c s = 1.29421.

c
To study further the differences between the Euler and Serre solitary waves, we consider three quantities pertinent to the propagation of the solitary waves: the amplitude, shape and phase. First, we define the normalized peak amplitude error as

AE[F ] ≐ F (x * (t), t) -F (0, 0) F (0, 0) , (3.1)
where x * (t) is the curve along which the computed solution F (x, t) achieves its maximum.

Monitoring AE as a function of time, we observe that although the Euler solitary waves do not propagate as traveling waves to the Serre system, their amplitude asymptotically tends towards a constant indicating that they evolve into a solitary wave solution of the Serre equation, see Figures 3 and4. We define the normalized shape error as

SE[F ] ≐ log 10 (min τ ζ(τ )) , ζ(τ ) ≐ F (x, t n ) -F exact (x, τ ) F exact (x, 0) . (3.2)
The minimum in (3.2) is attained at some critical τ = τ * (t n ). This, in turn, is used to define the phase error as

P E[F ] ≐ log 10 ( τ * -t n ). (3.3)
In order to find τ * , we use Newton's method to solve the equation ζ ′′ (τ ) = 0. The initial guess for Newton's method is chosen as τ 0 = t n -∆t. Figures 3 and4 contain plots of the shape and phase errors. We observe that the shape error is of O(10 -3 ) when c s = 1. for other model equations such as the classical Boussinesq system and the results are comparable [START_REF] Bona | Boussinesq equations and other systems for smallamplitude long waves in nonlinear dispersive media. I: Derivation and linear theory[END_REF][START_REF] Duran | On the Galilean Invariance of Some Nonlinear Dispersive Wave Equations[END_REF].

Head-on collision of solitary waves

The collision of two solitary waves of the Serre system has previously been studied theoretically and numerically in [START_REF] Choi | Exact Evolution Equations for Surface Waves[END_REF][START_REF] Craig | Solitary water wave interactions[END_REF][START_REF] Dutykh | Finite volume and pseudo-spectral schemes for the fully nonlinear 1D Serre equations[END_REF][START_REF] Li | A Numerical Study of the Exact Evolution Equations for Surface Waves in Water of Finite Depth[END_REF][START_REF] Marchant | The extended Korteweg-de Vries equation and the resonant flow of a fluid over topography[END_REF][START_REF] Mirie | Collision between two solitary waves. Part 2. A numerical study[END_REF][START_REF] Mitsotakis | On the Galerkin/Finite-Element Method for the Serre Equations[END_REF][START_REF] Benjamin | Solitary-Wave Solutions of Nonlinear Problems[END_REF]. While the phenomena related to these interactions have been understood quite well we summarize here the dynamics of the head-on collision of solitary waves and we focus on the related dynamics compared to experimental data and to numerical simulations of the full Euler equations. The interaction of solitary waves for the Serre equations is general more inelastic than in weakly nonlinear models such as the classical Boussinesq system [START_REF] Bona | Boussinesq equations and other systems for smallamplitude long waves in nonlinear dispersive media. I: Derivation and linear theory[END_REF]. Highly nonlinear interactions result in the generation of large amplitude dispersive tails.

In order to study the physical relevance of the head-on collision of two Serre solitary waves, we compare the Serre numerical solution with the experimental data of [START_REF] Craig | Solitary water wave interactions[END_REF]. In this experiment, the Serre system is written in dimensional and unscaled form with an initial condition that includes two counter-propagating solitary waves in the interval [-5, 5]. The speeds of these solitary waves are c s,1 = 0.7721 m s and c s,2 = 0.7796 m s. Their and the experimental data is impressive. The agreement in the generated dispersive tails in Figure 6 is even more impressive. Such agreement cannot be found in the case of head-on collisions of solitary waves of Boussinesq type models, [START_REF] Dutykh | Finite volume schemes for dispersive wave propagation and runup[END_REF], indicating that the high-order nonlinear terms are important in studying even these small-amplitude solutions. Finally, we mention that the maximum amplitude of the solution observed in Figure 5(c) during the collision is smaller than the real amplitude, possibly, because of a splash phenomenon that cannot be described by any model (see also [START_REF] Dutykh | Finite volume schemes for dispersive wave propagation and runup[END_REF]).

We now compare a head-on collision of two unequal solitary waves via numerical solutions of the Serre and Euler equations. For both models, we consider a right-traveling solitary wave with c s = 1.1 and a left-traveling solitary wave with c s = 1.2. These solitary waves are initially translated so that the maximum peak amplitudes are achieved at x = -100 and x = 100 respectively. Results from the numerical simulations are included in Figure 7. Both models show similar behavior, however the maximum amplitude observed during the collision using the Euler equations is larger than in the Serre system. Also the interaction in the Euler equations lasts longer and therefore a larger phase shift is observed. The leading waves of the dispersive tails are almost identical in the two models, but the amplitude of the tails in the Euler system decay to zero more slowly than the amplitude of the tails in the Serre system. These numerical simulations verify the ability of the Serre system to accurately model head-on collisions of solitary waves. They also show that the Serre system is consistent with the Euler equations during and after the head-on collision with almost identical solutions.

Stability of traveling waves

The previous experiment of the head-on collision of two solitary waves indicates that the solitary waves are robust. In this section, we present the behavior of a solitary wave under small perturbations. We explore the effects of modifying some of the high-order terms of the Serre system. We show that modifying one such term one can produce regularized shock waves, as opposed to classical dispersive shock waves. Finally, we examine the stability of the cnoidal wave solutions.

Stability of solitary waves

We consider perturbations of the amplitude, perturbations of the wavelength, and randomnoise perturbations of the shape. As we show below, all of the solitary waves we tested were stable to all of the perturbations we considered.

We chose a solitary wave with speed c s = 1.4 and amplitude A = 0.96 for all numerical simulations in this section. We perturb the amplitude by multiplying the pulse by a parameter p such that h p (x, 0) = 1 + p ⋅ a 1 sech 2 (K s x) , (4.1) while keeping the velocity component of the solution unperturbed as in (1.4). When p = 1.1 the initial condition sheds a small-amplitude dispersive tail and results in a new solitary wave with amplitude A = 1.02050. Figure 8 presents the initial condition and the resulting solution at t = 130. Similar observations resulted in all cases we tested.

We consider perturbations of the wavelength K s by taking the initial condition for h to be

h p (x, 0) = 1 + a 1 sech 2 (p ⋅ K s x) . (4.
2) The results in this case were very similar to the results we obtained when we perturbed the amplitude of the solitary waves and so we don't show the results here. Table 2 shows the amplitudes of the solitary waves that result from various amplitude and wavelength perturbations. p Amplitude perturbation Wavelength perturbation 0.8 0.83860 1.02710 0.9 0.89936 0.99225 1.1 1.02050 0.93017 1.2 1.08087 0.90260 Table 2. Amplitudes of the uniformly perturbed solitary waves.

Similar results were obtained when non-uniform perturbations were used. In order to consider non-uniform perturbations, we used pseudo-random noise distributed uniformly in [0, 1]. Denoting the noise function by N(x), the perturbed solitary wave is given by where the parameter p determines the magnitude of the noise. Figure 9 shows the perturbed solitary wave with p = 0.2. This type of perturbation is not only non-uniform, but is also non-smooth. Nevertheless this initial condition is the L 2 -projection of the actual solution which ensures the required by the FEM smoothness. Figure 10 shows the evolution of this perturbed solitary wave. The solution consists of a new solitary wave and a small-amplitude dispersive tail. It does not differ qualitatively from the solution shown in Figure 8. The values of the amplitudes of the emerging solitary waves for various values of p are presented in Table 3. These results suggest that the solitary waves of the Serre system are orbitally stable with respect to this class of perturbations.

h p (x, t) = 1 + 1 -p N(x) ⋅ a 1 sech 2 (K s x -c s t) , (4.3)
p Amplitude -0.2 0.90110 -0.1 0.93055 0.1 0.98942 0.2 1.01883 0.5 1.10694

Table 3. Amplitudes of the randomly perturbed solitary waves.

Persistence of the solitary waves

One other aspect related to the stability of the solitary waves is their ability to persist when some of the high-order terms in the PDE are perturbed. In this section, we examine if a Serre solitary wave retains its shape when some of the terms of the Serre system are perturbed. Introducing the parameters α, β, and γ, we rewrite the Serre equations in the following form

η t + u x + (ηu) x = 0 , u t + η x + uu x -1 3h [h 3 (α u xt + βuu xx -γ(u x ) 2 )] x = 0 . (4.4)
The unperturbed Serre equations correspond to α = β = γ = 1. We first study the persistence of the solitary waves when the system is perturbed by perturbing the parameters α, β, γ and considering a solitary wave of the unperturbed system as an initial condition.

In this section we use the solitary wave (1.4) with c s = 1.4 as an initial condition. If α = 0.9, or if β = 0.9, or if γ = 0.9 the solitary wave evolves in a manner similar to the amplitude perturbations in Section 4.1. The new solitary waves are very similar to the unperturbed solitary wave. This further indicates that the solitary waves of the Serre system are stable. More interesting phenomena is observed when the solitary waves are used as initial conditions to systems with small values of the parameters α, β and γ. When all the three parameters are very small, the solutions tend to break into dispersive shock waves or other forms of undular bores. In the first numerical simulation, we consider α = β = γ = 0.01. This is similar to the case of the small dispersion limit where the weakly nonlinear terms are dominant. Figure 11 demonstrates that the solution becomes a dispersive shock. This phenomenon has been previously observed in dispersive systems, cf. [START_REF] Dutykh | Finite volume schemes for dispersive wave propagation and runup[END_REF][START_REF] El | Unsteady undular bores in fully nonlinear shallow-water theory[END_REF]23]. Unexpectedly, taking the parameters β and γ to be very small, i.e. β = γ = 0.001 and keeping the parameter α = 1, the solitary wave persists and evolves into a new solitary wave which is similar to the unperturbed solitary wave qualitatively similar to those presinted in Section 4.1. This persistence is remarkable because the solitary wave remains almost the same even if two of the most important terms have been almost eliminated. If α = β = 1 and γ = 0.001 or if α = γ = 1 and β = 0.001 the behavior is similar.

The behavior changes dramatically if large perturbations α are considered. The results from the simulation with α = 0.001 and β = γ = 1 is shown in Figure 12. In this case, the initial condition breaks into different waves but instead of forming a dispersive shock wave, it forms a new kind of regularized shock wave. This suggests a new breaking mechanism by the elimination of the u xt term. Similar dissipative behavior has been observed in nonlinear KdV-type equations where some high-order nonlinear terms introduce dissipation to the system [6].

Stability of Cnoidal waves

We follow the work of Carter & Cienfuegos [START_REF] Carter | The kinematics and stability of solitary and cnoidal wave solutions of the Serre equations[END_REF] in order to study the linear stability of the solutions given in (1.5). We enter a coordinate frame moving with the speed of the solutions by defining χ = xc c t and τ = t. In this moving frame, the Serre equations are given by

h τ -c c h χ + hu χ = 0 , u τ -c c u χ + uu χ + h χ -1 3h h 3 u χτ -c c u χχ + uu χχ -(u χ ) 2 χ = 0 , (4.5)
and the solution given in (1.5) simplifies to the following time-independent solution

h = h 0 (χ) = a 0 + a 1 dn 2 K c χ, k , u = u 0 (χ) = c c 1 -h 0 h(χ) . (4.6)
We consider perturbed solutions of the form

h pert (χ, τ ) = h 0 (χ) + µh 1 (χ, τ ) + O(µ 2 ) , u pert (χ, τ ) = u 0 (χ) + µu 1 (χ, τ ) + O(µ 2 ) , (4.7)
where h 1 and u 1 are real-valued functions and µ is a small real parameter. Substituting (4.7) into (4.5) and linearizing leads to a pair of coupled, linear partial differential equations that are constant coefficient in τ . Without loss of generality, assume

h 1 (χ, τ ) = H(χ)e Ωτ + c.c. , u 1 (χ, τ ) = U(χ)e Ωτ + c.c. , (4.8)
where H(χ) and U(χ) are complex-valued functions, Ω is a complex constant, and c.c. denotes complex conjugate. If Ω has a positive real part, i.e. if Re(Ω) > 0, then the perturbations h 1 and u 1 grow exponentially in τ and the solution is said to be unstable. Substituting (4.8) into the linearized PDEs gives

L H U = Ω M H U , (4.9)
where L and M are the linear differential operators defined by where prime represents derivative with respect to χ and

L = -u ′ 0 + (c c -u 0 )∂ χ -η ′ 0 -η 0 ∂ χ L 21 L 22 , M = 1 0 0 1 -η 0 η ′ 0 ∂ χ -1 3 η 2 0 ∂ χχ , (4.10)
L 21 = -η ′ 0 (u ′ 0 ) 2 -c c η ′ 0 u ′′ 0 - 2 3 c c η 0 u ′′′ 0 + η ′ 0 u 0 u ′′ 0 - 2 3 η 0 u ′ 0 u ′′ 0 + 2 3 η 0 u 0 u ′′′ 0 + η 0 u 0 u ′′ 0 -g -η 0 (u ′ 0 ) 2 -c c η 0 u ′′ 0 ∂ χ , (4.11)
L 22 = -u ′ 0 + η 0 η ′ 0 u ′′ 0 + 1 3 η 2 0 u ′′′ 0 + c c -u 0 -2η 0 η ′ 0 u ′ 0 - 1 3 η 2 0 u ′′ 0 ∂ χ + η 0 η ′ 0 u 0 -c c η 0 η ′ 0 - 1 3 η 2 0 u ′ 0 ∂ χχ + 1 3 η 2 0 u 0 - 1 3 c c η 2 0 ∂ χχχ . (4.
12)

The Fourier-Floquet-Hill method described in Deconinck & Kutz [START_REF] Deconinck | Computing spectra of linear operators using the Floquet-Fourier-Hill method[END_REF] is then used to solve the differential eigenvalue problem given in (4.9). This method establishes that all bounded solutions of (4.9) have the form

H U = e iρχ H P U P , (4.13)
where H P and U P are periodic in χ with period 2K K c and ρ ∈ -πK c (4K), πK c (4K) . Using this method, Carter & Cienfuegos established that solutions of the form given in (4.6) with sufficiently small amplitude and steepness are spectrally stable and solutions with sufficiently large amplitude or steepness are spectrally unstable. For example, the solution with a 0 = 0.3, a 1 = 0.2 and k = 0.75 is unstable with respect to the perturbation shown in Figure 13. The period of this perturbation is twelve times the period of the exact solution (ρ = 1 12). The theory establishes that the magnitude of this perturbation will grow like e 0.00569t . We corroborated this result by using the following perturbation-seeded solution as the initial condition in our Serre solver h pert (x, 0) = h(x, 0) + 10 -7 h 1 (x) , u pert (x, 0) = u(x, 0) + 10 -7 u 1 (x) . Here h(x, 0), u(x, 0) is the solution given in equations (4.6) with a 0 = 0.3, a 1 = 0.2 and k = 0.75 and h 1 (x), u 1 (x) is the perturbation shown in Figure 13. Figure 14 contains a plot of the magnitude of the first Fourier mode of the solution versus t. This mode initially (up to t = 1500) grows exponentially with a rate of 0.00577, very close to the rate predicted by the linear theory. However after more time, the solution returns to a state close to the initial one. The first portion of this recurrence phenomenon is depicted in Figure 15. Note that the solution at t = 2870 has nearly returned to its initial state. Similar behavior has been observed to other shallow water models by Ruban [30] and it is referred to as the Fermi-Pasta-Ulam recurrence.

Dispersive shock waves

A simple DSW traveling to the right can be generated using the Riemann initial data, cf. [START_REF] El | Unsteady undular bores in fully nonlinear shallow-water theory[END_REF],

h(x, 0) = h -, for x < 0 h + , for x > 0 , u(x, 0) = u -, for x < 0 u + , for x > 0 , (5.1)
with the compatibility condition (Riemann invariant)

u - 2 - √ h -= u + 2 - √ h + . (5.2)
DSWs can also be generated during the dam-break problem simulation. In this case, the initial data for h(x, 0) are the same as in (5.1), but there is no flow at t = 0, i.e. , u(x, 0) = 0. As shown in [START_REF] El | Unsteady undular bores in fully nonlinear shallow-water theory[END_REF], this generates two counter-propagating DSWs, one on each side of the "dam", and two rarefaction waves that travel toward the center. We consider the initial condition for h to be a smooth step function that decays to zero as x → ∞. Specifically, we choose

η(x, 0) = 1 2 η 0 1 + tanh x 0 -x 2 , (5.3)
where η 0 = 0.1, x 0 = 350, and u(x, 0) = 0. A plot of this initial condition is included in Figure 16. Both the Euler and Serre equations generate two counter propagating DSWs and two rarefaction waves. Figure 16 demonstrates that the amplitude of the leading wave for both solutions is almost the same. For example, the amplitude of the Euler leading wave at t = 200 is A = 0.06372 while the amplitude of the Serre leading wave at t = 200 is A = 0.06356. Although the leading waves have almost the same amplitudes, the phase speeds are slightly different. The difference in phase speeds is demonstrated in Figure 16.

After verifying that the Serre system has dispersive shock waves that are comparable with the full Euler equations, we examine the interactions of simple DSWs starting with the head-on collision. For the head-on collision we again consider two initial waveforms similar to (5.1) but translated as is shown in Figure 17(a). These step functions generate two counter-propagating waves that begin to interact at approximately t = 27. The collision is inelastic. After the collision there are two DSWs propagating in different directions on the trailing edge of the DSWs. We now consider overtaking collisions of DSWs. For this situation we consider doublestep initial conditions as is shown in Figure 18. The first step has amplitude 0.1 while the shorter step has amplitude 0.05. This initial condition generate two DSWs that propagate to the right. Because shorter DSWs propagate with smaller phase speeds than taller DSWs, the taller DSW approaches the shorter one and they interact. The interaction is so strong that the symmetry of the leading wave of both DSWs is destroyed. The two waves appear to merge and propagate as one single-phase DSW. Similar behavior has been observed in NLS-type and KdV-type equations [START_REF] Ablowitz | Interactions and asymptotics of dispersive shock waves -Korteweg-de Vries equation[END_REF][START_REF] Hoefer | Interactions of dispersive shock waves[END_REF]. Finally, we mention that the solutions shown in Figures 17 and18 are magnifications of the actual solutions. The rest of the solution, not shown in these figures, consists of dispersive rarefaction waves that we do not study in this paper. For more information see [START_REF] Mitsotakis | On the Galerkin/Finite-Element Method for the Serre Equations[END_REF].

1 Figure 1 .

 11 Figure 1. The evolution of an Euler solitary wave with c s = 1.1 when used as initial condition in the Serre system.

2 Figure 2 .

 22 Figure 2. The evolution of an Euler solitary wave with c s = 1.2 when used as initial condition in the Serre system.

1 Figure 3 .

 13 Figure 3. The amplitude, shape and phase error of the Euler's solitary wave of c s = 1.1 propagating with Serre equations. See also Figure 1.

Figure 4 .

 4 Figure 4. The amplitude, shape and phase errors of the Euler's solitary wave of c s = 1.1 propagating with Serre equations. See also Figure 2.

5 Figure 5 .

 55 Figure 5. Comparison of the head-on collision of two solitary waves of the Serre system with experimental data.

Figure 6 .

 6 Figure 6. (Cont'd) Comparison of the head-on collision of two solitary waves of the Serre system with experimental data.

Figure 7 .

 7 Figure 7. Comparison between the head-on collisions of two solitary waves for the Serre and Euler systems.

Figure 8 .

 8 Figure 8. The perturbed solution generated by the perturbation of the amplitude of the solitary wave with c s = 1.4, A = 0.96, and perturbation parameter p = 1.1.

ηFigure 9 .

 9 Figure 9. The solitary wave perturbed by pseudo-random noise with p = 0.2.

Figure 10 .

 10 Figure 10. Evolution of the initial condition shown in Figure 9 with p = 0.2.

Figure 11 .

 11 Figure 11. The evolution of a solitary wave of a perturbed system (α = β = γ = 0.01).

Figure 12 .

 12 Figure 12. The evolution of a solitary wave of a perturbed system and the generation of a regularized shock (α = 0.001, β = γ = 1).

1 Figure 13 .

 113 Figure 13. Unstable perturbations corresponding to the solution given in (1.5) with a 0 = 0.3, a 1 = 0.2 and k = 0.75.

Figure 14 .

 14 Figure 14. The magnitude of the first Fourier mode of the perturbation versus time.

Figure 15 .

 15 Figure 15. The periodic instability of the perturbed cnoidal wave.

Figure 16 .

 16 Figure 16. The dam break problem.

Figure 17 .Figure 18 .

 1718 Figure 17. Head-on collision of two simple DSWs.

Table 1 .

 1 Amplitudes of the Euler and Serre solitary waves corresponding to different speeds.

	s	Euler Serre
	1.01 0.02012 0.0201
	1.05 0.10308 0.1025
	1.1 0.21276 0.2100
	1.15 0.33007 0.3225
	1.2 0.45715 0.4400
	1.28 0.70512 0.6384

Acknowledgments

D. Mitsotakis thanks Professor Boaz Ilan for suggestions, comments, and stimulating discussions related to dispersive waves. The authors acknowledge the invaluable help of Professor Paul Milewski for discussions related to the numerical schemes for the Euler equations and Professor Didier Clamond for discussions on pseudo-spectral methods. J. Carter was supported by the National Science Foundation under grant number DMS-1107476. D. Mitsotakis was supported by the Marsden Fund administered by the Royal Society of New Zealand.


Eric Lacot
email: eric.lacot@ujf-grenoble.fr

Bahram Houchmandzadeh

Vadim Girardeau

Olivier Hugon

Olivier Jacquin

Nonlinear modification of the laser noise power spectrum induced by frequency-shifted optical feedback

I. INTRODUCTION

Laser properties (power, polarization, coherence, dynamical behavior, etc.) can be significantly affected and modified by optical feedback [START_REF] Lang | External optical feedback effects on semiconductor injection laser properties[END_REF][START_REF] Erneux | Nonlinear dynamics of an injected quantum cascade laser[END_REF]. These properties have been widely studied in the particular case of weak optical feedback. For example, it has been shown that the steady state intensity of a laser subjected to coherent optical feedback from the reflection on an external surface depends on the reflectivity, distance, and motion of the surface. This led to the laser feedback interferometry technique [START_REF] Taimre | Laser feedback interferometry: A tutorial on the self-mixing effect for coherent sensing[END_REF]. It also has been demonstrated that the dynamical behavior of a laser can be several orders of magnitude more sensitive to optical feedback than its steady state properties. Since the pioneering work of Otsuka on self-mixing modulation effect in a class-B laser [START_REF] Otsuka | Effects of external perturbations on LiNdP 4 O 12 lasers[END_REF], the high dynamical sensitivity of lasers to frequencyshifted optical feedback has been used in metrology [START_REF] Otsuka | Self-mixing thin-slice solid-state laser metrology[END_REF], for example in self-mixing laser Doppler velocimetry [START_REF] Okamoto | Ultrahighly sensitive laser-Doppler velocity meter with a diode-pumped Nd:YVO 4 microchip laser[END_REF][START_REF] Kawai | Ultrahigh-sensitivity self-mixing laser Doppler velocimetry with laser-diode-pumped microchip LiNdP 4 O 12 lasers[END_REF][START_REF] Suddo | Determination of velocity of self-mobile phytoplankton using a self thin-slice solid-state laser[END_REF] and in laser optical feedback imaging (LOFI) [START_REF] Lacot | Laser optical feedback tomography[END_REF][START_REF] Hugon | Coherent microscopy by laser optical feedback imaging (LOFI) technique[END_REF][START_REF] Glastre | Demonstration of a plenoptic microscope based on laser optical feedback imaging[END_REF]. Compared to conventional optical heterodyne detection, frequency-shifted optical feedback allows for several orders of magnitude higher intensity modulation contrast [START_REF] Lacot | Coherent laser detection by frequency-shifted optical feedback[END_REF][START_REF] Lacot | Comparative study of autodyne and heterodyne laser interferometry for imaging[END_REF][START_REF] Jacquin | Experimental comparison of autodyne and heterodyne laser interferometry using an Nd:YVO 4 microchip laser[END_REF]. In the case of weak optical feedback, the laser dynamics is linear and the maximum of the modulation is obtained when the frequency shift is resonant with the laser relaxation oscillation frequency. In this condition, an optical feedback level as low as -170 dB (i.e., 10 17 times weaker than the laser intracavity power) has been detected [START_REF] Okamoto | Ultrahighly sensitive laser-Doppler velocity meter with a diode-pumped Nd:YVO 4 microchip laser[END_REF]. When the optical feedback becomes stronger, nonlinearities appear in the laser dynamics [START_REF] Tan | Power spectral characteristic of a microchip Nd:YAG laser subjected to frequency-shifted optical feedback[END_REF][START_REF] Tan | Response of microchip solid state laser to external frequency-shifted feedback and its applications[END_REF]. They can cause the apparition of chaos, bistability, and hysteresis phenomenon with the tuning (back and forth) of the frequency shift of the optical feedback [START_REF] Schwartz | Subharmonic hysteresis and period doubling bifurcations for periodically driven laser[END_REF].

The main objective of the present work is to study the modification of the noise power spectrum of the laser and more specifically its resonance frequency (i.e., the laser relaxation frequency) induced by the nonlinear-laser dynamics in the strong feedback situation. A better understanding of the laser nonlinear behavior could lead to a new generation of laser metrology techniques with improved performances. This article is organized as follows. In Secs. II and III, we recall the rate equations governing the dynamics of a laser submitted to frequency-shifted optical feedback. These equations are solved numerically to show characteristic examples of the modification of the laser noise power spectrum under different optical feedback conditions. In particular, the influence of (i) the amount of light returning to the laser cavity and (ii) the detuning between the frequency shift and intrinsic relaxation frequency are studied. Section III is devoted to analytical resolution of the rate equation: using a bifurcation analysis, we determine the amplitude of the laser output power modulation in the strong feedback situation. We show then how the nonlinear dynamical coupling links the value of the relaxation frequency to the strength of the output power modulation. Section IV is devoted to experimental results where we show their good agreement with the theoretical predictions of the preceding sections. The final section is devoted to the general discussion of these results and to their prospective applications.

For weak optical feedback (R e 1) and a short round-trip time delay (τ e 1/F e ), the dynamical behavior of a laser with frequency-shifted (F e ) optical feedback can be described by the following set of differential equations [START_REF] Lacot | Coherent laser detection by frequency-shifted optical feedback[END_REF]:

dI dt = BI N -γ c I + γ c 2 R e I cos ( e t + e ) + F I (t), (1a) dN dt = γ 1 [N 0 -N ] -BNI, ( 1b
)
F I (t) = 0 and F I (t)F I (t -τ ) = 2γ c I δ(τ ), (1c)
where I and N are, respectively, the laser intensity (photon unit) and the population inversion (atom unit). γ 1 is the decay rate of the population inversion, γ c is the laser cavity decay rate, γ 1 N 0 is the pumping rate, and B is related to the Einstein coefficient (i.e., the laser cross section). In Eq. (1a), the cosine function expresses the coherent interaction (i.e., the beating at the angular frequency: e = 2πF e ) between the lasing and the feedback electric field. The optical feedback is characterized by the effective power reflectivity R e and the optical phase shift e = ω c τ e induced by the round-trip time between the laser and the target (where ω c is the optical pulsation). Regarding the noise, the laser quantum fluctuations are described by the Langevin noise function F I (t), with a zero mean value and a white noise type correlation function [Eq. (1c)] [START_REF] Petermann | Laser diode modulation and noise[END_REF][START_REF] Kolobov | Role of pumping statistics and dynamics of atomic polarization in quantum fluctuations of laser sources[END_REF].

In the absence of optical feedback (R e = 0) and noise [F I (t) = 0], the steady state of Eqs. [1(a)-1(c)] is given by

N S = γ c /B, ( 2a
)
I S = I sat [η -1], ( 2b
)
where η = N 0 /N S is the normalized pumping parameter and I sat = γ 1 /B is related to the saturation intensity of the laser transition.

In this regime, the intrinsic dynamics of a class-B laser (γ c γ 1 η) is characterized by damped relaxation oscillations of the laser output power with a relaxation angular frequency R = √ γ 1 γ c (η -1) and a damping rate R = γ 1 η/2. Experimentally, this transient dynamics is constantly excited by the laser quantum noise described by the Langevin force F I (t).

In the presence of a strong optical feedback, the laser intrinsic dynamics, and in particular the relaxation oscillation frequency, is modified and will depend on the modulation conditions (R e , e ). In this article, "strong feedback" corresponds to the regime where the modulation frequency is nearly resonant ( e ≈ R ) and where the amount of optical feedback (R e ) is high enough to induce nonlinear dynamical behaviors (which give rise to the generation of harmonic and parametric dynamical frequencies) in the laser output power modulation [START_REF] Lacot | Comparative study of autodyne and heterodyne laser interferometry for imaging[END_REF][START_REF] Tan | Power spectral characteristic of a microchip Nd:YAG laser subjected to frequency-shifted optical feedback[END_REF][START_REF] Tan | Response of microchip solid state laser to external frequency-shifted feedback and its applications[END_REF]. In contrast, the "weak feedback" regime corresponds to the situation where the modulation frequency is far away from the resonance (| e -R | 0) and where the amount of optical feedback (R e ) is small enough to induce only linear dynamical behaviors in the laser output power modulation.

B. Numerical results

Using a standard Runge-Kutta method, we have solved a normalized form of Eqs. (1a)-1(c), where the laser intensity is divided by the saturation intensity (i = I /I sat ), the population inversion is divided by its stationary value (n = N/N S ), and the time is multiplied by the cavity damping rate (t = γ c t). The normalized Langevin force is a random variable with zero mean value and a standard deviation equal to 1/ √ I sat dt where the integration step is dt = 40. The number of integration steps is equal to 16 384.

Figure 1 shows the rf power spectra of the laser output power modulation for two different optical feedback conditions (R e , e ). When the feedback is weak [Fig. 1(a)], one can observe that the laser dynamics is principally composed of the superposition of the laser output power modulation at the modulation frequency F e and of the noise power spectrum related to the transient relaxation oscillations, with a resonant frequency F R = R /2π and a half width at half maximum R /2π . In the strong feedback regime [Fig. 1(b)], the power spectrum is composed of peaks at the modulation frequency (F e ) and its harmonics (2F e , . . . ), while the noise generates two sidebands beside each peak, located at F e ± R ,2F e ± R , . . . . In this study, we define the new value of the relaxation frequency FR as the position of the highest noise peak in the vicinity of F e . In Fig. 1(b), it corresponds to the left sideband: FR = F e -R [(2) in Fig. 1(b)]. One can notice that all other peaks can be obtained from a linear combination of F e and FR . For example, the right sideband of F e corresponds to F e + R = 2F e -FR [(4) in Fig. 1(b)]. The comparison of Figs. 1(a) and 1(b) shows how the noise power spectrum is modified by the strong optical feedback at the modulation frequency (F e ) and how the new relaxation frequency ( FR ) is moved down from its intrinsic position (F R ).

By comparison with the work made in [START_REF] Tan | Power spectral characteristic of a microchip Nd:YAG laser subjected to frequency-shifted optical feedback[END_REF][START_REF] Tan | Response of microchip solid state laser to external frequency-shifted feedback and its applications[END_REF], which is focused on the study of the signal power spectrum, the present work analyzes how the noise power spectrum of the laser is modified and how the relaxation frequency of the laser is shifted, when the laser is submitted to a strong optical feedback. For a given modulation frequency, the left column of Fig. 2 shows how FR (i.e., the noise left sideband) is shifted to lower frequencies when the amount of optical feedback increases. For a given amount of optical feedback, the right column of Fig. 2 shows how FR is pushed down to a lower frequency when the modulation frequency approaches the intrinsic laser relaxation frequency.

Also, one can observe that in the low feedback condition the noise power spectrum exhibits a single resonance [Figs. 1(a

C. Bifurcation analysis

To study the dynamical response of a laser subject to frequency-shifted optical feedback (R e = 0), we have made a bifurcation analysis using an asymptotic approximation similar to [START_REF] Schwartz | Subharmonic hysteresis and period doubling bifurcations for periodically driven laser[END_REF][START_REF] Erneux | Stabilizing or destabilizing lasers by continuous delayed feedback[END_REF][START_REF] Witomski | Parametric amplification of frequency-shifted feedback[END_REF]. As already mentioned, we focus our analysis on the case e ≈ R , where R is the intrinsic relaxation oscillation frequency of the laser. For a more convenient theoretical analysis, we have reformulated the set of Eqs. (1a)-(1c) by using the new variables (x and y), the new time (s), and the new modulation parameters (δ e ,σ e ) defined by [START_REF] Schwartz | Subharmonic hysteresis and period doubling bifurcations for periodically driven laser[END_REF][START_REF] Witomski | Parametric amplification of frequency-shifted feedback[END_REF]

I = I S (1 + y), N = N S 1 + R γ c
x , s = R t, (3a)

σ e = e R , δ e = 2 R e γ c R . ( 3b
)
The new variables x and y are, respectively, the deviation from the nonzero intensity laser steady state and the population inversion steady state. In the absence of noise, the set of Eqs. (1a) and (1b) can then be rewritten as

dx ds = -y -εx[1 + (η -1)(1 + y)], ( 4a
)
dy ds = (1 + y)[x + δ e cos (σ e s + e )]. (4b)
For our microchip laser, ε = γ 1 / R ≈ 10 -3 is a small quantity which motivates an asymptotic analysis of Eqs. (4a) and (4b). For σ e ≈ 1, the leading approximation of the solutions for the deviations x and y from the nonzero intensity state can be written as

x = A exp (iS) -i A 2 3 exp (i2S) + c.c., (5a)
y = -iA exp (iS) - 2 3 A 2 exp (i2S) + c.c., (5b)
where S = σ e s = e t. The complex amplitude A = 0(ε 1/2 ) of the small periodic oscillation is a function of the slow time εs and satisfies the following equations [START_REF] Schwartz | Subharmonic hysteresis and period doubling bifurcations for periodically driven laser[END_REF][START_REF] Witomski | Parametric amplification of frequency-shifted feedback[END_REF]:

dA ds = -i[σ e -1]A -i A 2 A * 6 - ηε 2 A + i δ e 4 exp (i e ). ( 6
)
Note that our asymptotic approximation is valid in the limit ε → 0 assuming that √ R e and |σ e -1| are at maximum 0(ε) quantities.

In the right-hand side of Eqs. (5a) and (5b), the first term corresponds to the laser output power modulation at the modulation frequency (i.e., at the cavity loss modulation frequency). The second term corresponds to the first harmonic modulation (i.e., nonlinear-laser dynamics) induced by the nonlinear coupling between the laser intensity and the laser population inversion through the stimulated emission.

III. ASYMPTOTIC SOLUTIONS

We now investigate the solutions of Eq. ( 6), and show how the intrinsic laser dynamics (and in particular the laser relaxation frequency) can be influenced by the optical feedback through the modulation conditions (δ e ,σ e ). In this study, the key is the third-order term (-iA 2 A * /6).

A. Linear case (

A 2 A * ≈ 0)
When the product A 2 A * can be neglected in Eq. ( 6), the solution of the amplitude equation can be written as small transient (with subscript T) relaxations of random excitations (due to quantum noise) around the stationary solution (with subscript S):

A(s) = A S,L + B T ,L exp (-βs), (7a)
where A S,L = iδ e exp (i e ) 4i(σ e -1) + 2ηε , (7b)

β = i(σ e -1) + ηε/2, ( 7c
)
and the B T ,L factor is given by the amplitude of the quantum noise.

By combining these results with Eqs. (3a) and (3b) and (5a) and (5b), we find that the stationary laser output power modulation is then given by

I S y S,L = 2 R e I S |G L ( e )| sin[ e t + arg(A S,L )], (8a) |G L ( e )| = γ c ( e + η γ 1 ) 2 R -2 e + i e η γ 1 ≈ γ c 4( e -R ) 2 + (η γ 1 ) 2 , ( 8b
)
where G L ( e ) is the amplification gain of the linear laser dynamics [START_REF] Lacot | Comparative study of autodyne and heterodyne laser interferometry for imaging[END_REF][START_REF] Jacquin | Experimental comparison of autodyne and heterodyne laser interferometry using an Nd:YVO 4 microchip laser[END_REF].

The above equation shows that the laser output power modulation exhibits a strong resonance (and therefore a high sensitivity to optical feedback) when e = R . The linear regime is valid when |A S,L | 1, which implies R e R e, lim where R e, lim = (γ 1 η/2γ c ) 2 .

The transient laser output power modulation can be obtained from the relaxing part of Eq. ( 7):

I S y T ,L = 2I S |B T ,L | exp - γ 1 η 2 t sin[ R t + arg(B T ,L )]. (9)
Experimentally, this transient dynamic is constantly excited by the laser quantum noise described by the Langevin forces [F I (t) in Eqs. (1a)-(1c)]. The rf power spectrum of the laser noise is obtained from the Fourier transform (F) of Eq. ( 9):

P S L ( ) ∝ |F [I S y T ,L (t)]| 2 ≈ [I S 2|B T ,L |/γ c ] 2 |G L ( )| 2 . ( 10
)
In previous studies [START_REF] Lacot | Coherent laser detection by frequency-shifted optical feedback[END_REF][START_REF] Lacot | Comparative study of autodyne and heterodyne laser interferometry for imaging[END_REF][START_REF] Jacquin | Experimental comparison of autodyne and heterodyne laser interferometry using an Nd:YVO 4 microchip laser[END_REF], we had already demonstrated that the detection of frequency-shifted optical feedback is shot noise limited and that the rf power spectrum is given by

P S L ( ) ≈ 2I out |G L ( )| 2 2 F, ( 11
)
where I out = γ c I S is the laser output power and F is the detection bandwidth which is supposed to be smaller than the resonance width ( F R /2π ). Equations [START_REF] Hugon | Coherent microscopy by laser optical feedback imaging (LOFI) technique[END_REF] and [START_REF] Glastre | Demonstration of a plenoptic microscope based on laser optical feedback imaging[END_REF] show that shape of the noise power spectrum is linked to the amplification gain of the laser dynamics [G L ( )] and that the amplitude of the intensity relaxation oscillations induced by the laser quantum noise is proportional to the square root of the laser output power:

I S |B T ,L | ∝ 2I out 2 F . ( 12
)

B. Nonlinear case (

A 2 A * = 0)
In this section, we study the coupling of the stationary modulation of the laser output power with the laser transient dynamics (i.e., the laser quantum noise) through the nonlinear cubic term (-A 2 A * /6) of Eq. ( 6). More precisely, the aim of this study is to determine how the intrinsic laser dynamics (and in particular the laser relaxation frequency) is modified by the optical modulation conditions induced by the optical feedback. When the product A 2 A * cannot be neglected in Eq. ( 6), the solution of the amplitude equation can be written as

A(s) = A S,NL + B T ,NL (s), (13)
where A S,NL is the stationary solution and B T ,NL (s) the transient dynamics of random excitations (due to quantum noise) around the stationary solution in the nonlinear regime (subscript NL).

In agreement with our numerical simulations [Figs. 1(b) and 2], we assume that the laser output power modulation is stronger than the laser quantum noise [|A S,NL | |B T ,NL (s)|]. Keeping only the first-order nonlinear terms in Eq. ( 13) gives

A 2 A * ≈ A S,NL A S,NL A * S,NL + 2A S,NL A * S,NL B T ,NL + A S,NL A S,NL B * T ,NL . ( 14
)
The insertion of Eqs. ( 13) and ( 14) into Eq. ( 6) gives the following equalities:

i δ e 4 exp (i e ) = i ηε 2 + σ e -1 A S,NL + i A 2 S,NL A * S,NL 6 , ( 15a
)
dB T ,NL ds = - ηε 2 + i(σ e -1) B T ,NL + i A S,NL A * S,NL 3 B T ,NL -i A 2 S,NL 6 B * T ,NL . (15b)
Equations (15a) and (15b) clearly show that the stationary and the transient dynamics (and therefore the laser noise) are coupled through the nonlinear-laser dynamics.

Nonlinear stationary solutions

Taking the complex conjugate of Eq. (15a) and multiplying them together, we find the square modulus R 2 S,NL = A S,NL A * S,NL of the modulation amplitude:

R 6 S,NL 36 + 2(σ e -1) 6 R 4 S,NL + (σ e -1) 2 + ηε 2 2 R 2 S,NL - δ e 4 2 = 0, (16)
which is a third-order equation in R 2 S,NL . For a given amount of optical feedback δ e , the above equation has three real positive roots when the frequency shift is in a range defined by (see [START_REF] Schwartz | Subharmonic hysteresis and period doubling bifurcations for periodically driven laser[END_REF] or the Appendix for details)

- δ e /(4 √ 6) ηε/2 2 < (σ e -1) < - 3 2 2 / 3 δ e 4 √ 6 2/3 . ( 17
)
The jump between the two stable solutions allows one to observe a hysteresis phenomenon which has already been studied in [START_REF] Schwartz | Subharmonic hysteresis and period doubling bifurcations for periodically driven laser[END_REF]. This kind of hysteresis has also been observed in the dynamics of a laser submitted to an optical injection [START_REF] Otsuka | Optical bistability and nonlinear resonance in resonant-type semiconductor laser amplifier[END_REF] or a modulation of the pumping power [START_REF] Phillps | Optical chaos and hystetresis in a laser-diode pumped Nd doped fibre laser[END_REF]. More specifically, the effect of noise on the size of the hysteresis zone has been studied in [START_REF] Celet | Slowly Passing through Resonance Strongly Depends on Noise[END_REF].

Let us recall that the goal of the present article is not to study the hysteresis phenomenon, but to determine how the noise power spectrum of the laser and more specifically the laser relaxation frequency are related to the strength of R 2 S,NL .

Nonlinear relaxation oscillations

To determine the noise power spectrum we have studied the transient dynamics of random excitations (B T ,NL ).

Writing B T ,NL (s)= exp[-(ηε/2)s][X(s)+iY (s)], Eq. (15b) gives

d ds X Y = β Y ( -β X ) -( + β X ) -β Y X Y , ( 18
)
where

β X + iβ Y = A 2 S,NL /6, = (σ e -1) + R 2 S,NL /3 and R 2
S,NL is the square modulus of the nonlinear-laser output power modulation given by Eq. ( 16).

The two eigenvalues (λ) of the above matrix obey the equation

λ 2 = R 2 S,NL 6 2 -(σ e -1) + R 2 S,NL 3 2 . ( 19
)
For the two stable solutions, one has λ 2 < 0, then λ ∈ C and λ = ±iδ R,NL , where δ R,NL allows one to measure the detuning of the transient frequency compared to the modulation frequency:

δ R,NL = σ e -σR = e R - ˜ R R , ( 20
)
where ˜ R is the new value of the laser relaxation frequency modified by the nonlinear coupling between the transient and the permanent dynamics. Note that knowing σ e (which is experimentally controlled) and measuring the shift of the lateral sidebands ±δ R,NL = ± R /F R [see Fig. 1(b)], one could determine the relative variation of the laser relaxation oscillation frequency:

˜ R -R R = σ e -1 -δ R,NL . ( 21
)
To make a simple physical analysis of this result, let us look at a particular situation where the optical feedback is not too strong [R 2

S,NL

(σ e -1)]. In this condition, Eq. ( 19) simply gives

δ R,NL ≈ (σ e -1) + R 2 S,NL 3 , ( 22
)
and using Eq. ( 19), one obtains the relative variation of the laser relaxation frequency:

˜ R -R R ≈ - R 2 S,NL 3 . ( 23
)
Equation ( 23) clearly shows that the frequency shift of the laser relaxation frequency is always negative and is proportional to the square of the amplitude of the laser output The lines correspond to results obtained from analytical equations [Eq. ( 16)], while the dots correspond to results obtained from numerical simulations. Solid lines and square dots (R e = 1 × 10 -7 ); dotted lines and circular dots (R e = 1 × 10 -8 ). Laser parameters are identical to Fig. 1. power modulation. Therefore the frequency shift is important for a strong optical feedback (R 2 S,NL ≈ 1), while it is roughly equal to zero for a weak optical feedback (R 2 S,NL 1). One can notice that due to the fact that R 2 S,NL (σ e ) can exhibit hysteresis [START_REF] Schwartz | Subharmonic hysteresis and period doubling bifurcations for periodically driven laser[END_REF], ˜ R (σ e ) can also exhibit hysteresis.

Figure 3 shows (i) the relative variation of the relaxation frequency [Fig. 3(a)] and (ii) the variation of the amplitude of the laser output power modulation [Fig. 3(b)] versus the detuning of the modulation frequency, for two different values of the optical feedback (R e = 1 × 10 -7 ,R e = 1 × 10 -8 ).

Figure 3 shows relatively good agreement between the analytical and the numerical results for both the amplitude of the laser output power modulation and the resonance frequency of the noise power spectrum (i.e., the laser relaxation frequency). In agreement with Eq. ( 23), one can observe that the decrease of the relaxation frequency is proportional to the increase of the laser amplitude modulation.

In Fig. 3, to avoid confusion, the hysteresis phenomenon is not shown (the frequency shift is not swept back and forth). The square and the circle correspond simply to two different values of the optical feedback. Note also that the results shown Finally, the rf power spectrum of the laser noise is obtained from the Fourier Transform (F) of B T ,NL :

P S NL ( ) ∝ |F[I S y T ,NL (t)]| 2 ≈ [I S 2|B T ,NL (0)|/γ c ] 2 |G NL ( )| 2 , ( 24a
)
|G NL ( )| ≈ γ c (i + γ 1 η) ( e -2π R ) 2 -2 + i γ 1 η + g NL γ c (i + γ 1 η) ( e + 2π R ) 2 -2 + i γ 1 η , (24b) g NL ≈ ⎛ ⎝ 2 R 2 S,NL 6 4 σ e -1 + 2 R 2 S,NL 6 ⎞ ⎠ 1 / 2 , ( 24c
)
where G NL ( ) is the amplification gain (i.e., the modulation transfer function) of the nonlinear-laser dynamics. This gain is composed of two resonance frequencies symmetrically located on both sides of the frequency shift: (F e ± R ), where g NL ≈ |G NL ( e -2π R )| / |G NL ( e + 2π R )| is the ratio between the two maxima which depends on the feedback conditions. In agreement with the numerical results shown in Fig. 2, Eqs. (24a)-(24c) show that the amplitude of the two maxima of the noise power spectrum are of the same order of magnitude for strong optical feedback (i.e., when R 2 S,NL ≈ 1). Inversely, the right maximum disappears for weak optical feedback (i.e., when R 2 S,NL 1) leading to the conventional linear gain with only one resonance frequency ( R ):

G NL ( ) ≈ R 2 S,NL →0 G L ( ).
One can notice an average • in Eq. (24a). This average is due to the fact that B T ,NL (s) is the transient dynamics of random excitations (due to quantum noise). With a uniform distribution of initial phase condition, this average leads to an incoherent interaction [modulus sum in Eq. (24b)] between the two resonance curves.

IV. EXPERIMENTAL RESULTS

A. Experimental setup

To study the nonlinear dynamics of a laser submitted to frequency-shifted optical feedback, we have used a LOFI setup [START_REF] Lacot | Laser optical feedback tomography[END_REF]. A schematic diagram of this setup is shown in Fig. 4. The laser is a diode pumped Nd:YAG microchip laser. The maximum available pump power is 400 mW at 810 nm and the total output power of the microchip laser is 80 mW with a central wavelength of λ = 1064 nm. This laser has a plane-parallel cavity which is stabilized by the thermal lens induced by the Gaussian pump beam. The two dielectric mirrors are directly coated on the laser material (full cavity). The input dichroic mirror transmits the pump power and totally reflects the infrared laser wavelength. On the other side, the dichroic output mirror allows one to totally reflect the pump power (to increase the pump power absorption and therefore the laser efficiency) and only partially reflects (95%) the laser wavelength. The microchip cavity is relatively short, L c ≈ 1 mm, which ensures a high damping rate of the cavity and therefore a good sensitivity to optical feedback. Part of the light diffracted and/or scattered by the target returns inside the laser cavity after a second pass through the frequency shifter. Therefore, the optical frequencies of the reinjected light are shifted by F e . This frequency can be adjusted and is typically of the order of the laser relaxation frequency F R , which is in the megahertz range for the microchip laser used in this study.

The optical feedback is characterized by the complex target reflectivity (r e = √ R e exp(j e ), where the phase e = (2π/λ)d e describes the optical phase shift induced by the round-trip time delay τ e (i.e., the distance d e = c τ e , where c is the velocity of light) between the laser and the target. The effective power reflectivity (R e = |r e | 2 ) takes into account the target albedo, the numerical aperture of the collection optics, the frequency shifters efficiencies, the transmission of all optical components, and the overlap of the retrodiffused field with the Gaussian cavity beam (confocal feature).

The coherent interaction (beating) between the lasing electric fields and the frequency-shifted reinjected fields leads to a modulation of the laser output power at the frequency F e . For detection purposes, a small part of the laser output beam is sent to a photodiode. The delivered voltage is analyzed by a numerical oscilloscope which allows fast Fourier transform calculations, and processed by a lock-in amplifier which gives the LOFI signal (i.e., the amplitude and the phase of the laser output power modulation) at the demodulation frequency F e . Experimentally, LOFI images can be obtained pixel by pixel by a full two-dimensional galvanometric scanning. In this study, the scanning device is not used because all the measures are realized on a single target point. One can notice that in contrast to a conventional heterodyne interferometer, the LOFI setup shown here does not require complex alignment. Indeed, the LOFI interferometer is always self-aligned because the laser simultaneously fulfills the functions of the source (i.e., photon emitter) and of the photodetector (i.e., photon receptor).

B. Experimental observations of the relaxation frequency shift

Figure 5 shows how the experimental noise power spectrum of the laser is modified when the amount of optical feedback (R e ) increases and/or when the modulation frequency approaches the intrinsic laser relaxation frequency (F e ≈ F R ).

In agreement with Eqs. (24a)-(24c), Fig. 5 shows that the noise power spectrum gain is composed of two resonance frequencies symmetrically located on both sides of the frequency shift: (F e ± R ), and that the amplitude of the two maxima are of the same order of magnitude for strong optical feedback (i.e., when R 2 S,NL -20 dB in arbitrary units). Inversely, the right maximum disappears for weak optical feedback (i.e., when R 2 S,NL -20 dB) leading to the conventional linear gain with only one resonance frequency (F R ).

For a given modulation frequency, the left column of Fig. 5 shows how FR (i.e., here the left sideband) is shifted down to lower frequencies when the amount of optical feedback increases.

For a given amount of optical feedback, the right column of Fig. 5 shows how FR is shifted to a lower frequency when the modulation frequency decreases. Note that the numerical results of Fig. 2 and the experimental results of Fig. 5 are qualitatively similar.

More quantitatively, Fig. 6 shows the relative variation of the relaxation frequency [Fig. 6(a)], and the variation of the amplitude of laser output power modulation [Fig. 6(b)] as a function of the modulation frequency. Good agreement between the analytical and the experimental results can be observed.

First, one can observe in Figs. 6(a) and 6(b) that the amplitude of the modulation (R 2 S,NL ) at the frequency shift (F e ) and the relaxation frequency shift FR as a function of the detuning vary in opposite direction. Indeed, in agreement with Eq. ( 23), the decrease of the relaxation frequency is proportional to the increase of the square of the laser amplitude modulation. Secondly, one can observe a small hysteresis induced by the direction (up or down) of the modulation frequency variation.

In agreement with Eq. ( 23), one can understand that the hysteresis observed in the relaxation frequency shift is induced by the hysteresis of the modulation amplitude previously studied in [START_REF] Schwartz | Subharmonic hysteresis and period doubling bifurcations for periodically driven laser[END_REF]. The experimental results shown in Fig. 6 have been analytically adjusted by using a feedback reflectivity of R e ≈ 1 × 10 -8 and the following values of the laser parameters: η = 1.2,γ 1 = 1 × 10 5 s -1 ,γ c = 1 × 10 9 s -1 , corresponding to an intrinsic laser relaxation frequency of F R ≈ 700 kHz.

V. CONCLUSION AND PERSPECTIVES

In this article we have demonstrated (analytically, numerically, and experimentally) how the stationary dynamics (i.e., the output power modulation) and the transient dynamics (i.e., the laser quantum noise) of a laser subjected to frequency-shifted optical feedback are coupled through the nonlinear-laser dynamics. Both the numerical results and the experimental ones show very good agreement with the analytical predictions. More precisely, this study shows how the noise power spectrum is related to the strength of the beating between the intracavity laser electric field and the frequency-shifted optical electric field (i.e., the amplitude of the laser output power modulation). The shape of the noise power spectrum gives an image of the transfer modulation function (i.e., of the amplification gain) of the nonlinear-laser dynamics. A better understanding of the laser nonlinear behavior could lead to a new generation of laser metrology techniques with improved performances.

More specifically, this study shows how the relaxation frequency is shifted and how this shift becomes significant in the strong feedback situation where the modulation frequency is nearly resonant ( e ≈ R ) and where the amount of optical feedback is high enough [R e > R e, lim = (γ 1 η/2γ c ) 2 ] to induce nonlinear dynamical behaviors in the laser output power modulation. Under these conditions, we have also observed that the relaxation frequency shift exhibits a small hysteresis, induced by the direction (up or down) of the variation of the modulation frequency. This hysteresis phenomenon appears simultaneously in the amplitude of the laser power modulation and in the shift of the relaxation oscillations frequency (i.e., resonant peak of the noise). This correlation demonstrates again the nonlinear dynamical coupling between these two physical quantities.

From the applied optics point of view, since the relaxation frequency, and consequently the noise power spectrum of the laser, can be controlled by the strength of the optical feedback, it opens interesting new perspectives in the field of interferometric phase measurement (and more particularly vibrometry) using a LOFI setup. For example, when the amount of optical feedback increases, the noise in the vicinity of the modulation frequency seems to be pushed far away, allowing lower phase noise and consequently the possibility to detect lower amplitudes of vibration [START_REF] Lacot | Phase-sensitive laser detection by frequency-shifted optical feedback[END_REF][START_REF] Otsuka | Real-time nanometer vibration measurement with self-mixing microchip solid-state laser[END_REF]. Also, the shape of the noise power spectrum gives an image of the transfer modulation function (i.e., of the amplification gain) of the nonlinear-laser dynamics. The knowledge of this gain seems to be of importance to allow measurement of transient vibrations with subnanometric amplitude and a broad vibration spectrum of several megahertz (around the carrier frequency F e ). This work is in progress.

APPENDIX: EXISTENCE OF THREE REAL ROOTS FOR THE AMPLITUDE OF THE LASER OUTPUT POWER MODULATION

The square modulus of the stationary modulation amplitude in the nonlinear case [Eq. ( 16)] obeys the equation where a = (σ e -1)/α and b = (ηε/2)/α. This equation is geometrically understood as the crossing of a parabola and a hyperbola (Fig. 7).

For a given value of b, there exists a range a ∈ [a 1 ,a 2 ] where the equation has three real roots. To find the boundary of this range, we have to solve simultaneously Eq. (A2) and 2(x + a) = -

1 x 2 , ( A3
)
where b is given and a and x are to be determined. The exact solution is obtained by solving first for x:

1 4x 4 - 1 x + b 2 = 0, (A4)
and then use Eq. (A3) to get the value of a i . For b 1, the computation is simplified and the boundaries are approximately given by (Fig. 7)

a 1 = -3/2 2/3 , a 2 = -1/b 2 . (A5)
In terms of real parameters, the above condition is written

FIG. 1 .

 1 FIG. 1. Numerical simulation. Power spectra of the laser intensity dynamics |F[I (t)/I sat ]| 2 . (a) Weak feedback: R e = 10 -11 ,F e /F R = 1.3; peaks: (1) F R , (2) F e . (b) Strong feedback: R e = 10 -7 ,F e /F R = 1.01; peaks: (1) F e -FR , (2) FR , (3) F e , (4) 2F e -FR , (5) F e + FR , (6) 2F e , (7) 3F e -FR . Laser parameters: I sat = 3 × 10 9 photons,η = 1.2,γ c = 1 × 10 9 s -1 ,γ 1 /γ c = 1 × 10 -5 ,F R /γ c = 2.251 × 10 -4 . The dashed curve is a fit of the laser noise power spectrum without feedback.

FIG. 2 .

 2 FIG. 2. Numerical simulation. Power spectra of the laser intensity dynamics. Left column: fixed modulation frequency (F e /F R = 1.01) and an increasing amount of optical feedback (R e ); (a) R e = 10 -9 ; (b) R e = 10 -8 , (c) R e = 10 -7 . Right column: fixed amount of optical feedback (R e = 10 -7 ) and a modulation frequency going closer to the intrinsic relaxation frequency; (d) F e /F R = 1.3; (e) F e /F R = 1.1; (f) F e /F R = 1.01. Laser parameters are identical to Fig. 1. The dashed curve is a fit of the laser noise power spectrum without feedback.

Figure 2

 2 Figure 2 captures the general feature of laser noise power spectrum modification when the amount of optical feedback (R e ) increases and/or when the modulation frequency approaches the intrinsic laser relaxation frequency (F e → F R ).For a given modulation frequency, the left column of Fig.2shows how FR (i.e., the noise left sideband) is shifted to lower frequencies when the amount of optical feedback increases. For a given amount of optical feedback, the right column of Fig.2shows how FR is pushed down to a lower frequency when the modulation frequency approaches the intrinsic laser relaxation frequency.Also, one can observe that in the low feedback condition the noise power spectrum exhibits a single resonance [Figs.1(a) and 2(d)], while in the in strong feedback condition it exhibits a double resonance [see Figs. 2(a)-2(c) and 2(e)-2(f)] due to

  Figure 2 captures the general feature of laser noise power spectrum modification when the amount of optical feedback (R e ) increases and/or when the modulation frequency approaches the intrinsic laser relaxation frequency (F e → F R ).For a given modulation frequency, the left column of Fig.2shows how FR (i.e., the noise left sideband) is shifted to lower frequencies when the amount of optical feedback increases. For a given amount of optical feedback, the right column of Fig.2shows how FR is pushed down to a lower frequency when the modulation frequency approaches the intrinsic laser relaxation frequency.Also, one can observe that in the low feedback condition the noise power spectrum exhibits a single resonance [Figs.1(a) and 2(d)], while in the in strong feedback condition it exhibits a double resonance [see Figs. 2(a)-2(c) and 2(e)-2(f)] due to

FIG. 3 .

 3 FIG. 3. (a) Relative variation of the relaxation frequency and (b) square of the amplitude of the laser output power modulation. The lines correspond to results obtained from analytical equations [Eq. (16)], while the dots correspond to results obtained from numerical simulations. Solid lines and square dots (R e = 1 × 10 -7 ); dotted lines and circular dots (R e = 1 × 10 -8 ). Laser parameters are identical to Fig. 1.

FIG. 4 .

 4 FIG. 4. Schematic diagram of the LOFI setup. L 1 , L 2 , and L 3 : lenses; BS: beam splitter; GS: galvanometric scanner; FS: frequency shifter with a round-trip frequency shift F e ; PD: photodiode.

FIG. 5 .

 5 FIG. 5. Experimental power spectra of the laser intensity dynamics. Left column: fixed modulation frequency (F e /F R = 1.075) and an increasing amount of optical feedback R e characterized by an increase of the laser output power modulation; (a) R 2 S,NL = -23.1 dB; (b) R 2 S,NL = -14.3 dB; (c) R 2 S,NL = -9.27 dB. Right column: fixed amount of optical feedback and a decreasing modulation frequency around the intrinsic relaxation frequency; (d) F e /F R = 1.075,R 2 S,NL = -24.72 dB; (e) F e /F R ≈ 1.0,R 2 S,NL = -19.18 dB; (f) F e /F R = 0.925,R 2 S,NL = -12.15 dB.

FIG. 6 .

 6 FIG. 6. (a) Relative variation of the relaxation frequency and (b) square of the amplitude of the laser output power modulation versus the normalized scanning of the modulation frequency. The lines correspond to results obtained from analytical equations, while the dots correspond to results obtained from experimental measurement. Dashed lines and circular dots: the scanning of the modulation frequency increases. Solid lines and square dots: the scanning of the modulation frequency decreases. Estimated laser parameters: η = 1.2,γ 1 = 1 × 10 5 s -1 ,γ c = 1 × 10 9 s -1 ,F R ≈ 700 kHz. Estimated feedback parameter: R e = 1 × 10 -8 .

4 √ 6 FIG. 7 .

 467 FIG. 7. Log-log plot of the hyperbola g(x) = 1/x (solid line) and three parabola f (x) = (x + a) 2 + b 2 (dashed, dotted, and dashdotted lines) for a fixed value of b = 0.23 and three values of a: (a 1 = -3/2 2/3 ≈ -1.89, a = -5, a 2 ≈ -1/b 2 = -18.85).

  ) is identical to Eq. (17).