Frédéric Marazzato
email: frederic.marazzato@enpc.fr

Alexandre Ern
email: alexandre.ern@enpc.fr

Christian Mariotti
email: christian.mariotti@cea.fr

Laurent Monasse
email: laurent.monasse@inria.fr

An explicit pseudo-energy conserving time-integration scheme for Hamiltonian dynamics

HAL is a

Introduction

Energy and momentum conservation is an important property of numerical schemes for a large number of physical problems. For instance, in statistical physics, accurately conserving first integrals constitutes a fundamental requirement to capture the correct behaviour of the system. In mechanics, conservation of the mechanical energy (together with momentum) is an important feature for systems such as the acoustics in a piano [START_REF] Chabassier | Energy preserving schemes for nonlinear Hamiltonian systems of wave equations: Application to the vibrating piano string[END_REF] or nonlinear contact dynamics [START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF][START_REF] Fetecau | Nonsmooth Lagrangian mechanics and variational collision integrators[END_REF]. In this work, we consider Hamiltonian systems consisting of N particles in dimension d (typically, d = 1, 2 or 3) where q i , p i ∈ R d are the position and momentum of the particle i ∈ {1, . . . , N }. We assume that the Hamiltonian has the following split form:

H(q, p) = 1 2 p T M -1 p + V (q), (1)
where q = (q 1 , . . . , q N ) ∈ R dN is the position vector of the particles, p = (p 1 , . . . , p N ) ∈ R dN is the momentum vector of the particles, M is the symmetric positive definite mass matrix and V is the potential energy. The system is thus driven by the equations q = M -1 p, ṗ = -∇V (q). (

Several approaches have been proposed to tackle the issue of conservation when integrating numerically [START_REF] Chabassier | Introduction and study of fourth order theta schemes for linear wave equations[END_REF]. A first possibility consists in the use of symplectic schemes [START_REF] Hairer | Geometric numerical integration: structurepreserving algorithms for ordinary differential equations[END_REF], which integrate a modified (not explicitly exhibited, except in certain simple cases) Hamiltonian and thus preserve the first integrals of the dynamics over exponentially long times (with respect to the time-step), up to fluctuations whose amplitudes grow with the timestep. However, in the case of variable time-steps, symplectic schemes lose their conservation properties since the modified Hamiltonian changes with the time-step [START_REF] Calvo | The development of variable-step symplectic integrators, with application to the two-body problem[END_REF]. When the time-step size is driven by the shape of the Hamiltonian (e.g. in Kepler's problem with high eccentricity), a workaround consists in adding a perturbation accounting for the time-step variation in order for a rescaled dynamic to remain Hamiltonian [START_REF] Hairer | Variable time step integration with symplectic methods[END_REF]. In practice, for mechanical problems, such a condition on the time-step can become impractical, since the time-step could be imposed due to coupling or stiffness phenomena not accounted for in the Hamiltonian part. For an extended review of variational integrators in mechanics, we refer the reader to [START_REF] Marsden | Discrete mechanics and variational integrators[END_REF]. Another approach consists in imposing the exact conservation of energy and momentum at each step of the numerical scheme. Integrating on the constant energy manifold can be carried out using projection [START_REF] Hughes | Transient finite element formulations that preserve energy[END_REF] or Lie group integration [START_REF] Iserles | Lie-group methods[END_REF], but these methods are computationally expensive as soon as the manifold of constant energy and momentum has a complex shape. Another class of methods, energy-momentum conserving schemes, have been proposed in [START_REF] Simo | A new approach to the analysis and simulation of strain softening in solids[END_REF][START_REF] Gonzalez | On the stability of symplectic and energy-momentum algorithms for non-linear Hamiltonian systems with symmetry[END_REF][START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF][START_REF] Chabassier | Energy preserving schemes for nonlinear Hamiltonian systems of wave equations: Application to the vibrating piano string[END_REF] for nonlinear mechanics, contact mechanics and nonlinear wave equations, among others. The general principle is to integrate the nonlinear forces at a special time during the time-step, which is determined through a nonlinear implicit procedure. A higher-order version of these implicit schemes has been derived for linear wave propagation in [START_REF] Chabassier | Introduction and study of fourth order theta schemes for linear wave equations[END_REF]. Variational integrators combining features of symplectic and energy-momentum schemes have been developed for variable time-step strategies [START_REF] Kane | Symplectic-energy-momentum preserving variational integrators[END_REF] and nonlinear mechanical problems in [START_REF] Groß | Conservation properties of a time FE method. part IV: Higher order energy and momentum conserving schemes[END_REF].

To the best of our knowledge, no explicit pseudo-energy conserving scheme has been proposed to date for nonlinear problems. With the motivation that explicit schemes often result in greater computational efficiency, the goal of the present work is to develop such an explicit scheme for nonlinear mechanics, where pseudo-energy conservation holds exactly for exact force integration and up to second-order accuracy in the presence of quadratures. The present scheme hinges on two key ideas. The first one, already considered in [START_REF] Mariotti | A new leapfrog scheme for rotational motion in 3d[END_REF], is to approximate the dynamics of the particles by free-flight trajectories during each time-step. The second one is to use momentum jumps at the discrete time nodes to approximate the acceleration. In doing so, we circumvent the negative result on the existence of explicit schemes in [START_REF] Chabassier | Energy preserving schemes for nonlinear Hamiltonian systems of wave equations: Application to the vibrating piano string[END_REF]Lemma 3.3] through the use of a two-step strategy. This idea has some links with the implicit energy-conserving average vector field method [START_REF] Quispel | A new class of energy-preserving numerical integration methods[END_REF] where the conservation of the Hamiltonian is formulated using an implicit integral of the forces derived from the potential V over the time-step. A high-order generalization of the average vector field method using collocation has been developed in [START_REF] Hairer | Energy-preserving variant of collocation methods[END_REF]. The present numerical scheme shares with average vector field methods the salient feature of average force integration over each time-step. However, the two schemes differ on the discretisation of the acceleration, which is based here on momentum jumps.

A further development of the present work is to devise an asynchronous version of our scheme that lends itself to slow-fast decompositions as presented in [START_REF] Hairer | Geometric numerical integration: structurepreserving algorithms for ordinary differential equations[END_REF], with the goal to further reduce the computational cost of the simulation. In the case of mechanical systems with local stiffness, the conditional stability of an explicit time-integration scheme typically involves small time-steps for the whole system. A promising direction to mitigate this drawback consists in using a local time-stepping strategy. In the linear case, explicit high-order energy-momentum conserving methods with local time-stepping have been proposed in [START_REF] Diaz | Energy conserving explicit local time stepping for secondorder wave equations[END_REF]. In the nonlinear case, a modified Störmer-Verlet method for Hamiltonian systems containing slow and fast components is developed in [START_REF] Hairer | Geometric numerical integration: structurepreserving algorithms for ordinary differential equations[END_REF]. It is proved that this time-integrator remains symplectic, but the ratio of the fast and slow time-steps strongly influences the error on the total energy and, in general, a good balance has to be found experimentally. This phenomenon is called resonance since it is encountered for certain slow/fast ratios. Similarly, asynchronous variational integrators generally exhibit resonances when the local time-steps are close to certain rational ratios, so that ensuring stability requires adequate fitting of the local time-steps [START_REF] Fong | Stability of asynchronous variational integrators[END_REF]. In contrast, the asynchronous version of the present scheme allows one to make slow-fast time-integration while conserving a pseudo-energy (in the absence of quadrature errors). Our numerical tests show that the asynchronous scheme still exhibits second-order accuracy; a mathematical proof of this property is postponed to future work.

This paper is organized as follows. In Section 2, we present the scheme for a Hamiltonian system of interacting particles with a synchronous time-integration and establish the main properties of the scheme including second-order accuracy, time-reversibility, linear stability under a CFL condition, and pseudo-energy conservation under exact force integration. In Section 3, we test the synchronous scheme on various benchmarks from the literature including a nonlinear wave propagation problem. In Section 4, we present the slow-fast time-stepping capabilities of the asynchronous version of the scheme, together with numerical results on model particle systems connected by springs and on an inhomogeneous wave equation. These results demonstrate the efficiency gains of the asynchronous scheme with respect to the synchronous scheme.

Synchronous scheme

In this section, we present our scheme in its synchronous version and establish its main properties.

Definition of the scheme

We consider a sequence of discrete time nodes t n , n = 0, 1, . . ., with time-steps h n = t n+1 -t n and time intervals I n = [t n , t n+1 ]. The scheme is written at step n as follows: knowing p n-1/2 , q n , and [p] n , one computes

p n+1/2 = p n-1/2 + [p] n , ( 3a
)
q n+1 = q n + h n M -1 p n+1/2 , (3b) 1 2 [p] n+1 + [p] n = - In ∇V ( qn (t))dt, (3c)
with the free-flight trajectory over the time interval I n defined by

qn (t) = q n + (t -t n )M -1 p n+1/2 ∀t ∈ I n . ( 4
)
Here, [p] n represents the jump of the momentum vector at time t n , q n the position vector at time t n , and p n+1/2 is the momentum vector between t n and t n+1 . We observe that q n+1 = qn (t n+1 ). We initialize the scheme as follows:

p -1/2 = p(t 0 ), q 0 = q(t 0 ), [p] 0 = 0, (5)
where q(t 0 ), p(t 0 ) are the given position and momentum vectors at the initial time t 0 . The scheme (3) can alternatively be written as the following 2-step scheme without jumps: knowing p n-1/2 , q n , and p n+1/2 , one computes

q n+1 = q n + h n M -1 p n+1/2 , (6a) 1 2 p n+3/2 -p n-1/2 = - In ∇V ( qn (t))dt, (6b)
with the free-flight trajectory defined by [START_REF] Diaz | Energy conserving explicit local time stepping for secondorder wave equations[END_REF]. The initialization of the scheme, equivalent to [START_REF] Fetecau | Nonsmooth Lagrangian mechanics and variational collision integrators[END_REF], is as follows:

p -1/2 = p(t 0 ), q 0 = q(t 0 ), p 1/2 = p(t 0 ). ( 7
)
This initialization is tailored to achieve exact pseudo-energy conservation under exact force integration, as shown in Theorem 1 below. Other choices for the initialization are possible, for instance using a one-step method.

In the numerical implementation of the scheme, the integral in (3c) (or in (6b)) is usually not computed exactly but with a quadrature of the form

Q n (f (t); t n ; t n+1 ) = h n I i=0 ω i f (λ i t n + (1 -λ i )t n+1 ) ≈ In f (t)dt, (8)
where the real numbers ω i are the weights and the real numbers λ i ∈ [0, 1] define the quadrature points. Applying the quadrature componentwise for the calculation of the forces and exploiting that the position of the particles varies linearly in time during the free flight, we obtain

Q n ∇V ( qn (t)); t n ; t n+1 = h n I i=0 ω i ∇V (λ i q n + (1 -λ i )q n+1 )) ≈ In ∇V ( qn (t))dt, (9)
and we replace (3c) with

