Vincent Acary
email: vincent.acary@inria.fr

Energy conservation and dissipation properties of time-integration methods for the nonsmooth elastodynamics with contact

Keywords: contact dynamics, impact, computational contact mechanics, numerical time integration, energy conservation, dissipation properties Dynamique du contact, impact, mécanique numérique du contact, intégration numérique en temps, conservation d'énergie, propriétés de dissipation

This research report is devoted to the study of the conservation and the dissipation properties of the mechanical energy of several time-integration methods dedicated to the elastodynamics with unilateral contact. Given that the direct application of the standard schemes as the Newmark schemes or the generalized-α schemes leads to energy blow-up, we study two schemes dedicated to the time-integration of nonsmooth systems with contact: the Moreau-Jean scheme and the nonsmooth generalized-α scheme. The energy conservation and dissipation properties of the Moreau-Jean is firstly shown. In a second step, the nonsmooth generalized-α scheme is studied by adapting the previous works of Krenk and Høgsberg in the context of unilateral contact. Finally, the known properties of the Newmark and the Hilber-Hughes-Taylor (HHT) scheme in the unconstrained case are extended without any further assumptions to the case with contact.

Introduction and motivations

The numerical time integration of mechanical systems with unilateral contact is known to be a difficult task, mainly due to the nonsmoothness of the dynamic response when a contact occurs. For two recent reviews of the existing methods in the literature and the associated issues, we refer to [Doyen et al., 2011, Krause and[START_REF] Krause | Presentation and comparison of selected algorithms for dynamic contact based on the Newmark scheme[END_REF] and the standard textbooks [START_REF] Laursen | Computational Contact and Impact Mechanics -Fundamentals of Modeling Interfacial Phenomena in Nonlinear Finite Element Analysis[END_REF][START_REF] Wriggers | Computational Contact Mechanics[END_REF]. One of the main conclusions is that standard schemes widely used in computational contact Mechanics, such as the Newmark scheme, the Hilber-Hughes-Taylor scheme (HHT) or the generalized-α scheme cannot be directly applied to the simulation of systems with unilateral contact and impact.

In the most favorable cases, these schemes exhibit artificial oscillations of the contact velocities and forces, that blur the whole stresses in the structure. The source of these artificial oscillations is the nonsmoothness of the contact conditions that yields a jump in the velocities when a contact is closing. In finite-freedom Mechanics, when one deals with space-discretized structures after a semi-discretization, the velocity jump of a finite mass is associated with an impact. Hence, an impact law has to be specified to close the equations of the system. The low regularity of the velocity and the distributional character of the reaction prevents the use of schemes with an high order of accuracy. A partial remedy for theses problems is to use a fully implicit first-order approximation of the contact forces as it was suggested in [Jean and Moreau, 1987, Moreau, 1988a[START_REF] Carpenter | Lagrange constraints for transient finite element surface contact[END_REF][START_REF] Kane | Finite element analysis of nonsmooth contact[END_REF] and/or a treatment of the constraints at the velocity level together with an (possibly perfectly inelastic or implicitly defined) impact law [Moreau, 1988a, Laursen and[START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF]. This work has been extended to nonlinear elastodynamics by [START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impacts[END_REF]. In [START_REF] Ayyad | A frictionless viscoelastodynamic contact problem with energy consistent properties: numerical analysis and computational aspects[END_REF], the nonlinear elastodynamics with a constraint at the position level is also considered. This yields a similar numerical scheme for the energetic properties but a additional step, in the same vein as in [START_REF] Laursen | Improved implicit integrators for transient impact problems -geometric admissibility within the conserving framework[END_REF], is added to correct the constraint at the position level.

In the worst cases, the standard schemes exhibit numerical instabilities, and possibly, energy blow-ups. An alternative solution is to design numerical schemes that conserve or dissipate energy. In [START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF], the energy-conserving scheme based on the pioneering works Simo and co-workers [Simo andTarnow, 1992, Simo et al., 1995] for nonlinear elastodynamics is extended to the elasto-dynamical with unilateral contact. This work results in the use of the mid-point rule together with a velocity-level constraints and an implicit treatment of the contact forces.

The velocity-level formulation has two main advantages: it controls the dissipation of energy at contacts and it stabilizes the contact velocity. One of the main drawback is the violation of the constraints on the position level which is proportional to the time-step. In [START_REF] Laursen | Improved implicit integrators for transient impact problems -geometric admissibility within the conserving framework[END_REF], the authors propose to satisfy the constraints at the position level together the energy conservation by introducing an artificial velocity variable at the price to have oscillations of the contact velocity. In the latter case, the velocity-level constraints is not satisfied. Alternatively, the constraints at the velocity level and at the position level can be both satisfied by adding an artificial multiplier to perform a projection on the position constraint [START_REF] Acary | Projected event-capturing time-stepping schemes for nonsmooth mechanical systems with unilateral contact and Coulomb's friction[END_REF]. As remarked in [START_REF] Krause | Presentation and comparison of selected algorithms for dynamic contact based on the Newmark scheme[END_REF], there exists no algorithm satisfying the constraints at both position and velocity levels and ensuring the energy conservation simultaneously.

In the previous attempts to adapt the Newmark-based schemes in the context of computational contact mechanics, very few results are available on the energy conservation or dissipation when we deal with the unilateral contacts at the position level. In the unconstrained case, the classical Newmark scheme, with the special choice of parameters γ = 2β = 1/2, leads to an algorithm conserving the total energy of the system. With position based unilateral constraints, it is shown in [START_REF] Krause | Presentation and comparison of selected algorithms for dynamic contact based on the Newmark scheme[END_REF] that the scheme cannot conserve the energy even if γ = 2β = 1/2. In [START_REF] Khenous | Problèmes de contact unilatéral avec forttement sec de Coilomb en élastostatique et élastodynamique[END_REF], the author can only conclude to the energy dissipation of the Newmark scheme is the special case of a very dissipative order one scheme with γ = 2β = 1. With a full implicit treatment of the constraint as in [START_REF] Kane | Finite element analysis of nonsmooth contact[END_REF], as far as we know there is no general study of the energy properties for all admissible values of γ and β. When we consider the Newmark scheme with a velocity-level formulation of the constraints, the only available results are those in [START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF] that can be adapted to the special case of the Newmark scheme with a fully implicit treatment of the Lagrange multiplier. Indeed, the midpoint rule together with a velocity-level formulation is very similar to the Moreau-Jean scheme [Jean and Moreau, 1987, Moreau, 1988a, 1999[START_REF] Jean | Dynamics in the presence of unilateral contacts and dry friction: a numerical approach[END_REF] based on the θ-method when θ is equal to 1/2. This latter scheme is in turn equivalent to the Newmark scheme with γ = 2β = 1/2. Besides these special cases, the general case of the HHT scheme and generalized-α scheme are not treated from the energy properties point of view. This report attempts to bridge this gap.

Very recently, a new class of schemes has been proposed in [START_REF] Chen | A Newmark-Type Integrator for Flexible Systems Considering Nonsmooth Unilateral Constraints[END_REF][START_REF] Chen | A nonsmooth generalized-α scheme for flexible multibody systems with unilateral constraints[END_REF] which takes advantage of the Moreau-Jean scheme in terms of robustness and stability while adding some key properties of the Newmark-based schemes, that are the second-order approximation of the smooth force terms and the controlled damping of the high-frequency dynamics. This work yields the so-called nonsmooth Newmark, nonsmooth HHT and nonsmooth generalized-α schemes that deal with the contact forces trough their associated impulses in a fully implicit way, and treat the constraints at the velocity level together with Newton's impact law. In this report, the main goal is to show that the Moreau-Jean scheme and the nonsmooth schemes have the same energetic behavior as their counterparts in the unconstrained case. To this aim, the detailed list of the objectives is as follows:

• to show the energy conserving and decaying properties of the Moreau-Jean scheme.

• to provide results on the algorithmic energy conservation and dissipation of the nonsmooth Newmark scheme.

• to apply and to extend the techniques developed in [START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF]Høgsberg, 2005, Krenk, 2006] to study the dissipation properties of the nonsmooth generalized-α schemes.

• to show by means of the previous method that the nonsmooth HHT scheme dissipates a kind of algorithmic energy.

• to propose an alternative α-scheme as in [START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF]Høgsberg, 2005, Krenk, 2006] that dissipates a kind of algorithmic energy.

Inria

The report is organized as follows. Section 2 recalls the basic ingredients of the nonsmooth modeling of finite-dimensional mechanical systems subjected to unilateral contact and impact. In Section 3, the energy balance of a mechanical system with jumps in the velocity is formulated. The main schemes, studied in this report, are detailed in Section 4. The energy analysis of the Moreau-Jean scheme is done in Section 5. Section 6 starts by the presentation of the Krenk-Høgsberg method for the analysis of the discrete energy balance over a time-step for the α-schemes. After a first general result on the nonsmooth generalized-α schemes and the alternative nonsmooth Krenk-Høsberg generalized-α scheme, the nonsmooth HHT case and the nonsmooth Newmark case are fully developed. Section 7 concludes the report.

Notation

The following notation is used throughout the paper. The Euclidean norm for a vector x ∈ IR n is denoted by x . For positive definite (respectively positive semi-definite) matrix M ∈ IR n×n , x M denotes the norm (respectively the semi-norm) in the metric defined by M . Let I denote a real time interval of any sort. For a function f : I → IR n of Bounded Variation (BV), we denote the right-limit function by f + (t) = lim s→t,s>t f (s), and respectively the left-limit by f -(t) = lim s→t,s<t f (s). We denote by 0 = t 0 < t 1 < . . . < t k < . . . < t N = T a finite partition (or a subdivision) of the time interval [0, T ] (T > 0). The integer N stands for the number of time intervals in the subdivision. The length of a time step is denoted by h k = t k+1 -t k . For simplicity's sake, the schemes are presented in the sequel with a time step denoted by h for short. The value of a real function x(t) at the time t k , is approximated by x k . In the same way, the notation

x k+θ = (1 -θ)x k + θx k+1 is used for θ ∈ [0, 1].
The following notation is introduced to analyze the energetic behavior of the nonsmooth generalized-α scheme

x k,γ = γx k+1 + (1 -γ)x k , x k-1,γ = γx k + (1 -γ)x k-1 . (1)
This notation generalizes the previous notation x k+θ = (1 -θ)x k + θx k+1 to avoid the ambiguity when a multi-step method is studied. For a function f : IR → IR n and h > 0, one writes f (h) = O(h) if and only if there exist positive numbers δ and M such that f (h) M h for h < δ. The notation dt defines the Lebesgue measure on IR.

Nonsmooth mechanical systems with unilateral contact

Let us consider the equations of motion of a mechanical system subjected to unilateral constraints in the linear case:

           q(t 0 ) = q 0 , v(t 0 ) = v 0 , (2a)
q(t) = v(t), (2b)
M v(t) + Kq(t) + Cv(t) = Gλ(t), (2c)
g(q(t)) = G ⊤ q(t) + w 0, λ(t) 0, g ⊤ (q(t))λ(t) = 0, (2d) where
• q(t) ∈ IR n is the vector of generalized coordinates and v(t) = q(t) the associated vector of generalized velocities,

• the initial conditions are q 0 ∈ IR n and v 0 ∈ IR n ,

• M ∈ IR n×n is the symmetric inertia matrix, K ∈ IR n×n and C ∈ IR n×n are respectively the stiffness and the damping matrices,

• the function g(q(t)) ∈ IR m , called the gap function is used to define the unilateral constraints; with a slight abuse of notation we will also write the gap function as g(t),

• the Jacobian matrix of g with respect to q is G ⊤ = ∇ ⊤ q g(q(t)) and is assumed to be constant in the linear setting, w ∈ IR m is a constant vector,

• λ(t) ∈ IR m is the Lagrange multiplier vector associated with the constraints.

The constitutive law (2d) for the perfect unilateral constraints is also termed the Signorini condition and can be written equivalently as

0 g(t) ⊥ λ(t) 0, (3)
where the inequalities involving vectors are understood to hold component-wise and the ⊥ symbol means that y ⊤ λ = 0. Finally, let us define the following variables relative to the constraints, called local variables: the local velocity U (t) and the (local) Lagrange multiplier λ(t) which is associated with the generalized reaction forces r(t) such that

U (t) = ġ(q(t)) = G ⊤ v(t), r(t) = Gλ(t). (4)
For finite-freedom mechanical systems with unilateral constraints as in (2), it is well-known that impacts may occur if the relative velocity at contact U (t) is not compatible with the constraints. In other terms, if the contact indexed by α is closing with a negative relative velocity at t ⋆ , that is U α,-(t ⋆ ) < 0, the velocity has to jump in order to satisfy the constraints after the time of impact. However, the right velocity at the impact U α,+ (t ⋆ ) is not determined by the system (2d). This is the reason why an impact law must be added to close the system of equations. In particle or rigid body dynamics, the most simple impact law is Newton's impact law

U α,+ (t) = -e α U α,-(t), if g α (t) = 0 and U α,-(t) 0, α ∈ I (5)
where e α is the coefficient of restitution at contact α and I = {1 . . . m} is the set of indices of constraints.

Usually, the condition in (5) only involves the condition g α (t) = 0. Indeed, if the condition g α (t) 0 is always satisfied on the time interval of study, the relative pre-impact velocity U α,-(t) cannot be positive when g α (t) = 0, except at the initial time. We will see further that this is not the case when the condition (5) is only prescribed at discrete time instants. We will have therefore to adapt the condition (5) in a suitable way in Section 4 when a time-discretization is performed.

The fact that the velocity v(t) may encounter jumps yields some difficulties to define the acceleration everywhere in time. It is usual to assume that the velocity is a function of bounded variation that admits a right and left limit everywhere and which can be associated with a differential measure dv (see [Moreau, 1988b] for details). Almost everywhere with respect to the Lebesgue measure dt, the velocity is differentiable in the classical way with respect to time and we get dv dt = q(t). When a jump occurs, the standard differentiation cannot be applied since the acceleration is given by a Dirac distribution. Let us write the equation of motion in terms of differential measure in the linear case

               M dv + Kq(t) dt + Cv(t) dt = F (t) dt + G dI, q(t) = v + (t), U (t) = G ⊤ v(t), g(t) = G ⊤ q(t) + w, if g α (t) 0 and U -(t) < 0, then 0 U + (t) + eU -(t) ⊥ dI 0, ( 6
)
where di is the impulse reaction measure. Its local variant is defined by di = G dI. When the evolution is smooth, r(t) is considered as the density of di with respect to the Lebesgue measure, that is

r(t) = di dt (t), or equivalently λ(t) = dI dt (t). (7)
The last line of ( 6) defines the second-order Moreau sweeping process [Moreau, 1988a]. It can be interpreted as a reformulation of the unilateral constraint (3) at the velocity level together with an impact law. It can also be equivalently viewed as an index-reduction technique, standard in the Inria Differential Algebraic Equations (DAE) Theory when a constraint is differentiated to decrease the index of the DAE. It also contains to the so-called persistency conditions, defined in [START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF] when the contact is closed, but adds the Newton impact law when a new contact occurs.

It can also be shown that the system (6) contains in a single formulation the nonimpulsive dynamics and the impact dynamics. If we omit the singular part in the decomposition of the measure, we can split the measure dv and di as follows splitting of measures

dv = q(t) dt + i (v + (t i ) -v -(t i ))δ ti , dI = λ dt + i P i δ ti . ( 8
)
Inserting this decomposition in the measure differential system (6), we get the standard nonimpulsive equation of motion (2) almost everywhere and the impact equations at the time of impacts:

         M (v + (t i ) -v -(t i )) = G P i , U ± (t i ) = G ⊤ v ± (t i ), g(t) = G ⊤ q(t) + w, if U -(t i ) < 0, then 0 U + (t i ) + eU -(t i ) ⊥ P i 0. (9)
For more details on the modeling of multibody systems with unilateral constraints, we refer to [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics[END_REF][START_REF] Pfeiffer | Multibody Dynamics with Unilateral Contacts[END_REF], Moreau, 1988a] and for the mathematical analysis, we refer to [START_REF] Schatzman | A class of nonlinear differential equations of second order in time[END_REF][START_REF] Monteiro Marques | Differential Inclusions in Nonsmooth Mechanical Problems. Shocks and Dry Friction[END_REF][START_REF] Stewart | Rigid body dynamics with friction and impact[END_REF], Ballard, 2000].

Energy balance analysis

In the case on nonsmooth motion with impact, we recall that the equation of motion in terms of measures are given by

M dv + (Kq(t) + Cv(t)) dt = F dt + di, dq = v(t)dt. (10)
A detailed analysis of the energy balance for nonsmooth systems can be found in [START_REF] Brogliato | Dissipative Systems Analysis and Control[END_REF] and [START_REF] Leine | Stability and Convergence of Mechanical Systems with Unilateral Constraints[END_REF]. Let us recall in this paper the basic formulae. The energy balance can be obtained by multiplying the equation of motion by

v + + v - (v + + v -) ⊤ M dv + (v + + v -) ⊤ (Kq + Cv) dt = (v + + v -) ⊤ F dt + (v + + v -) ⊤ di, (11)
that is equivalently

d(v ⊤ M v) + (v + + v -)(Kq + Cv) dt = (v + + v -)F dt + (v + + v -) di. ( 12
)
With the standard definition of the total mechanical energy of the system,

E := 1 2 v ⊤ M v + 1 2 q ⊤ Kq, (13)
we get the following energy balance

2dE := d(v ⊤ M v) + 2q ⊤ Kdq = 2v ⊤ F dt -2v ⊤ Cv dt + (v + + v -) ⊤ di. ( 14
)
If we split the differential measure in di = r(t) dt + i p i δ ti , we get

2dE := d(v ⊤ M v) + 2q ⊤ Kdq = 2v ⊤ (F + r) dt -2v ⊤ Cv dt + i (v + + v -) ⊤ p i δ ti . (15)
By integration over a time interval [t 0 , t 0 ] such that t i ∈ (t 0 , t 1 ), we obtain an energy balance equation (also named the dissipation equality [START_REF] Brogliato | Dissipative Systems Analysis and Control[END_REF]) as

∆E = E(t 1 ) -E(t 0 ) = t1 t0 v ⊤ (F + r) dt - t1 t0 v ⊤ Cv dt + i 1 2 (v + (t i ) + v -(t i )) ⊤ p i . ( 16
)
The right hand side of the energy balance equation represents the work done in the time interval [t 0 , t 1 ] that can be decomposed as follows:

• the term

W ext = t1 t0 v ⊤ F dt, (17)
is the work done by the external forces,

• the term

W damping = - t1 t0 v ⊤ Cv dt, (18)
is the work done by the damping internal forces. If the damping matrix C is a positive semi-definite matrix, we then conclude on the sign of the work, W damping 0.

• the term

W contact = t1 t0 v ⊤ r dt, (19)
is the work done by the contact forces. If we consider perfect unilateral constraints, we have

v ⊤ r = v ⊤ Gλ = U ⊤ λ = 0 (20)
and then W contact = 0,

• the term

W impact = 1 2 (v + (t i ) + v -(t i )) ⊤ p, (21)
represents the work done by the contact impulse at the time of impact t i . Since p = GP and U = G ⊤ v , using the Newton impact law, we get

W impact = 1 2 (v + (t i ) + v -(t i )) ⊤ p = 1 2 (v + (t i ) + v -(t i )) ⊤ GP = 1 2 (U + (t i ) + U -(t i )) ⊤ P = 1 2 ((1 -e)U -(t i )) ⊤ P 0 for 0 e 1. (22)
The formula (21) of the work done by the impulse during an impact is also known as Thomson and Tait's formula [Brogliato, 1999, Section 4.2.12]. The fact that the work has to be negative such that an impact dissipates some energy is also related the Clausius-Duhem inequality applied to an impulse motion [START_REF] Frémond | Non-Smooth Thermo-mechanics[END_REF].

Background on the time-integration methods

Leaving aside the time-integration methods based on an accurate event detection procedure (event-tracking schemes or event-driven schemes [Acary and Brogliato, 2008, Chap. 8]), the time-integration of nonsmooth mechanical systems is performed by means of event-capturing time-stepping schemes. In these schemes, the impact are not accurately located but captured by the refinement of the time-step size. One of the most well-proven method is the Moreau-Jean scheme [Jean and Moreau, 1987, Moreau, 1988a[START_REF] Jean | Dynamics in the presence of unilateral contacts and dry friction: a numerical approach[END_REF]. This method which is sound from Inria the mathematical analysis point of view (convergence analysis can be found in [START_REF] Monteiro Marques | Differential Inclusions in Nonsmooth Mechanical Problems. Shocks and Dry Friction[END_REF][START_REF] Dzonou | A sweeping process approach to inelastic contact problems with general inertia operators[END_REF][START_REF] Dzonou | A convergence result for a vibro-impact problem with a general inertia operator[END_REF]) and which takes advantage of a strong practical experience, is of low order of accuracy, but very robust in many practical applications. Very recently, an attempt has been made in [START_REF] Chen | A Newmark-Type Integrator for Flexible Systems Considering Nonsmooth Unilateral Constraints[END_REF][START_REF] Chen | A nonsmooth generalized-α scheme for flexible multibody systems with unilateral constraints[END_REF] to improve its accuracy at least on the smooth terms of the equations of motion. This work has led to the nonsmooth Newmark and generalized-α scheme as it extends the standard schemes for computational mechanics to the nonsmooth dynamical case. Both schemes will be briefly presented in the following sections.

Moreau-Jean's scheme

The Moreau-Jean scheme [Jean and Moreau, 1987, Moreau, 1988a[START_REF] Jean | Dynamics in the presence of unilateral contacts and dry friction: a numerical approach[END_REF] is based on the Moreau sweeping process (6). The numerical time integration of ( 6) is performed on an interval (t k , t k+1 ] of length h as follows

                 M (v k+1 -v k ) + hKq k+θ + hCv k+θ -hF k+θ = p k+1 = GP k+1 , (23a)
q k+1 = q k + hv k+θ , (23b)
U k+1 = G ⊤ v k+1 , ( 23c
) 0 U α k+1 + eU α k ⊥ P α k+1 0, α ∈ I 1 , P α k+1 = 0, α ∈ I \ I 1 , (23d)
with θ ∈ [0, 1]. The following approximations are considered:

v k+1 ≈ v + (t k+1 ); U k+1 ≈ U + (t k+1 ); p k+1 ≈ di(]t k , t k+1 ]), P k+1 ≈ dI(]t k , t k+1 ]). ( 24
)
Note that the unknown variable P k+1 is equivalent to an impulse. This characteristic feature of the Moreau-Jean scheme renders the numerical integration consistent when an impact occurs. Indeed, when the time-step vanishes, a choice of a variable equivalent to a contact force would lead to some unbounded values.

