Denys Dutykh

Didier Clamond
email: diderc@unice.fr

AND Paul Milewski
email: p.a.milewski@bath.ac.uk

Dimitrios Mitsotakis

FINITE VOLUME AND PSEUDO-SPECTRAL SCHEMES FOR THE FULLY NONLINEAR 1D SERRE EQUATIONS

Keywords: Serre equations, finite volumes, UNO scheme, IMEX scheme, spectral methods, Euler equations, free surface flows 2010 Mathematics Subject Classification. 76B15 (primary), 76B25, 65M08 (secondary)

After we derive the Serre system of equations of water wave theory from a generalized variational principle, we present some of its structural properties. We also propose a robust and accurate finite volume scheme to solve these equations in one horizontal dimension. The numerical discretization is validated by comparisons with analytical, experimental data or other numerical solutions obtained by a highly accurate pseudo-spectral method.

Introduction

The full water wave problem consisting of the Euler equations with a free surface still is a very difficult to study theoretically and even numerically. Consequently, water wave theory has always been developed through the derivation, analysis and comprehension of various approximate models (see the historical review of Craik [24] for more information). For this reason, a plethora of approximate models have been derived under various physical assumptions. In this family, the Serre equations have a particular place and they are the subject of the present study. The Serre equations can be derived from the Euler equations, contrary to Boussinesq systems or the shallow water system, without the small amplitude or the hydrostatic assumptions respectively.

The Serre equations are named after François Serre, an engineer at École Nationale des Ponts et Chaussées, who derived this model for the first time in 1953 in his prominent paper entitled "Contribution à l'étude des écoulements permanents et variables dans les canaux" (see [59]). Later, these equations were independently rediscovered by Su and Gardner [START_REF] Su | Korteweg-de Vries equation and generalizations. III. Derivation of the Korteweg-de Vries equation and Burgers equation[END_REF] and by Green, Laws and Naghdi [START_REF] Green | On the theory of water waves[END_REF]. The extension of Serre equations for general uneven bathymetries was derived by Seabra-Santos et al. [START_REF] Seabra-Santos | Numerical and Experimental study of the transformation of a Solitary Wave over a Shelf or Isolated Obstacle[END_REF]. In the Soviet literature these equations were known as the Zheleznyak-Pelinovsky model [START_REF] Zheleznyak | Physical and mathematical models of the tsunami climbing a beach[END_REF]. For some generalizations and new results we refer to recent studies by Barthélémy [START_REF] Barthélémy | Nonlinear shallow water theories for coastal waves[END_REF], Dias & Milewski [START_REF] Barth | Aspects of unstructured grids and finite-volume solvers for the Euler and Navier-Stokes equations[END_REF] and Carter & Cienfuegos [12].

A variety of numerical methods have been applied to discretize dispersive wave models and, more specifically, the Serre equations. A pseudo-spectral method was applied in [START_REF] Barth | Aspects of unstructured grids and finite-volume solvers for the Euler and Navier-Stokes equations[END_REF], an implicit finite difference scheme in [START_REF] Mirie | Collision between two solitary waves. Part 2. A numerical study[END_REF][START_REF] Barthélémy | Nonlinear shallow water theories for coastal waves[END_REF] and a compact higher-order scheme in [16,17]. Some Galerkin and Finite Element type methods have been successfully applied to Boussinesq-type equations [START_REF] Dougalis | Theory and numerical analysis of Boussinesq systems: A review[END_REF]54,4,[START_REF] Antonopoulos | Numerical solution of Boussinesq systems of the Bona-Smith family[END_REF]. A finite difference discretization based on an integral formulation was proposed by Bona & Chen [START_REF] Bona | A Boussinesq system for two-way propagation of nonlinear dispersive waves[END_REF].

Recently, efficient high-order explicit or implicit-explicit finite volume schemes for dispersive wave equations have been developed [START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF][START_REF] Dutykh | Finite volume methods for unidirectional dispersive wave models[END_REF][START_REF] Dutykh | Finite volume methods for unidirectional dispersive wave models[END_REF]. The robustness of the proposed numerical schemes also allowed simulating the run-up of long waves on a beach with high accuracy [START_REF] Dutykh | Finite volume methods for unidirectional dispersive wave models[END_REF]. The present study is a further extension of the finite volume method to the practically important case of the Serre equations. We develop also a pseudo-spectral Fourier-type method to validate the proposed finite volume scheme. In all cases where the spectral method is applicable, it outperforms the finite volumes. However, the former is applicable only to smooth solutions in periodic domains, while the area of applicability of the latter is much broader including dispersive shocks (or undular bores) [START_REF] El | Unsteady undular bores in fully nonlinear shallowwater theory[END_REF], non-periodic domains, etc.

The present paper is organized as follows. In Section 2 we provide a derivation of the Serre equations from a relaxed Lagrangian principle and discuss some structural properties of the governing equations. The rationale on the employed finite volume scheme are given in Section 3. A very accurate pseudo-spectral method for the numerical solution of the Serre equations is presented in Section 4. In Section 5, we present convergence tests and numerical experiments validating the model and the numerical schemes. Finally, Section 6 contains the main conclusions.

Mathematical model

Consider an ideal incompressible fluid of constant density ρ. The vertical projection of the fluid domain Ω is a subset of R 2 . The horizontal independent variables are denoted by x = (x 1 , x 2 ) and the upward vertical one by y. The origin of the Cartesian coordinate system is chosen such that the surface y = 0 corresponds to the still water level. The fluid is bounded below by an impermeable bottom at y = -d and above by the free surface located at y = η(x, t). We assume that the total depth h(x, t) ≡ d + η(x, t) remains positive h(x, t) h 0 > 0 at all times t. The sketch of the physical domain is shown in Figure 1.

Remark 1. We make the classical assumption that the free surface is a graph y = η(x, t) of a single-valued function. This means in practice that we exclude some interesting phenomena, (e.g., wave breaking) which are out of the scope of this modeling paradigm.

Assuming that the flow is incompressible and irrotational, the governing equations of the classical water wave problem are the following [START_REF] Lamb | Hydrodynamics[END_REF][START_REF] Stoker | Water waves, the mathematical theory with applications[END_REF][START_REF] Mei | The applied dynamics of ocean surface waves[END_REF]71]