1 2 [p] n+1 + [p] n = -Q n ∇V ( qn (t)); t n ; t n+1 , ( 10
)
and a similar modification for (6b). In what follows, we assume that the quadrature is symmetric: ∀i ∈ {0, . . . , I}, ω i = ω I-i and λ i = 1 -λ I-i , [START_REF] Hairer | Geometric numerical integration: structurepreserving algorithms for ordinary differential equations[END_REF] and at least of order two (i.e., that the quadrature integrates exactly affine polynomials). We also assume that V is of class C 2 , i.e., V ∈ C 2 (R dN ; R). This implies that

Q n ∇V ( qn (t)); t n ; t n+1 = In ∇V ( qn (t))dt + O(h 3 n ). ( 12
)

Properties of the scheme

We now establish various properties of the scheme: pseudo-energy conservation (in the absence of quadrature errors), symmetry (or time-reversibility), second-order accuracy, and linear stability (with constant time-step).

Theorem 1 (Pseudo-energy conservation). Assume that the quadrature is exact. Then, the scheme (3) exactly conserves the following pseudo-energy:

Hn := V (q n ) + 1 2 p n-1/2 T M -1 p n+1/2 . ( 13
)
Moreover, denoting H 0 := H(q(t 0 ), p(t 0 )) the value of the exact Hamiltonian at the initial time time t 0 , we have Hn = H 0 for all n ≥ 0 if the scheme is initialized using [START_REF] Fetecau | Nonsmooth Lagrangian mechanics and variational collision integrators[END_REF].

Proof. Using (3.a), (4), and the chain rule, we obtain

d dt (V ( qn (t))) = ∇V ( qn (t)) • ( qn ) (t) = M -1 p n+1/2 T ∇V ( qn (t)) .
Integrating in time and using (3c) and the symmetry of M , we infer that

V (q n+1 ) -V (q n ) = p n+1/2 T M -1 t n+1 t n ∇V ( qn (t))dt = -p n+1/2 T M -1 1 2 [p] n+1 + [p] n = -p n+1/2 T M -1 1 2 p n+3/2 -p n-1/2 .
This leads to

V (q n+1 ) + 1 2 p n+1/2 T M -1 p n+3/2 = V (q n ) + 1 2 p n-1/2 T M -1 p n+1/2 ,
showing that Hn+1 = Hn , thereby proving the first assertion. Finally, using the initialization (5), we obtain H0 = H 0 , and this concludes the proof.

Remark 2 (Quadratures). In practice, the integral in Equation (3c) can be computed exactly only for polynomial potentials V . For instance, using the n-point Gauss-Lobatto (resp., Gauss-Legendre) quadrature, polynomials of degree up to 2n -3 (resp., 2n -1) are integrated exactly. The use of quadratures instead of exact integration for general nonlinear potentials entails only an approximate conservation of the pseudo-energy. Since the scheme is second-order accurate (see Theorem 5 below), we expect that pseudo-energy conservation is second-order accurate at best:

Hn = H0 + O(h 2 ). ( 14
)
where h := sup n h n . The order of the quadrature has an influence on the multiplicative constant in O(h 2 ), with the constant (swiftly) decreasing when increasing the quadrature order. Numerical results are presented on a nonlinear wave propagation problem in Section 3.3.

Remark 3 (Momentum conservation). Let (1, ..., 1) be the vector of size dN filled with ones. Assume that the system is isolated, i.e., (1, .., 1) T • ∇V (q) = 0 for all q ∈ R dN . Then, the total momentum, defined as P n+1/2 := (1, .., 1) T • p n+1/2 for all n ≥ 0, is conserved. This follows by taking the product of (1, ..., 1) T with Equation (6b) and using the null initialisation of the momentum jump which follows from Equation [START_REF] Fetecau | Nonsmooth Lagrangian mechanics and variational collision integrators[END_REF].

Proposition 4 (Symmetry). If the quadrature ( 9) is exact or symmetric, then the scheme [START_REF] Fong | Stability of asynchronous variational integrators[END_REF] is symmetric (or time-reversible).

Proof.

Let Y n = (q n , p n-1/2 +p n+1/2 2
, p n+1/2 -p n-1/2 ) T . Since we are going to consider positive and negative time-steps in this proof, we denote by sign(h n ) the sign of the timestep. The numerical scheme can be written as

Y n+1 = Φ hn (Y n ), where for a generic column vector Y = (Y 1 , Y 2 , Y 3 ) T , we have Φ hn (Y ) =             Y 1 + h n M -1 Y 2 + sign(h n ) 2 Y 3 Y 2 -Q n ∇V Y 1 + tM -1 Y 2 + sign(h n ) 2 Y 3 ; 0; h n -Y 3 -2 sign(h n )Q n ∇V Y 1 + tM -1 Y 2 + sign(h n ) 2 Y 3 ; 0; h n            
, where we used the invariance by translation of the quadrature Q n . Therefore, we have

Φ -hn (Y ) =             Y 1 -h n M -1 Y 2 - sign(h n ) 2 Y 3 Y 2 -Q n ∇V Y 1 + tM -1 Y 2 - sign(h n ) 2 Y 3 ; 0; -h n -Y 3 + 2 sign(h n )Q n ∇V Y 1 + tM -1 Y 2 - sign(h n ) 2 Y 3 ; 0; -h n            
.

It remains to verify that Φ hn • Φ -hn (Y ) = Y using that the quadrature is symmetric or exact. To fix the ideas, we assume that h n > 0. Let us write Y = Φ -hn (Y ) and

Y = Φ hn (Y ). Since Y 1 = Y 1 -h n M -1 (Y 2 -1 2 Y 3 ) and Y 2 + 1 2 Y 3 = Y 2 -1 2 Y 3 , we infer that Y 1 = Y 1 + h n M -1 Y 2 + 1 2 Y 3 = Y 1 . ( 15
)
For the second component, we obtain

Y 2 = Y 2 -Q n ∇V Y 1 + (t -h n )M -1 Y 2 - 1 2 Y 3 ; 0; h n = Y 2 -Q n ∇V Y 1 + tM -1 Y 2 - 1 2 Y 3 ; -h n ; 0 = Y 2 + Q n ∇V Y 1 + tM -1 Y 2 - 1 2 Y 3 ; 0; -h n = Y 2 ,
where we used [START_REF] Kane | Symplectic-energy-momentum preserving variational integrators[END_REF] in the first line and invariance by translation and symmetry of the quadrature Q n in the second and third lines respectively. The proof that Y 3 = Y 3 follows along similar lines.

Theorem 5 (Consistency error). Assume that V ∈ C 2 (R dN ; R). If the quadrature ( 9) is exact or at least of order two, the scheme (3) has second-order accuracy in time.

Proof. Let q(t), p(t) be the exact solution to [START_REF] Chabassier | Introduction and study of fourth order theta schemes for linear wave equations[END_REF]. Let us consider the column vector

Y (t n ) = (q(t n ), p(t n-1/2 )+p(t n+1/2 ) 2 , p(t n+1/2 ) -p(t n-1/2 )) T .
The consistency error is defined as

η n+1 := Y (t n+1 ) -Φ hn (Y (t n )) h n ,
where Φ hn is defined in the previous proof. Let us denote by η n+1

1 , η n+1 2 , η n+1 3
the three components of the consistency error. We have

h n η n+1 1 = q(t n+1 ) -q(t n ) -h n M -1 p(t n+1/2 ), h n η n+1 2 = p(t n+3/2 ) -p(t n-1/2 ) 2 + Q n ∇V ( qn (t)); t n ; t n+1 , h n η n+1 3 = 2h n η n+1 2 ,
where qn

(t) = q(t n ) + M -1 p(t n+1/2 )(t -t n ).
Using a Taylor expansion and the equation q(t) = M -1 p(t), we infer that

h n η n+1 1 = h n q(t n+1/2 ) -h n M -1 p(t n+1/2 ) + O(h 3 n ) = O(h 3 n ).
Moreover, since the quadrature is of second-order (at least) so that it can be replaced by the mid-point quadrature up to O(h 3 n ), and using the equations q(t) = M -1 p(t) and ṗ(t) = -∇V (q)(t), we obtain

h n η n+1 2 = h n ṗ(t n+1/2 ) -h n ∇V ( q(t n+1/2 )) + O(h 3 n ) = h n ṗ(t n+1/2 ) -h n ∇V (q(t n ) + 1 2 h n q(t n+1/2 )) + O(h 3 n ) = h n ṗ(t n+1/2 ) -h n ∇V (q(t n+1/2 )) + O(h 3 n ) = O(h 3 n ).
We conclude that η n+1 = O(h 2 n ), i.e., the scheme is second-order accurate. Proposition 6 (Linear stability). Assume that the potential V is quadratic with a positive definite Hessian H := D 2 V . Let λ be the largest eigenvalue of H. Let µ > 0 be the smallest eigenvalue of M . Then the scheme (3) is conditionally stable for a constant time-step h under the following CFL condition:

h < 2 µ λ . ( 16
)
Proof. Since the potential V is quadratic, the dynamical system (2) is linear. Let Z n be the column vector such that Z n = q n , p n-1/2 , p n+1/2 T . Adding a linear functional to V does not change the nature of the Hamiltonian system. We thus consider ∇V (0) = 0 and V (0) = 0. Since ∇V (q) is by assumption linear in q, we have ∇V ( qn

(t)) = ∇V (q n ) + (t -t n )HM -1 p n+1/2 , so that In ∇V ( qn (t))dt = hHq n + 1 2 h 2 HM -1 p n+1/2 .
Therefore, the scheme (3) can be written as

Z n+1 = AZ n with A =    I dN 0 dN hM -1 0 dN 0 dN I dN -2hH I dN -h 2 HM -1    .
The matrix M being symmetric definite positive, its square root M 1/2 is well-defined. We then observe that

à =    M 1/2 0 dN 0 dN 0 dN M -1/2 0 dN 0 dN 0 dN M -1/2    A    M -1/2 0 dN 0 dN 0 dN M 1/2 0 dN 0 dN 0 dN M 1/2    =    I dN 0 dN hI dN 0 dN 0 dN I dN -2hS I dN -h 2 S    ,
where we introduced the symmetric positive definite matrix S = M -1/2 HM -1/2 . Up to an adequate change of variable for each of the coordinates, it is possible to assume that S is diagonal. Denoting (σ i ) 1≤i≤dN the eigenvalues of S and scaling the momenta in Z n by the factors ( √ 2σ i ) 1≤i≤dN , -Ã is block diagonal in the following matrices of order 3, for all i ∈ {1, . . . , dN }:

a i =    -1 0 -h √ 2σ i 0 0 -1 h √ 2σ i -1 h 2 σ i    . The characteristic polynomial χ a i of a i is χ a i (X) = (X -1)(X 2 -X(h 2 σ i -2) + 1), which
shows that 1 is an eigenvalue of a i . Moreover, the polynomial X 2 -X(h 2 σ i -2) + 1 is positive as long as h < 2 √ σ i , and the two complex conjugate eigenvalues, written b i and b i , have a modulus equal to 1 and a nonzero imaginary part. Thus, the three eigenvalues are distinct, and the 3×3 matrix a i can be diagonalized for all i ∈ {1, . . . , dN }. Let P i the matrix such that P -1 i a