The index set I 1 results from the time-discretization of the conditional statement in (6) if g α (t) 0 and U -(t) < 0 (25) that allows us to apply the Signorini condition at the velocity level. In the numerical practice, we choose to define this set by

I 1 = {α ∈ I | G ⊤ (q k + hv k ) + w 0 and U k 0}. ( 26
)
Other strategies for defining I 1 can be found in [START_REF] Acary | Projected event-capturing time-stepping schemes for nonsmooth mechanical systems with unilateral contact and Coulomb's friction[END_REF]. The numerical scheme which solves (6) enforces in discrete time the Newton law at each time step. Conversely, the constraints g(t) 0 are not satisfied. A violation of the constraints, proportional to the time-step, can occur at the activation of the contact, that is, when the contact is closing. The violation may be corrected by a projection technique onto the constraints extending the Gear-Gupta-Leimkuhler approach for DAEs as it has recently been done in [START_REF] Acary | Projected event-capturing time-stepping schemes for nonsmooth mechanical systems with unilateral contact and Coulomb's friction[END_REF].

Nonsmooth Newmark and generalized-α scheme

In [START_REF] Chen | A Newmark-Type Integrator for Flexible Systems Considering Nonsmooth Unilateral Constraints[END_REF][START_REF] Chen | A nonsmooth generalized-α scheme for flexible multibody systems with unilateral constraints[END_REF], a new family of time-integration schemes has been developed that keep the advantages of the Moreau-Jean in terms of robustness and efficiency while adding the some key properties of the Newmark [Newmark, 1959], the Hilber-Hughes-Taylor (HHT) [START_REF] Hughes | The Finite Element Method, Linear Static and Dynamic Finite Element Analysis[END_REF] and the generalized-α [START_REF] Chung | A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-α method[END_REF] schemes. The most well-known property of the latter schemes is the controlled damping of the high-frequency contents of the dynamics.