∇ 2 φ + ∂ 2 y φ = 0 -d(x, t) y η(x, t), (2.1)
∂ t η + (∇φ) • (∇η) -∂ y φ = 0 y = η(x, t), (2.2)
∂ t φ + 1 2 |∇φ| 2 + 1 2 (∂ y φ) 2 + g η = 0 y = η(x, t), (2.3)
d t + (∇d) • (∇φ) + ∂ y φ = 0 y = -d(x, t), (2.4)
with φ being the velocity potential (by definition, the irrotational velocity field (u, v) = (∇φ, ∂ y φ), g the acceleration due to the gravity force and ∇ = (∂ x 1 , ∂ x 2 ) denotes the gradient operator in horizontal Cartesian coordinates and |∇φ| 2 ≡ (∇φ) • (∇φ).

The incompressibility condition leads to the Laplace equation for φ. The main difficulty of the water wave problem lies on the nonlinear free surface boundary conditions and that the free surface shape is unknown. Equations (2.2) and (2.4) express the free-surface kinematic condition and bottom impermeability respectively, while the dynamic condition (2.3) expresses the free surface isobarity.

The water wave problem possesses several variational structures [55, [START_REF] Whitham | A general approach to linear and non-linear dispersive waves using a Lagrangian[END_REF][START_REF] Luke | A variational principle for a fluid with a free surface[END_REF][START_REF] Zakharov | Stability of periodic waves of finite amplitude on the surface of a deep fluid[END_REF][START_REF] Broer | On the Hamiltonian theory of surface waves[END_REF]. In the present study, we will focus mainly on the Lagrangian variational formalism but not exclusively. The surface gravity wave equations (2.1)-(2.4) can be derived by minimizing the following functional proposed by Luke [START_REF] Luke | A variational principle for a fluid with a free surface[END_REF]:

L = ˆt2 t 1 ˆΩ L ρ d 2 x dt, L = - ˆη -d g y + ∂ t φ + 1 2 (∇ φ) 2 + 1 2 (∂ y φ) 2 dy. (2.5)
In a recent study, Clamond and Dutykh [START_REF] Clamond | Practical use of variational principles for modeling water waves[END_REF] proposed using Luke's Lagrangian (2.5) in the following relaxed form

L = (η t + μ • ∇η -ν) φ + (d t + μ • ∇d + ν) φ -1 2 g η 2 + ˆη -d µ • u -1 2 u 2 + νv -1 2 v 2 + (∇ • µ + ν y ) φ dy, (2.6) x y O d h(x, t) η(x, t)
Figure 1. Sketch of the physical domain.

where {u, v, µ, ν} are the horizontal, vertical velocities and associated Lagrange multipliers, respectively. The additional variables {µ, ν} (Lagrange multipliers) are called pseudovelocities. The 'tildes' and 'wedges' denote, respectively, a quantity computed at the free surface y = η(x, t) and at the bottom y = -d(x, t). We shall also denote below with 'bars' the quantities averaged over the water depth. While the original Lagrangian (2.5) incorporates only two variables (η and φ), the relaxed Lagrangian density (2.6) involves six variables {η, φ, u, v, µ, ν}. These additional degrees of freedom provide us with more flexibility in constructing various approximations. For more details, explanations and examples we refer to [START_REF] Clamond | Practical use of variational principles for modeling water waves[END_REF].

Derivation of the Serre equations

Now, we illustrate the practical use of the variational principle (2.6) on an example borrowed from [START_REF] Clamond | Practical use of variational principles for modeling water waves[END_REF]. First of all, we choose a simple shallow water ansatz, which is a zerothorder polynomial in y for φ and for u, and a first-order one for v, i.e., we approximate flows that are nearly uniform along the vertical direction

φ ≈ φ(x, t), u ≈ ū(x, t), v ≈ (y + d) (η + d) -1 ṽ(x, t).
(2.7)

We have also to introduce suitable ansatz for the Lagrange multiplier µ and ν

µ ≈ μ(x, t), ν ≈ (y + d) (η + d) -1 ν(x, t).
In the remainder of this paper, we will assume for simplicity the bottom to be flat d(x, t) = d = Cst (the application of this method to uneven bottoms can be found in [START_REF] Dutykh | Shallow water equations for large bathymetry variations[END_REF]31], for example). With this ansatz the Lagrangian density (2.6) becomes

L = (η t + μ • ∇η) φ -1 2 g η 2 + (η + d) μ • ū -1 2 ū2 + 1 3 ν ṽ -1 6 ṽ2 + φ ∇ • μ . (2.8)
Finally, we impose a constraint of the free surface impermeability, i.e. ν = η t + μ • ∇η.

After substituting the last relation into the Lagrangian density (2.8), the Euler-Lagrange equations and some algebra lead to the following equations:

h t + ∇ • [ h ū ] = 0, (2.9) ūt + 1 2 ∇| ū| 2 + g ∇h + 1 3 h -1 ∇[ h 2 γ ] = ( ū • ∇h) ∇(h∇ • ū) -[ ū • ∇(h∇ • ū) ] ∇h, (2.10)
where we eliminated φ, μ and ṽ and where

γ ≡ ṽt + ū • ∇ṽ = h (∇ • ū) 2 -∇ • ūt -ū • ∇ [ ∇ • ū ] , (2.11)
is the fluid vertical acceleration at the free surface. The vertical velocity at the free surface ṽ can be expressed in terms of other variables as well, i.e., ṽ =

η t + (∇ φ) • (∇η) 1 + 1 3 |∇η| 2 .
In two dimensions (one horizontal dimension), the sum of two terms on the right hand side of (2.10) vanish and the system (2.9), (2.10) reduces to the classical Serre equations [59].

Remark 2. In [START_REF] Clamond | Practical use of variational principles for modeling water waves[END_REF] it is explained why equations (2.9), (2.10) cannot be obtained from the classical Luke's Lagrangian. One of the main reasons is that the horizontal velocity ū does not derive from the potential φ using a simple gradient operation. Thus, a relaxed form of the Lagrangian density (2.6) is necessary for the variational derivation of the Serre equations (2.9), (2.10) (see also [START_REF] Kim | A derivation of the Green-Naghdi equations for irrotational flows[END_REF] & [START_REF] Miles | Weakly dispersive nonlinear gravity waves[END_REF]).

Remark 3. In some applications in coastal engineering it is required to estimate the loading exerted by water waves onto vertical structures [START_REF] Clauss | The New Year Wave in a sea keeping basin: Generation, propagation, kinematics and dynamics[END_REF]. The pressure can be computed in the framework of the Serre equations as well. For the first time these quantities were computed in the pioneering paper by M. [START_REF] Zheleznyak | Influence of long waves on vertical obstacles[END_REF] [74]. Here for simplicity we provide the expressions in two space dimensions which were derived in [74]. The pressure distribution inside the fluid column being given by

P(x, y, t) ρgd = η -y d + 1 2 h d 2 -1 + y d 2 γ d g h ,
one can compute the force F exerted on a vertical wall:

F(x, t) ρgd 2 = ˆη -d P ρgd 2 dy = 1 2 + γ 3 g h d 2 .
Finally, the tilting moment M relative to the sea bed is given by the following formula:

M(x, t) ρgd 3 = ˆη -d P ρgd 3 (y + d) dy = 1 6 + γ 8 g h d 3 .

Generalized Serre equations

A further generalization of the Serre equations can be obtained if we modify slightly the shallow water ansatz (2.7) following again the ideas from [START_REF] Clamond | Practical use of variational principles for modeling water waves[END_REF]:

φ ≈ φ(x, t), u ≈ ū(x, t), v ≈ y + d η + d λ ṽ(x, t).
In the following we consider for simplicity the two-dimensional case and put µ = u and ν = v together with the constraint ṽ = η t + ũη x (free surface impermeability). Thus, the Lagrangian density (2.6) becomes

L = ( h t + [ h ū ] x ) φ -1 2 g η 2 + 1 2 h ū2 + 1 2 β h ( η t + ū η x ) 2 , (2.12)
where β = (2λ + 1) -1 . After some algebra, the Euler-Lagrange equations lead to the following equations

h t + [ h ū ] x = 0, (2.13) ūt + ū ūx + g h x + β h -1 [ h 2 γ ] x = 0, (2.14)
where γ is defined as above (2.11). If β = 1 3 (or, equivalently, λ = 1) the classical Serre equations (2.9), (2.10) are recovered.

Using equations (2.13) and (2.14) one can show that the following relations hold

[ h ū ] t + h ū2 + 1 2 g h 2 + β h 2 γ x = 0, [ ū -β h -1 (h 3 ūx ) x ] t + 1 2 ū2 + g h -1 2 h 2 ū2 x -β ūh -1 (h 3 ūx ) x x = 0, [ h ū -β (h 3 ūx ) x ] t + h ū2 + 1 2 g h 2 -2 β h 3 ū2
x -β h 3 ū ūxx -h 2 h x ū ūx x = 0, (2.15)

1 2 h ū2 + 1 2 β h 3 ū2 x + 1 2 g h 2 t + 1 2 ū2 + 1 2 β h 2 ū2 x + g h + β h γ h ū x = 0.
Physically, these relations represent conservations of the momentum, quantity q = ūβ h -1 (h 3 ūx ) x , its flux q := h ū -β (h 3 ūx ) x and the total energy, respectively. Moreover, the Serre equations are invariant under the Galilean transformation. This property is naturally inherited from the full water wave problem since our ansatz does not destroy this symmetry [START_REF] Benjamin | Hamiltonian structure, symmetries and conservation laws for water waves[END_REF] and the derivation is made according to variational principles.

Equations (2.13)-(2.14) admit a (2π/k)-periodic cnoidal traveling wave solution

ū = c η d + η , (2.16) η = a dn 2 1 2 κ(x -ct)|m -E/K 1 -E/K = a -H sn 2 1 2 κ(x -ct)|m , (2.17)
where dn and sn are the Jacobian elliptic functions with parameter m (0 m 1), and where K = K(m) and E = E(m) are the complete elliptic integrals of the first and second kind, respectively [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF]. The wave parameters are given by the relations

k = π κ 2 K , H = m a K K -E , (κd) 2 = g H m β c 2 , (2.18) m = g H (d + a) (d + a -H) g (d + a) 2 (d + a -H) -d 2 c 2 .
(2.19) However, in the present study, we are interested in the classical solitary wave solution which is recovered in the limiting case m → 1

η = a sech 2 1 2 κ(x -ct), ū = c η d + η , c 2 = g(d + a), (κd) 2 = a β(d + a)
.

(2.20)

For illustrative purposes, a solitary wave along with a cnoidal wave of the same amplitude a = 0.05 are depicted in Figure 2. Using the exact solitary wave solution (2.20) we can assess the accuracy of the Serre equations (with β = 1

3 ) by making comparisons with corresponding solutions to the original full Euler equations. The procedure we use to construct traveling wave solutions to Euler equations is described in [18]. The Matlab script used to generate these profiles (up to machine precision) can be freely downloaded from the File Exchange server [19]. The results of comparison for several values of the speed parameter c are presented on Figure 3. We can see that solitary waves to the Serre equations approximate fairly well the full Euler solutions approximately up to the amplitude a/d = 1 2 . We note that similar conclusions were obtained in a previous study by [START_REF] Li | A Numerical Study of the Exact Evolution Equations for Surface Waves in Water of Finite Depth[END_REF] [START_REF] Li | A Numerical Study of the Exact Evolution Equations for Surface Waves in Water of Finite Depth[END_REF].

Invariants of the Serre equations

Henceforth we consider only the two-dimensional case. As pointed out by Y. [START_REF] Li | Hamiltonian structure and linear stability of solitary waves of the Green-Naghdi equations[END_REF] [45], the classical Serre equations possess a non-canonical Hamiltonian structure which can be easily generalized for the model (2.13), (2.14) where the Hamiltonian functional H and the symplectic operator J are defined as

h t qt = J • δ H / δ q δ H / δh ,
H = 1 2 ˆR h ū2 + β h 3 ū2 x + g η 2 dx, J = - h x 0 qx + q∂ x h∂ x .
The variable q is defined by

q ≡ h ū -β [ h 3 ūx ] x .
The conservation of the quantity q was established in equation (2.15).

According to [45], one-parameter symmetry groups of Serre's equations include the space translation (x + ε, t, h, u), the time translation (x, t + ε, h, u), the Galilean boost (x + εt, t, h, u + ε) and the scaling e ε (e ε x, t, e ε h, u). Using the first three symmetry groups and the symplectic operator J, one may recover the following invariants:

Q = ˆR η q d + η dx, H, ˆR [ t q -x η ] dx. (2.21)
Obviously, the equation (2.13) leads to an invariant closely related to the mass conservation property ´R η dx. The scaling does not yield any conserved quantity with respect to the symplectic operator J. Below, we are going to use extensively the generalized energy H and the generalized momentum Q conservation to assess the accuracy of the numerical schemes in addition to the exact analytical solution (2.20).

Finite volume scheme and numerical results

In the present study, we propose a finite volume discretization procedure [5,[START_REF] Barth | Finite volume methods: foundation and analysis[END_REF] for the Serre equations (2.13), (2.14) that we rewrite here as

h t + [ h u ] x = 0, (3.1)
u t + 1 2 u 2 + g h x = β h -1 h 3 (u xt + u u xx -u 2 x ) x , (3.2)
where the over-bars have been omitted for brevity. (In this section, over-bars denote quantities averaged over a cell, as explained below.) We begin our presentation by a discretization of the hyperbolic part of the equations (which are simply the classical Saint-Venant equations) and then discuss the treatment of dispersive terms. The Serre equations can be formally put under the quasilinear form

V t + [ F (V ) ] x = S(V ), (3.3)
where V , F (V ) are the conservative variables and the advective flux function, respectively

V ≡ h u , F (V ) ≡ h u 1 2 u 2 + g h .
The source term S(V ) denotes the right-hand side of (3.1), (3.2) and thus, depends also on space and time derivatives of V . The Jacobian of the advective flux F (V ) can be easily computed

A(V ) = ∂ F (V ) ∂V = u h g u .
The Jacobian A(V ) has two distinctive eigenvalues

λ ± = u ± c s , c s ≡ gh.
The corresponding right and left eigenvectors are provided here

R = h -h c s c s , L = R -1 = 1 2 h -1 c -1 s -h -1 c -1 s .
We consider a partition of the real line R into cells (or finite volumes)

C i = [x i-1 2 , x i+ 1 2 ] with cell centers x i = 1 2 (x i-1 2 + x i+ 1 2 ) (i ∈ Z).
Let ∆x i denotes the length of the cell C i . In the sequel we will consider only uniform partitions with ∆x i = ∆x, ∀i ∈ Z. We would like to approximate the solution V (x, t) by discrete values. In order to do so, we introduce the cell average of V on the cell C i (denoted with an overbar), i.e.,

V i (t) ≡ hi (t) , ūi (t) = 1 ∆x ˆCi V (x, t) dx.
A simple integration of (3.3) over the cell C i leads the following exact relation:

d V dt + 1 ∆x F (V (x i+ 1 2 , t)) -F (V (x i-1 2 , t)) = 1 ∆x ˆCi S(V ) dx ≡ Si .
Since the discrete solution is discontinuous at cell interfaces x i+ 1 2 (i ∈ Z), we replace the flux at the cell faces by the so-called numerical flux function

F (V (x i± 1 2 , t)) ≈ F i± 1 2 ( V L i± 1 2 , V R i± 1 2 ), where V L,R i± 1 2
denotes the reconstructions of the conservative variables V from left and right sides of each cell interface (the reconstruction procedure employed in the present study will be described below). Consequently, the semi-discrete scheme takes the form

d V i dt + 1 ∆x F i+ 1 2 -F i-1 2 = Si . (3.4)
In order to discretize the advective flux F (V ), we use the FVCF scheme [START_REF] Ghidaglia | Une méthode volumes-finis à flux caractéristiques pour la résolution numérique des systèmes hyperboliques de lois de conservation[END_REF][START_REF] Ghidaglia | On the numerical solution to two fluid models via cell centered finite volume method[END_REF]:

F(V , W ) = F (V ) + F (W ) 2 -U(V , W ) • F (W ) -F (V ) 2 .
The first part of the numerical flux is centered, the second part is the upwinding introduced through the Jacobian sign-matrix U(V , W ) defined as

U(V , W ) = sign A 1 2 (V + W ) , sign(A) = R • diag(s + , s -) • L,
where s ± ≡ sign(λ ± ). After some simple algebraic computations, one can find

U = 1 2 s + + s - (h/c s ) (s + -s -) (g/c s ) (s + -s -) s + + s - ,
the sign-matrix U being evaluated at the average state of left and right values.

High order reconstruction

In order to obtain a higher-order scheme in space, we need to replace the piecewise constant data by a piecewise polynomial representation. This goal is achieved by various so-called reconstruction procedures such as MUSCL TVD [43, [START_REF] Van Leer | Towards the ultimate conservative difference scheme V: a second order sequel to Godunov' method[END_REF][START_REF] Van Leer | Upwind and High-Resolution Methods for Compressible Flow: From Donor Cell to Residual-Distribution Schemes[END_REF], UNO [START_REF] Harten | Uniformly high-order accurate nonscillatory schemes[END_REF], ENO [39], WENO [START_REF] Xing | High order finite difference WENO schemes with the exact conservation property for the shallow water equations[END_REF] and many others. In our previous study on Boussinesq-type equations [START_REF] Dutykh | Finite volume schemes for dispersive wave propagation and runup[END_REF], the UNO2 scheme showed a good performance with small dissipation in realistic propagation and run-up simulations. Consequently, we retain this scheme for the discretization of the advective flux in Serre equations.

Remark 4. In TVD schemes, the numerical operator is required (by definition) not to increase the total variation of the numerical solution at each time step. It follows that the value of an isolated maximum may only decrease in time which is not a good property for the simulation of coherent structures such as solitary waves. The non-oscillatory UNO2 scheme, employed in our study, is only required to diminish the number of local extrema in the numerical solution. Unlike TVD schemes, UNO schemes are not constrained to damp the values of each local extremum at every time step.

The main idea of the UNO2 scheme is to construct a non-oscillatory piecewise-parabolic interpolant Q(x) to a piecewise smooth function V (x) (see [START_REF] Harten | Uniformly high-order accurate nonscillatory schemes[END_REF] for more details). On each segment containing the face x i+ 1 2 ∈ [x i , x i+1 ], the function Q(x) = q i+ 1 2 (x) is locally a quadratic polynomial and wherever v(x) is smooth we have

Q(x) -V (x) = 0 + O(∆x 3 ), d Q dx (x ± 0) - d V dx = 0 + O(∆x 2 ).
Also, Q(x) should be non-oscillatory in the sense that the number of its local extrema does not exceed that of V (x). Since q i+ 1 2

(x i ) = V i and q i+ 1 2 (x i+1 ) = V i+1
, it can be written in the form

q i+ 1 2 (x) = V i + d i+ 1 2 {V } × x -x i ∆x + 1 2 D i+ 1 2 {V } × (x -x i )(x -x i+1 ) ∆x 2 ,
where

d i+ 1 2 {V } ≡ V i+1 -V i and D i+ 1 2
V is closely related to the second derivative of the interpolant since D i+ 1 2 {V } = ∆x 2 q ′′ i+ 1 2 (x). The polynomial q i+ 1 2 (x) is chosen to be the least oscillatory between two candidates interpolating V (x) at (x i-1 , x i , x i+1 ) and (x i , x i+1 , x i+2 ). This requirement leads to the following choice of

D i+ 1 2 {V } ≡ minmod (D i {V }, D i+1 {V }) with D i {V } = V i+1 -2 V i + V i-1 , D i+1 {V } = V i+2 -2 V i+1 + V i ,
and where minmod(x, y) is the usual minmod function defined as minmod(x, y) ≡ 1 2 [ sign(x) + sign(y) ] × min(|x|, |y|). To achieve the second order O(∆x 2 ) accuracy, it is sufficient to consider piecewise linear reconstructions in each cell. Let L(x) denote this approximately reconstructed function which can be written in this form

L(x) = V i + S i × x -x i ∆x , x ∈ [x i-1 2 , x i+ 1 2 ].
In order to L(x) be a non-oscillatory approximation, we use the parabolic interpolation Q(x) constructed below to estimate the slopes S i within each cell

S i = ∆x × minmod d Q dx (x i -0), d Q dx (x i + 0) .
In other words, the solution is reconstructed on the cells while the solution gradient is estimated on the dual mesh as it is often performed in more modern schemes [5,[START_REF] Barth | Finite volume methods: foundation and analysis[END_REF]. A brief summary of the UNO2 reconstruction can be also found in [START_REF] Dutykh | Finite volume schemes for dispersive wave propagation and runup[END_REF][START_REF] Dutykh | Finite volume methods for unidirectional dispersive wave models[END_REF].

Treatment of the dispersive terms

In this section, we explain how we treat the dispersive terms of Serre equations (3.1), (3.2). We begin the exposition by discussing the space discretization and then, we propose a way to remove the intrinsic stiffness of the dispersion by partial implicitation.

For the sake of simplicity, we split the dispersive terms into three parts:

M(V ) ≡ β h -1 h 3 u xt x , D 1 (V ) ≡ β h -1 h 3 u u xx x , D 2 (V ) ≡ β h -1 h 3 u 2
x x . We propose the following approximations in space (which are all of the second order O(∆x 2 ) to be consistent with UNO2 advective flux discretization presented above)

M i ( V ) = β h -1 i h 3 i+1 (ū xt ) i+1 - h 3 i-1 (ū xt ) i-1 2 ∆x = β h -1 i 2 ∆x h 3 i+1 (ū t ) i+2 -(ū t ) i 2 ∆x - h 3 i-1 (ū t ) i -(ū t ) i-2 2 ∆x = β h -1 i 4 ∆x 2 h 3 i+1 (ū t ) i+2 -( h 3 i+1 + h 3 i-1 ) (ū t ) i + h 3 i-1 (ū t ) i-2 .
The last relation can be rewritten in a short-hand form if we introduce the matrix M( V ) such that the i-th component of the product M( V )• V t gives exactly the expression M i ( V ).

In a similar way, we discretize the other dispersive terms without giving here the intermediate steps

D 1i ( V ) = β h-1 i 2 ∆x 3 h 3 i+1 ūi+1 (ū i+2 -2ū i+1 + ūi ) - h 3 i-1 ūi-1 (ū i -2ū i-1 + ūi-2 ) , D 2i ( V ) = β h -1 i 8∆x 3 h 3 i+1 (ū i+2 -ūi ) 2 - h 3 i-1 (ū i -ūi-2 ) 2 .
In a more general nonperiodic case asymmetric finite differences should be used near the boundaries. If we denote by I the identity matrix, we can rewrite the semi-discrete scheme (3.4) by expanding the right-hand side S i

d h dt + 1 ∆x F (1)
+ ( V ) -F (1)
-( V ) = 0, (3.5) (I -M) • d ū dt + 1 ∆x F (2) + ( V ) -F (2)
-( V ) = D( V ) • ū, (3.6)
where F

(1,2) ± ( V ) are the two components of the advective numerical flux vector F at the right (+) and left (-) faces correspondingly and

D( V ) ≡ D 1 ( V ) -D 2 ( V ).
Finally, in order to obtain the semidiscrete scheme, one has to solve a linear system to find explicitly the time derivative dū/dt. A mathematical study of the resulting matrix I -M is not straightforward to perform. However, in our numerical tests we have never experienced any difficulties to invert it.

Temporal scheme

We rewrite the inverted semi-discrete scheme (3.5)-(3.6) as a system of ODEs:

∂ t w = L(w, t), w(0) = w 0 .
In order to solve numerically the last system of equations, we apply the Bogacki-Shampine method [9]. It is a third-order Runge-Kutta scheme with four stages. It has an embedded second-order method which is used to estimate the local error and, thus, to adapt the time step size. Moreover, the Bogacki-Shampine method enjoys the First Same As Last (FSAL) property so that it needs three function evaluations per step. This method is also implemented in the ode23 function in Matlab [START_REF] Shampine | The MATLAB ODE Suite[END_REF]. A step of the Bogacki-Shampine method is given by

k 1 = L(w (n) , t n ), k 2 = L(w (n) + 1 2 ∆t n k 1 , t n + 1 2 ∆t), k 3 = L(w (n) ) + 3 4 ∆t n k 2 , t n + 3 4 ∆t), w (n+1) = w (n) + ∆t n × 2 9 k 1 + 1 3 k 2 + 4 9 k 3 , k 4 = L(w (n+1) , t n + ∆t n ), w (n+1) 2 = w (n) + ∆t n × 4 24 k 1 + 1 4 k 2 + 1 3 k 3 + 1 8 k 4 .
Here w (n) ≈ w(t n ), ∆t is the time step and w

(n+1) 2
is a second order approximation to the solution w(t n+1 ), so the difference between w (n+1) and w (n+1) 2

gives an estimation of the local error. The FSAL property consists in the fact that k 4 is equal to k 1 in the next time step, thus saving one function evaluation.

If the new time step ∆t n+1 is given by ∆t n+1 = ρ n ∆t n , then according to H211b digital filter approach [START_REF] Söderlind | Digital filters in adaptive time-stepping[END_REF]62], the proportionality factor ρ n is given by:

ρ n = δ ε n β 1 δ ε n-1 β 2 ρ -α n-1 , (3.7)
where ε n is a local error estimation at time step t n , δ is the desired tolerance and the constants β 1 , β 2 and α are defined as

α = 1 4 , β 1 = β 2 = 1 4 p .
The parameter p is the order of the scheme (p = 3 in our case).

Remark 5. The adaptive strategy (3.7) can be further improved if we smooth the factor ρ n before computing the next time step ∆t n+1

∆t n+1 = ρn ∆t n , ρn = ω(ρ n ).
The function ω(ρ) is called the time step limiter and should be smooth, monotonically increasing and should satisfy the following conditions

ω(0) < 1, ω(+∞) > 1, ω(1) = 1, ω ′ (1) = 1.
One possible choice is suggested in [62]:

ω(ρ) = 1 + κ arctan ρ -1 κ .
In our computations the parameter κ is set to 1.

Pseudo-spectral Fourier-type method for the Serre equations

In this Section we describe a pseudo-spectral solver to integrate numerically the Serre equations in periodic domains. In spectral methods, it is more convenient to take as variables the free surface elevation η(x, t) and the conserved quantity q(x, t)

η t + [ (d + η) ū ] x = 0, (4.1)
q t + q u -1 2 ū2 + g η -1 2 (d + η) 2 ū2 x x = 0, (4.2) q -ū + 1 3 (d + η) 2 ūxx + (d + η)η x ūx = 0. (4.3)
The first two equations (4.1), (4.2) are of evolution type, while the third one (4.3) relates the conserved variable q to the primitive variables: the free surface elevation η and the velocity ū. In order to solve relation (4.3) with respect to the velocity ū, we extract the linear part as

ū -1 3 d 2 ūxx -q = 1 3 (2dη + η 2 ) ūxx + (d + η) η x ūx N (η,ū)
.

Then, we apply to the last relation the following fixed point type iteration in Fourier space

ûj+1 = q 1 + 1 3 (kd) 2 + F {N(η, ūj )} 1 + 1 3 (kd) 2 j = 0, 1, 2, • • • , (4.4)
where ψ ≡ F{ψ} denotes the Fourier transform of a quantity ψ. The last iteration is repeated until the desired convergence. For example, for moderate amplitude solitary waves (≈ 0.2), the accuracy 10 -16 is attained in approximatively 20 iterations if the velocity ū0 is initialized from the previous time step. We note that the usual 3/2 rule is applied to the nonlinear terms for anti-aliasing [65, [START_REF] Clamond | A fast method for fully nonlinear water-wave computations[END_REF]35].

Remark 6. One can improve the fixed point iteration (4.4) by employing the so-called relaxation approach [START_REF] Isaacson | Analysis of Numerical Methods[END_REF]. The relaxed scheme takes the following form

ûj+1 = q 1 + 1 3 (kd) 2 + F {N(η, ūj )} 1 + 1 3 (kd) 2 θ + (1 -θ) ûj j = 0, 1, 2, • • • ,
where θ ∈ [0, 1] is a free parameter. We obtained the best convergence rate for θ = 1 2 . In order to improve the numerical stability of the time stepping method, we will integrate exactly the linear terms in evolution equations

η t + d ūx = -[ η ū ] x , q t + g η x = 1 2 ū2 + 1 2 (d + η) 2 ū2
x -q u x . Taking the Fourier transform and using the relation (4.3) between ū and q, we obtain the following system of ODEs:

ηt + ikd 1 + 1 3 (kd) 2 q = -ik F{ηū} - ikd F {N(η, ūj )} 1 + 1 3 (kd) 2 , qt + ikg η = ik F 1 2 ū2 + 1 2 (d + η) 2 ū2 x -qu .
The next step consists in introducing the vector of dimensionless variables in Fourier space V ≡ (ik η, iω q/g), where ω 2 = gk 2 d/[1 + 1 3 (kd) 2 ] is the dispersion relation of the linearized Serre equations. With unscaled variables in vectorial form, the last system becomes

V t + L • V = N( V ), L ≡ 0 iω iω 0 .
On the right-hand side, we put all the nonlinear terms

N( V ) = k 2 F{ηū} + dk 2 F {N(η, ūj )} /(1 + 1 3 (kd) 2 ) -(kω/g) F 1 2 ū2 + 1 2 (d + η) 2 ū2 x -qu .
In order to integrate the linear terms, we make a last change of variables [51, 35]:

Ŵ t = e (t-t 0 )L • N e -(t-t 0 )L • Ŵ , Ŵ (t) ≡ e (t-t 0 )L • V (t), Ŵ (t 0 ) = V (t 0 ).
Finally, the last system of ODEs is discretized in time by Verner's embedded adaptive 9(8) Runge-Kutta scheme [START_REF] Verner | Explicit Runge-Kutta methods with estimates of the local truncation error[END_REF]. The time step is chosen adaptively using the so-called H211b digital filter [START_REF] Söderlind | Digital filters in adaptive time-stepping[END_REF]62] to meet some prescribed error tolerance (generally of the same order of the fixed point iteration (4.4) precision). Since the numerical scheme is implicit in the velocity variable ū, the resulting time step ∆t is generally of the order of the spatial discretization O(∆x).

Numerical results

In this section we present some numerical results using the finite volume scheme described hereinabove. First, we validate the discretization and check the convergence of the scheme using an analytical solution. Then we demonstrate the ability of the scheme to simulate the practically important solitary wave interaction problem. Throughout this section we consider the initial value problem with periodic boundary conditions unless a special remark is made.

Convergence test and invariants preservation

Consider the Serre equations (3.1), (3.2) posed in the periodic domain [-40, 40]. We solve numerically the initial-periodic boundary value problem with an exact solitary wave solution (2.20) posed as an initial condition. Then, this specific initial disturbance will be translated in space with known celerity under the system dynamics. This particular class of solutions plays an important role in water wave theory [START_REF] Drazin | Solitons: An introduction[END_REF][START_REF] Dougalis | On some Boussinesq systems in two space dimensions: Theory and numerical analysis[END_REF] and it will allow us to assess the accuracy of the proposed scheme. The values of the various physical parameters used in the simulation are given in Table 1.

The error is measured using the discrete L ∞ norm for various successively refined discretizations. The result is shown on Figure 4. As anticipated, the finite volume scheme (black solid line with circles) shows a fairly good second order convergence (with estimated slope ≈ 1.99). During all the numerical tests, the mass conservation was satisfied with accuracy of the order ≈ 10 -14 . This impressive result is due to excellent local conservative properties of the finite volume method. We also investigate the numerical behavior of the scheme with respect to the less obvious invariants H and Q defined in (2.21). These Undisturbed water depth: d 1 Gravity acceleration: g 1 Solitary wave amplitude: a 0.05 Final simulation time: T 2 Free parameter:

β 1/3 Table 1.
Values of various parameters used in convergence tests. invariants can be computed exactly for solitary waves. However, we do not provide them to avoid cumbersome expressions. For the solitary wave with parameters given in Table 1, the generalized energy and momentum are given by the following expressions:

H 0 = 21 √ 7 100 + 7 √ 3 10 log √ 21 -1 √ 21 + 1 ≈ 0.0178098463, Q 0 = 62 √ 15 225 + 2 √ 35 5 log √ 21 -1 √ 21 + 1 ≈ 0.017548002.
These values are used to measure the error on these quantities at the end of the simulation. Convergence of this error under the mesh refinement is shown on Figure 5. One can observe a slight super-convergence phenomenon of the finite volume scheme. This effect is due to the special nature of the solution we use to measure the convergence. This solution is only translated under the system dynamics. For more general initial conditions we expect a fair theoretical 2 nd order convergence for the finite volume scheme. As anticipated, the pseudo-spectral scheme shows the exponential error decay.

Solitary wave interaction

Solitary wave interactions are an important phenomenon in nonlinear dispersive waves which have been studied by numerical and analytical methods and results have been compared to experimental evidence. They also often serve as one of the most robust nonlinear benchmark test cases for numerical methods. We mention only a few works among the existing literature. For example, in [48,56,23] solitary wave interactions were studied experimentally. The head-on collision of solitary waves was studied in the framework of full Euler equations in [23,[START_REF] Chambarel | Head-on collision of two solitary waves and residual falling jet formation[END_REF]. Studies of solitary waves in various approximate models can be found in [START_REF] Li | A Numerical Study of the Exact Evolution Equations for Surface Waves in Water of Finite Depth[END_REF][START_REF] Dougalis | Solitary waves of the Bona-Smith system[END_REF]2,[START_REF] Dutykh | Finite volume schemes for dispersive wave propagation and runup[END_REF][START_REF] Dutykh | Finite volume methods for unidirectional dispersive wave models[END_REF]. To our knowledge, solitary wave collisions for the Serre equations were studied numerically for the first time in the PhD thesis of Seabra-Santos [START_REF] Seabra-Santos | Contribution à l'étude des ondes de gravité bidimensionnelles en eau peu profonde[END_REF]. Finally, there are also a few studies devoted to simulations with full Euler equations [START_REF] Li | A Numerical Study of the Exact Evolution Equations for Surface Waves in Water of Finite Depth[END_REF]35,23].

Head-on collision

Consider the Serre equations posed in the domain [-40, 40] with periodic boundary conditions. In the present section, we study the head-on collision (weak interaction) of two solitary waves of equal amplitude moving in opposite directions. Initially, two solitary waves of amplitude a = 0.15 are located at x 0 = ±20 (other parameters can be found in Table 1). The computational domain is divided into N = 1000 intervals (finite volumes in 1D) of the uniform length ∆x = 0.08. The time step is chosen to be ∆t ≈ 10 -3 . The process is simulated up to time T = 36. The numerical results are presented in Figure 6. As expected, the solitary waves collide quasi-elastically and continue to propagate in opposite directions after the interaction. The value of importance is the maximum amplitude during the interaction process, sometimes referred to as the run-up. Usually, it is larger than the sum of the amplitudes of the two initial solitary waves. In this case, we obtain a run-up of 0.3130 > 2a = 0.3. In order to validate the finite volume simulation, we performed the same computation with the pseudo-spectral method presented briefly in Section 4. We used a fine grid of 1024 nodes and adaptive time stepping. The overall interaction process is visually identical to the finite volume result shown in Figure 6. The run-up value according to the spectral method is 0.3127439 showing again the accuracy of our simulation. The small inelasticity is evident from the small dispersive wave train emerging after the interaction (for an example in a slightly different setting described below, see Figure 16, as first found numerically and experimentally by Seabra-Santos [START_REF] Seabra-Santos | Contribution à l'étude des ondes de gravité bidimensionnelles en eau peu profonde[END_REF].

Overtaking collision

A second type of solitary wave interaction is the overtaking collision (or strong interaction) of two solitary waves of different amplitudes moving in the same direction. Sometimes this situation is also referred to as the following collision or strong interaction. For this case we consider a physical domain [-75, 75] divided into N = 1000 equal control volumes.The initial data consists of two separated solitary waves of different amplitudes moving in the same direction. The solitary wave with larger amplitude moves faster and will overtake the smaller wave. This situation was simulated with the finite volume scheme and the numerical results are presented in Figure 7. The parameters used in this simulation are given in Table 2. The strong interaction is also inelastic with a small dispersive tail emerging after the over-taking (see Figure 15 for a zoom).

-

Experimental validation

In this Section we present a comparison between the classical Serre model solved with our finite volume scheme and one head-on collision experiment from [23]. This specific experiment was already considered in the context of Boussinesq-type systems [START_REF] Dutykh | Finite volume schemes for dispersive wave propagation and runup[END_REF].

We simulate a portion of the wave tank [-0.9, 2.7] (divided into N = 1000 equal control volumes) where the interaction process takes place. The initial data consists of two solitary waves (of different amplitudes in this case) moving in opposite directions. The exact parameters are given in Table 3. Simulation snapshots are presented in Figures 8910111213141516. The general agreement is very good, validating the Serre equations in water wave theory, along with our numerical developments. Figure 16 shows visible dispersive oscillations after the interaction process, numerical evidence of the inelastic character of solitary waves interactions in the framework of the Serre equations.

- x η(x, t)

Figure 16. Dispersive tail after head-on collision of two solitary waves (weak interaction). Small wavelets between two solitary waves clearly indicate that the collision is inelastic.

Conclusions

The current study is devoted to the Serre equations stemming from water wave modeling [59, [START_REF] Barthélémy | Nonlinear shallow water theories for coastal waves[END_REF][START_REF] Barth | Aspects of unstructured grids and finite-volume solvers for the Euler and Navier-Stokes equations[END_REF]. First, we presented a derivation of this model using a relaxed variational principle [START_REF] Clamond | Practical use of variational principles for modeling water waves[END_REF]. We then described an Implicit-Explicit finite volume scheme to discretize the equations. The overall theoretical accuracy of the discretization scheme is of secondorder. This conclusion is confirmed by comparisons with an exact solitary wave solution. The energy conservation properties of our scheme are also discussed and quantified. In order to validate further our numerical scheme, we present a Fourier-type pseudo-spectral method. Both numerical methods are compared on solitary wave interaction problems. The proposed discretization procedure was successfully validated with several numerical tests along with experimental data. In contrast with the highly accurate spectral method, the finite volume method has the advantage of being robust and generalizable to realistic complex situations with variable bathymetry, very steep fronts, dry areas, etc. The present study should be considered as the first step to further generalisations to 2D cartesian meshes [START_REF] Mingham | High-Resolution Finite-Volume Method for Shallow Water Flows[END_REF]13,69].

Figure 2 .

 2 Figure 2. Two exact solutions to the Serre equations. The solitary wave amplitude is equal to a = 0.05. For the cnoidal wave parameters m and a are equal to 0.99 and 0.05 respectively. Other cnoidal wave parameters are deduced from relations (2.18), (2.19).

25 Figure 3 .

 253 Figure 3. Comparison of solitary wave solutions to the Serre and the full Euler equations.

Figure 4 .

 4 Figure 4. Convergence of the numerical solution in the L ∞ norm computed using the finite volume method.

Figure 5 .

 5 Figure 5. Hamiltonian and generalized momentum conservation convergence computed using the finite volume and spectral methods under the mesh refinement. The conserved quantities are measured at the final simulation time.

Figure 6 .

 6 Figure 6. Head-on collision of two equal solitary waves simulated with the finite volume scheme.

Figure 7 .

 7 Figure 7. Overtaking (or following) collision of two solitary waves simulated with the finite volume scheme.

Figure 14 .

 14 Figure 14. Head-on collision of two solitary waves of different amplitudes. Comparison with experimental data [23].

Figure 15 . 5 Table 3 .

 1553 Figure 15. Dispersive tail after overtaking collision of two solitary waves (strong interaction) at T = 120.0.

Table 2 .

 2 Values of various parameters used to simulate the overtaking collision.

	Undisturbed water depth: d	1
	Gravity acceleration: g	1
	Large solitary wave amplitude: a 1 0.6
	Initial position: x 1	-60
	Small solitary wave amplitude: a 2 0.1
	Initial position: x 2	-45
	Final simulation time: T	96
	Free parameter: β	1/3

  Figure 8. Head-on collision of two solitary waves of different amplitudes. Comparison with experimental data [23].

	x 10 -3											x 10 -3
									Simulation										Simulation
									Experimental data	18								Experimental data
												16
												14
												12
	η(x, t)										η(x, t)	8 10
												6
												4
												2
												0
	-2											-2
	0.2	0	0.2	0.4	0.6	0.8 x	1	1.2	1.4	1.6	1.8	-0.2	0	0.2	0.4	0.6	0.8 x	1	1.2	1.4	1.6	1.8
				(a) t = 18.5 s								(b) t = 18.6 s
	x 10 -3											x 10 -3
									Simulation										Simulation
									Experimental data	18								Experimental data
												16
												14
												12
	η(x, t)										η(x, t)	8 10
												6
												4
												2
												0
	-2											-2
	-0.2	0	0.2	0.4	0.6	0.8 x	1	1.2	1.4	1.6	1.8	-0.2	0	0.2	0.4	0.6	0.8 x	1	1.2	1.4	1.6	1.8
				(a) t = 18.7 s								(b) t = 18.8 s
	Figure 9. Head-on collision of two solitary waves of different amplitudes.
					Comparison with experimental data [23].

  Figure 10. Head-on collision of two solitary waves of different amplitudes. Comparison with experimental data [23]. Note the difference in vertical scales on the left and right images.Figure 11. Head-on collision of two solitary waves of different amplitudes. Comparison with experimental data [23].

		x 10 -3
		18								Simulation Experimental data	0.025								Simulation Experimental data
		16
		14											0.02
		12
	η(x, t)	6 8 10										η(x, t)	0.01 0.015
		4
													0.005
		2
		0
													0
		-2
		-0.2	0	0.2	0.4	0.6	0.8 x	1	1.2	1.4	1.6	1.8	-0.2	0	0.2	0.4	0.6	0.8 x	1	1.2	1.4	1.6	1.8
					(a) t = 18.92 s								(b) t = 19.0 s
										Simulation										Simulation
		0.025								Experimental data	0.025								Experimental data
		0.02											0.02
	η(x, t)	0.01 0.015										η(x, t)	0.01 0.015
		0.005											0.005
		0											0
		-0.2	0	0.2	0.4	0.6	0.8 x	1	1.2	1.4	1.6	1.8	-0.2	0	0.2	0.4	0.6	0.8 x	1	1.2	1.4	1.6	1.8
					(a) t = 19.05 s								(b) t = 19.1 s

  Figure 12. Head-on collision of two solitary waves of different amplitudes. Comparison with experimental data [23].Figure 13. Head-on collision of two solitary waves of different amplitudes. Comparison with experimental data [23].

		x 10 -3 x 10 -3											x 10 -3 x 10
										Simulation Simulation										Simulation
		18								Experimental data Experimental data	18 18								Experimental data
		16											16 16
		14											14 14
		12											12 12
	η(x, t) η(x, t)	8 10										η(x, t)	8 10 8 10
		6											6 6
		4											4 4
		2											2 2
		0											0 0
		-2 -2											-2 -2
		0.2 -0.2	0 0	0.2 0.2	0.4 0.4	0.6 0.6	0.8 x 0.8 x	1 1	1.2 1.2	1.4 1.4	1.6 1.6	1.8 1.8	-0.2 -0.2	0 0	0.2 0.2	0.4 0.4	0.6 0.6	0.8 x 0.8	1 1	1.2 1.2	1.4 1.4	1.6 1.6	1.8 1.8
					(a) t = 19.15 s (a) t = 19.85 s								(b) t = 19.19 s
		x 10 -3											x 10 -3
										Simulation										Simulation
		18								Experimental data	18								Experimental data
		16											16
		14											14
		12											12
	η(x, t)	8 10										η(x, t)	8 10
		6											6
		4											4
		2											2
		0											0
		-2											-2
		-0.2	0	0.2	0.4	0.6	0.8 x	1	1.2	1.4	1.6	1.8	-0.2	0	0.2	0.4	0.6	0.8 x	1	1.2	1.4	1.6	1.8
					(a) t = 19.33 s								(b) t = 19.5 s

Acknowledgments References Acknowledgments D. Dutykh acknowledges the support from French "Agence Nationale de la Recherche", project "MathOcéan" (Grant ANR-08-BLAN-0301-01) along with the support from ERC under the research project ERC-2011-AdG 290562-MULTIWAVE. P. Milewski acknowledges the support of the University of Savoie during his visits in 2011.


Julien Veytizou
email: julien.veytizou@grenoble-inp.fr

Hugo Xuereb

Guillaume Thomann

Design of a clip product based on customer needs for playing acoustic music

Keywords: user environment, design, disabled people, acoustic music, robotic musical

Inclusion of uses and users in product design remains a challenge to take up; especially when their characterizations are very specific (it's the case with disabled persons). In the musical domain, a lot of adapted interfaces are manufactured to enable users with disabilities to play music from digital audio. But few of them allow the music practice on acoustic instruments, which is one of the goal identified by the AE2M non-profit association (Ergonomic Adaptation of the Musical Material). In this specific context, this paper presents the design of a universal product which transforms any user environment object to a personalized interface, to play percussion instruments.

Introduction

Traditionally, designers create products or interfaces for able-bodied users, assumed to have normal motor skills [START_REF] Keates | Towards a practical inclusive design approach[END_REF]. Atypical users, including people with disabilities, remain outside of these skills and capabilities standards. They must adapt to with existing systems, using technical aids. Unfortunately, these technologies are unusable for people with temporary impairments. They are also not adapted to users whose abilities change over time. The data, detailed in a research report conducted by INSEE (France's National Institute of Statistics and Economic Studies), show a significant need of assistance in handling (only 35.3% satisfaction for remote handling) [START_REF] Brouard | Le handicap en chiffres[END_REF]. Nevertheless, these means are often abandoned because of their price and complex use. The main reasons for giving up are due to a minimum consideration of future users in design activities.

There are two main groups among system designers for disabled people [START_REF] Sperandio | Designing technological devices for a normal population, namely also including disabled people and the elderly[END_REF]:  The first which go to certain well-defined needs and seek a technological solution (social pull). These proponents generally have a good knowledge of disabilities. They are looking primarily practical systems, efficient and robust, sometimes ignoring the more advanced technologies. The risk here might be to focus too much on obsolete technology, unusable or too expensive.

 The second which go to a technology they have mastered and are looking for applications in the field of disabilities (technology push). The risk here might be to propose solutions for imaginary or non-priority needs. Integrate multidisciplinary teams in the design process can avoid these risks.

Our research work consists to the development and the manufacture of physical interfaces facilitating access to instrumental music for children with heavy physical disabilities. Currently, disabled user environment are not adapted to their needs and human machine interface. One solution is a custom design of products for each user motor skills. But this process can be lengthy and difficult because of the myriad of user profiles. Moreover, one of the most objectives is to reach low cost production and the large diffusion for disabled customers.

The research question can be written as below: How to design a universal product which would transform any user environment to a personalized interface?

To answer this question, our methodology is the following: Section 2 presents an identification of needs in our context thanks to a literature review. Three main parts are exposed: future user characteristics, existing adapted system to play music, human machine interaction in musical practice. Thanks to the needs list defined previously, section 3 offers a first prototype based on piezoelectric sensors technology. Its structure is defined with a representation of human machine interaction in our specific context. Its architecture is described in functional component decomposition. After a presentation of test methods, section 4 shows the first results when using this interface for a musical playing.

Identification of needs

As the NF EN ISO 9241-210 UCD Process 1 recommends, designers take into account the contextual factors as well as the requirements of the future users before producing design solutions.

User with disability

User characteristics define the context in which the product is used. The futures users are children with disability. We are focused on the disability ICF-CY definition (International Classification of Functioning, Disability and Health for Children and Youth), proposed by the WHO (World Health Organization) to understand the terminology of disability [START_REF] Icf | International Classification of Functioning, Disability and Health (ICF) Disability and Health[END_REF]. This definition provides a common language that can record the problems occurring during infancy, childhood and adolescence. It is based on body functions and body structures, activity limitations, restrictions participation and environmental factors significant for children and youth (Figure 1). Disability therefore results from the interaction between a person with health problems and contextual factors (personal and environmental). Impairments refer to modifications of body functions or anatomical structures, such as paralysis. Activity Limitations are difficulties an individual may experience in these activities (walking, eating, playing etc...). Participation restrictions are problems a person may meet to participate in a situation of daily life, such as transportation inaccessibility.

Existing adapted system

This section analyzes existing systems for the practice of music. This analysis can reveal needs, problems and constraints that have not been taken into account.

To play music with digital audio

Interfaces designed for disabled people allow musical practice from digital audio:  The BAO-PAO fits perfectly with the computer-assisted music because it needs to operate a computer and special software. It consists of four rigid steel arches, each terminated by two spheres at the ends. A laser beam passes between the upper sphere and the lower sphere. It generates a sound when the musician crossed it with a drumstick or hand [START_REF][END_REF].  TouchTone is an electronic musical instrument. Its goal is to develop musical ability, to develop bimanual coordination and to increase social participation of children with hemiplegic cerebral palsy [START_REF] Bhat | TouchTone: an electronic musical instrument for children with hemiplegic cerebral palsy[END_REF].  The GUI//RO is a stylus controlled instrument, augmented with resistive force feedback. The system is manufactured through an electromagnetic coil controlled by Arduino board. Sounds and haptic effects are generated by the position on the stylus in touch-sensitive screen [START_REF] Müller | Reflective Haptics : Resistive Force Feedback for Musical Performances with Stylus-Controlled Instruments[END_REF]. With these solutions, users can play any instrument thanks to synthesize sound. Nerveless, the musical instrument is dematerialized. Indeed, the user doesn't visualize a real musical instrument for the music game. Moreover, transportation and equipment installation can be difficult and not adapted with the environment of disabled people.

2.2.2

To play acoustic music with AE2M project Few systems allow the music practice on acoustic instruments. This is the goal of the AE2M non-profit association (Ergonomic Adaptation of the Musical Material) [START_REF] Ae2m | AE2M Project[END_REF]. It provides access to instrumental music for children with disabilities, with the same autonomy level as non-disabled. It's a multidisciplinary team which is surrounded by three main applications specialties (Figure 2): Musicians (by specialized professors in instrumental music and musical manufactory), Paramedical specialists (by doctors, specialist therapists) and Engineers (principally by teachers and engineering students).

Fig. 2. Competencies triangle of AE2M project

For a better understanding of this competencies triangle, it is necessary to develop the roles and the work of these specialties in the project. Engineers should be familiar with the properties of the musical instrument which children would try to play. So they consult the music specialists of the project. Moreover, engineers have to discuss with the paramedical specialists who spend all their times with disabled children. They know some physical capacities of these children, those needed for the project development.

Generally, interfaces manufactured by the association AE2M allow children to play percussion instruments A specific solution realized is an electromechanical system to allow a strike of the drumstick on percussion instruments (Figure 3). The user activates the system by the means of accessories (push button, pressure effort, etc.) positioned around him according to his ability to perform body motions. With this system, the musician couldn't perform different strike velocity on the instrument and the preparation of musical activities could be long because of the myriad of user profiles. Depending on motor skills (movement, strength ...) and the user environment (wheelchair, tablet ...), this system may have difficulty to be integrated naturally.

Human machine interaction for musical practice

To define the user tasks, it is important to design the user-system interaction. Metois proposes a generic definition that considers the "musical gesture" as anything that fills the space between the musical intentions (cognition, psychology, and musicology) and musical sounds (waveforms physical) [START_REF] Métois | Musical Sound Information -Musical Gestures and Embedding Systems[END_REF]. "Musical gesture" includes the "instrumental gesture" concept. Cadoz defines "Gestural Canal" as a means of acting on the physical world and as a means of two way communication: sending and receiving information. Then, his definition considers the Instrumental Gesture as a communication modality specific at "Gestural Canal" as follow [START_REF] Cadoz | Gestural channel of the man machine communication: instrumental communication[END_REF]:  It applies to a material object and there is physical interaction with it,  In this context of interaction, it will occur physical phenomena,  These phenomena can then become communicational messages.

Synthesis

These different parts allow describing the needs as well as requirements of our interface realization. In our musical context, we must offer to disabled child with particular health condition, a product which would allow them to enjoy musical play (activities) and even to participate in concert (participation) with the same autonomy level as non-disabled (body structure and function). For musicians and paramedical team, the system must be easy to install and must be integrated naturally to the environment of the user. In addition, the delay between the musical intention of the user and the system response must be transparent. The musician recommends a delay time less than ten milliseconds (between the action of the children and the action of mechanical system). The system must also perform different strike velocity on the instrument. For engineers, the interface must contain a minimum of electronic equipment. If maintenance must be performed, it should be easy. Our proposition is to transform any object in the user environment as an action unit on electromechanical systems.

3

Concept solution proposition

This section offers a concept solution using the previous synthesis. The product structure is defined with a representation of human machine interaction in our specific context. Its architecture is described in functional component decomposition.

Product structure

The product structure is based on our proper representation of human machine interaction using the definition of "musical and instrumental gesture". In our context, the user has a musical intention (Figure 4). He will perform an action on an object environment. This interaction will produce physical phenomena (force, pressure or vibration) that can be controlled by the user (duration, amplitude). These phenomena can thus become communicational messages to the electromechanical system for the strike of the drumstick on the percussion instrument. The essential elements in designing of our future product are:  A clip sensor for the physical phenomena detection.  An electronic board to transform these phenomena in action message for actuator positioned on the musical instruments.

Product architecture

From this previous definition as well as the identification of needs, this part describes the product architecture. It's divided principally in two parts (Figure 5):  Control part with sensor clip, electronic board and microcontroller input  Actuator part with microcontroller output, electronic board and mechanical system

The sensor emits a signal. It is interpreted by a microcontroller to control a mechanical system placed on the instrument. These control actions are based on the musical intentions of the user. In our case, the Arduino board corresponds to our expectations. This material is open source, free and easy for everyone to use. After this interpretation, Arduino can generate PWM signals to control a solenoid. But the output current is too weak to command this system. A power signal processing is therefore necessary. All this process must be realized in less than 10 milliseconds. This section shows an evaluation of the design solution against user needs and requirements. Test methods are presented for the sensor choice and the study of the strike velocity.

Tests methods

The sensor choice will be based on the following main criteria: accuracy, sensitivity, size and of price. Three different sensors performances are studied. They are different in size and form: two circular (Multicomp and Murata) and one rectangular (Prowave). They are arranged in a circular arc at an equal distance (d) from a disturbance area where it will be throw, at the same height, different weights. Sensitivity is measure by an oscilloscope (Figure 6). To study the strike velocity, Musician asks to the user to perform three different strikes (Piano, Mezzo, and Forte) on a plastic box where the clip sensor is located (n°1 on Figure 7). After signal processing (n°2 on Figure 7), the drumstick movement (n°3 on Figure 7) on the instrument (n°4 on Figure 7) is filmed. Speed strike is analyzed with the Kinovea software, a free and open source solution for video analysis. It is mostly used by sports coaches, animation artists and ergonomics engineers.

Sensor choice

Figure 8 shows the sensitivity analyzes of piezoelectric sensors Murata, Multicomp and Prowave. The signal amplitude is measured when the weight is throwing in the disturbance area. For a 500 grams weight, Murata gives a signal of 30.2 volts, Multicomp 9.8 volts and Prowave 2.9 volts. Multicomp exceeds 5 volts for weight greater than 120 grams, Murata for few grams. However, Prowave never exceed this value. For the next experimentation, the strike velocity is analyzed with Murata.

Strike velocity

Figure 9 shows the strike speed of the drumstick controlled by the mechanical system. Three strike speeds are saved: the first at 0.4 m/s, the second at 0.7 m/s and the third at 1.9 m/s.

Discussion

Our interface allows transforming any user environment object to a percussion musical instrument. Its architecture is simple with low-cost and easy system Arduino, little electronic equipment to facilitate maintenance. Our system is integrated more naturally to the environment of user. For the sensor choice, Murata sensitivity is more efficient. Indeed, a child with low motion disabilities could control musical instrument. The response time of the system is less than two milliseconds (inferior at ten milliseconds recommended by the musician). It's possible to adjust the strike velocity on the instrument with different actions on an object.

Conclusion

Identification of needs and collaboration with the AE2M multidisciplinary team has allows to generate uses specifications of a product. These specifications have been destined to a team of Research and Development, expert in new technology. Thus, this methodology has led to the design of a universal product which would transform any user environment to a personalized interface. The product structure has defined by a representation of human machine interaction in our context. The product archi-tecture is described in functional component decomposition. After these first encouraging results, this system facilities access to instrumental music for users with heavy physical disabilities, but also to increase the uses performance by musicians. The perspectives are focused by the clip design, the increase of the performance of the mechanical system and the "intelligence" of this system (the actuator part is able to adapt to user motor skills).

Fig. 1 .

 1 Fig. 1. ICF-CY definition: International Classification of Functioning, Disability and Health for Child and Youth

Fig. 3 .

 3 Fig. 3. Use of electromechanical magnet for a playing game

  The development is inspired from the different work realize by the LEMUR (League of Electronic Musical Urban Robots) which is the producer of the famous Pat Metheny's Orchestrion composed of forty robotic musical instruments. It uses principally stepper motor, hobby servo and solenoids as actuator. They are controlled by a Pulse Width Modulation Signal (PWM)[START_REF] Kapur | A comparison of solenoid-based strategies for robotic drumming[END_REF][START_REF] Kapur | A Pedagogical Paradigm for Musical Robotics[END_REF].

Fig. 4 .

 4 Fig. 4. Human Machine interaction of our interface

Fig. 8 .

 8 Fig. 8. Architecture of the smart clip sensor

Fig. 6 . 7 .

 67 Fig. 6. Test for sensor choice Fig. 7. Test for signal processing and velocity

Fig. 8 .

 8 Fig. 8. Sensitivity analyze of piezoelectric sensors

Fig. 9 .

 9 Fig. 9. Strike velocity study with Murata sensor


Pierre Bosch

Thomas Simon

On the self-decomposability of the Fréchet distribution

Keywords: 2000 Mathematics Subject Classification. 60E07, 60G51, 60G70, 62E10 Exponential functional -Fréchet distribution -Gamma subordinator -Infinite divisibility -Self-decomposability. 1

. We also review different analytical and probabilistic interpretations of the infinite divisibility of Γ -α t for t, α > 0.

Introduction

The extreme value theorem -see e.g. Theorem 8.13.1 in [START_REF] Bingham | Regular variation[END_REF] -states that non-degenerate distribution functions arising as limits of properly renormalized running maxima of i.i.d. random variables belong to one of the families

F 0 (x) = e -e -x , x ∈ R, or F ξ (x) = 1 -e -x 1/ξ if ξ > 0 e -x 1/ξ if ξ < 0 , x > 0.
The distribution F 0 is known as the Gumbel distribution, whereas F ξ is called a Weibull distribution for ξ > 0 and a Fréchet distribution for ξ < 0. In the following, we denote by X ξ the random variable with distribution function F ξ . Observe that 1 -X ξ ξ d -→ X 0 as ξ → 0, so that the above parametrization is continuous in ξ. In the present paper we are interested in the self-decomposability (SD) of X ξ , referring e.g. to Section 15 in [START_REF] Sato | Lévy processes and infinitely divisible distributions[END_REF] for an account on self-decomposability. The Gumbel distribution is SD because of the identities

X 0 d = -log L d = -α log L + α log S α
for every α ∈ (0, 1), where here and throughout L stands for the standard exponential variable and S α for the standard positive α-stable variable -see e.g. Exercise 29.16 in [START_REF] Sato | Lévy processes and infinitely divisible distributions[END_REF] for a proof of the second identity. If ξ ∈ (0, 1) then the variable X ξ is not infinitely divisible (ID) and hence not SD, because of its superexponential distribution tails -see e.g. Theorem 26.1 in [START_REF] Sato | Lévy processes and infinitely divisible distributions[END_REF]. When ξ ≥ 1, the variable X ξ has a completely monotone density and is ID by Goldie's criterion -see e.g. Theorem 4.2 in [START_REF] Steutel | Some recent results in infinite divisibility[END_REF], or by the ME property which makes it the first-passage time of some continuous time Markov chain -see e.g. Chapter 9 in [START_REF] Bondesson | Generalized Gamma convolutions and related classes of distributions and densities[END_REF] for an account. When |ξ| ≥ 1, the identity in law

X ξ d = L ξ
and the HCM theory of Thorin and Bondesson [START_REF] Bondesson | Generalized Gamma convolutions and related classes of distributions and densities[END_REF] show that the distribution of X ξ is a generalized Gamma convolution (GGC) and is hence SD -see Example 4.3.4 in [START_REF] Bondesson | Generalized Gamma convolutions and related classes of distributions and densities[END_REF]. The natural question whether X ξ is SD or even ID for ξ ∈ (-1, 0) was first raised by Steutel in 1973 -see Section 3.4 in [START_REF] Steutel | Some recent results in infinite divisibility[END_REF], and has remained open ever since. In section 4.5 of [START_REF] Bondesson | Generalized Gamma convolutions and related classes of distributions and densities[END_REF] -see also the Appendix B.3 of [START_REF] Steutel | Infinite divisibility of probability distributions on the real line[END_REF], this problem is rephrased in the broader context of generalized Gamma distributions. The latter are power transformations of Γ t where {Γ t , t ≥ 0} is the Gamma subordinator, and can be thought of as size-biased samplings of X ξ when ξ < 0, in view of the formulae

E[f (Γ ξ t )] = E[f (X ξ )X u ξ ] E[X u
ξ ] valid for every f bounded continuous and t > 0, with u = (t -1)/ξ. Recall in passing that Steutel's equation -see e.g. Theorem 51.1 in [START_REF] Sato | Lévy processes and infinitely divisible distributions[END_REF] -establishes a precise link between sizebiased sampling of order one and infinite divisibility for integrable positive random variables. In this note, we provide an answer to the above questions of [START_REF] Steutel | Some recent results in infinite divisibility[END_REF][START_REF] Bondesson | Generalized Gamma convolutions and related classes of distributions and densities[END_REF].

Theorem. For every ξ ∈ (-1, 0) and t > 0, the random variable Γ ξ t is SD. As a direct consequence of this result, all Fréchet distributions are SD and the extreme value distribution F ξ is ID if and only if ξ ∈ (0, 1). Contrary to the case |ξ| ≥ 1, our argument is probabilistic and consists in showing that Γ ξ t is distributed as the exponential functional of some spectrally negative Lévy process. This extends a classical result of Dufresne [START_REF] Dufresne | The distribution of a perpetuity, with applications to risk theory and pension funding[END_REF] for the case ξ = -1. The identification is made possible thanks to a entire moment method due to Bertoin-Yor [START_REF] Bertoin | On the entire moments of self-similar Markov processes and exponential functionals[END_REF], which applies in our context as a case study. The proof is given in the next section.

In Section 3, we review the possible interpretations of the infinite divisibility of Γ ξ t for ξ < 0. The classical case ξ = -1 allows at least four different formulations in terms of processes, and also an explicit computation of the Lévy density which shows the GGC property without the HCM argument. For ξ < -1 the ID property is only known by analytical means and there is no direct probabilistic explanation, save for the case t = 1 by subordination or, tentatively, the spectral theory of a certain spectrally positive Markov processes. The situation for ξ ∈ (-1, 0) is exactly the opposite since in addition to the exponential functional argument, the ID property can also be obtained rigorously by a first-passage time argument for a spectrally positive Markov processes. On the other hand there is no analytic proof of the ID property for ξ ∈ (-1, 0). In this situation the GGC character of Γ ξ t remains in particular an open question, which we plan to tackle in some further research.

Proof of the Theorem

We begin with a computation on the Gamma function.

Lemma. For every α ∈ (0, 1) and u, t > 0 one has

uΓ(t + α(u + 1)) Γ(t + αu) = Γ(t + α) Γ(t) u + 0 -∞ (e ux -1 -ux)f α,t (x)dx,
where

f α,t (x) = e (1+t/α)x (α + e x/α + t(1 -e x/α )) αΓ(1 -α)(1 -e x/α
) α+2 is the density of a Lévy measure on (-∞, 0).

Proof. We set λ = t + αu > 0 and compute Γ(λ + α) Γ(λ) = λβ(λ + α, 1 -α) Γ(1 -α) = λ Γ(1 -α) +∞ 0 e -(α+λ)x (1 -e -x ) α dx = α Γ(1 -α) +∞ 0 (1 -e -λx ) e -αx (1 -e -x ) α+1 dx
where the second equality comes from a change of variable and the third from an integration by parts. This yields

uΓ(t + α(u + 1)) Γ(t + αu) = αu Γ(1 -α) +∞ 0 (1 -e -(t+αu)x ) e -αx (1 -e -x ) α+1 dx = Γ(t + α) Γ(t) u + αu Γ(1 -α) 0 -∞ (1 -e αux ) e (α+t)x (1 -e x ) α+1 dx = Γ(t + α) Γ(t) u + 0 -∞ (e ux -1 -ux)f α,t (x)dx
where again, the second equality comes from a change of variable and the third from an integration by parts.

Remarks. (a) The above proof follows [START_REF] Bertoin | On subordinators, self-similar Markov processes and some factorizations of the exponential variable[END_REF] p. 102. Notice in passing that some computations performed in [START_REF] Bertoin | On subordinators, self-similar Markov processes and some factorizations of the exponential variable[END_REF] are slightly erroneous. For example the subordinator whose exponential functional is distributed as τ -α α (with the notation of [START_REF] Bertoin | On subordinators, self-similar Markov processes and some factorizations of the exponential variable[END_REF]) has no drift, but it is also killed at rate 1/Γ(1 -α).

(b) The above decomposition extends to α = 1 since

uΓ(t + (u + 1)) Γ(t + u) = u(t + u)
is the Lévy-Khintchine exponent of a drifted Brownian motion (the latter was already noticed in [START_REF] Bertoin | On the entire moments of self-similar Markov processes and exponential functionals[END_REF] -see Example 3 therein -in order to recover Dufresne's identity). However, such a formula does not seem to exist for α > 1.

End of the proof. Fix ξ ∈ (-1, 0), t > 0, and set α = -ξ ∈ (0, 1) for simplicity. The entire moments of Γ α t are given for every n ≥ 1 by

E[Γ αn t ] = Γ(t + αn) Γ(t) = Γ(t + α) Γ(t) × • • • × Γ(t + αn) Γ(t + α(n -1)) = m ψ(1) . . . ψ(n -1) (n -1)!
with the notation

ψ(u) = uΓ(t + α(u + 1)) Γ(t + αu) = Γ(t + α) Γ(t) u + 0 -∞ (e ux -1 -ux)f α,t (x)dx
by the Lemma, and

m = Γ(t + α) Γ(t) = ψ ′ (0+).
It is clear that ψ is the Lévy-Khintchine exponent of a spectrally negative Lévy process Z with infinite variation and mean m > 0. By Proposition 2 in [START_REF] Bertoin | On the entire moments of self-similar Markov processes and exponential functionals[END_REF], this entails

E[Γ αn t ] = E[I -n ] for every n ≥ 1, where I is the exponential functional of Z: I = ∞ 0 e -Zs ds.
Since Z has no positive jumps, Proposition 2 in [START_REF] Bertoin | On the entire moments of self-similar Markov processes and exponential functionals[END_REF] shows also that the random variable 1/I is moment-determinate, whence Γ ξ where Z ′ is an independent copy of Z and the second equality follows from the Markov property at T y . This shows that for every c ∈ (0, 1) there is an independent factorization

I = cI + I c
for some random variable I c , in other words that

I d = Γ ξ t is self-decomposable.
Remarks. (a) By the above Remark 1 (b), it does not seem that Γ ξ t is distributed as the exponential functional of a Lévy process for ξ < -1. It would be interesting to have an explanation of the infinite divisibility of Γ ξ t in terms of processes when ξ < -1. See next section for a more precise conjecture in the case t = 1.

(b) The self-decomposability of S α α for every α ∈ (0, 1) has been shown by Patie [START_REF] Patie | A refined factorization of the exponential law[END_REF] in using the same kind of argument. Specifically, one can write

E[S nα α ] = Γ(1 + n) Γ(1 + αn) = m ψ(1) . . . ψ(n -1) (n -1)!
where we use the same notation as above and, correcting small mistakes made in Paragraph 3.2 of [START_REF] Patie | A refined factorization of the exponential law[END_REF],

ψ(u) = u Γ(1 + α) + 0 -∞ (e ux -1 -ux) (1 -α)e x/α ((2 -α)e x/α + (1 -e x/α )) α 2 Γ(1 + α)(1 -e x/α ) 3-α dx
is the Lévy-Khintchine exponent of some spectrally negative Lévy process with positive mean. Setting α = t = 1/2 and comparing the above expression to the one in the Lemma, one can check the well-known identity in law (2.1)

S 1/2 = 1 2 Γ 1/2 •
The present paper shows that all positive powers of S 1/2 are SD and one may wonder if the same is true for S α with any α ∈ (0, 1). See [START_REF] Jedidi | Further examples of GGC and HCM densities[END_REF] for related results and also for a characterization of the SD property of negative powers of S α when α ≤ 1/2.

Further remarks and open questions

In this section we would like to review several existing or tentative approaches for the ID, SD and GGC properties of the distribution of Γ ξ t or X ξ when ξ ≤ 0.

3.1. The case ξ = 0. This is a rather specific situation but we include it here for completeness. As mentioned in Section 3.4 of [START_REF] Steutel | Some recent results in infinite divisibility[END_REF], the SD property of the two-sided X 0 is a direct consequence of the extreme value theory because

(3.1) L 1 + L 2 2 + • • • + L n n -log n d = max(L 1 , . . . , L n ) -log n d -→ X 0
as n → +∞, where L 1 , . . . , L n are independent copies of L ∼ Exp (1). The above identity and convergence in law, readily obtained in comparing Laplace transforms and distribution functions, yield after some further computations the following closed expression for the Laplace transform of X 0 :

E[e -λX 0 ] = Γ(1 + λ) = exp -γλ + ∞ 0 (e -λx -1 + λx) dx x(e x -1)
,

where γ is Euler's constant. The complete monotonicity of 1/(e x -1) shows then that X 0 is an extended GGC in the sense of Chapter 7 in [START_REF] Bondesson | Generalized Gamma convolutions and related classes of distributions and densities[END_REF]. See also Exercise 18.19 in [START_REF] Sato | Lévy processes and infinitely divisible distributions[END_REF] and Example 7.2.3 in [START_REF] Bondesson | Generalized Gamma convolutions and related classes of distributions and densities[END_REF] for another argument based on Pick functions, recovering (3.1).

3.2.

The case ξ = -1. This is the classical situation, very well-known, but we give some details for comparison purposes. The ID property of X -1 can first be understood by the sole fact that

lim n→+∞ nx 1 + nx n = e -1/x
because the left-hand side is the first-passage time distribution function of a certain birth and death process -see Theorem 3.1 and (3.3) in [START_REF] Steutel | Some recent results in infinite divisibility[END_REF]. The random variable Γ -1 t is also a GIG and is hence distributed as the unilateral first-passage time of a diffusion [START_REF] Barndorff-Nielsen | First hitting time models for the generalized inverse Gaussian distribution[END_REF], which explains its infinite divisibility by continuity and the Markov property. More precisely one has

(3.2) 1 4Γ t d = inf{u > 0, X t u = 0}
where {X t u , u ≥ 0} is a Bessel process of dimension 2(1 -t) starting from one. The SD property follows as above from Dufresne's identity [START_REF] Dufresne | The distribution of a perpetuity, with applications to risk theory and pension funding[END_REF], which reads

(3.3) 2 Γ t d = ∞ 0 e Bu-tu/2 du
where {B u , u ≥ 0} is a standard linear Brownian motion. Also, Exercise 16.4 in [START_REF] Sato | Lévy processes and infinitely divisible distributions[END_REF] shows that Γ -1 t is the one-dimensional marginal of a certain self-similar additive process, whence its self-decomposability by Theorem 16.1 in [START_REF] Sato | Lévy processes and infinitely divisible distributions[END_REF]. The link between this latter interpretation and (3.2) and (3.3) has been explained in depth in [START_REF] Yor | On certain exponential functionals of real-valued Brownian motion[END_REF].

It does not seem that these four interpretations can provide any explicit information on the Lévy-Khintchine exponent of Γ -1 t . But in this case analytical computations can also be performed. More precisely, taking for simplicity the same normalization as in (3.2) and setting ϕ t (λ) = -log E[e -λ/4Γt ], formulae (7.12.23), (7.11.25) and (7.11.26) in [START_REF] Erdelyi | Higher transcendental function[END_REF] entail

(3.4) ϕ ′ t (λ) = K t-1 ( √ λ) 2 √ λK t ( √ λ)
where K t is the Macdonald function. This shows ϕ ′ 1/2 (λ) = 1/2 √ λ viz. ϕ 1/2 (λ) = √ λ when t = 1/2, and one recovers the identity (2.1). For t = 3/2, one obtains

2ϕ ′ 3/2 (λ) = 1 1 + √ λ = E[e -λ(L 2 ×S 1/2 ) ] = ∞ 0 1 λ + x √ x dx π(1 + x)
where the first equality follows from Formula (7.2.40) in [START_REF] Erdelyi | Higher transcendental function[END_REF], and the third equality from Exercise 29.16 in [START_REF] Sato | Lévy processes and infinitely divisible distributions[END_REF] and (2.2.5) in [START_REF] Bondesson | Generalized Gamma convolutions and related classes of distributions and densities[END_REF]. This means precisely -see (3.1.1) in [START_REF] Bondesson | Generalized Gamma convolutions and related classes of distributions and densities[END_REF] -that the distribution of 1/4Γ 3/2 is a GGC with zero drift and Thorin measure

U 3/2 (dx) = √ x dx 2π(1 + x) •
The latter property can be extended to all values of t thanks to a result originally due to Grosswald [START_REF] Grosswald | The Student t-distribution of any degree of freedom is infinitely divisible[END_REF] on Student's distribution. Together with (3.4), the main theorem in [START_REF] Grosswald | The Student t-distribution of any degree of freedom is infinitely divisible[END_REF] entails namely that the distribution of 1/4Γ t is a GGC with zero drift and Thorin measure

U t (dx) = 1 π 2 x(J 2 t ( √ x) + Y 2 t ( √ x))
where J t and Y t are the usual Bessel functions of the first kind -see [START_REF] Erdelyi | Higher transcendental function[END_REF] p. 4.

3.3.

The case ξ ∈ (-1, 0). In this situation, the present paper yields an interpretation of the self-decomposability of Γ ξ t by the identification

Γ ξ t d = ∞ 0 e -Zu du,
where Z is a spectrally negative Lévy process. Another explanation, similar to (3.2), can then be obtained by the Lamperti transformation -see e.g. the introduction of [START_REF] Bertoin | On the entire moments of self-similar Markov processes and exponential functionals[END_REF] for an account and references. More precisely, setting

Y u = exp[-Z τu ]
with the notation τ u = inf{s > 0, s 0 e -Zv dv > u}, then Y = {Y u , 0 ≤ u < Γ ξ t } is a spectrally positive Markov process (which is also self-similar) starting from one and we have

Γ ξ t d = inf{u > 0, Y u = 0},
so that the infinite divisibility of Γ ξ t (but not, or at least not directly, its self-decomposability) is a consequence of the Markov property and the absence of negative jumps for Y. It would be interesting to see if Γ ξ t could be embedded in some self-similar additive process analogous to the Brownian escape process of the case ξ = 1, described in Exercise 16.4 of [START_REF] Sato | Lévy processes and infinitely divisible distributions[END_REF].

Our main result can also be interpreted analytically in terms of generalized Bessel functions. Setting α = -ξ and writing down

(3.5) E[e -λΓ ξ t ] = 1 αΓ(t) ∞ 0 x -tα-1 e -λx+x -1/α dx = Z t/α 1/α (λ) αΓ(t)
with the notation (1.7.42) of [START_REF] Kilbas | Theory and applications of fractional differential equations[END_REF], the infinite divisibility of Γ ξ t entails that the function

(3.6) λ → - Z ν ′ ρ (λ) Z ν ρ (λ)
is completely monotone for any ρ > 1 and ν > 0. One might ask if the latter function is also a Stieltjes transform, which is equivalent to the GGC property for the distribution of Γ ξ tsee Chapter 3 in [START_REF] Bondesson | Generalized Gamma convolutions and related classes of distributions and densities[END_REF]. Indeed, it is very natural to conjecture such a property for ξ ∈ (-1, 0) in view of the above cases ξ = 0 and ξ = -1. Compared to classical Bessel functions, the theory of generalized Bessel functions is however rather incomplete, and proving like in [START_REF] Grosswald | The Student t-distribution of any degree of freedom is infinitely divisible[END_REF] that the function (3.6) is a Stieltjes transform is believed to be challenging.

3.4. The case ξ < -1. In this situation the GGC property of the distribution of Γ ξ t is most quickly obtained from the HCM character of the density function -see Chapter 5 and especially Example 5.5.2 in [START_REF] Bondesson | Generalized Gamma convolutions and related classes of distributions and densities[END_REF]. Notice that this analytical argument extends to ξ = -1 but not to ξ ∈ (-1, 0) since otherwise Γ -ξ t would also have a HCM density and hence be ID, which is false. This entails that the function in (3.6) is indeed a Stieltjes transform for any ρ ∈ (0, 1) and ν > 0, and it would be interesting to identify the underlying Thorin measure as in Grosswald's theorem.

A probabilistic interpretation of the self-decomposability of Γ ξ 1 = L ξ can also be given by Bochner's subordination. Setting α = -1/ξ ∈ (0, 1), one has indeed

L ξ d = L -1 × S α d = S α L -α .
where {S α u , u ≥ 0} stands for the α-stable subordinator with marginal S α 1 d = S α . Since L -α is SD by our result, this means that L ξ is the marginal of some subordinator which is itself subordinated to the α-stable one, and Proposition 4.1. in [START_REF] Sato | Subordination and self-decomposability[END_REF] shows that L ξ is SD. Besides, setting ϕ α resp. ϕ ξ for the Lévy-Khintchine exponent of L -α resp. L ξ , one deduces from Theorem 30.4 in [START_REF] Sato | Lévy processes and infinitely divisible distributions[END_REF] the following relationship

ϕ ξ (λ) = ϕ α (λ α ).
Another, tentative, probabilistic interpretation of the self-decomposability of L ξ could be given in terms of a certain spectrally positive Markov process. Setting α = -1/ξ and y α (λ) = E[e -λL ξ ], Theorem 4.17 p. 258 in [START_REF] Kilbas | Theory and applications of fractional differential equations[END_REF] and (3.5) above show that y α is a solution to the fractional differential equation xD α+1

-y α -αy α = 0, where D α+1 -is a fractional Riemann-Liouville derivative -see Section 2.1 in [START_REF] Kilbas | Theory and applications of fractional differential equations[END_REF]. When α = 1 viz. ξ = -1 the above amounts to a Bessel equation and Feller's theory applies, making L -1 the first-passage time of a Bessel process of index 0 -see [START_REF] Kent | Some probabilistic properties of Bessel functions[END_REF]. When α ∈ (0, 1) the operator D α+1 -is the infinitesimal generator of a spectrally positive (1 + α)-stable Lévy process reflected at its minimum, which is a spectrally positive Markov process -see Section 3 in [START_REF] Patie | Intertwining certain fractional derivatives[END_REF] and the references therein. By downward continuity, one may wonder if L ξ cannot be viewed as the first-passage time of some scale-transformation of the latter, eventhough no Feller's theory is available for fractional operators whose order lies in [START_REF] Barndorff-Nielsen | First hitting time models for the generalized inverse Gaussian distribution[END_REF][START_REF] Bertoin | On subordinators, self-similar Markov processes and some factorizations of the exponential variable[END_REF].

The above probabilistic interpretations do not seem to hold for t = 1. On the one hand, Theorem 4.17 in [START_REF] Kilbas | Theory and applications of fractional differential equations[END_REF] yields then an equation with two fractional derivatives of different order for the Laplace transform of Γ ξ t . On the other hand, keeping the notation α = -1/ξ, Theorem 1 in [START_REF] Shanbhag | An extension of Goldie's result and further results in infinite divisibility[END_REF] shows the factorization

Γ ξ t = Γ -1 αt × S (t) α where S (t)
α is the size-biased sampling of S α of order -αt, viz.

E[f (S (t) α )] = E[f (S α )S -αt α ] E[S -αt α ] •
The GGC character of S (t) α follows from that of S α and Theorem 6.2.4. in [START_REF] Bondesson | Generalized Gamma convolutions and related classes of distributions and densities[END_REF], which shows by the above case ξ = -1 that Γ ξ t is the independent product of two SD random variables. By Theorem 16.1 in [START_REF] Sato | Lévy processes and infinitely divisible distributions[END_REF], this entails that there exist two independent 1-self-similar additive increasing processes Y and Z such that

Γ ξ t d = Y Z 1 .
Unfortunately, contrary to Bochner's subordination the independent composition of two additive processes is not necessarily an additive process anymore, so that the above identity does not provide a probabilistic proof of the self-decomposability of Γ ξ t .

  The self-decomposability of I is a direct consequence of the Markov property. More precisely, introducing the stopping-time T y = inf{s > 0, Z s = y} for every y > 0, the fact that Z s → +∞ a.s. as s → +∞ and the absence of positive jumps entail that T y < +∞ a.s. Decomposing, we get I = Ty 0 e -Zs ds + ∞ Ty e -Zs ds d = Ty 0 e -Zs ds + e -y ∞ 0 e -Z ′ s ds

Acknowledgements. The authors are grateful to L. Bondesson for the interest he took in this work, and several useful comments. Ce travail a bénéficié d'une aide de l'Agence Nationale de la Recherche portant la référence ANR-09-BLAN-0084-01.