i P i = Diag(1, b i , b i ).
Then, writing P -1 the block-diagonal matrix composed of the elementary matrices (P -1 i ) 1≤i≤dN and D the block-diagonal matrix composed of the diagonal matrices of eigenvalues of (a i ) 1≤i≤dN , we infer that, for all n ∈ N:

Ãn ≤ P -1 P D n ≤ P -1 P ,
because the diagonal matrix D has diagonal entries of modulus 1, and thus D = 1. Hence, for all n ∈ N, we obtain Z n ≤ C Z 0 , for a constant C independent of n. Since the eigenvalues σ i of S are positive and smaller than λ µ , we conclude that linear stability holds true under the CFL condition [START_REF] Kane | Variational integrators and the newmark algorithm for conservative and dissipative mechanical systems[END_REF].

Remark 7 (Comparison with Störmer-Verlet).

A possible writing of the Störmer-Verlet method is the following:

q n+1 = q n + h n M -1 p n+1/2 , p n+3/2 = p n+1/2 -h n+1 ∇V (q n+1 ).
Both the present scheme and the Störmer-Verlet scheme are of leapfrog-type, have a similar CFL condition for linear stability, and are second-order accurate. The main difference is that, using a mid-point quadrature, the forces used to update the momenta at t n+1 are computed at t n+1/2 with the present scheme, i.e., 1 2 (p n+3/2 -p n-1/2 ) = -h n ∇V ( q(t n+1/2 )), whereas the momentum update can be rewritten as

1 2 (p n+3/2 -p n-1/2 ) = - 1 2 h n+1 ∇V (q n+1 ) + h n ∇V (q n ) ,
in the Störmer-Verlet scheme. Moreover, the origin of energy conservation is different for the two schemes. The Störmer-Verlet scheme is energy-conserving only for constant time-steps due to its symplecticity. The present scheme enjoys an algebraic pseudo-energy preservation property for every time-step (constant or not) up to quadrature errors.

Remark 8 (Adaptive time-stepping for discrete energy control). The conservation of the pseudo-energy does not imply stability since p n-1/2 T M -1 p n+1/2 does not have a sign a priori. However, defining the discrete energy

H n := V (q n ) + 1 8 p n-1/2 + p n+1/2 T M -1 p n-1/2 + p n+1/2 , ( 17
)
a straightforward calculation shows that

H n = Hn + 1 8 ([p] n ) T M -1 [p] n = Hn + 1 8 |M -1/2 [p] n | 2 . ( 18
)
This implies that 0 ≤ H n -Hn ≤ O(h 2 n ), where we used the identity [START_REF] Krysl | Explicit Newmark/Verlet algorithm for time integration of the rotational dynamics of rigid bodies[END_REF] for the lower bound and we invoked Theorem 5 for the upper bound. One can use the identity [START_REF] Krysl | Explicit Newmark/Verlet algorithm for time integration of the rotational dynamics of rigid bodies[END_REF] during the computations for an on-the-fly monitoring of possible departures of the conserved pseudo-energy Hn from the discrete energy H n . The idea is to check whether 1 8

|M -1/2 [p] n | 2 ≤ fly
Hn after every momenta computation and to halve the time-step if this bound is not met (note that the momentum jumps converge to zero with the time-step). The benefits of such an adaptive time-stepping strategy are illustrated in Section 3.3.

Numerical results

In this section, we present numerical results for the scheme (3). We consider classical benchmarks from the literature and a nonlinear wave equation from [START_REF] Chabassier | Energy preserving schemes for nonlinear Hamiltonian systems of wave equations: Application to the vibrating piano string[END_REF].

Convergence study

We perform a convergence study with a single particle in dimension d = 1. The reference solution is q(t) = sin(t) 4 + 1, and the corresponding potential energy is

V (q(t)) = 8 (q(t) -1) 3/2 -(q(t) -1) 2 .
We apply the scheme (3) to this Hamiltonian system over 10 3 seconds using the mid-point quadrature as well as the three-and five-point Gauss-Lobatto quadratures of order 3 and 7, respectively, for the integration of the forces. We report the 1 -error with respect to the reference solution (the sum of the errors at the discrete time nodes divided by the number of time-steps) in Figure 1 as a function of the number of force evaluations. We observe that for the three quadratures, the convergence is of second order as expected. The quadrature order does not impact the convergence rate but has an influence on the computational efficiency. We note that in this case, the mid-point quadrature is more efficient than the three-and five-point Gauss-Lobatto quadratures of order 3 and 7. The Hamiltonian is

H(q, p) = 1 2 m i=1 (p 2 2i-1 + p 2 2i ) + ω 2 4 m i=1 (q 2i -q 2i-1 ) 2 + m i=0 (q 2i+1 -q 2i ) 4 , q 1 q 2 • • • q 2m-1 q 2m stiff harmonic soft nonlinear
x i = (q 2i + q 2i-1 )/ √ 2, y i = (p 2i + p 2i-1 )/ √ 2, x m+i = (q 2i -q 2i-1 )/ √ 2, y m+i = (p 2i -p 2i-1 )/ √ 2,
the Hamiltonian can be rewritten as

H(x, y) = 1 2 2m i=1 y 2 i + ω 2 2 m i=1 x 2 m+i + 1 4 (x 1 -x m+1 ) 4 + m-1 i=1 (x i+1 -x m+i+1 -x i -x m+i ) 4 + (x m + x 2m ) 4 . ( 19
)
As the system is Hamiltonian, the total energy of the system should be conserved by the numerical scheme. The Fermi-Pasta-Ulam system has yet another quasi-invariant. Letting I j (x m+j , y m+j ) = 1 2 y 2 m+j + ω 2 x 2 m+j be the oscillatory energy of the jth stiff spring, the total oscillatory energy I = I 1 + I 2 + • • • + I m is close to a constant value as proved in [11, p.22]: I(x(t), y(t)) = I(x(0), y(0)) + O(ω -1 ).

In our numerical experiment, we set m = 3 and ω = 50. Figure 3a (left panel) shows the variation of the oscillating energies and of the pesudo-energy Hn over time for a constant time-step h = 10 -3 . The energy exchange between the oscillatory modes is remarkably similar to the reference solution given in [START_REF] Hairer | Geometric numerical integration: structurepreserving algorithms for ordinary differential equations[END_REF]Chap. I.4] and represented in Figure 3b (right panel). The reference solution was computed with high accuracy using a Runge-Kutta 4 integrator with a time-step of h = 10 -4 . In particular, the total oscillatory energy I displays fast oscillations around a fixed constant. The conservation of energy is verified up to machine precision, even with a mid-point quadrature. The results being already very satisfactory, the results computed with higher order quadratures are omitted for brevity. A more detailed study of the influence of the order of quadrature on pseudo-energy conservation is presented in the next section.

Nonlinear wave equation

The setting comes from [START_REF] Chabassier | Energy preserving schemes for nonlinear Hamiltonian systems of wave equations: Application to the vibrating piano string[END_REF]. The interval Ω = [0, 1] represents a one-dimensional string. Let V : R d → R be the potential energy, with dimension d = 2. It is assumed that V verifies the following conditions:

• Smoothness: V is of class C 2 ;
• Convexity: V is strictly convex; • Coercivity: ∃K > 0 so that V (u) ≥ K|u| 2 for all u ∈ R 2 ;

• Boundedness:

∃M > 0 so that |∇V (u)| 2 ≤ M min(V (u), (1 + |u| 2 )) for all u ∈ R 2 .
The problem of interest is to find

u : Ω × R + → R 2 such that        ∂ 2 tt u -∂ x (∇V (∂ x u)) = 0, u(0, t) = 0, u(1, t) = 0, u(x, 0) = u 0 (x), ∂ t u(x, 0) = v 0 (x), (20)
with given initial conditions u 0 : Ω → R 2 and v 0 : Ω → R 2 . For a pair (u 1 , u 2 ) ∈ R 2 , the functional V takes the following value:

V (u 1 , u 2 ) = u 2 1 + u 2 2 2 -α (1 + u 1 ) 2 + u 2 2 -(1 + u 1 ) ,
where the parameter α ∈ [0, 1) is related to the tension of the string, such that the string behaviour is nonlinear when α > 0 and the strength of the nonlinearity increases with α.

The following variational formulation in

V := H 1 0 (Ω; R 2 ) is considered: d 2 dt 2 Ω u • v + Ω ∇V (∂ x u) • ∂ x v = 0, ∀v ∈ V, ∀t > 0.
We use H 1 -conforming P 1 Lagrange finite elements for the space discretization. Let N be the number of nodes discretizing the string and (ϕ i ) 1≤i≤2N be the nodal basis functions associated with the degrees of freedom of the string in the two directions. These basis functions span the finite-dimensional subspace V N V. The space semi-discrete function approximating the exact solution is u N (t) = 2N i=1 q i (t)ϕ i (x) ∈ V N and solves the following space semi-discrete problem:

d 2 dt 2 Ω u N • v N + Ω ∇V (∂ x u N ) • ∂ x v N = 0, ∀v N ∈ V N , ∀t > 0.
Introducing the vector q = (q 1 , . . . , q 2N ) ∈ R 2N , the following Hamiltonian system has to be integrated in time:

H(q, p) = 1 2 p T M -1 p + V (q), V (q) = Ω V 2N i=1 q i ∂ x ϕ i ,
where M is the classical P 1 Lagrange finite element mass matrix. Assuming that all the components of q associated with the first direction are enumerated first and then those associated with the second direction, the matrix M is a 2 × 2 block-diagonal matrix and each diagonal block is a tridiagonal matrix of size N ×N equal to ∆x tridiag(1/6, 2/3, 1/6).