Acary

In the linear time invariant dynamics with unilateral constraints, the new scheme can be written as follows

                                                       M vk+1 + Kq k+1 + Cv k+1 = F k+1 , (27a)
M w k+1 = GP k+1 , (27b)
U k+1 = G ⊤ v k+1 , (27c)
0 U α k+1 + eU α k ⊥ P α k+1 0, α ∈ I 1 , P α k+1 = 0, α ∈ I \ I 1 , (27d)
with

(1 -α m )a k+1 + α m a k = (1 -α f ) vk+1 + α f vk , (27e)
qk+1 = q k + hv k + h 2 (1/2 -β) a k + h 2 βa k+1 , (27f)
ṽk+1 = v k + h(1 -γ)a k + hγa k+1 , (27g)
v k+1 = ṽk+1 + w k+1 , (27h)
q k+1 = qk+1 + h 2 w k+1 . (27i)
The numerical scheme ( 27) has been designed such that it deals, as the Moreau-Jean scheme with the contact forces and impact through their associated impulses P k+1 in a fully implicit way. In this manner, we ensure that the scheme will be consistent when the time-step vanishes if an impact occurs. Furthermore, it also includes a treatment of the unilateral constraint together with the Newton-impact law at the velocity level. This aspect is crucial for the conservation and dissipation properties as we will see in Section 6. Finally, the last important property is the second order approximation of the smooth terms given by the generalized-α schemes that allows us to take advantage of the controlled damping of the high-frequency dynamics. It is indeed wellknown that one of the main advantages of the α-scheme with respect to the Newmark scheme is the possibility to introduce some damping of the high-frequency dynamics without altering the order. In the smooth case when the contact is not taken into account, the condition of second order accuracy reads as

γ = 1 2 + α f -α m . (28)
The optimal parameters are usually chosen according to the spectral radius at infinity

ρ ∞ ∈ [0, 1] such that α m = 2ρ ∞ -1 ρ ∞ + 1 , α f = ρ ∞ ρ ∞ + 1 , β = 1 4 (γ + 1 2 ) 2 . ( 29
)
The nonsmooth Newmark algorithm is obtained with α m = α f = 0 and the nonsmooth HHT scheme in the form published in [START_REF] Hughes | The Finite Element Method, Linear Static and Dynamic Finite Element Analysis[END_REF] is obtained with α m = 0 and

α f ∈ [0, 1/3].
Variant of the Moreau-Jean scheme for 2β = γ = θ. Let us remark that the nonsmooth Newmark algorithm can be reformulated as

   M (v k+1 -v k ) = h(F k+γ -Kq k+γ -Cv k+γ ) + GP k+1 (30a)
q k+1 = q k + hv k + 1 2 hM -1 [h(F k+2β -Kq k+2β -Cv k+2β ) + GP k+1 )] (30b)
with (27d). With straightforward manipulations, one obtains if

2β = γ = θ                  M (v k+1 -v k ) + hKq k+θ + hCv k+θ -hF k+θ = p k+1 = GP k+1 , (31a)
q k+1 = q k + hv k+1/2 , (31b)
U k+1 = G ⊤ v k+1 , (31c)
0 U α k+1 + eU α k ⊥ P α k+1 0, α ∈ I 1 , P α k+1 = 0, α ∈ I \ I 1 . (31d) Inria
This scheme appears as a variant of the Moreau-Jean scheme. We will see in Section 6.4 that these scheme has better dissipation properties than the original one. Note that the original Moreau-Jean scheme cannot be viewed as a special case of the nonsmooth Newmark scheme when θ = 1/2.

4.3 Nonsmooth Krenk-Høgsberg (KH) generalized-α scheme

In [START_REF] Krenk | Properties of time integration with first order filter damping[END_REF]] and [START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF], an alternative collocation method is proposed for the generalized-α scheme

(1 -α m )[M a k+1 + Cv k+1 -F k+1 ] + α m [M a k + Cv k + F k ] = (1 -α f )[-Kq k+1 ] + α f [-Kq k ].
(32) This scheme, that will be termed in the sequel as the Krenk-Høgsberg (KH) generalized-α scheme, differs only from the original generalized-α scheme by the fact that the damping terms and the load terms have the same weight as the inertial term. In the original generalized-α scheme as it is given in (61), the weighting of the damping and the load terms follow the stiffness term. The nonsmooth KH generalized-α scheme is obtained by replacing (27d) by ( 32) in ( 27). Although the order of the method decreases, we will see in the sequel that this scheme has better dissipation properties than the standard generalized-α scheme.

The analysis of the KH scheme in [START_REF] Krenk | Properties of time integration with first order filter damping[END_REF] shows that the scheme introduces a slight increase of the frequency response of the mechanical system. In other terms, the original generalized-α scheme has a slightly improved frequency response. This is mainly related to the alternative choice of weighting of the structural damping term. On the contrary, the choice of weighting of the load term in the KH scheme is superior.

Most of the original α-schemes are contained in the KH generalized-α scheme. Note that in the original paper on the HHT scheme [START_REF] Hilber | Improved numerical dissipation for the time integration algorithms in structural dynamics[END_REF], the weighting of the load term follows the inertial term. This is the reason why the original HHT scheme in [START_REF] Hilber | Improved numerical dissipation for the time integration algorithms in structural dynamics[END_REF] can be obtained from (32) with α m = 0 and α f = α. In the same manner, the α-method of [START_REF] Wood | An alpha modification of Newmark's method[END_REF] can be obtained from (56) by choosing α f = 0.

Energy conservation and dissipation properties of Moreau-Jean scheme

In this section, we give a first result on the energy conservation and dissipation of the Moreau-Jean scheme. This result gives a criteria on the parameter θ which depends on the coefficient of restitution. Let us define the discrete approximation of the work done by the external forces within the step by

W ext k+1 = hv ⊤ k+θ F k+θ ≈ t k+1 t k F v dt, (33)
and the discrete approximation of the work done by the damping term by

W damping k+1 = -hv ⊤ k+θ Cv k+θ ≈ - t k+1 t k v T Cv dt. ( 34
)
We have the following estimate for the variation of the total mechanical energy.

Lemma 1 The discrete-time dissipation equality of the Moreau-Jean scheme (23) over a timestep [t k , t k+1 ] is given by

∆E -W ext k+1 -W damping k+1 = ( 1 2 -θ) v k+1 -v k 2 M + (q k+1 -q k ) 2 K + U ⊤ k+θ P k+1 . ( 35
) RR n°8602

Acary

Proof Let us first compute the variation of the energy E over the time-step

∆E = E(q k+1 , v k+1 ) -E(q k , v k ) = 1 2 (v k+1 + v k ) ⊤ M (v k+1 -v k ) + (q k+1 + q k ) ⊤ K(q k+1 -q k ) , (36)
since we assume that M = M ⊤ and K = K ⊤ . Let us also remark that for the θ-method, we have

1 2 (v k+1 + v k ) = 1 h (q k+1 -q k ) + ( 1 2 -θ)(v k+1 -v k ), (37)
and

v k+θ = 1 2 (v k+1 + v k ) -( 1 2 -θ)(v k+1 -v k ). (38)
Using (37) and then (23a), the energy balance (36) becomes

∆E = ( 1 2 -θ)(v k+1 -v k ) ⊤ M (v k+1 -v k ) + 1 h (q k+1 -q k )M (v k+1 -v k ) + 1 2 (q k+1 + q k ) ⊤ K(q k+1 -q k ) = ( 1 2 -θ)(v k+1 -v k ) ⊤ M (v k+1 -v k ) + 1 h (q k+1 -q k ) [-hKq k+θ -hCv k+θ + hF k+θ + GP k+1 ] + 1 2 (q k+1 + q k ) ⊤ K(q k+1 -q k ) = ( 1 2 -θ)(v k+1 -v k ) ⊤ M (v k+1 -v k ) + ( 1 2 -θ)(q k+1 -q k ) ⊤ K(q k+1 -q k ) -hv ⊤ k+θ Cv k+θ + hv ⊤ k+θ F k+θ + v ⊤ k+θ GP k+1 .
(39) Using the expression of the norms • M and the semi-norm • K , the expression (39) can be easily simplified in

∆E = ( 1 2 -θ) v k+1 -v k 2 M + (q k+1 -q k ) 2 K + hv ⊤ k+θ F k+θ -hv ⊤ k+θ Cv k+θ + v ⊤ k+θ GP k+1 .
(40) Using the definition of the discrete works in (33) and ( 34) and the fact that v ⊤ k+θ GP k+1 = U ⊤ k+θ P k+1 yields the result.

Remark 1 The previous result may be specified at the order h. Using the following approximation for function of bounded variations [START_REF] Acary | Higher order event capturing time-stepping schemes for nonsmooth multibody systems with unilateral constraints and impacts[END_REF] hv ⊤ k+θ F k+θ -

t k+1 t k F (t)v(t) dt = O(h), (41)
and

hv ⊤ k+θ Cv k+θ - t k+1 t k v ⊤ (t)Cv(t) dt = O(h), (42)
we get for h small enough ∆E -

t k+1 t k F v dt + t k+1 t k Cv dt = ( 1 2 -θ) v k+1 -v k 2 M + (q k+1 -q k ) 2 K + U ⊤ k+θ P k+1 + O(h). (43)
Let us give now a first result concerning the dissipation of the Moreau-Jean scheme

Proposition 1 The Moreau-Jean scheme dissipates energy in the sense that

E(t k+1 ) -E(t k ) W ext k+1 + W damping k+1 , ( 44
) Inria if 1 2 θ 1 1 + e 1. ( 45
)
where e = max e α , α ∈ I. In particular, for e = 0, we get 1 2 θ 1 and for e = 1, we get θ = 1 2 .

In other words, providing that ( 45) is satisfied, the variation of the total mechanical energy of the system is always less than the energy supplied by the external and damping forces.

Proof Obviously, we have

( 1 2 -θ) v k+1 -v k 2 M + (q k+1 -q k ) 2 K 0, if and only if θ 1/2. ( 46
)
I remains to prove that U ⊤ k+θ P k+1 0. Let us define the following index set of contacts

I 0 1 = {α ∈ I 1 | P α k+1 = 0}, (47)
and its complement in I 1 is denoted by I 0 1 = I 1 \ I 0 1 . We can therefore develop U ⊤ k+θ P k+1 as

U ⊤ k+θ P k+1 = α∈I U α k+θ P α k+1 = α∈I1 U α k+θ P α k+1 since for α ∈ I 1 , P α k+1 = 0 = α∈I 0 1 (1 -θ(1 + e))U α k P α k+1 (48)
since for α ∈ I 0 1 , P α k+1 = 0 and α ∈ I 0 1 , U α k+1 = -eU α k . Since we have P α k+1 0 and U α k 0 for all α ∈ I 1 , the constraint on θ is therefore

θ 1 1 + e α
1, for all α ∈ I.

By combining the constraints on θ in ( 46) and ( 49), the result is proved.

The following comments can be made on Proposition 1:

1. The variation of energy (35) may be formulated in another form as

∆E -W ext k+1 -W damping k+1 = ( 1 2 -θ) v k+1 -v k ) 2 M + q k+1 -q k 2 K + 1 2 (v k+1 + v k ) ⊤ GP k+1 -( 1 2 -θ)(v k+1 -v k ) ⊤ GP k+1 = ( 1 2 -θ) v k+1 -v k ) 2 M + q k+1 -q k 2 K -P ⊤ k+1 (U k+1 -U k ) +P ⊤ k+1 U k+1/2 . ( 50
)
The term P ⊤ k+1 U k+1/2 appears as the discrete dissipated energy at impact. We have also

-P ⊤ k+1 (U k+1 -U k ) = α∈I1 (1 + e)P α k+1 U α k 0. ( 51
)
This alternative form (50) shows that the scheme is always dissipative for θ = 1 2 . 2. The bounds (45) are not sharp since a part of the energy potentially generated at impact

(1 -θ(1 + e))U ⊤ k P k+1 is dissipated by the inertial term ( 1 2 -θ) v k+1 -v k ) 2 M .
3. For e = 0, the scheme is dissipative for the whole range θ ∈ [ 1 2 , 1]. We can also observe that the trend is somehow opposed the standard property of the θ-method. If the system is more dissipative at contact when e α = 1, we have to use the most conservative case with θ = 1/2.

Dissipation properties of nonsmooth generalized-α scheme

In this section, the behavior of the nonsmooth generalized-α scheme concerning the energyconserving or dissipating properties is studied. The proposed method of study is an extension of the pioneering works of [START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF][START_REF] Krenk | Properties of time integration with first order filter damping[END_REF]Høgsberg, 2005, Krenk, 2006] on the α-schemes.

Principle of the Krenk-Høgsberg method and its extension

One of the fundamental properties of the generalized-α scheme is the introduction of a controllable damping of the high-frequency dynamics without altering the second-order accuracy. In [START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF]Høgsberg, 2005, Krenk, 2006], this property is studied by explicitly exhibiting and adding a first-order filter and an associated additional state variable in the time-continuous dynamics. Once the augmented time-continuous dynamics is defined, it is shown that the generalized-α scheme applied to the original dynamics is equivalent to the application of the standard Newmark scheme to the augmented dynamics over two consecutive time-steps and performing a weighting procedure. Hence, the study that have been done for the Newmark scheme in the previous section can be adapted to the generalized-α scheme.

The original Krenk-Høgsberg method in [START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF]Høgsberg, 2005, Krenk, 2006] is based on the introduction one additional filter and one additional variable to the original dynamics. The augmented dynamics M a(t)

+ Cv(t) + Kq(t) = F (t) + η ν [Kz(t)], (52)
is introduced, together with the following auxiliary dynamics that filter the previous one

νh ż(t) + z(t) = νh q(t), (53)
where the time scale of the filter is given by νh. The parameter η is a non-dimensional parameter that permits to tune the effect of the filter on the original dynamics. The dynamics of the first order filter in ( 53) is discretized by means of a mid-point rule:

ν(z k+1 -z k ) + 1 2 (z k+1 + z k ) = ν(q k+1 -q k ), (54)
Rearranging the terms, we equivalently write (54) as

( 1 2 + ν)z k+1 + ( 1 2 -ν)z k = ν(q k+1 -q k ). ( 55
)
Let us consider now a linear combination of the augmented equation of motion ( 52) with the weight (1/2 + ν) at time t k+1 and the weight (1/2 -ν) at time t k :

( 1 2 + ν)[M a k+1 + Cv k+1 -F k+1 ] + ( 1 2 -ν)[M a k + Cv k -F k ] = ( 1 2 + ν -η)[-Kq k+1 ] + ( 1 2 -ν + η)[-Kq k ]. (56)
By choosing the values of ν and η such that

ν = 1 2 -α m η = ν - 1 2 + α f = α f -α m , (57)
we obtain the KH generalized-α scheme as in (32). The standard energetic analysis of the Newmark scheme can be then extended to the KH generalized-α scheme by adding the following damping force in the energetic analysis of the Newmark scheme

f A = η ν [Kz] (58)
Inria as it has been shown in [START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF]. This result will be extended to the nonsmooth case in Section 6.2. Unfortunately, the previous approach does not longer directly apply to the study of the standard generalized-α scheme since the weighting of the damping and the load terms follow the stiffness term in the method presented by [START_REF] Chung | A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-α method[END_REF]. In the following, we use three additional filters and three associated variables x, y, z to the original dynamics. Let us introduce the augmented dynamics

M a(t) + Cv(t) + Kq(t) = F (t) + η ν [Kz(t) + Cx(t) -y(t)], (59)
and the following auxiliary dynamics that filter the previous one

νh ż(t) + z(t) = νh q(t), νh ẋ(t) + x(t) = νh v(t), νh ẏ(t) + y(t) = νh Ḟ (t). ( 60
)
As previously, let us consider now a linear combination of the augmented equation of motion ( 59) with the weight (1/2 + ν) at time t k+1 and the weight (1/2 -ν) at time t k . By choosing the values of ν and η as in ( 57), the following discretization is obtained

(1 -α m )M a k+1 + α m M a k = (1 -α f )[-Kq k+1 -Cv k+1 + F k+1 ] + α f [-Kq k -Cv k + F k ]. (61)
The relation ( 61) is the characteristic relation that defines the generalized-α scheme. Contrary to the work in [START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF], the physical meaning of the filters are more difficult to justify, but it enables to retrieve the second order accurate generalized-α scheme developed in [START_REF] Chung | A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-α method[END_REF].

Remark 2 In the sequel, we will assume the ratio η/ν is finite which is not necessarily the case if ν = 0. For instance, if ρ ∞ = 1, we get from ( 29)

α m = 1 2 , α f = 1 2 , (62)
that yields ν = 0 and η = 0.

We known that the case ρ ∞ = 1 corresponds to the case with the minimal damping. The filters whose time-scale vanishes do not act as a filter since we get from (60) that z(t) = 0. However, from (29), we obtain also

η ν = 1/2 -α m α f -α m = 1 2 - 2ρ ∞ -1 ρ ∞ + 1 ρ ∞ ρ ∞ + 1 - 2ρ ∞ -1 ρ ∞ + 1 = 2 3 , (64)
so that the subsequent analysis remains valid.

The standard energetic analysis of the Newmark scheme can be then extended to the generalizedα scheme by adding the following damping force in the energetic analysis of the Newmark scheme

f A = η ν [Kz + Cx -y]. (65)
Let us define a discrete "algorithmic energy" of the form

H(q, v, a, z) = E(q, v) + h 2 4 (2β -γ)a ⊤ M a + η 2ν 2 (ν -(γ - 1 2 ))z ⊤ Kz. ( 66
) RR n°8602

Acary

Let us define the discrete approximation of the work done by the external forces within the step by

W ext k+1 = (q k+1 -q k ) ⊤ F k,γ ≈ t k+1 t k F v dt (67)
and the discrete approximation of the work done by the damping term by

W damping k+1 = -(q k+1 -q k ) ⊤ Cv k,γ ≈ - t k+1 t k v T Cv dt. (68)
The following result can be obtained.

Lemma 2 The variation of the "algorithmic" energy ∆H over a time-step performed by the nonsmooth generalized-α scheme ( 27) is

∆H -W ext k+1 -W damping k+1 + (q k+1 -q k ) ⊤ η ν [y k+γ -Cx k+γ ] = U ⊤ k+1/2 P k+1 + 1 2 h 2 ( 1 2 -γ)(2β -γ) (a k+1 -a k ) 2 M +(η + 1 2 -γ) q k+1 -q k 2 K + η ν (γ -ν - 1 2 ) z k+1 -z k 2 K . (69)
The proof of this Lemma is given in Appendix A.

Let us remark that the analysis of the dissipation properties of the nonsmooth generalized-α scheme only differs from the non-impulsive case by the term U ⊤ k+1/2 P k+1 . This is mainly the result of the design of the nonsmooth generalized-α scheme which deals with the nonsmooth terms with a low-order approximation scheme. By the way, a direct analysis of the nonsmooth Newmark scheme can be carried out by directly extending the work of [START_REF] Hughes | The Finite Element Method, Linear Static and Dynamic Finite Element Analysis[END_REF]. For the sake of space, the nonsmooth Newmark scheme will be treated as a special case in Section 6.4. Let us give a result of the sign of U ⊤ k+1/2 P k+1 that appears as the additional term due to the nonsmooth terms in the dynamics.

Lemma 3 Let us consider that the local velocities U k+1 and impulses P k+1 satisfies (27d). Then the discrete work of the contact forces is negative

U ⊤ k+1/2 P k+1 0 (70)
Proof By introducing the sets of indices as in (47), we have

U ⊤ k+1/2 P k+1 = α∈I U α k+1/2 P α k+1 = α∈I1 U α k+1/2 P α k+1 since for α ∈ I 1 , P α k+1 = 0 = α∈I 0 1 1 2 (1 -e)U α k P α k+1 (71)
since for α ∈ I 0 1 , P α k+1 = 0 and α ∈ I 0 1 , U α k+1 = -eU α k . We conclude that U ⊤ k+1/2 P k+1 0 since P α k+1 0 and U α k 0 for all α ∈ I 1 . Lemma 2 and Lemma 3 do not permit to conclude in the general case to the dissipation of the scheme. This is mainly due to the presence of the terms related to y k+γ and x k+γ in the left-hand side of (69). Although these terms are only related to the external forces and the damping terms, it seems difficult to expressed them in terms of the original dynamical system. One can only conclude on special cases when the external forces are constant and the damping matrix vanishes. In the next sections, we prefer to focus our effort on the nonsmooth KH generalized-α scheme in Section 6.2, on the HHT scheme in Section 6.3 and on the nonsmooth Newmark scheme in Section 6.4. Inria 6.2 The nonsmooth KH generalized-α scheme case.

The following proposition is a direct application of Lemma 2 and Lemma 3 Proposition 2 The variation of the "algorithmic" energy ∆H over a time-step performed by the nonsmooth KH generalized-α scheme (32) is

∆H -W ext k+1 -W damping k+1 = U ⊤ k+1/2 P k+1 - 1 2 h 2 (γ - 1 2 )(2β -γ) (a k+1 -a k ) 2 M -(γ - 1 2 -η) q k+1 -q k 2 K - η ν (ν -γ + 1 2 ) z k+1 -z k 2 K . (72)
Moreover, the nonsmooth KH generalized α-scheme is stable in the following sense

∆H -W ext k+1 -W damping k+1 0, ( 73
) if 2β γ 1 2 and 0 η γ - 1 2 ν. ( 74
)
In terms of α m and α f , the condition ( 74) is equivalent to

2β γ 1 2 and 0 α f -α m γ - 1 2 1 2 -α m . (75)
Proof: The proof of the equation ( 72) follows exactly the same lines as the proof of Lemma 2 by cancelling the term y k+γ and x k+γ . The inequality ( 73) is directly obtained with the conditions (74) and the fact that U ⊤ k+1/2 P k+1 0 comes from Lemma 3. The equivalent form of the condition in ( 75) is obtained with the help of (57).

Note that with the second order accuracy condition (28) γ = 1/2 + α f -α m , the condition (75) simplifies in

2β γ 1 2 and 0 α f -α m 1 2 -α m . (76)
The nonsmooth KH generalized-α scheme appears as an interesting alternative for the computation of the linear elastodynamics of mechanical system with unilateral contact and impact.