In our numerical experiments, we consider the values α = 0 (which corresponds to the linear case), α = 0.8 (which corresponds to a mildly nonlinear behavior), and α = 0.99 (which corresponds to a strongly nonlinear behavior). The space discretisation is such that ∆x = 0.01 and thus N = 99 basis functions are used in each direction. The time-step is ∆t = 0.0033. Using the same space discretisation, the greatest stable constant timestep has been found to be ∆t max = 0.0055. Three numerical simulations are performed in every case by letting the amplitude of the initial condition u 0 at time t 0 be 0.01, 0.1, or 0.3. The initial velocity at time t 0 is always taken to be zero. The results are reported in Figure 4 where in all cases, a mid-point quadrature is used. Six snapshots of the R 2valued deformation vector u 0 (x) + u(x, t) of the string over one second are represented horizontally in various colors; specifically, at each snapshot in time, the deformation vector is plotted in the corresponding vertical plane. The role played by the nonlinearity can be observed in the fact that the amplitude of u 0 influences the vibration of the string. The tension which causes nonlinearity also changes the wave celerity. We observe an excellent agreement between the present results and the results reported in [START_REF] Chabassier | Energy preserving schemes for nonlinear Hamiltonian systems of wave equations: Application to the vibrating piano string[END_REF]. The time-variation of the discrete pseudo-energy Hn defined by ( 13) and of the discrete energy H n defined by [START_REF] Kane | Finite element analysis of nonsmooth contact[END_REF] are shown in Figure 5a in the most challenging case where α = 0.99 and u 0 = 0.3. We first observe that the variations the discrete energy H n are very moderate. In Figure 5b, we illustrate the adaptive time-stepping strategy discussed in Remark 8, where we take fly = 0.03% to control the departure of Hn from H n at each iteration. In this situation, the adaptive time steps take values in the range [0.0024, 0.005]. Concerning the discrete pseudo-energy Hn , we observe conservation up to machine precision when employing a five-point Gauss-Legendre quadrature of order 9. To illustrate the impact of quadratures on pseudo-energy conservation, we perform two numerical experiments. First, Table 1 reports the maximal variation of the discrete pseudo-energy Hn depending on the quadrature used with a constant time-step ∆t = 0.0033. The mid-point quadrature is precise enough when used on the linear equation (α = 0). The three-point Gauss-Legendre quadrature of order 5 is found to give very satisfactory results for the two nonlinear cases (α = 0.8 and α = 0.99). As announced in Remark 2, the maximal error on pseudo-energy conservation is observed to decrease when increasing the quadrature order. The total number of force evaluations is 300, 900, and 1500 when using the mid-point quadrature and the three-and five-point Gauss-Legendre quadratures of order 5 and 9, respectively. In the second experiment, we illustrate the second-order accuracy of pseudo-energy conservation when using the mid-point quadrature. We consider again the most challenging case where α = 0.99 and u 0 = 0.3. Figure 6 shows the maximal variation of the discrete pseudo-energy Hn with respect to the value of the (fixed) time-step used in the simulation, confirming the second-order accuracy. To conclude this section, we present some comments on the relative costs of the present scheme with respect to an implicit scheme, e.g., the one devised in [START_REF] Chabassier | Energy preserving schemes for nonlinear Hamiltonian systems of wave equations: Application to the vibrating piano string[END_REF]. In the most challenging case where α = 0.99 and u 0 = 0.3, the energy conservation in [START_REF] Chabassier | Energy preserving schemes for nonlinear Hamiltonian systems of wave equations: Application to the vibrating piano string[END_REF] is machine error. For a comparable error, we can consider the present scheme with the five-point Gauss-Legendre quadrature of order 9. The convergence criterion of the Newton's method in the implicit scheme can be estimated to require at least a couple of iterations per timestep and the same number of Hessian computations and global matrix inversions per time-step, whereas the explicit method with the five-point Gauss-Legendre quadrature of order 9 requires only 5 force evaluations per time-step. Thus, although no general conclusions can be drawn, the present explicit method stands good chances to be quite competitive with respect to an implicit method.

Asynchronous scheme

Owing to the CFL condition ( 16), the time-step can be required to be small in regions with stiff or nonsmooth dynamics. The overall efficiency of the computation would be compromised by the large number of integral calculations in the whole domain, while most of these would be redundant in smooth regions. We therefore propose an asynchronous version of the scheme which preserves the general properties of the synchronous version. In this section, we first present the idea behind slow-fast decomposition of the particles and we devise an asynchronous scheme for which we prove pseudo-energy conservation at the slow time nodes under exact force integration. Second-order accuracy is expected and is illustrated numerically on two test cases including an inhomogeneous wave equation.

Slow-fast splitting

In order to simplify the presentation of the asynchronous scheme, we limit ourselves here to the integration of a slow-fast dynamics, i.e., we consider a system with essentially two distinct time scales. The forces between the particles are supposed to be split into a "fast" set with an associated time-step h F and a "slow" set with an associated time-step h S > h F . For example, the splitting can result from the relative stiffness of the forces in the system. Consequently, the particles are split into three sets: the slow particles are subjected only to slow forces, the fast particles are subjected only to fast forces, and the remaining particles, which are called mixed particles, are subjected to both slow and fast forces. This definition means that the slow particles do not interact with the fast particles directly, so that the potential V can be decomposed as follows:

V (q) = V S (q S ) + V M (q M , q S ) + V F (q F , q M ),
where q F , q S , and q M denote respectively the positions of the fast, slow, and mixed particles, the potential V S describes the interactions between slow particles, V M the interactions between slow and mixed particles, and V F the interactions between mixed and fast particles (or between themselves). For instance, the purple particle in Figure 7 and the Particle 3 in Figure 8 are mixed particles. The mixed particle in Figure 7 is subjected to a "fast" force by the stiff spring on its right and a "slow" force by the soft spring on its left. The fast particle (in red) is only subjected to a "fast" force by the stiff spring. The slow particles (in blue) are only subjected to "slow" forces by the two soft springs. Finding a slow-fast decomposition is not possible for every Hamiltonian system. For example, in the case where all the particles interact with each other, no slow-fast splitting is available. The most favorable configuration is the one where the slow and the fast particles interact essentially among themselves and have very few interactions with mixed particles. This configuration is encountered in inhomogeneous problems where an interface separates two zones where the properties are different; the slow and the fast particles are then located in the two zones, whereas the mixed particles are located at the interface.

In what follows, we abuse the notation by denoting F , M and S the sets collecting the indices in {1, . . . , N } of the fast, mixed and slow particles, respectively. For simplicity, we assume that the mass matrix M is diagonal and denote M F , M M and M S the restriction of M to the F , M and S particles respectively. Still for simplicity, we assume that both time-steps h S and h F are kept constant.

Presentation of the asynchronous scheme

Without much loss of generality, we can suppose that the slow and fast time-steps are commensurate so that h S = Kh F with K ∈ N * . We then define the coarse time nodes t n = nh S and the fine time nodes t n,m = t n + mh F for all m ∈ {0, . . . , K}. The asynchronous scheme consists in integrating K times the dynamics of the F and M particles with the "fast" forces computed at each time-step of length h F and in updating the S particles with the "slow" forces computed once at the end of each time-step of length h S . The general procedure is depicted in Figure 8 for four particles in the same configuration as in Figure 7. The efficiency of the asynchronous scheme hinges on the fact that each particle has a free-flight movement during each time-step, with the neighbouring particle forces acting only at the end of the time-step.

Soft Soft Stiff Particle 1 Slow • Particle 2 Slow • Particle 3 Mixed • Particle 4 Fast • h S h F t n t n+1 t n,0 t n,1 t n,2 t n,3 t n,4 h S 0 ∇V S ( qn S ) h F 0 ∇V M ( qn S , qn,0 M ) 2h F h F ∇V M ( qn S , qn,1 M ) 3h F 2h F ∇V M ( qn S , qn,2 M ) 4h F 3h F ∇V M ( qn S , qn,3 M ) h F 0 ∇V F ( qn,0 F , qn,0 M ) 2h F h F ∇V F ( qn,1 F , qn,1 M ) 3h F 2h F ∇V F ( qn,2 F , qn,2 M ) 4h F 3h F ∇V F ( qn,3 F , qn,3 M )
Figure 8: Asynchronous integration of four particles with a slow-fast synamics, h S = 4h F Let us now describe in more detail the asynchronous scheme over the coarse time interval I n = [t n , t n+1 ]. At the beginning, we have at our disposal the triple (p

n-1/2 i , q n i , p n+1/2 i
) for the slow particles (i ∈ S) and the triple (p

n,-1/2 i = p n-1,K-1/2 i , q n,0 i , p n,1/2 i
) for the fast and the mixed particles (i ∈ F ∪ M ). The asynchronous scheme then proceeds as follows (we use here the two-step formulation which reduces to [START_REF] Fong | Stability of asynchronous variational integrators[END_REF] in the synchronous case):

• For the fast particles (i ∈ F ), one computes for all m ∈ {0, . . . , K -1},

q n,m+1 i = q n,m i + h F 1 m i p n,m+1/2 i , ( 21a
) p n,m+3/2 i = p n,m-1/2 i -2 t n,m+1 t n,m ∂V F ∂q i ( qn,m F (t), qn,m M (t))dt, (21b)
with the free-flight trajectories for the fast and the mixed particles defined as qn,m

j (t) = q n,m j + 1 m j p n,m+1/2 j (t -t n,m ), ∀t ∈ [t n,m , t n,m+1 ], ∀j ∈ F ∪ M. ( 22
)
• For the mixed particles (i ∈ M ), one computes for all m ∈ {0, . . . , K -1}, the position q n,m+1 i as in (21a), whereas equation (21b) is replaced by

p n,m+3/2 i = p n,m-1/2 i -2 t n,m+1 t n,m ∂V F ∂q i ( qn,m F (t), qn,m M (t)) + ∂V M ∂q i ( qn,m M (t), qn S (t)) dt, ( 23
)
where the free-flight trajectories of the slow particles are computed over the coarse time interval as follows:

qn j (t) = q n j + 1 m j p n+1/2 j (t -t n ), ∀t ∈ [t n , t n+1 ], ∀j ∈ S. (24)
• For the slow particles (i ∈ S), one computes

q n+1 i = q n i + h S 1 m i p n+1/2 i , ( 25a
)
p n+3/2 i = p n-1/2 i -2 K-1 m=0 t n,m+1 t n,m ∂V M ∂q i ( qn,m M (t), qn S (t))dt -2 t n+1 t n ∂V S ∂q i ( qn S (t))dt, (25b)
with the free-flight trajectories defined above.

Note that the slow forces between slow and mixed particles need to be evaluated at every fine time-step. In the worst case scenario, every slow force links a slow particle with a mixed particle, which results in the asynchronous scheme reverting to the synchronous scheme. Such a case typically occurs when the particles all interact or when the system alternates fast and slow forces. On the other hand, the efficiency of the asynchronous scheme compared to the synchronous scheme is maximal in the case where the mixed particles constitute a small fraction of the particles and their interaction is limited to a small fraction of the slow particles. A typical case is a nearest-neighbour interaction with slow and fast particles located in distinct regions, the mixed particles being confined in a lower dimensional delimiting interface. In the limit of a large number of particles, the computational cost per large time-step h S reduces to K integrals of the fast forces and one integral of the slow forces.

Proposition 9 (Synchronization of particles). Assume that the numerical integration is exact. Then the numerical scheme ( 21)- [START_REF] Wohlmuth | Variationally consistent discretization schemes and numerical algorithms for contact problems[END_REF] exactly conserves the following pseudo-energy at the coarse time nodes t n :