The nonsmooth HHT case

With the special choice α m = 0, the generalized-α scheme reduces to the HHT scheme in the form presented in [START_REF] Hughes | The Finite Element Method, Linear Static and Dynamic Finite Element Analysis[END_REF]. The equivalent filter parameters are ν = 1/2 and η = α f for the HHT scheme. For the sake of simplicity, the parameter α f will be denoted as α := α f . The HHT scheme is given by

M a k+1 + (1 -α)[Kq k+1 + Cv k+1 ] + α[Kq k + Cv k ] = (1 -α)F k+1 + αF k . ( 77
)
The application of Lemma 2 in this context leads to the following definition of the approximation of works as follows:

W ext k+1 = (q k+1 -q k ) ⊤ [(1 -α)F k,γ + αF k-1,γ ] ≈ t k+1 t k F v dt (78)
and

W damping k+1 = -(q k+1 -q k ) ⊤ C [(1 -α)v k,γ + αv k-1,γ ] ≈ - t k+1 t k v T Cv dt. ( 79
)
The following result is straightforwardly derived as a consequence of Lemma 2.

RR n°8602

Acary Proposition 3 The variation of the "algorithmic" energy ∆H over a time-step performed by the nonsmooth HHT scheme (scheme ( 27) with α m = 0 and

α f = α) is ∆H -W ext k+1 -W damping k+1 = U ⊤ k+1/2 P k+1 - 1 2 h 2 (γ - 1 2 )(2β -γ) (a k+1 -a k ) 2 M -(γ - 1 2 -α) q k+1 -q k 2 K -2α(1 -γ) z k+1 -z k 2 K . (80)
Moreover, the nonsmooth HHT scheme dissipates the "algorithmic" energy H in the following sense

∆H -W ext k+1 -W damping k+1 0, ( 81
) if 2β γ 1 2 and 0 α γ - 1 2 1 2 . ( 82
)
Proof: The right-hand side of ( 80) is directly obtained from Lemma 2 by writing ν = 1/2 and η = α. We have still to simplify the term related to the external and damping forces. With ν = 1/2, the discretization of the filters amounts to solving

1 2 (x k+1 -x k ) + 1 2 (x k+1 + x k ) = 1 2 (v k+1 -v k ), 1 2 (y k+1 -y k ) + 1 2 (y k+1 + y k ) = 1 2 (F k+1 -F k ). (83)
Simplification yields

2x k+1 = v k+1 -v k , 2y k+1 = F k+1 -F k . (84)
Therefore, we obtained for the additional terms

F k,γ - η ν y k,γ = F k,γ -2αy k,γ = γF k+1 -(1 -γ)F k -α [γ(F k+1 -F k ) + (1 -γ)(F k -F k-1 )] = (1 -α) [γF k+1 + (1 -γ)F k ] + α [γF k + (1 -γ)F k-1 ] = (1 -α)F k,γ + αF k-1,γ . (85)
By applying the same manipulations to the term involving Cx k+γ the result is obtained. The inequality ( 81) is straightforwardly obtained thanks to the positiveness of the quadratic terms when the conditions (82) are applied to (80). Since the remaining term U ⊤ k+1/2 P k+1 from Lemma 3, the proof is completed.

Note that with the second order accuracy condition (28) γ = 1/2 + α, the condition simplifies in

2β γ 1 2 and 0 α 1 2 . ( 86
)

The nonsmooth Newmark case

With the Newmark scheme (α m = α f = 0), the value of the parameters are ν = 0, η = 1/2. The algorithmic energy reduces to

H(q, v, a) = E(q, v) + h 2 4 (2β -γ)a ⊤ M a. ( 87
)
Although there is no direct mechanical interpretation of this quantity, it allows one to conclude on the dissipation property of the scheme since K is a semi-norm for 2β γ. Let us remark that we retrieve the algorithmic energy introduced by [START_REF] Hughes | A note on the stability of Newmark's algorihtm in non-linear structural dynamics[END_REF]. The following result can be obtained.

Inria

Proposition 4 The variation of the "algorithmic" energy ∆H over a time-step performed by the nonsmooth Newmark scheme (30) is

∆H -W ext k+1 -W damping k+1 = ( 1 2 -γ) q k+1 -q k 2 K + h 2 (2β -γ) (a k+1 -a k ) 2 M + U ⊤ k+1/2 P k+1 .
(88) Moreover, the nonsmooth Newmark scheme dissipates the "algorithmic" energy H in the following sense

∆H -W ext k+1 -W damping k+1 0, (89)
for 2β γ 1 2 . ( 90
)
Proof: The relation ( 88) is direct application of Lemma 2 with ν = 0, η = 1/2. The inequality ( 89) is straightforward to obtain. For 2β γ 1 2 , the first of the right-hand side of ( 88) is non positive and the fact that U ⊤ k+1/2 P k+1 0 comes from Lemma 3.

The very interesting fact in the nonsmooth Newmark scheme is that the standard dissipation properties of the original Newmark scheme are conserved when unilateral contacts and impacts are included. It is well-known that considering the unilateral contact at the position level and a semi-implicit rule for the reaction forces leads to blow-up in energy [START_REF] Chawla | Energy consistent algorithms for frictional contact problem[END_REF]] and an ad-hoc restitution rule at contact (see Figure 4 in [START_REF] Chen | A nonsmooth generalized-α scheme for flexible multibody systems with unilateral constraints[END_REF]). With the nonsmooth Newmark scheme, the energy balance is equivalent to the unconstrained case. The following result is a corollary of Proposition 4 that specifies the properties of the nonsmooth Newmark scheme for some particular choices of the parameters γ and β.

Corollary 1 For the following specific values of the parameters γ and β, the variation of the "algorithmic" energy ∆H over a time-step performed by the nonsmooth Newmark scheme (30) can be specified as follows

• For γ = 1 2 , we get ∆H -W ext k+1 -W damping k+1 = U ⊤ k+1/2 P k+1 . (91)
The algorithmic energy K only changes at each time-step by the amount of work that is done by the dissipative effects in the system and the work supplied by the external forces.

• For γ = 2β 1/2, E(q, v) = H(q, v, a) and we have

∆E -W ext k+1 -W damping k+1 = ( 1 2 -γ) (q k+1 -q k ) 2 K + P ⊤ k+1 U k+1/2 (92)
Since P ⊤ k+1 U k+1/2 0, the scheme always dissipates the mechanical energy of the system.

• For γ = 2β = 1 2 , we get

∆E -W ext k+1 -W damping k+1 = P ⊤ k+1 U k+1/2 (93)
The total mechanical energy E only changes at each time-step by the amount of work that is done by the dissipative effects in the system and the work supplied by the external forces.

The results stated in ( 92) and ( 93) are better than for the standard Moreau-Jean scheme with the θ-method since there is no condition on γ. The results apply to the variant of the Moreau-Jean scheme (31).

Conclusion

The main results of this paper and some perspectives may be summarized now:

• A study of the energy conservation and dissipation of the Moreau-Jean scheme has been carried out. Under suitable assumptions on the numerical parameter θ with respect to the coefficients of restitution, we show that the Moreau-Jean scheme dissipates the total mechanical energy. For θ = 1/2, the dissipated energy is only due to the discrete works of the damping, external forces and the impact impulses. Moreover, a simple variant of the Moreau-Jean scheme that always dissipates energy has been proposed. This latter method removes the condition on θ with respect to the coefficients of restitution.

• As we said in the introduction, there is a lot of difficulties to establish energy conservation or dissipation results in the context of Newmark schemes with unilateral constraints on the generalized coordinates. One of the difficulties raises from the fact that for discrete (or spacediscretized) mechanical systems, we have to introduce the notion of coefficient of restitution, or at least to give a law that defines the velocity after the impact. The second difficulty is related to higher order approximation of nonsmooth terms like impulses. With the new schemes developed in [START_REF] Chen | A Newmark-Type Integrator for Flexible Systems Considering Nonsmooth Unilateral Constraints[END_REF][START_REF] Chen | A nonsmooth generalized-α scheme for flexible multibody systems with unilateral constraints[END_REF], these two difficulties are overcome. This allows us to extend the known results on the energy conservation and dissipation properties to the nonsmooth case with unilateral contacts and impacts. In particular, we prove that the Newmark scheme is able to conserve or dissipates an algorithmic energy H(q, v, a) that depends on the coordinates, the velocities and the acceleration. For HHT scheme, we are also to prove the same kind of results with an extended algorithmic energy H(q, v, a, z) based on the introduction of additional filters. For the generalized-α scheme, we obtain equivalent results as in the smooth case. However, this cannot allow to easily conclude. Nevertheless, the contribution of the contact terms in the discrete energy balance (69), that is U ⊤ k+1/2 P k+1 , is identical for the Newmark and the HHT scheme. We can infer that this result combined with the classical result of [START_REF] Arnold | Convergence of the generalized-α scheme for constrained mechanical systems[END_REF], [START_REF] Jay | Extensions of the HHT-method to differential-algebraic equations in mechanics[END_REF], [START_REF] Ch | Solving constrained mechanical systems by the family of Newmark and α-methods[END_REF] on stability of the generalized-α when unilateral constraints remain closed should imply the stability of the scheme in practice. Nevertheless, we extend a variant of the generalized-α scheme due to [START_REF] Krenk | Properties of time integration with first order filter damping[END_REF] to the contact case. For this latter case, we obtain the same results as in the smoooth case.

• The dissipation properties of the schemes studied in this paper allows us to conclude to the boundedness of the total mechanical energy of the system E(q, v) for the Moreau-Jean scheme and to the boundedness the algorithmic energy H(q, v, a, z) for the Newmark and the α-schemes by adding some standard assumption on the works of the external forces. Let us remark that the total mechanical energy E(q, v) and the algorithmic energy H(q, v, a, z) are positive semi-definite functions if we assume that M is positive definite, K is positive semi-definite and the conditions of dissipation of the schemes are satisfied. The boundedness of theses energies implies the boundedness of the discrete velocities v k and acceleration a k which in some sense guarantee the stability of the scheme in the numerical practice. The boundedness of the discrete generalized coordinates can also easily conclude if the stiffness matrix is positive definite, or by simply inspecting the relations that relate the discrete generalized coordinates to the discrete generalized velocities and accelerations. In the smooth case (Lipschitz ordinary differential equation), boundedness properties are equivalent to the stability of the linear multi-step time integration method thanks to the theory developed by Dalhquist, Lax and Richtmyer [START_REF] Dahlquist | 33 years of numerical instability, part i[END_REF][START_REF] Lax | Survey of the stability of linear finite difference equations[END_REF][START_REF] Richtmyer | Difference Methods For Initial-Value Problems[END_REF]. Classical approaches to study the stability of numerical scheme for ordinary differential equations trough its amplification matrix or the roots of the stability function were also extended with success to the case of differential algebraic equations [START_REF] Cardona | Time integration of the equations of motion in mechanism anslysis[END_REF][START_REF] Ascher | Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems[END_REF]. Unfortunately, in our case, we cannot directly conclude to the stability of the schemes by any of these methods since Inria we are in a nonsmooth case and the numerical scheme does not reduce to a linear-multistep formula. Nevertheless, it can be a first step to prove the convergence of the scheme.

• The first perspective that can be drawn for this work is the possible adaptation of the results to the case with Coulomb's friction. One of the main difficulties is that even in the continuous time-case the modeling of Coulomb's friction at the impulse level together with an impact law does not lead necessarily to a dissipative system (see [START_REF] Leine | Stability and Convergence of Mechanical Systems with Unilateral Constraints[END_REF] for details). The second difficulty should be to prove that the numerical scheme correctly discretize the friction such that it still dissipates energy. The second perspective is the extension towards to the nonlinear case as it has been done some special case of the mid-point rule in [START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impacts[END_REF]Le Tallec, 2006, Ayyad et al., 2009].