Hn = V S (q n S ) + V M (q n,0 M , q n S ) + V F (q n,0 F , q n,0 M ) + i∈S 1 2m i p n-1/2 i T p n+1/2 i + i∈F ∪M 1 2m i p n,-1/2 i T p n,1/2 i . ( 26
) Proof. Let us set Hn S =V S (q n S ) + i∈S 1 2m i p n-1/2 i T p n+1/2 i , Hn,m F M =V F (q n,m F , q n,m M ) + V M (q n,m M , qn S (t n,m )) + i∈F ∪M 1 2m i p n,m-1/2 i T p n,m+1/2 i
, for all m ∈ {0, . . . , K}, so that Hn = Hn S + Hn,0 F M . Following the same calculations as in the proof of Theorem 1 for equation [START_REF] Wohlmuth | Variationally consistent discretization schemes and numerical algorithms for contact problems[END_REF], we infer that

Hn+1 S = Hn S - K-1 m=0 t n,m+1 t n,m ∂V M ∂q S ( qn,m M (t), qn S (t)) • M -1 S p n+1/2 S dt.
Similarly, for all m ∈ {0, . . . , K -1}, using ( 21) and ( 23), we have Hn,m+1 Remark 10 (Asynchronous pseudo-energy conservation). The pseudo-energy Hn of Theorem 1 is not conserved after every integration over a fast time-step h F in the asynchronous setting. This results from the fact that during a "slow" time-step h S , the effect of forces has been taken into account for the "fast" particles but not for the "slow" particles.

F M = Hn,m F M + t n,m+1

Numerical results

In this section, we present numerical results on the asynchronous scheme. We first consider a variant of the Fermi-Pasta-Ulam system with a slow-fast dynamics and then an inhomogeneous wave propagation problem.

Fermi-Pasta-Ulam system with slow-fast dynamics

We propose a slight variation of the Fermi-Pasta-Ulam test case in order to assess the efficiency of the asynchronous scheme. Contrary to the usual setting where stiff and soft springs alternate, we suppose here that the system is composed of one stiff region and one soft region, delimited by an interface in the middle of the domain. Figure 9 illustrates the setting. There are (m -1) fast particles, 1 mixed particle, and m slow particles.

q 1 q 2 q m • • • q 2m-1 q 2m
stiff harmonic soft nonlinear Figure 9: Setting for the Fermi-Pasta-Ulam system with slow-fast dynamics We consider a problem in dimension d = 1. The Hamiltonian is given by

H(p, q) = 1 2 2m i=1 p 2 i + ω 2 4 m i=1 (q i -q i-1 ) 2 + 2m i=m (q i+1 -q i ) 4 .
In the present experiment, we take m = 3 and ω 2 = 10. The fast forces being generated by stiff linear springs, the fast time-step h F should respect the CFL condition from Equation [START_REF] Kane | Variational integrators and the newmark algorithm for conservative and dissipative mechanical systems[END_REF], which here leads to h F < 3 • 10 -3 . The small time-step h S being controlled by soft nonlinear springs, the CFL condition ( 16) is not applicable. A constant stable time-step has been found empirically to be h S ≤ 10 -1 using the five-point Gauss-Lobatto quadrature of order 7. The dynamics of the particles is presented in Figure 10 for h S = 0.01 and h F = 2 • 10 -4 , so that 50 iterations of the fine time-step are carried out for each iteration of the coarse time-step. Observe that, as expected, the fast particles (1 ≤ i ≤ m) exhibit oscillations with a typical frequency ω, whereas the slow particles (m + 1 ≤ i ≤ 2m) have tame nonlinear oscillations with a frequency smaller than 1. Figure 11 shows that the conservation of the discrete pseudo-energy Hn defined by ( 13) is as perfect for the asynchronous scheme as for the synchronous scheme with the five-point Gauss-Lobatto quadrature of order 7.

The computational cost of the scheme is proportional to the number N of force evaluations. With an n-point quadrature and a total integration time T , the numbers of force evaluations N s and N a for the synchronous and asynchronous schemes respectively on the present slow-fast problem are given by:

N s = (n -1)T 2m + 1 h F , N a = (n -1)T m + 1 h F + m h S .
Recalling that K = h S h F ≥ 1 is the number of fast steps per slow step, the cost reduction η of the asynchronous scheme with respect to the synchronous scheme is given by

η = N a N s = 1 + m (m+1)K 1 + m m+1 .
For h S = 0.01, h F = 2 • 10 -4 and T = 100, N a = 1.015 • 10 8 , to be compared with N s = 1.75 • 10 8 . As m increases,

η m→+∞ ----→ 1 + 1 K 2 .
When the number of fast subiterations K increases, η tends to 0.5, which means that the computational cost reduction of the asynchronous scheme compared to the synchronous scheme approaches 50%. This is the best-case scenario, since the computational cost is concentrated on the fast dynamics where frequent evaluations are required, whereas the slow dynamics is almost costless. In order to assess the accuracy of the asynchronous scheme, we consider the L ∞ -error of the position of the asynchronous solution with respect to the synchronous solution using the small time-step h F . Figure 12a (left panel) shows the evolution of the error as the coarse time-step h S is refined, with fixed fine time-step h F = 10 -4 . We observe a second-order convergence of the error. Figure 12b (right panel) displays the evolution of the error as the fine time-step h F is further refined, with fixed coarse time-step h S = 10 -2 . We observe that the error decreases until it reaches a plateau, which is due to the error on the slow particles. These observations confirm that reducing the fine time-step beyond h F = h S /50 does not significantly improve the error since the error is dominated by the error on the slow particles. Conversely, the error reduction due to the coarse time-step reduction is not compromised by the asynchronous scheme.

Finally, a convergence test is carried out with a constant ratio h S h F = 25 and using the five-point Gauss-Lobatto quadrature of order 7. The error is measured as previously by the L ∞ -error on the positions between the synchronous and asynchronous schemes. The results are presented in Figure 13. The synchronous method is used with a constant time-step h F . We observe second-order convergence as both time-steps are refined simultaneously.

Inhomogeneous wave propagation

As a physically relevant example of the slow-fast test case, we consider the propagation of a wave in a linear elastic material in dimension d = 1, with an inhomogeneous speed of sound. Denote the domain Ω, u 0 : Ω → R and v 0 : Ω → R initial conditions for Figure 13: Fermi-Pasta-Ulam system with slow-fast dynamics: L ∞ -error on the positions between the synchronous scheme (with time-step h F ) and the asynchronous scheme (with time-steps h S and h F having fixed ratio h S h F = 25).

reflections and transmissions occur, which result in the final state in Figures 14f and15f.

The exact solution can be expressed as follows for all t > 0:

∀x ∈ (0, 0.5), u(x, t) = k≥0 c 2 -c 1 c 1 + c 2 k (u 0 (x + k -c 1 t) -u 0 (k -x -c 1 t)), ∀x ∈ (0.5, 1), u(x, t) = 2c 1 c 1 + c 2 k≥0 c 2 -c 1 c 1 + c 2 k u 0 c 1 c 2 (x -0.5) + k + 0.5 -c 1 t .
The numerical solution matches very well the exact solution. We can observe slight overshoots near the extrema and at the tail of the peaks, especially in the slow domain. This can be explained by the fact that the space-discretization (28) is slightly dispersive so that steep variations tend to generate oscillations (similar to the Gibbs phenomenon). Figure 16 presents the behavior of the error for the asynchronous and synchronous schemes with respect to the number of force evaluations. We consider the maximal error at the end of the simulation, i.e., at t = 0.5. The results are obtained by letting h S , h F → 0 while keeping the ratio h S /h F fixed. The convergence is of order 1 with respect to the number of force evaluations. Since the number of force evaluations scales like (∆x • ∆t) -1 , this is compatible with a second-order convergence in space and in time after taking into account the CFL condition. In conclusion, the asynchronous scheme displays similar errors to the synchronous scheme, with roughly half of the number of force evaluations involved as noted in Section 4.3.1. This confirms the efficiency of the asynchronous scheme.

Conclusion

In this paper, a new explicit pseudo-energy conserving time-integration scheme has been proposed. It is capable of handling general nonlinear Hamiltonian systems and has been tested on classical numerical benchmarks and on a nonlinear wave propagation problem. The present scheme enables the use of local time-stepping strategies to circumvent stiff CFL conditions on the time-step and to enhance computational efficiency in the context of slow-fast dynamics.

Various perspectives of the present work can be considered. We believe that the timeintegration of dissipative systems should be a straightforward extension of the present scheme. Variational integrators have been proposed for dissipative systems and have proven to be able to accurately track the physical dissipation of energy [START_REF] Kane | Variational integrators and the newmark algorithm for conservative and dissipative mechanical systems[END_REF]. Other possible developments lie in the adaptation of the scheme to constrained Hamiltonian systems [START_REF] Leyendecker | Variational integrators for constrained dynamical systems[END_REF], such as mechanical contact problems [START_REF] Kane | Finite element analysis of nonsmooth contact[END_REF][START_REF] Wohlmuth | Variationally consistent discretization schemes and numerical algorithms for contact problems[END_REF] and rigid body rotations [START_REF] Krysl | Explicit Newmark/Verlet algorithm for time integration of the rotational dynamics of rigid bodies[END_REF][START_REF] Salomon | Energy-conserving algorithms for a corotational formulation[END_REF][START_REF] Mariotti | A new leapfrog scheme for rotational motion in 3d[END_REF]. Another perspective is the high-order extension of the present scheme.

Figure 1 :

 1 Figure 1: Convergence test: 1 -convergence for a single particle

Figure 2 :

 2 Figure 2: Fermi-Pasta-Ulam test case

Figure 3 :

 3 Figure 3: Fermi-Pasta-Ulam test case: (a) Energy variation, present scheme, h = 10 -3 ; (b) reference RK4 solution, h = 10 -4

Figure 4 :

 4 Figure 4: Nonlinear wave equation: Deformations of the string over time with nonlinearity parameter α = 0 (top), α = 0.8 (middle), and α = 0.99 (bottom); the amplitude of u 0 is u 0 = 0.01 (left), u 0 = 0.1 (middle), and u 0 = 0.3 (right)

Figure 6 :

 6 Figure 6: Nonlinear wave equation: Maximal error on the conservation of the pseudoenergy Hn as a function of the time-step in the case α = 0.99, u 0 = 0.3, and a mid-point quadrature

Figure 7 :

 7 Figure 7: Example of system of particles with a slow-fast splitting

Figure 10 : 4 )- 2 × 10 - 14 - 1 . 5 × 10 - 14 - 1 × 10 - 14 - 5 × 10 - 15 0 5 × 10 - 15 1 × 10 - 14 1. 5 × 10 Figure 11 :

 10421014151014110145101551015101451011 Figure 10: Fermi-Pasta-Ulam system with slow-fast dynamics: Position dynamics for the asynchronous scheme (h S = 0.01, h F = 2 • 10 -4 )

Figure 12 :

 12 Figure 12: Fermi-Pasta-Ulam system with slow-fast dynamics: Convergence of the asynchronous scheme (a) with respect to the coarse time-step h S , with fixed fine time-step h F = 10 -4 , and (b) with respect to the fine time-step h F , with fixed coarse time-step h S = 10 -2

Figure 14 :

 14 Figure 14: Inhomogeneous wave propagation: Displacement u for ∆x = 5 × 10 -4 at times (a) t = 0, (d) t = 0.1, (b) t = 0.2, (e) t = 0.3, (c) t = 0.4, (f) t = 0.5