A Proof of Lemma 2

To prove the result, we start from the dissipation analysis of the Newmark scheme. Following the method introduced in [START_REF] Hughes | A note on the stability of Newmark's algorihtm in non-linear structural dynamics[END_REF] and fully developed in [START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF], we define an intermediate discrete "algorithmic energy" of the form

K(q, v, a) = E(q, v) + h 2 4 (2β -γ)a ⊤ M a. ( 94
)
Although there is no direct mechanical interpretation of this quantity, it allows one to conclude on the dissipation property of the scheme since K is a semi-norm for 2β γ. From the definition of the nonsmooth Newmark scheme, we get

q k+1 -q k = h 2 (v k+1 + v k ) + h 2 2 (2β -γ)(a k+1 -a k ). (95)
and

v k+1 -v k = h 2 (a k+1 + a k ) + h(γ - 1 2 )(a k+1 -a k ) + w k+1 . (96)
Using ( 95) and ( 96), another formulation of ∆E in ( 36) is given by

2∆E = (v k+1 + v k ) ⊤ M h 2 (a k+1 + a k ) + h(γ - 1 2 )(a k+1 -a k ) + w k+1 +(q k+1 + q k ) ⊤ K h 2 (v k+1 + v k ) + h 2 2 (2β -γ)(a k+1 -a k ) (97)
Inria

RR n°8602

Using the definition of the nonsmooth Newmark scheme, the relation (97) can be developed as in

(98) Using again ( 96) in ( 98), we get with some manipulations

(99) Since the term (a k+1 -a k ) ⊤ M (a k+1 + a k ) appears in the last line of (99), the computation of ∆K can be developed as follows

Substituting the definition of W ext k+1 and W damping k+1 in (100) yields

(101) From ( 67), the additional damping force (65) generates an additional term in (101) given by

Acary

For the sake of simplicity, we consider for a moment that C = 0. The damping term will be discussed later. Developing the additional term (102) using the rule

we obtain

Let us now focus on the second term in the right hand side of (104). Using

we obtain

(106) Let us restart from (100) with the additional term (q k+1 -q k ) ⊤ f A k+γ , we get

With the definition of H in (66), we can simplify (107) in

The term depending the damping matrix can be included by adding to the external forces a damping force -Cv and the additional term η ν x k+γ . Using the definition of the discrete works in ( 67) and (68) the result is obtained.

Inria

RESEARCH CENTRE GRENOBLE -RHÔNE-ALPES


Julien Letemplier

Thomas Simon

THE AREA OF A SPECTRALLY POSITIVE STABLE PROCESS STOPPED AT ZERO

Keywords: 2010 Mathematics Subject Classification. 60E07, 60G51, 60G52. Key words and phrases. Exponential functional -Hitting time -Integrated process -Moments of Gamma type -Self-decomposability -Series representation -Stable Lévy process

An identity in law for the area of a spectrally positive Lévy α-stable process stopped at zero is established. Extending that of Lefebvre [14] for Brownian motion, it involves an inverse Beta random variable and the square of a positive stable random variable. This identity entails that the stopped area is distributed as the perpetuity of a spectrally negative Lévy process, and is hence self-decomposable. We also derive a convergent series representation for the density, whose behaviour at zero is shown to be Fréchet-like.

Introduction and statement of the results

Let {B t , t ≥ 0} be a standard linear Brownian motion, starting from one, and let T = inf{t > 0, B t = 0} be its first hitting time of zero. The random variable A = T 0 B s ds has been investigated by Lefebvre, who obtained in Theorem 2 of [START_REF] Lefebvre | First-passage densities of a two-dimensional process[END_REF] the simple identity in law (1.1)

A d = 2 9Γ 1/3
where, here and throughout, Γ a stands for the Gamma random variable with density x a-1 Γ(a) e -x 1 {x>0} .

The identity (1.1) is obtained as a consequence of the closed expression for the Laplace transform of the bivariate random variable (T, A) in terms of the Airy function -see Theorem 1 in [START_REF] Lefebvre | First-passage densities of a two-dimensional process[END_REF]. As observed in [START_REF]Sur le premier instant de passage de l'intégrale du mouvement brownien[END_REF] p. 402, this latter expression can be easily derived thanks to the Feynman-Kac formula. Notice that Airy functions also appear in the expression of the Laplace transform of many other Brownian areas -see [START_REF] Janson | Brownian excursion area, Wright's constants in graph enumeration, and other Brownian areas[END_REF], whose laws are more complicated than (1.1).

In this paper, our concern is to generalize (1.1) to the random variables

A α = T 0 L (α) s ds
where {L (α) t , t ≥ 0} is a strictly α-stable Lévy process without negative jumps, starting from one, and T = inf{t > 0, L (α) t = 0} is its first hitting time of zero. Without loss of generality we choose the normalization E[e -tL (α) 1 ] = e t α , t ≥ 0, where α ∈ [START_REF] Bingham | Regular Variation[END_REF][START_REF] Bondesson | Generalized Gamma convolutions and related classes of distributions and densities[END_REF] is the self-similarity parameter. We refer e.g. to Chapter 3 in [START_REF] Sato | Lévy processes and infinitely divisible distributions[END_REF] for more information on stable Lévy processes and the above normalization. The boundary cases α = 1, 2 correspond to the unit drift resp. the Brownian motion with variance √ 2, so that we have (1.2)

A 1 = 1 2 and A 2 d = 1 9Γ 1/3 •
The above second identity, which is actually the precise statement of Theorem 2 in [START_REF] Lefebvre | First-passage densities of a two-dimensional process[END_REF], follows from (1.1) and the self-similarity of Brownian motion. In the case α ∈ (1, 2) the generator of the process L (α) is non-local, and it seems unappropriate to appeal to Feynman-Kac formulae in order to obtain a tractable expression for the Laplace transform for A α .

Moreover, the absence of transition densities written in closed form prevents from using explicit computations as in the Gaussian case -see [START_REF] Lefebvre | First-passage densities of a two-dimensional process[END_REF] and the references therein, to handle the random variable A α . Instead, we will compute the fractional moments of A α and exhibit a multiplicative identity in law. Introduce the Beta random variable B a,b with density

Γ(a + b) Γ(a)Γ(b) x a-1 (1 -x) b-1 1 (0,1) (x),
and the positive a-stable random variable Z a with Laplace transform

E[e -λZa ] = e -λ a , λ ≥ 0.
Our main observation is the following.

Theorem. With the above notation, one has the independent factorization

(1.3) A α d = α + 1 4 × Z 2 2 α+1 × B -1 1 2 , α-1 2(α+1)
for every α ∈ (1, 2).

Observe that (1.3) is in accordance with the two boundary cases: when α = 1 the two random variables on the right-hand side boil down to one, whereas when α = 2 the following identity obtained in Theorem 1 of [START_REF] Simon | Comparing Fréchet and positive stable laws[END_REF]:

Z 2 2 3 d = 4 27 Γ -1 2 3 × B -1 1 3 , 1 6
, combined with the elementary factorization Γ a d = Γ a+b × B a,b , allows to recover the second identity in (1.2). The proof of (1.3) relies on an identification of the fractional moments of A α . The explicit computation of the latter in terms of Gamma functions -see (2.1) belowis made possible by the strong Markov property and some exact results on the stable Kolmogorov process recently obtained in [START_REF] Profeta | Persistence of integrated stable processes[END_REF].

The inverse Gamma random variable involved in Lefebvre's identity shares a number of distributional properties related to infinite divisibility. Recall that a non-negative random variable X is said to be self-decomposable if the following identities hold

X d = cX + X c
for every c ∈ (0, 1) with X c independent of X, or equivalently if its log-Laplace transform is of the type

-log E[e -λX ] = aλ + ∞ 0 (1 -e -λx ) k(x) x dx, λ ≥ 0
for some a ≥ 0 and a non-negative, non-increasing function k integrating 1 ∧ x -1 . See again Chapter 3 in [START_REF] Sato | Lévy processes and infinitely divisible distributions[END_REF] for details. The fact that all inverse Gamma random variables are selfdecomposable can be observed either by a direct and non-trivial computation on the Laplace transform, or by Dufresne's celebrated identity for the perpetuity of a Brownian motion with drift -see Section 3.2 in [START_REF] Bosch | On the self-decomposability of the Fréchet distribution[END_REF] for details and references. In this paper we will use the second approach and show the same property for A α .

Corollary 1. For every α ∈ (1, 2), the random variable A α is self-decomposable.

More precisely, it will be proved that A α is distributed as the perpetuity of a spectrally negative Lévy process which drifts towards +∞. The latter background integrated Lévy process turns out to be tightly connected with the dual process -L (α) conditioned to stay positive, and also to the Fréchet distribution which is hidden in the factorization (1.3) and can be viewed as another perpetuity -see Remark 3 below.

It is known that inverse Gamma distributions also satisfy a property which is more stringent than self-decomposability. The law of a non-negative random variable X is called a generalized Gamma convolution (X ∈ G for short) if it is the weak limit of an independent sum of Gamma random variables, or equivalently if its log-Laplace transform reads

-log E[e -λX ] = aλ + ∞ 0 log x x + λ U(dx), λ ≥ 0,
for some a ≥ 0 and some non-negative measure U satisfying certain integrability conditions, which is called the Thorin measure of X. We refer to Chapter 3 in [START_REF] Bondesson | Generalized Gamma convolutions and related classes of distributions and densities[END_REF] for more details on this notion refining that of self-decomposability. The fact that inverse Gamma random variables are all in G, with explicit Thorin measure, is a consequence of a computation on Bessel functions -see again Section 3.2 in [START_REF] Bosch | On the self-decomposability of the Fréchet distribution[END_REF] for details and references. Although we cannot write the Laplace transform of A α in any tractable way, we can prove the following result.

Corollary 2. For every α ∈ (1, 2), the law of √ A α is a generalized Gamma convolution, with infinite Thorin measure.

A conjecture by Bondesson -see [START_REF] Bondesson | Generalized Gamma convolutions and related classes of distributions and densities[END_REF] p. 97 and also Section 7 in [START_REF] Bondesson | A class of probability distributions that is closed with respect to addition as well as multiplication of independent random variables[END_REF] -states that the G-property is stable by power transformation of order greater than one. If this conjecture is true, then Corollary 2 entails that A α ∈ G as well, a reinforcement of Corollary 1. It will be shown in Section 3.1. below that A 5/3 ∈ G, but we do not know as yet how to handle the other values of α.

The closed expression for the fractional moments of A α can be inverted in order to give a convergent series representation for its density f Aα . Throughout this paper, we will set f X for the density of an absolutely continuous random variable X.

Corollary 3. The density of A α has a convergent series representation:

f Aα (x) = Γ α α + 1 × ∞ n=0 (-1) n (α + 1) n+1 α+1 -1 x -n+1 α+1 -1 n! Γ(1 -n+1 α+1 )Γ(1 -n+2 α+1 ) , x > 0.
Observe that when α = 2 the above summation is made over n = 3p only, and that further simplifications lead to

(1.4) f A 2 (x) = x -4/3 3 2/3 Γ(1/3) ∞ p=0 (-1) p (9x) -p p! = Γ(2/3) x -4/3 e -1 9x 2π 3 1/6 ,
which is the expression to be found in Theorem 2 of [START_REF] Lefebvre | First-passage densities of a two-dimensional process[END_REF]. The above corollary also entails the first order asymptotics

f Aα (x) ∼ (α + 1) 1 α+1 -1 x -1 α+1 -1 Γ( α-1 α+1 ) as x → +∞,
which has, up to the multiplicative constant, the same speed as that of the density of the factor Z 2 2/(α+1) at infinity -see Formula (14.31) in [START_REF] Sato | Lévy processes and infinitely divisible distributions[END_REF]. On the other hand, it does not seem possible to deduce from the above series representation the exact behaviour of f Aα at zero. Nevertheless, using the identity (1.3) we can show the following estimate.