Figure 15 :Figure 16 :

 1516 Figure 15: Inhomogeneous wave propagation: Velocity du dt for ∆x = 5 × 10 -4 at times (a) t = 0, (d) t = 0.1, (b) t = 0.2, (e) t = 0.3, (c) t = 0.4, (f) t = 0.5

Table 1 :

 1

	0.07 0.08		Relative variation of energy (in %)	H n Hn	0.07 0.08	Relative variation of energy (in %)	H n Hn
	0.06							0.06
	0.05							0.05
	0.04							0.04
	E						E
	0.03							0.03
	0.02							0.02
	0.01							0.01
	0.00							0.00
	0.0 0.01	0.2	0.4	t	0.6	0.8	1.0	0.0 0.01	0.2	0.4	t	0.6	0.8	1.0
				(a)							(b)
	α		0				0.8				0.99
	u 0	0.01 0.1 0.3	0.01	0.1	0.3	0.01		0.1	0.3
	MP	me me me 5.5e-08 5.6e-05 2.4e-05 1.5e-08 1.2e-05 1.1e-04
	GL5 me me me 1.7e-13 1.2e-14 4.4e-14 3.3e-12 6.4e-13 1.8e-13
	GL9 me me me 1.6e-13 1.2e-14	me	3.1e-12		me	me

Figure 5: Nonlinear wave equation: time-variation of the discrete energy H n and pseudoenergy Hn over a unit time interval for α = 0.99 and an amplitude of 0.3 for u 0 : (a) Fixed time-step ∆t = 0.0033 ; (b) Adaptive time-step as in Remark 8 Nonlinear wave equation: Maximal errors in the conservation of the pseudoenergy Hn for the mid-point (MP), three-point Gauss-Legendre of order 5 (GL5) and five-point Gauss-Legendre of order 9 (GL9) quadratures; 'me' means machine error

Acknowledgements

The authors would like to thank F. Legoll (University Paris-Est, Navier Laboratory) for stimulating discussions on the integration of Hamiltonian dynamics. The authors are also thankful to the anonymous referees for their insightful remarks.

displacement and velocity respectively, and u : Ω × R + → R the displacement, u follows the equations:

in Ω,

We take Ω = (0, 1) and we set

Setting N ∈ N * , ∆x = 1 N and x i = i∆x for all i ∈ {0, . . . , N }, the partial differential equation ( 27) can be semi-discretized in space with the following centered finite difference scheme (which is equivalent to a discretization using H 1 -conforming P 1 Lagrange finite elements after lumping the mass matrix):

, the ordinary differential equation in ( 28) is derived from the following Hamiltonian:

The CFL condition [START_REF] Kane | Variational integrators and the newmark algorithm for conservative and dissipative mechanical systems[END_REF] becomes

For the indices i such that x i ≤ 0.5, one must then take h ≤ 0.2∆x, while for x i > 0.5, it suffices that h ≤ 2∆x. In what follows, we therefore set the slow (resp. fast) particules as the elements i such that x i > 0.5 (resp. x i < 0.5) and define h S = ∆x and h F = 0.1∆x. The mixed particle is the particle at the interface between the fast and slow particles.

The numerical solution and the exact solution for the displacement and the velocity computed with ∆x = 5×10 -4 are presented in Figures 14 and15 respectively. The system is initialized with the functions u 0 (x) = 10 -2 e -(20(x-0.2)) 2 1 (0,0.5) (x), v 0 (x) = 8(x -0.2)e -(20(x-0.2)) 2 1 (0,0.5) (x).

The initial condition propagates to the right with the speed of sound c 1 = 10, until it reaches x = 0.5. At the boundary between the slow and fast domain, it is partly transmitted to the right with speed of sound c 2 = 1 and partly reflected with speed -c 1 . The reflected wave reflects again on the left boundary x = 0 of the domain. Successive


A Taoutioui

A Dubois

N Sisourat

A Makhoute

Classical and semiclassical non-perturbative treatments of the electron transfer process in the collisional system proton-hydrogen involving initial excited states

Keywords: electron transfer, ion-atom collision, non-perturbative treatments

We investigate the total cross sections of the electron transfer process in the protonhydrogen collisional system involving initial excited states at impact energies in the 2.5 eV /u-100 keV /u range. The calculations are based on two non-perturbative treatments, namely the semiclassical Atomic Orbitals Close-coupling (SC-AOCC) approach and the Classical Trajectory Monte Carlo (CTMC) method. The results are presented and discussed for the collisions H + + H(n = 1, 2, 3) taking into account the degeneracy of each level in orbital momentum . We find that the SC-AOCC results show empirical n 3 and n 4 scaling laws at low and high impact velocities, respectively.

Introduction

The capture process in ion-atom collisions has attracted special interest over many decades [1], due to its utmost importance in various domains, from fusion to interstellar space research [2]. For example the corresponding cross sections are key input data for plasma fusion since they are required in collisional radiative modelling and diagnostics tools, as neutral beam spectroscopy [3], as stored e.g. in the Atomic Data and Analysis Structure (OPEN-ADAS) [START_REF]Atomic data and analysis structure[END_REF].

As a benchmark system, the capture in collisions between proton and ground-state hydrogen H(1s) has been extensively investigated both experimentally [START_REF] Hvelplund | [END_REF][6][7][8][9][10][11][12][13] and with quantal, classical and semiclassical approaches developed to understand the general behavior of electron transfer processes in ion-atom collisions [14][15][START_REF] Bransden | Charge exchange and the theory of ion-atom collisions[END_REF][START_REF] Belkić | Principles of quantum scattering theory[END_REF][START_REF] Belkić | Quantum theory of high-energy ion-atom collisions[END_REF].

On the other hand the electron capture from excited hydrogen H(n = 2) by various ions has been treated using, among others, two non-perturbative approaches: the Classical Trajectory Monte Carlo (CTMC) method [START_REF] Jorge | [END_REF] and the semiclassical Atomic Orbitals Close-coupling (SC-AOCC) approach [20]. In the very low velocity region where the resonant capture is extremely dominant, the electron transfer process was investigated for H + + H(nl) within a pure quantal treatment, in which the total wave function is expanded in the basis set of the eigenstates of the transient molecule formed during the collision [21][22][23]. With increasing collision velocity, a pure quantal study becomes more computationally demanding and semiclassical treatments are justified and tractable [24]: in the latter, the nuclear motion is treated classically and a quantum approach is used to describe the electronic part. In particular, the well-know semiclassical non-perturbative SC-AOCC is a good approximation in the intermediate velocity range [25][26][27].

However, for collisional systems involving very high excited initial hydrogen states, the numerical implementation of the SC-AOCC approach requires extremely large basis sets and therefore very large computation time, even for proton projectile which corresponds to the simplest collision system in ion-atom collisions. Alternatively, for highly excited targets, the CTMC method presents several advantages : (i) fundamentally, through Bohr's equivalence principle when n 1, (ii) its fairly simple numerical implementation and (iii) much lower computational costs compared to the SC-AOCC approaches. In addition, the CTMC method is a good approximation at intermediate velocities to investigate electron transfer and ionization processes induced by impact of bare ions [28]. Note that CTMC has been recently used even for collisions between ions and complex targets, e.g., N a(3s), N a * (3p), Li(2s), Ar and He, where the active electron is subject to non-coulomb interactions and requires the use of model potentials [29,30] Nowadays the impressive progress of highlyparallel computers allows the use of the SC-AOCC approach with very large basis sets, which were unusable before due to prohibiting CPU execution time. Therefore, our purpose in this paper is to present a study of the capture processes for proton-hydrogen collisions involving initial excited states, up to H(n ), using non-perturbative SC-AOCC and, for comparison, CTMC approaches. Indeed the CTMC gives accurate results in the intermediate velocity range [START_REF] Olson | Physics of electronic and atomic collisions[END_REF] and it is expected to be a good approximation for highly excited targets even in the low velocity range. We also present comparisons with experimental data, when available, and with the analytical fit proposed by Janev et al. [START_REF] Janev | Collision processes in low-temperature hydrogen plasmas[END_REF], based on many theoretical and experimental data depending on the collision velocity range (see [START_REF] Janev | Collision processes in low-temperature hydrogen plasmas[END_REF] and references within). In addition, we compare our results with those used in plasma physics, available in the ADAS database [START_REF]Atomic data and analysis structure[END_REF].

The paper is organised as follows. In Sec. 2, we briefly review the CTMC method and the SC-AOCC approach that are used to obtain the results presented in Sec. 3. This is followed by their comparisons and a discussion which ends up by the writing of two empirical scaling laws and a general fitting function obtained from the SC-AOCC results for the shell selected capture cross sections. Finally, conclusions are presented in Sec. [START_REF]Atomic data and analysis structure[END_REF].

Atomic units (a.u.) are used throughout the paper unless stated otherwise.

Theoretical approaches

The proton-hydrogen collision system is a three-body problem formed by the electron e of mass m e , targetproton T of mass M p and the projectile-proton P of mass M p . The r e , r T and r P denote the position vectors for the electron, target proton and for the proton projectile, respectively, in the laboratory frame, see Fig. 1. In order to reduce the three-body problem to a two-body problem, it is convenient to use the Jacobi coordinates ( r, R, R CM ) where r is the coordinate of the electron with respect to the proton target, R the relative coordinate of the projectile with respect to the electron-target-proton center of mass and R CM the centre-of-mass coordinate. In the center-of-mass frame, the three-body system reduces to a two-body problem with the corresponding Hamiltonian

H = p 2 2µ + P 2 2M + V T e (r) + V P e (x) + V P T (s), (1)
where p and P are the momenta related to r and R, respectively, and the reduced masses µ and M are µ = m e M p /(m e +M p ) and M = M p (M p +m e )/(m e +2M p ). In Eq. 1, V T e , V P e and V P T are respectively, the electron-target, the electron-projectile and the protonproton Coulomb potentials with the corresponding separation distances r, x and s defined in Fig. 1.

In the following, we review the applied theoretical approaches used to study the capture processes for proton-hydrogen collisions.

Classical Trajectory Monte Carlo

In order to calculate the cross sections for electron transfer processes classically in H + + H(nl) collisions, we use the Classical Trajectory Monte Carlo model proposed by Abrines and Percival [14]. The dynamics of the collision system is determined by the classical Hamilton equations

ṙi = p i µ , ṗi = - ∂H ∂r i , i = 1, 2, 3 (2a)
and