Corollary 4. The asymptotic behaviour of the density of

A α when x → 0+ is (1.5) f Aα (x) ∼ κ α x α 2 1-α 2 e -cα x 1 1-α , with κ α = Γ( α α+1 ) α+1 α-1 2π (α + 1) α α 2 -1 and c α = (α -1)(α + 1) α 1-α .
This shows that the behaviour of f Aα at zero is that of the generalized Fréchet density

f c α-1 α Γ 1-α 1/(α+1) (x) = κα x α 2 1-α 2 e -cα x 1 1-α , up to the normalizing constant κα = (α -1) -α α+1 Γ( 1 α+1 ) (α + 1) α α 2 -1
which does not coincide with κ α except for α = 2. Observe also that making α = 2 on the right-hand side of (1.5) yields the density in (1.4). It should be possible to obtain a full asymptotic expansion of f Aα at zero with our method -see Remark 5 below. But we have not adressed this issue, which is believed to be awfully technical, in the present paper.

Proofs

Proof of the Theorem.

To simplify the notation we will set L = L (α) . Introducing the area process

A t = t 0 L s ds, t ≥ 0,
recall that the bivariate process X = {(A t , L t ), t ≥ 0} is strongly Markovian and denote by P (x,y) its law starting from (x, y). Consider the stopping time

S = inf{t > 0, A t = 0}
and observe that under P (0,1) one has a.s. S > T, A T > 0, and L S < 0. Setting {F t , t ≥ 0} for the natural completed filtration of X and applying the strong Markov property at T entails that for every s ∈ R, one has

E (0,1) [|L S | s-1 ] = E (0,1) [ E[|L S | s-1 |F T ]] = E (0,1) [ E (A T ,0) [|L S | s-1 ]] = E [A s-1 α+1 α ] × E (1,0) [|L S | s-1 ]
, possibly with infinite values on both sides, where the second equality follows from the absence of negative jumps for L, and the third equality from the self-similarity of L and A with respective indices 1/α and 1 + 1/α. Applying Theorem B in [START_REF] Profeta | Persistence of integrated stable processes[END_REF] in the particular case ρ = 1/α (beware that we consider here the dual process, with no positive jumps), we get

E (0,1) [|L S | s-1 ] = sin( πs α+1 ) sin( παs α+1
) and

E (1,0) [|L S | s-1 ] = (α + 1) 1-s α+1 Γ( α+2 α+1 )Γ( 1-s α+1 ) sin( π α+1 ) Γ( s α+1 )Γ(1 -s) sin( παs α+1
) for all |s| < 1 + 1/α. Dividing and simplifying with the help of the complement formula for the Gamma function, we deduce

(2.1) E [A s α ] = (α + 1) s × Γ( α α+1 )Γ(1 -(α + 1)s) Γ( α α+1 -s)Γ(1 -
s) for all s < 1/(α + 1). Applying the Legendre-Gauss multiplication formula for the Gamma function entails

E [A s α ] = α + 1 4 s × Γ(1 -(α + 1)s) Γ(1 -2s) × Γ( α α+1 )Γ( 1 2 -s) Γ( 1 2 )Γ( α α+1 -s) ,
and we can conclude by a fractional moment identification, recalling (see e.g. Formula (25.5) in [START_REF] Sato | Lévy processes and infinitely divisible distributions[END_REF] for the second expression) that

E [B s a,b ] = Γ(a + s)Γ(a + b) Γ(a)Γ(a + b + s) and E [Z s a ] = Γ(1 -s a ) Γ(1 -s) • Remark 1. (a)
It is well-known and easy to see -see e.g. Theorem 46.3 in [START_REF] Sato | Lévy processes and infinitely divisible distributions[END_REF] -that under P (0,1) , the random variable T is distributed as Z 1/α . The above theorem provides hence a connection between A α = A T and the random variable Z 2/(α+1) . Notice that one can also derive from (2.1) another factorization:

Z 1 α+1 d = (α + 1) -1 × Γ -1 α α+1 × A α .
However, it seems difficult with our method to obtain some valuable information on the Mellin transform of the bivariate random variable (T, A T ).

(b) With the notation of our above proof, it is possible to derive the law of A T under P (x,y) for any x ∈ R and y > 0, by the self-similarity of L (α) . One finds

A α d = x + (α + 1) y α+1 4 × Z 2 2 α+1 × B -1 1 2 , α-1 2(α+1)
.

Proof of Corollary 1.

Let us first observe that A α is infinitely divisible, by a simple pathwise argument not relying on (1.3). Setting T x = inf{t > 0, L t = x} for all x > 0 and using the fact that L has no negative jumps, it is easy to see from the Markov property that under P (0,1) , for every n ≥ 2, there is an independent decomposition

A α = X (n) 1 + • • • + X (n) n where X (n) i d = A T n-i n under P (0, n+1-i n )
for every i = 1, . . . , n. Moreover, one has T 1-1/n → 0 a.s. under P (0,1) as n → +∞ (by the well-known fact -see e.g. Theorem 47.1 in [START_REF] Sato | Lévy processes and infinitely divisible distributions[END_REF] -that L visits immediately the negative half-line when starting from 0), so that

X (n) 1
→ 0 a.s. under P (0,1) when n → +∞ as well. Last, it is straightforward that

P[|X (n) i | > ε] ≤ P[|X (n) 1 | > ε]
for every ε > 0 and i = 1, . . . , n. Putting everything together and applying Khintchine's theorem on triangular arrays -see e.g. Theorem 9.3 in [START_REF] Sato | Lévy processes and infinitely divisible distributions[END_REF], entails that A α is infinitely divisible.

Remark 2. The above argument does not make use of the self-similarity of L, and hence applies to any spectrally positive Lévy process which is not a subordinator.

We now proceed to the proof of the self-decomposability of A α . We will use the same argument as in [START_REF] Bosch | On the self-decomposability of the Fréchet distribution[END_REF][START_REF] Bosch | On the infinite divisibility of inverse Beta distributions[END_REF], expressing A α as the perpetuity of a certain spectrally negative Lévy process. Setting

Ψ α (u) = u E[A -(u+1) α ] E[A -u α ]
for every u > 0, we first deduce from (1.3) the formula Ψ α (u) = Φ α (s), with the notation s = (α + 1)u and

Φ α (s) = Γ(α + s) Γ(s) •
Applying the Lemma in [START_REF] Bosch | On the self-decomposability of the Fréchet distribution[END_REF] with t = 1 shows, after some simplifications, that

Φ α (s) = Γ(α) s + 0 -∞ (e sx -1 -sx)f α (x) dx with f α (x) = e αx Γ(-α)(1 -e x ) α+1 • Since f α integrates x 2 ∧ 1 on (-∞, 0
), this entails that Ψ α is the Laplace exponent of a spectrally negative Lévy process. Applying Bertoin-Yor's criterion for perpetuities, as stated in [START_REF] Bosch | On the infinite divisibility of inverse Beta distributions[END_REF] pp. 8-9, shows that

(2.2) A α d = ∞ 0 e -Z (α) t dt
where {Z (α) t , t ≥ 0} is the spectrally negative Lévy process with Laplace exponent

E[e λZ (α) 1 ] = e Ψα(λ) = e Φα((α+1)λ) , λ ≥ 0.
It is then easy to see from the representation (2.2) and the spectral negativity of Z (α) that A α is self-decomposable -see the end of the proof of the Theorem in [START_REF] Bosch | On the self-decomposability of the Fréchet distribution[END_REF].

Remark 3. (a) The above expression (2.2) extends to the boundary cases α = 1, 2. When α = 1, the Lévy process Z (1) has Laplace exponent 2λ, so that (2.2) boils down to

A 1 d = ∞ 0 e -2t dt = 1 2 •
When α = 2, the Lévy process Z (2) has Laplace exponent 3λ + 9λ 2 , and (2.2) reads

A 2 d = ∞ 0 e 3 √ 2Bt-3t dt d = 1 18 ∞ 0 e Bt-t/6 dt d = 1 9Γ 1/3
, the third identity in law being a particular case of Dufresne's identity.

(b) It follows from Corollary 2 in [START_REF] Caballero | Conditioned stable Lévy processes and the Lamperti representation[END_REF] that the spectrally negative Lévy process {ξ ↑,n t , t ≥ 0} appearing in the Lamperti transform of the dual process -L (α) conditioned to stay positive, has log-Laplace exponent

-log E[e λξ ↑,n 1 ] = Φ α (λ), λ ≥ 0,
with the above notation for Φ α . This can be shown from the second formula in Corollary 2 of [START_REF] Caballero | Conditioned stable Lévy processes and the Lamperti representation[END_REF] written in an appropriate way, bewaring the unusual notation (7) therein for the negativity parameter and correcting a misprint (the + before c -should be a -) in Formula ( 17) therein. This entails the curious identity in law

(2.3) A T d = ∞ 0 e -(α+1)ξ ↑,n t dt.
Recall that Brownian motion conditioned to stay positive is distributed as a three-dimensional Bessel process, whose Lamperti process ξ ↑,n is the drifted Brownian motion {B t + t/2, t ≥ 0} -see [START_REF] Caballero | Conditioned stable Lévy processes and the Lamperti representation[END_REF] p. 969 and the references therein, so that by Dufresne's identity we obtain

∞ 0 e -3ξ ↑,n t dt d = 2 9Γ1 3 d = A.
This can be viewed as a particular case of (2.3), with the proper normalization. It is quite interesting to compare (2.3) with the identity

(2.4) T d = ∞ 0 e -αξ ↑,n t dt,
which follows from the above Remark 1 and Theorem 7 of [START_REF] Chaumont | Some explicit identities associated with positive selfsimilar Markov processes[END_REF]. We do not have any sensible interpretation for the structural relationship between (2.3) and (2.4). Notice last that the perpetuities of certain Lévy processes with positive mean and jumping density

K e bx (1 -e x ) α+1 1 {x<0}
for some K, b > 0 and α ∈ (1, 2) have been studied in [START_REF] Patie | Exponential functional of a new family of Lévy processes and self-similar continuous state branching processes with immigration[END_REF]. Observe in particular that the factorization obtained in Theorem 4.6 (2) therein shares some similarities with (1.3).

(c) Combining the Kanter factorization -see Corollary 4.1 in [START_REF] Kanter | Stable densities under change of scale and total variation inequalities[END_REF] -and (1.3) shows the identity

A α d = α + 1 4 × Γ 1-α 1 × B 1 2 , α-1 2(α+1) × K α+1 2 α+1 -1
, where K a is the so-called Kanter random variable of index a ∈ (0, 1) -see Section 3 in [START_REF] Simon | Comparing Fréchet and positive stable laws[END_REF] for more details about this random variable. Consider now the Fréchet random variable Γ 1-α 1 appearing in the above factorization of A α . In [START_REF] Bosch | On the self-decomposability of the Fréchet distribution[END_REF], it was shown that Γ 1-α 1 is also distributed as the perpetuity of a spectrally negative Lévy process, and the latter turns out to be quite close to the above process Z (α) . More precisely, it follows from the Lemma in [START_REF] Bosch | On the self-decomposability of the Fréchet distribution[END_REF] with t = 1, the proof of the Theorem in [START_REF] Bosch | On the self-decomposability of the Fréchet distribution[END_REF] and a change of variable that (2.5)

Γ 1-α 1 d = ∞ 0 e -Z(α) t dt
where { Z(α) t , t ≥ 0} is the spectrally negative Lévy process having log-Laplace exponent