Ṙi = P i M , Ṗi = - ∂H ∂R i , i = 1, 2, 3 (2b)
where the subscript i denotes the three Cartesian coordinates. The system of equations ( 2) is solved numerically by the means of a Runge-Kutta method with an adaptive step size, starting from a given configuration of the collision partners which defines the initial conditions ( r 0 , p 0 ) and ( R 0 , P 0 ). In order to represent classically an initial state H(n ) of principal quantum number n and angular momentum , we have chosen to use the Becker and MacKellar approach [START_REF] Becker | [END_REF] which generalises the procedure of Abrines and Percival by subdividing the initial microcanonical ensemble (for a complete n shell) into sub-ensembles corresponding to the quantal state H(n ). Each initial electronic trajectory is then defined by five random parameters: the ellipse eccentricity ε, where ε 2 is sampled uniformly in the interval 1-( +1) 2 n 2 , 1-2 n 2 , the mean anomaly α, where 0 ≤ α ≤ 2π, and the three Euler angles (ψ, θ, φ), which determine the random orientation of the ellipse in space. The eccentric angle ξ is defined from the parameters α and ε by solve the Kepler's equation

α = ξ -ε sin(ξ). (3)
Then the vectors r 0 and p 0 can be obtained in terms of the energy E n = -1/n 2 , ξ, ε, and the three Euler angles.

In order to initialize the incident projectile, we consider the collision frame illustrated in Fig. 2 for an azimuthal angle ϕ b = 0, i.e., the collision plane is the (XOZ) plane. The projectile is specified by its momentum P 0 = M v where the relative velocity v is chosen along Z and its initial distance R 0 from the electron target-proton center-of-mass is set sufficiently large to be able to consider that the two collision partners are isolated. The uniform beam of proton projectile is constructed by choosing randomly the square of the impact parameter b in the interval [0, b 2 max ] where b max is the value of the impact parameter beyond which the studied processes do not take place. Its value varies depending on the process, the collisional target and the collision velocity. For a given set of initial conditions ( r 0 , p 0 ) and ( R 0 , P 0 ), the set of the twelve coupled equations, resulting from (2), is solved numerically. This procedure is repeated for a large number of trajectories, N tot , in order to reach the reasonable statistical convergence of the cross sections.

During the interaction between the projectile and the hydrogen target, the electron moves under the influence of the potentials, V T e and V P e so that, at the end of the collision, three different processes can be distinguished. The first kind of process corresponds to the electron which remains bound to the target, i.e., an elastic or excitation process has occurred. The second process is when the electron is captured by the proton projectile (electron transfer). In this case, the electron passes from a bound state on the hydrogen target to a bound state on the proton projectile so that its energy relative to the projectile ion becomes negative while its energy with respect to the target gets positive. Thirdly, if none of these processes happen, ionization takes place. For a given collision velocity the total cross section (TCS) of a process τ (capture, ionization or excitation) in the H + + H(n ) collision is given by

σ n τ (v) = N n τ N tot πb 2 max , (4)
and for a reaction

H + + H(n) σ n τ (v) = n-1 =0 (2 + 1)σ n τ (v) n 2 , ( 5
)
where N n τ is the number of events of the process τ . We have run many tests to determine the appropriate b max for each H + + H(n ) collision and for each considered velocity v. For example, the capture process from excited hydrogen target may be likely for large values of b max in low velocity range: b max = 20 a.u. for H(2 ) targets, and b max = 48 a.u. for H(3s), H(3p) and H(3d) targets.

The semiclassical Atomic Orbitals Close-Coupling

The scattering problem is treated by the semi-classical atomic orbital close-coupling approach within the impact parameter approximation [24]. The relative target-projectile motion is described by straightline, constant velocity, trajectories where the position vector R reads,

R = b + vt. ( 6
)
while the electronic motion is described quantally through the time-dependent Schrödinger equation

H -i ∂ ∂t ψ( r, t) = 0, (7)
In our approach, see [34,35], the equation ( 7) is solved as follows: the electronic wave function ψ is expanded into target and projectile-centered traveling atomic orbitals ξ T k and ξ P k , where k denotes the quantum numbers n, and m [25][26][27]. The electronic wave function ψ is thus written as

ψ( r, t) = N T k=1 a T k (t)ξ T k ( r, t) + N P k=1 a P k (t)F ( r, t)ξ P k ( x, t), (8) where
F ( r, t) = e i v. r-i 1 2 v 2 t . (9)
denotes the Electron Translation Factor (ETF) [24] which ensures the Galilean invariance of the results. The time-dependent complex amplitudes a T k and a P k related to the states ξ T k and ξ P k , respectively, have to be obtained numerically and become probability amplitudes in the asymptotic regions (t → ∞) where all states in the expansion (8) are orthogonal. The calculations are performed in the collision frame as considered in the CTMC method, i.e., the impact parameter b is parallel to X and the velocity v is parallel to the quantization axis Z. The first step is to obtain numerically the eigenstates (and pseudo-states) ξ T,P k and their associated eigenvalues (positive energies) by diagonalizing the Hamiltonians h T = H -V P e -V P T and h P = H -V T e -V P T , respectively. In our calculations, we used a basis set consisting of 160 states of the target and projectile hydrogen expressed in terms of Gaussian type orbitals s, p, d, f and g, on each center. Energetically, this basis set describes well the bound states up to the level n = 7 with negligible differences with the exact energies. The complex amplitudes a T k and a P k are obtained asymptotically by solving numerically the equation ( 7) along a certain nuclear trajectory (b,v). The transition probability from the initial state |n m to another state |n m leading to the electronic process τ under consideration is defined as follows

P n m→n m τ (b, v) = |a K n m (t → +∞)| 2 , ( 10
)
where a K n m denotes the complex amplitudes a T n m or a P n m . The partial cross section of an electronic process τ for a given velocity corresponding to the transition |n m → |n m reads

σ n m→n m τ (v) = 2π bmax 0 P n m→n m τ (b, v)bdb. (11)
Then, the total cross section of a process τ is given by

σ n m τ (v) = n m σ n m→n m τ (v). ( 12
)
Note that in the case of the capture processes, the values of b max established in the SC-AOCC approach are larger than those found in the CTMC approach for a given collision velocity. This is mainly due the quantum description of the electronic motion which does not possess the cut-off in the probability density inherent to the classical description of bound trajectories around an attraction center. Convergence of the cross sections with respect to the basis set was carefully checked by using a larger basis set composed of 191 states on each center: it was found that the cross sections evaluated with the large basis do not differ by more than 3% from those obtained with the first basis for all impact energies considered in this work, for the initial targets H(n = 1, 2). However for target H(n=3), the difference between both basis sets at low and high velocities does not exceed 4% while at intermediate velocities the difference can reach about 10%. TCS of capture (cm²)

Results and discussion

Hvelplund [START_REF] Hvelplund | [END_REF] Newman [6] W ittkover [7] Bayfield [8] Schwab [9] McClure [10] Gealy [11] Fite [12] Fite [13] Present CTMC Fig. 3a shows the comparison between the total cross sections of the capture process obtained with our CTMC and SC-AOCC approaches and the available experimental data. In the intermediate velocity range (v ≈ 1 a.u.), the present CTMC results are in agreement with those of the SC-AOCC and with the experimental data reported in the references [START_REF] Hvelplund | [END_REF][6][7][8][9][10][11][12][13].

In the low velocity range (v < 1 a.u.), the CTMC cross sections show a serious discrepancy, showing a nearly constant value, thus underestimating the experimental and semiclassical data. For initial ground state target, a quantum or semiclassical treatment is therefore needed in this latter region where molecular mechanisms play a dominant role, not described by classical mechanics.

The present SC-AOCC results are in excellent agreement with those found experimentally in the all velocity region considered, 0.01 < v < 2 a.u. In the low velocity range, the resonant capture, for which the energy transferred to the electron is almost zero, is dominant and well taken into account by the SC-AOCC approach. The observed very good agreement with the experimental data suggests also that even for the lowest velocities considered, the nuclear relative motion can be treated classically, see Eq. ( 6). This was also demonstrated by Hunter et al. [36] by comparing their purely quantum results (Perturbed Stationary States in Partial Wave treatment) with the semiclassical results of Bates et al. [37].

In Fig. 3b our results are compared with two recent close-coupling calculations [38,39]. In [38] the authors have used a fully quantum-mechanical approach based on the atomic orbital close coupling method (QM-AOCC). Their TCS are in agreement with our present SC-AOCC results and with the experiments except in the intermediate velocity range (0.2 < v < 2 a.u). Results in [39], obtained with the quantum-mechanical convergent close-coupling method (QM-CCC) using Laguerre basis functions, are in excellent agreement with our semiclassical results in the intermediate velocity range.

It should be mentioned that the SC-AOCC approaches as well as, intrinsically, CMTC method do not formally satisfy the boundary conditions, as defined in the most refined continuum disturbed waves theories, see for example [40][41][42]. To evaluate the effects of the boundary conditions, we compare in Fig. 3c our results with those of da Silva et al. [43]. These authors have calculated the cross section of the resonant capture 1s-1s using two-state close-coupling approximation with the continuum distorted wave basis set (TSCC-CDW) proposed by Cheshire [40], which satisfies the boundary conditions. Fig. 3c shows our cross sections for total capture with those reported in [43] and deduced from their resonant capture cross sections by using the scaling rule σ 1s→all s states capture ≈ 1.202σ 1s→1s capture proposed by Oppenheimer [44]. These TCS are in excellent agreement with the present SC-AOCC results in the overlapping velocity range, showing negligible effects of the boundary conditions in this energy range.

Total cross sections of capture in H + + H(n = 2) collisions

The electron capture processes from excited states have been much less investigated. To compare our results, we use the available results for capture from H(2s) and H(2p) targets from the ADAS database which uses two series of theoretical data: at high velocities, the results of Eichler [45], based on the perturbative eikonal theory, and, for velocities in the low and intermediate ranges, the SC-AOCC results of Reinhold and Falcon [46]. ADAS DATA [START_REF]Atomic data and analysis structure[END_REF] CTMC results In Figs. 4 and5, we present the capture TCS for the excited targets H(2s) and H(2p), respectively. In the intermediate velocity range, down to v ≈ 0.5 a.u., the present SC-AOCC results are in a good agreement with those obtained by the CTMC approach, but there is a significant disagreement with those stored in ADAS [START_REF]Atomic data and analysis structure[END_REF], especially for H(2p) initial state, certainly related to the limited basis set used in [46] . At low velocities (v < 0.5 a.u.), there is a pronounced discrepancy between the SC-AOCC results and those of CTMC corresponding to H(2s) and H(2p) targets. This disagreement between the quantum and classical results becomes larger for very low velocities. Note however that the discrepancy is less pronounced that the one shown in Fig. 3 for H(1s) initial state. In the high velocity region (v > 1 a.u.), in the case of H(2s) excited target, the TCS of capture stored in ADAS database decay faster than those obtained with the CTMC and SC-AOCC methods whereas for the H(2p) target the ADAS data are closed to those obtained with ADAS DATA [START_REF]Atomic data and analysis structure[END_REF] CTMC results the present two approaches. In order to compare the semiclassical and the CTMC results with the analytical fit of Janev [START_REF] Janev | Collision processes in low-temperature hydrogen plasmas[END_REF], we calculate the average capture cross sections (see Eq. 5) from all possible quantum states of the energy level n = 2. The results are shown in Fig. 6. We note that in the regime of low and intermediate velocities, our SC-AOCC results are in agreement with those recommended by the analytical fit while the CTMC results and the ADAS data underestimate the TCS of capture. At high velocities the analytical fit decreases faster than those calculated by CTMC and SC-AOCC approaches.