-log E[e λ Z(α) 1 ] = (α -1) -1 Φ α ((α -1)λ)), λ ≥ 0.
We do not have any explanation for this proximity between the two Laplace exponents of Z (α) and Z(α) . Notice that the identity (2.5) can also been derived, in a complicated way, from Theorem 4.1 in [START_REF] Patie | Exponential functional of a new family of Lévy processes and self-similar continuous state branching processes with immigration[END_REF].

Proof of Corollary 2.

From (1.3) and a general result by Bondesson -see Theorem 1 in [START_REF] Bondesson | A class of probability distributions that is closed with respect to addition as well as multiplication of independent random variables[END_REF], in order to show

√ A α ∈ G it is enough to prove that B -1/2 1 2 , α-1 2(α+1) ∈ G and Z 2 α+1 ∈ G.
The second fact is well-known -see Example 3.2.1 in [START_REF] Bondesson | Generalized Gamma convolutions and related classes of distributions and densities[END_REF], whereas the first one follows at once from Theorem 2 in [START_REF] Bosch | On the infinite divisibility of inverse Beta distributions[END_REF]. Finally, the fact that √ A α has infinite Thorin measure follows from Theorem 4.1.4 in [START_REF] Bondesson | Generalized Gamma convolutions and related classes of distributions and densities[END_REF] and the above formula (2.1), which entails that √ A α has negative moments of all orders. Remark 4. A combination of Corollaries 1 and 2 shows the following identities in law between three random integrals:

∞ 0 e -Z (α) t dt d = T 0 L (α) t dt d = ∞ 0 a (α) t dΓ t 2 ,
where {Γ t , t ≥ 0} is the Gamma subordinator and {a 2.4. Proof of Corollary 3. We will reason along the same lines as in Proposition 2 in [START_REF] Simon | Hitting densities for spectrally positive stable processes[END_REF], and omit some details. Applying the Mellin inversion formula yields first

f Aα (x) = 1 2πx R M α (s) x -is ds,
with the notation

M α (s) = E [A is α ] = (α + 1) is × Γ( α α+1 )Γ(1 -(α + 1)is) Γ( α α+1 -is)Γ(1 -
is) for every s ∈ R. Suppose first x > 1. We evaluate the above integral with the help of the residue theorem applied to the contour joining -R to R along the real axis, and R to -R along the half-circle plotted in the lower half-plane. It is easy to see that the integral along this half-circle vanishes as R → +∞, so that it remains to consider the singularities inside the big contour, which are located at t n = -i(n + 1)/(α + 1), n ≥ 0. Computing

Res tn (M α (s)x -is ) = -i Γ α α + 1 × (-1) n (α + 1) n+1 α+1 -1 x -n+1 α+1 n! Γ(1 -n+1 α+1 )Γ(1 -n+2 α+1 ) we deduce f Aα (x) = Γ α α + 1 × ∞ n=0 (-1) n (α + 1) n+1 α+1 -1 x -n+1 α+1 -1 n! Γ(1 -n+1 α+1 )Γ(1 -n+2 α+1
) for every x > 1, and hence for every x > 0 by analytic continuation (Stirling's formula shows indeed that the series on the right-hand side converges absolutely for every x > 0).

Proof of Corollary 4. We will work on the random variable X

α = A 1 1-α
α , in order to simplify the notation. Changing the variable, the required estimate is tantamount to

(2.6) f Xα (x) ∼ Γ( α α+1 ) √ α 2 -1 2π (α + 1) α α 2 -1 x -α α+1 e -cαx , x → 0+.
Evaluating with (2.1) the positive entire moments

E[X n α ] = (α + 1) n 1-α × Γ( α α+1 ) Γ( α α+1 + n α-1 ) × Γ(1 + (α+1)n α-1 ) Γ(1 + n α-1
) for every n ≥ 0, and applying Stirling's formula shows that

(E[X n α ]) 1 n n → (α + 1) α α-1 e(α -1) as n → ∞.
By a theorem of Davies-Kasahara (see Corollary 4.12.5 in [START_REF] Bingham | Regular Variation[END_REF], or Lemma 3.2 in [START_REF] Csörgö | Some asymptotic properties of the local time of the uniform empirical process[END_REF] for a more appropriate formulation), we deduce that

x -1 log P[X α > x] → -c α as x → +∞.
This yields the required asymptotic behaviour, at the logarithmic scale, for the survival function of X α . Moreover, writing down via Fubini's theorem the moment generating function

E[e xXα ] = ∞ n=0 a n x n , x > 0,
with

a n = (α + 1) n 1-α × Γ( α α+1 )Γ(1 + (α+1)n α-1 ) n! Γ( α α+1 + n α-1 )Γ(1 + n α-1 ) ∼ Γ( α α+1 ) √ α 2 -1 2π (α -1) 1 α+1 c -n α n -α α+1
as n → ∞, and applying Karamata's theorem for power series -see Corollary 1.7.3 in [START_REF] Bingham | Regular Variation[END_REF] -shows that (2.7)

E[e cαxXα ] ∼ Γ( 1 α+1 ) Γ( α α+1 ) √ α 2 -1 2π (α -1) 1 α+1 (1 -x) 1 α+1 as x → 1-.
At this stage, it is worth mentioning that (2.7) can be obtained from (2.6) by integration. However, it does not seem that we can infer the reverse inclusion without any further assumption, such as the existence of a meromorphic extension in the neighbourhood of c α for the moment generating function -see Theorem 4 in [START_REF] Flajolet | Mellin transforms and asymptotics: Harmonic sums[END_REF], or a monotonicity condition on f Xα at infinity -see Theorem 4.12.11 in [START_REF] Bingham | Regular Variation[END_REF], which we both could not prove a priori.

In order to show (2.6) rigorously and finish the proof, we will use the following power transformation of (1.3):

X α d = 4 α + 1 1 α-1 × B 1 α-1 1 2 , α-1 2(α+1) × Z 2 1-α 2 α+1
.

The multiplicative convolution formula and a change of variable entails

(2.8) f Xα (x) = ∞ 0 1 1 + y f Uα 1 1 + y f Vα (x + xy) dy,
where we have set

U α = B 1 α-1 1 2 , α-1 2(α+1) and V α = 4 α + 1 1 α-1 × Z 2 1-α 2 α+1
.

On the one hand, we have

1 1 + y f Uα 1 1 + y = (α -1) α-1 2(α+1) Γ( α α+1 ) √ π Γ( α-1 2(α+1) ) y α-1 2(α+1) -1 (1 + O(y))
as y → 0+.

On the other hand, Formula (14.35) in [START_REF] Sato | Lévy processes and infinitely divisible distributions[END_REF] entail after a change of variable and several simplifications

f Vα (z) = √ α -1 (α + 1) α 2 2(α 2 -1) -1 α-1 2 √ π z -1/2 e -cαz (1 + O(z -1/2 )) as z → +∞.
Plugging these two first order expansions in the integral (2.8), and making further simplifications, yields finally the required asymptotic behaviour (2.6).

Remark 5. The two asymptotic expansions for the above f Uα and f Vα can be continued at every order -see again Formula (14.35) in [START_REF] Sato | Lévy processes and infinitely divisible distributions[END_REF] for the second function. This could be used to obtain a refined expansion for f Xα at infinity, or equivalently for f Aα at zero. Displaying the full asymptotic expansion of f Aα at zero seems however to be a very painful task.

Final remarks

3.1. On the G-property. As mentioned in the introduction, the fact that √ A α ∈ G gives some credit to the G-property for A α itself. This refinement of Corollary 1 would also shed some analytic light on the self-decomposability of A α . We can show this property for α = 5/3, where a combination of (1.3) and Theorem 1 in [START_REF] Simon | Comparing Fréchet and positive stable laws[END_REF] entails the identity in law

A 5/3 d = 3 2 13/3 × Γ -2 3 3 4 × B -2 3 1 4 , 1 12 × B -2 3 1 2 , 1 6 × B -1 1 2 , 1 8 .
Indeed, all random variable on the right-hand side are in G by Theorems 2 and 4 in [START_REF] Bosch | On the infinite divisibility of inverse Beta distributions[END_REF], and we can conclude by Theorem 1 in [START_REF] Bondesson | A class of probability distributions that is closed with respect to addition as well as multiplication of independent random variables[END_REF]. Notice that the same kind of argument would show the property for all α ∈ (1, 2) if we could prove that Z 2 β ∈ G for all β ∈ (2/3, 1) and that by the stability of the G-property with respect to weak convergence -see Theorem 3.1.5 in [START_REF] Bondesson | Generalized Gamma convolutions and related classes of distributions and densities[END_REF], it is enough to consider β rational. Again, Theorem 1 in [START_REF] Simon | Comparing Fréchet and positive stable laws[END_REF] shows a certain factorisation of Z 2 β with β rational into the product of one Fréchet and several Beta random variables at a unique negative power, which unfortunately becomes too small when α = 5/3, so that cannot use Theorem 5.2 in [START_REF] Bosch | On the infinite divisibility of inverse Beta distributions[END_REF] as above. Observe also from Section 3.1 in [START_REF] Bosch | On the infinite divisibility of inverse Beta distributions[END_REF] that small negative powers of Beta random variables may not be in G. To show the plausible fact that Z 2 β ∈ G for all β ∈ (0, 1) rational, it could be useful to write down the density of the Beta-Gamma product given in Formula (2.3) of [START_REF] Simon | Comparing Fréchet and positive stable laws[END_REF] as a Meijer G-function -see Theorem 9 in [START_REF] Springer | The distribution of products of Beta, Gamma and Gaussian random variables[END_REF], although it does not seem easy at first sight to express the Laplace transform in a tractable way.

3.2. On the bell-shape property. It follows from Corollaries 3 and 4 that the density function of A α is real analytic on (0, +∞), and that all its derivatives vanish at zero and at infinity. Moreover, a consequence of Corollary 2, Wolfe-Yamazato's theorem -see e.g. Theorem 51.3 in [START_REF] Sato | Lévy processes and infinitely divisible distributions[END_REF], and the principle of isolated zeroes, is the strict unimodality of this density function, that is its first derivative vanishes only once on (0, +∞). By Rolle's theorem, we deduce that its n-th derivative vanishes at least n times on (0, +∞), and one can ask whether it vanishes exactly n times for all n ≥ 1. Such a property, which is called the bell-shape in the literature, has been conjectured in [START_REF] Simon | Positive stable densities and the bell-shape[END_REF] for all positive self-decomposable distributions having an infinite spectral function at zero -see Conjecture 1 therein. Observe that A α has such a self-decomposable law by Corollary 2, Corollary 4 and a theorem by Zolotarev -see Remark 28.6 in [START_REF] Sato | Lévy processes and infinitely divisible distributions[END_REF]. Observe also that drawing the density with the help of the series representation of Corollary 3 and some plotting software exhibits the visual bell-shape property for f Aα , whose second derivative does seem to vanish only twice on (0, +∞). where {L (α) t , t ≥ 0} is a spectrally positive stable process. It is natural to ask if (1.3) could be extended in the presence of negative jumps, with T defined as the first passage time below zero. In this case, it is known that the stable process, when starting positive, crosses zero by a negative jump with an explicit expression for the law of the undershoot -see e.g. Remark 42.18 in [START_REF] Sato | Lévy processes and infinitely divisible distributions[END_REF] and the references therein. In order to apply our method, it would be hence necessary to find a closed expression of E (x,y) [|L S | s-1 ] with x > 0 and y < 0, with the notation of our above proof. When y = 0, Theorem B in [START_REF] Profeta | Persistence of integrated stable processes[END_REF] provides a formula for these fractional moments in terms of the Gamma function, in full generality on the stable process. However, it seems hard to get a tractable formula when y < 0.

  (α) t , t ≥ 0} some deterministic function which is related to the Thorin measure of √ A α -see Proposition 1.1 in [10].

3. 3 . 0 L

 30 Stable processes with negative jumps. The present paper studies the law of the random variableT (α)s ds