SC-AOCC results

SC-AOCC results

ADAS DATA [4]

Janev Fit [START_REF] Janev | Collision processes in low-temperature hydrogen plasmas[END_REF] SC-AOCC results

Total cross sections of capture in H + + H(n = 3) collisions

To our knowledge, there is no non-perturbative calculations published for the capture TCS from the excited states H(3s), H(3p) and H(3d). In absence of theoretical and experimental data, we compare in Figs. 7,8, and 9 the present SC-AOCC results only with the CTMC ones.

At low velocities, ie. v < v n=3 e ≈ 0.333 a.u. which corresponds to the matching velocity regime where the projectile velocity is close to the average electronic velocity in the initial state, we note that the present CTMC results for the capture TCS from H(3 ) are again lower than the semiclassical results. In contrast, a good agreement between both approaches is found in the intermediate velocity range

(v ≈ v n=3 e ). At high velocities (v > v n=3 e
), for capture from the excited H(3s) and H(3p) sub-shells, a good agreement is found between the SC-AOCC and the CTMC results. For the capture from H(3d), the CTMC results decrease faster than those obtained by the SC-AOCC approach.

In order to compare the results from these three sub-shells, we report in Fig. 10 the capture cross sections for H(3 ) initial states obtained with the SC-AOCC approach. First, in the low velocity range, the three series of TCS are of the same order of magnitude for all quantum states. At intermediate velocities, the capture from the state H(3d) is more likely than the capture from the excited state H(3p), which is larger than the capture from the state H(3s). This latter behaviour is opposite and more pronounced in the case of high velocities. Therefore we can write,

     σ 3d c (v) ≈ σ 3p c (v) ≈ σ 3s c (v) (v 0.333 a.u.) σ 3d c (v) > σ 3p c (v) > σ 3s c (v) (v ≈ 0.333 a.u.) σ 3d c (v) < σ 3p c (v) < σ 3s c (v) (v 0.333 a.u.) (13)
where σ 3s c , σ 3p c and σ 3d c are the capture TCS from the excited states H(3s), H(3p) and H(3d), respectively.

In Fig. 11, we present the classical and semiclassical results for the capture averaged TCS from the H(n = 3) shell, together with the fit derived in Janev et al. [START_REF] Janev | Collision processes in low-temperature hydrogen plasmas[END_REF] and the ADAS data, ie. data obtained by a perturbative eikonal theory [45]. Our semiclassical results agree well at intermediate and low velocities with the fit [START_REF] Janev | Collision processes in low-temperature hydrogen plasmas[END_REF] and with those stored in ADAS [START_REF]Atomic data and analysis structure[END_REF]. As noted before, the CTMC approach underestimates the TCS at low velocities, but performs well at higher velocities. On the other hand, in the high velocity range, the capture TCS from the fit decrease faster than the CTMC and semiclassical results.

It is important to note that at low velocities, the discrepancy between classical and semiclassical TCS of capture becomes less and less important when considering higher and higher excited initial targets. It is thus expected that this disagreement becomes negligible for capture from highly excited states. In these cases, the CTMC can therefore be used even at low velocities. On the other hand the semiclassical calculations of capture cross section from higher excited initial states H(n > 3) need larger basis sets that the one we employ here. In particular, it should include higher angular momenta ( > 4) so that the SC-AOCC approach cannot currently be used in a reasonable CPU time. These are the reasons why we only used the CTMC approach to investigate capture from high excited states, H(n > 3). In Fig. 12, we present the corresponding cross sections for capture from the H(n = 1, ..., 7) targets. We can obviously remark that the TCS of capture from any state H(n ≤ 7) has a similar behaviour which will be analysed in the next section. n 3 at low velocities (ṽ < 1). This behaviour was first predicted by Bates and Reid in an earlier work [22] using a two-state molecular quantal treatment. The authors have performed the calculations for resonant capture from initial states up to the H(n = 5) shells; they report that the cross sections of the resonant electron transfer can be written as

σ n→n c (v) = n 3 A(v) (n ≥ 2) ( 14
)
where A is a function of the velocity v.

On the other hand, as shown in Fig. 14, the TCS of capture follow a scaling law in n 4 in the region of intermediate and high velocities, ie. (ṽ ≥ 1). The capture TCS from any excited state are then proportional to the geometrical cross section S n (S n = πn 4 ) corresponding to a shell H(n). Thus, the TCS of electron transfer in the intermediate and high velocity can be described by an analytical function as

σ n c (v) = n 4 B(v) (n ≥ 1) (15)
where B is a function of v. Note that this scaling works perfectly for the CTMC results and explains the similar behaviour of the TCS observed in Fig. 12.

For capture from excited hydrogen atoms (n ≥ 2) by protons, we combine these two observations to provide an universal analytical function that can be used to predict the TCS of electron capture in a wide range of impact velocities including low, intermediate and high velocity regimes. The proposed analytical TCS reads,

σ n c (v) = n 3 (1 -g n (v))A n (v) + n 4 g n (v)B n (v) [10 -16 cm 2 ], (16)
where g n (v) = 1/(1 + e -α(v-1/n) ) and v is in atomic unit. In order to find the appropriate functions A n and B n , we employed the analytical functions used in [START_REF] Janev | Collision processes in low-temperature hydrogen plasmas[END_REF]: As a final note we would like to mention that we found that the TCS from [START_REF] Janev | Collision processes in low-temperature hydrogen plasmas[END_REF] decrease in the high velocity range faster than our semiclassical results (see Figs. 6 to 11), failure which is corrected in the present fit. However in this latter region, the analytical fit of Janev et al. [START_REF] Janev | Collision processes in low-temperature hydrogen plasmas[END_REF] is only based on the theoretical data of Eichler [45] which were obtained by the eikonal perturbative treatment. Alternatively, the proposed analytical fitting (Eq. 16) shows a good agreement with the exact SC-AOCC results for all velocities, as shown in the Figs. 15a and15b, which can validate this latter fit: the difference of the TCS evaluated by our fit with respect to the calculated SC-AOCC results for the initial states H(n = 2) and H(n = 3) does not exceed 15% in the high and intermediate velocities and about 10% at low velocities.

A n (v) = 4.8 log(2 n 2 v 2 + 10 n-2 ) , v < v e B n (v) = A 1 + B(nv) 2 + C(nv) 7.2 + D(nv) 11 , v ≥ v e , (17)
As a final task to check the validity of the proposed fit, we compare in Figs. 15c and15d the TCS obtained with our fit and those calculated with the CTMC method, for the excited states H(n = 4) and H(n = 7). As expected, in the intermediate and high velocities, the fit is in a good agreement with the classical results for these highly excited states. Compared to the excited states H(n = 2 -4), the TCS of capture from the excited state H(n = 7) by protons are in excellent agreement with those established by CTMC at low velocities where the discrepancy between the approaches manifests only towards very low velocities, below to 0.04 a.u. (40 eV/u).

Conclusions

In this paper we have presented total cross sections corresponding to the capture processes from excited atomic hydrogen targets H(n = 2, 3) by impact of protons, using two non-perturbative approaches (CTMC and SC-AOCC) at impact velocities in the (0.01 -2)a.u. range. The non-perturbative SC-AOCC provides accurate results in the low, intermediate and high velocity regions. For capture from higher excited targets H(n > 3) the SC-AOCC cannot currently be employed due to prohibiting CPU execution time. Alternatively, we searched for an analytical fit according to the SC-AOCC results in order to overcome this numerical limitation. The proposed analytical fit for capture TCS is relevant in the whole range of impact energies considered in this work and gives approximatively the prediction of SC-AOCC for any H(n) exited target.

At intermediate and high velocities for the excited targets H(n = 2, 3) we show an excellent agreement between the present non-perturbative classical and semiclassical treatments.

Furthermore, at low velocities, the discrepancy between classical and semiclassical capture TCS becomes less important when considering higher and higher excited targets. Therefore, in the case of capture from excited states H(n 1), the CTMC can be considered as a solid approach, not only for the intermediate and high velocities [47] but also for low velocities, following the corresponding shift of the matching velocity v n e . Finally, we should mention that the cross sections for initial H(n = 2, 3) target stored in ADAS should be updated especially for the low and high velocity ranges.

Figure 1 .

 1 Figure 1. (color online ) Schematic illustration of the relation between the position vectors re, r T and r P in the laboratory frame and the Jacobi coordinates r and R, where the e, T and P denote the electron, the proton target and the proton projectile, respectively.

Figure 2 .

 2 Figure 2. (color online )The collision frame and geometry: the collision plane is defined by the vectors R, v and b. In our case, the chosen collision plane is the (XOZ) one, i.e., the relative velocity v is chosen along the Z axis and the impact parameter b along X

3. 1 .

 1 Total cross sections of capture in H + + H(1s) collisions

PresentFigure 3 .

 3 Figure 3. (color online) Total cross section of capture from H(1s) by protons as function of the relative velocity v.

Figure 4 .

 4 Figure 4. (color online) TCS of capture from H(2s) by protons as function of the relative velocity v.

Figure 5 .

 5 Figure 5. (color online) TCS of capture from H(2p) by protons as function of the relative velocity v.

Figure 6 .

 6 Figure 6. (color online) TCS of capture from H * (n = 2) by protons as function of relative velocity v.

Figure 7 .Figure 8 .

 78 Figure 7. (color online) TCS of capture from the excited state H(3s) by protons as function of relative velocity v.

Figure 9 .Figure 10 .

 910 Figure 9. (color online) TCS of capture from the excited state H(3d) by protons as function of relative velocity v.

3. 4 .

 4 n 3 and n 4 scaling laws for total cross sections of capture from H(n) by protons In Fig.13, we present the semiclassical TCS scaled in n 3 as a function of the scaled velocity ṽ = v/v n e , where v n e = 1/n is the velocity of the electron target on the Bohr orbit corresponding to the energy level n. The present semiclassical non-perturbative treatment shows that the TCS of capture from the excited states H(n = 2, 3) follow fairly well an empirical scaling in

Figure 11 .Figure 12 .

 1112 Figure 11. (color online) TCS of capture from H(n = 3) by protons as function of relative velocity v.

Figure 13 .Figure 14 .

 1314 Figure 13. (color online) The scaled TCS of capture from H(n = 1, 2, 3) by protons as a function of the scaled velocity ṽ. The results stem from the SC-AOCC calculations.

  where the constants A, B, C and D are all positive, with A = 7, B = 0.75, C = 0.4 and D = 9.77 10 -4 .

Figure 15 .

 15 Figure 15. (color online) The comparison between the proposed fit and CTMC results.

Acknowledgments

AT was supported by a grant of Institut de Recherche pour le Développement IRD Bondy-France. This work was partly supported by LABEXs Plas@par and MiChem under Grant No. ANR-11-IDEX-0004-02.