CONSTRUCTION AND CONVERGENCE ANALYSIS OF CONSERVATIVE SECOND ORDER LOCAL TIME DISCRETISATION FOR LINEAR WAVE EQUATIONS

. In this work we present and analyse a time discretisation strategy for linear wave equations that aims at using locally in space the most adapted time discretisation among a family of implicit or explicit centered second order schemes. The proposed family of schemes is adapted to domain decomposition methods such as the mortar element method. They correspond in that case to local implicit schemes and to local time stepping. We show that, if some regularity properties of the solution are satisfied and if the time step verifies a stability condition, then the family of proposed time discretisa-tions provides, in a strong norm, second order space-time convergence. Finally, we provide 1D and 2D numerical illustrations that confirm the obtained theoretical results and we compare our approach on 1D test cases to other existing local time stepping strategies for wave equations

-Local implicit time discretisation, see for instance [7,18,21,30,31,39]. The strategy is to treat by an implicit time integration scheme the ODEs acting on the degrees of freedom corresponding to the region where the pertubations occur. By doing so, the time step restriction (CFL) is decorrelated from the perturbations. The price to pay is that a (hopefully small) linear problem must be solved at each iteration. -Local Time Stepping (LTS), see for instance [17,19,22,26,38]. The strategy is to use a first time marching scheme in the whole domain and a second one in the perturbed region. The chosen type of time discretisation used in both regions is often the same but time steps differ: a smaller time step is used locally. One can distinguish non-conservative strategies (see for instance [22]) from conservative strategies. The latter are based upon Leap-Frog schemes and can be separated into two categories depending on how sub-domains are coupled.
• Implicit LTS. A domain decomposition strategy is introduced at the continuous level together with some coupling conditions at the interface of the subdomains (typically by introducing a Lagrange multiplier to enforce in a weak sense those conditions as in the mortar element method, see [42] or [28]). This idea can be traced back to the work of Collino et al. [13,14] and has been pursued and improved in [2,17,33,37,38]. Such strategy is referred as implicit since the treatment of the transmission conditions is done implicitly at the fully discrete level. • Fully explicit LTS. The decomposition of the domain is done at the discrete level through the use of a discrete restriction operator on the region (and its surrounding) where perturbations occur. The resulting scheme does not introduce transmission conditions in the classical sense but is fully explicit. It has first been proposed in [19] and various extensions exist: Maxwell's equations (see [25]) and multi-level LTS (see [20]). Recently, in [26,27], a proof of space-time convergence is given. It shows that, for the scalar wave equation, a second order space-time convergence holds in the 2 norm in space.
In general, the space-time convergence analysis for these numerical methods is not always available, Grote et al. [26,27] being one exception, and although the methods have been implemented in realistic situations as 2D and 3D frameworks [18,21,22], their theoretical background does not always rely on a robust analysis. The question of convergence of the method as both the space and time steps vanish together is of crucial importance when dealing with PDEs as linear wave equations, and to this purpose, the energy based analysis has proved very efficient [17] and will be followed in the present work. In this work we construct and analyse local time discretisations that gather in an original framework both local implicit time discretisation and conservative implicit LTS. Moreover: -We show that the proposed time discretisations provide, under some regularity and stability conditions, second order space-time convergence, in a strong norm (for scalar wave equations, it provides convergence for the 1 -norm in space). -We provide extensive numerical convergence experiments for a 1D scalar wave propagation problem. The results show that our approach provides better space-time convergence properties, in the 1 -norm, than existing LTS approaches. In particular we study some situations where the LTS of [19] converges in ∆ 3/2 in ∞ (0, ; 1 (Ω)) whereas our approach always provides second order convergence.
The outline of the article is the following. In Section 2 we give all the necessary notations and assumptions related to the discretisation in space of linear conservative wave type problems. Section 3 is devoted to the introduction of a class of time discretisations -parameterised by two polynomial functions and -for which we show stability and second order space-time convergence results under some assumptions on the parameters (i.e. the coefficients of the polynomials and ) and some CFL conditions. In Section 4 we first present two preliminary applications of our discretisation framework. By adequately choosing the polynomial functions and we construct a local implicit time discretisation (Sect. 4.1) and, in Section 4.2, a first local timestepping scheme (with a ratio 2, see Sect. 5.2 for an accurate definition of the term ratio). Finally, in Section 4.3 we propose a strategy to construct general local time-stepping schemes. This strategy is based on the use of Chebychev polynomials (more precisely on Leap-Frog Chebychev method as introduced in [24]). Space-time numerical convergence results illustrate the developed theory for 1D test cases in Section 5 and for 2D test cases in Section 6. Finally, in Section 7 we compare our approach algorithmically to other local time stepping strategies: we first compare our approach numerically with the Fully explicit LTS of [19]. Moreover we explain why the proposed schemes can be seen as a generalisation of the ones proposed in [17].
The source code used to obtain convergence curves of Sections 5 and 7 are available as supplementary materials at the web link [44].

Semi-discrete wave propagation problem
We are interested in the simulation of coupled linear wave propagation problems. The most simple example one could think of is given by the following problem: being given a bounded connected open domain Ω partitioned as two disjoint connected domains Ω and Ω , find ( ) ∈ 1 (Ω ) and ( ) ∈ 1 (Ω ), for all ∈ [0, ] such that with for instance homogeneous Neumann boundary conditions on the domain's boundary ∇ · = 0 on Ω ∩ Ω, ∇ · = 0 on Ω ∩ Ω, and some transmission conditions on the complementary boundary Σ, that satisfies Σ ∩ Ω = ∅, where is the outward unitary normal of Ω . That is completed with initial conditions for and and their time derivative. The scalar functions ∈ ∞ (Ω ) and ∈ ∞ (Ω ) are positive and bounded from below. Such problems find applications in the wave scattering by obstacles and are of interest for modeling non destructive experiments for instance.

Continuous abstract formulation and main assumptions
In the following stands for either or . In this section we formulate the coupled wave propagation in a more abstract framework. To do so we use notations from [16], chapter XVIII, and [6]. We assume given Hilbert spaces ( , ). The space is equipped with the scalar product (·, ·) , the norm in is denoted |·| whereas the norms on is denoted ‖·‖ . Moreover we assume that is dense and continuously embedded in . We assume given a continuous hermitian bilinear form : × → R that satisfies for some real positive scalars and , We assume also being given another Hilbert space equipped with the norm ‖ · ‖ as well as a continuous bilinear form ( , ) on × . We consider the following abstract wave propagation problem: Let ∈ 0 ([0, ], ) and ∈ 0 ([0, ], ) be given, find ( ( ), ( ), ( )) ∈ × × solution, for all ∈ [0, ], to the coupled system of equations that is completed with initial conditions for and and their time derivative. Note that the scalar wave equation problem (2.1)-(2.5) enters the abstract framework presented above by choosing where is equipped with the standard 2 scalar product, and for all and in and for all in −1/2 (Σ) However, the setting we consider is rather general and for instance elastodynamics equations could also enter the abstract framework by writing transmission problems (continuity of displacements and stresses) and using vectorial forms of all the spaces and scalar products introduced. Moreover the abstract structure is also adapted to the domain decomposition method with overlapping introduced in [1] to deal with scattering problems in transient acoustics. System (2.3) can be rewritten in a more compact form using the following notation: bold letters are used to define unknowns in = × , e.g., = ( , ) and we introduce the bilinear forms We complete (2.4) with initial conditions Existence and uniqueness results for this abstract problem rely on the assumption that a so-called inf-sup condition holds. More precisely we assume that there exists > 0 such that where ‖ ‖ 2 = ‖ ‖ 2 + ‖ ‖ 2 (similarly we denote by | · | the composite norm in = × ). Notice that this condition is the same encountered in steady domain decomposition problems. We refer the reader to Chap. 7 of [5] for a proof. We do not provide here proofs of existence and uniqueness for such problems -they rely on energy analysis and/or Laplace transform -and refer, for instance to the work of [1,3] for related analysis. See also the Appendix A of this manuscript for more details on the matter. We assume that the solution has the following properties
Inequality (2.2) implies that the operator ,ℎ is self-adjoint and positive semi-definite. Its spectrum, denoted Sp( ,ℎ ), is composed by a finite number of non-negative eigenvalues. The spectral radius is defined as the maximum eigenvalue in the set Sp( ,ℎ ), i.e., ,ℎ := max Sp ( ,ℎ ) .
We also introduce the operators ,ℎ : ,ℎ ↦ → ℎ and ,ℎ : ℎ ↦ → ,ℎ as As done previously we define ℎ = ,ℎ × ,ℎ and represent by bold letters unknowns in ℎ . The semi-discrete equation we consider reads: Let ℎ ∈ 0 ([0, ], ,ℎ ) be given: find ( ℎ ( ), ℎ ( )) ∈ ℎ × ℎ and solution, for all ∈ [0, ], of together with the initial conditions In the rest of the work we assume that the following discrete inf-sup condition holds where is independent of ℎ. Since discretisation by finite elements of wave equations is now a well understood subject, most of the difficulty in constructing System (2.8) is to choose ℎ such that (2.10) holds. In fact, (2.10) is not a consequence of (2.6) and this question requires dedicated analysis. On this specific topic, we refer the reader to [4,28,41,42] for reference work concerning the mortar finite element method, that is a domain decomposition method without overlapping, and to [1], for a work concerning a domain decomposition method with overlapping.

Introduction of local time discretisations
We define the sequences { ℎ = ( ,ℎ , ,ℎ )} and { ℎ } as the approximations of ℎ ( ) and ℎ ( ) at time = ∆ for a time step ∆ > 0, and ∈ {1, 2, . . . , }. We define the final time of computation as = ∆ . These sequences are constructed by solving the following problem: where with the initial conditions We give more detail in Remark 3.5 on how the term 1 ℎ can be computed. The scheme (2.3) is consistent only if some conditions are satisfied on the polynomials ( ) and ( ). Since we want to construct perturbations, for small ∆ , of the standard centered scheme it seems natural to do the the following hypothesis. For stability reasons the time step ∆ can not be chosen arbitrarily. A so called CFL-condition has to be satisfied to obtain a stable scheme. In our case it corresponds to the assumption that follows.
Assumption 3.2. The following CFL-condition holds: there exists 0 < ≤ 1 such that Note that the spectral radius ,ℎ of equation (3.3) is (ℎ −2 ) for regular discretization when standard finite element methods are used, leading to a classical hyperbolic-type CFL condition. Note also that since and because of Assumption 3.1, we know that there exists ∆ small enough or equivalently small enough, such that (3.4) is satisfied for any fixed ℎ. As shown later, these conditions ensure the positivity of a preserved discrete energy.
We describe now more in detail an algorithm that computes the solution to (3.1). At each iteration, one needs to compute the Lagrange multiplier ℎ first, then compute +1 ,ℎ and +1 ,ℎ . More precisely, using the property that we can re-write equation (3.1b) in the following form Note that ,ℎ is a positive symmetric operator -hence invertible -if equation (3.4) holds. We now use a Schur complement technique: applying the operator ,ℎ to equation (3.7), applying the operator ,ℎ to (3.1a), we obtain by subtraction and thanks to (3.1c) the following system for the Lagrange multiplier ℎ The well-posedness of the above problem is a consequence of the surjectivity of either ,ℎ or ,ℎ which is a consequence of the inf-sup condition (2.10). Assuming that ,ℎ and ,ℎ are known, then ℎ can be computed using (3.8), it follows that +1 ,ℎ and +1 ,ℎ can be computed using respectively (3.7) and (3.1a).

Remark 3.3. A drastic simplification occurs when
In that case ,ℎ is the identity operator on ,ℎ and the volumic unknown +1 ,ℎ can be explicitly updated.
Remark 3.4. With the choice ( ) = 1 and ( ) = 1 − /4 we obtain a standard coupled explicit leap-frog schemes. It is not difficult to prove that the corresponding stability condition reads Condition (3.10) is penalizing since it depends in the same way in ,ℎ and ,ℎ but the latter can be large compared to ,ℎ .
Remark 3.5. The computation of the initial data using formula (3.2) involves the computation of the term 2 ℎ (0). This term is obtained by evaluating (2.8) at time = 0. More precisely, we have Using a Schur complement, the Lagrange multiplier is given As already mentioned the well-posedness of the above problem is a consequence of the surjectivity of either ,ℎ or ,ℎ which is a consequence of the inf-sup condition (2.10). Once ℎ (0) is computed, the value of 2 ℎ (0) is obtained easily from (3.11). Finally note that the initial data satisfy by construction the constraints

Space-time convergence analysis
We define the error terms ℎ = ( ,ℎ , ,ℎ ) and ℓ ℎ as In this section we show that the terms ℎ tends to 0 as ℎ and ∆ go to 0. More precisely we show that under Assumptions 3.1, 3.2 and 3.7 (given below) we obtain a uniform estimation with respect to ∆ and ℎ in the norm ∞ (0, ; ) of the error in (∆ 2 ) + ( ℎ ). The section is organized as follows -Definition of the consistency errors: we write a system of equations for the sequence ,ℎ , ,ℎ and ℓ ℎ that is similar to (3.1) with source terms that correspond to consistency errors that we will then specify. -Energy identity for the error equation: we proceed by energy analysis and write an energy identity satisfied by the error terms ,ℎ and ,ℎ . The introduced energy is positive under the CFL-condition of Assumption 3.2. -Stability result: we prove that the energy associated to the error { ℎ } is stable. To do so we use a discrete by-parts integration and a discrete energy analysis including the use of a discrete Gronwall's lemma. -Space-time convergence results: using Theorem 2.3, we deduce space-time convergence results in the norm ∞ (0, ; ),

Definition of the consistency errors.
We always assume that Assumption 2.2 holds and therefore that the solution is sufficiently smooth, in particular, With this assumption, all the manipulations and expression used below make sense in a standard continuous setting. Using equations (2.8) and (3.1) we obtain, for ∈ {1, 2, . . . , − 1}, with the consistency errors given by (3.14) Standard Taylor expansions allow us to simplify equations (3.13) and (3.14). There exist intermediate times such that, using equation (2.8a), Then using equation (2.8b) one can further simplify ,ℎ as If Assumptions 3.1 and 3.2 hold then there exists a rational function such that where ( ) is given by The consistency error ,ℎ has then the final expression

Energy identity for the error equation
To obtain an energy identity on the error equations (3.12)-(3.15)-(3.17) we use a standard discrete energy technique. The main ingredients of the strategy is to observe that, if Assumption 3.2 holds then is a non-negative symmetric operator. Moreover, with the same assumption, if we introduce the following notation, )︀ is well defined and is a non-negative symmetric operator. Note that from (3.18) and (3.16) we deduce that After standard algebraic manipulations (similar to the computations done in [7]) one can show the following lemma.
Lemma 3.6. Let Assumption 3.1 and 3.2 hold. Then, for ∈ {1, 2, . . . , − 1}, where ,ℎ is the identity operator in , and with Proof. For the sake of conciseness, we only list here the main steps of the proof: -compute the scalar product (·, ·) of the first equation of (3.12) with -compute the scalar product (·, ·) of the second equation of (3.12) with -sum the two obtained equations and use the third equation of (3.12) to get rid of the term involving ℓ .
are positive quadratic energy functionals if Assumption 3.2 holds.

Stability results.
To obtain meaningful results we need more assumptions on how the spectral radius of ,ℎ behaves with respect to ℎ compare to the spectral radius of ,ℎ . More precisely we assume the following property (3.23) Let us now suppose that Assumption 3.2 holds. We introduce the positive scalar ℛ , independent of ℎ, as |ℛ( )|, (3.24) where ℛ( ) is given by (3.18). Since ∆ 2 ,ℎ ≤ 4 2 2 one can show that for all ℎ in ,ℎ the following inequality holds where ( ) is given by (3.19). Again, since ∆ 2 ,ℎ ≤ 4 2 2 one can show that for all ℎ in ,ℎ the following inequality holds Theorem 3.8. Let Assumptions 2.2, 3.1, 3.2 and 3.7 hold. Then, there exists a scalar independent on , , ∆ and ℎ such that we have for ∈ {1, . . . , }, Proof. In what follows the scalar -independent on , , ∆ and ℎ -is allowed to change from one line to the other. After summing equation (3.20) over = 1 to = − 1 and taking into account equations (3.15) and (3.17), we find The proof then proceeds in five steps. One step for the estimation of each of the four terms above and a final step that collects all the obtained estimations in order to obtain (3.26) using a discrete Gronwall's lemma.
Step 1. Estimation of Ξ . Following the proof given in [8] (proof 2.4 of Lem. 2.3 and appendix) it is possible to show that It has to be noted that this inequality holds uniformly with respect to the time step (in the limit given by Assumption 3.2) and in particular it is valid if ∆ = 2/ √ ,ℎ . This result is not trivial: it is proven using a decomposition into low and high frequency components of the solution ,ℎ . Then using Cauchy-Schwarz inequality and the estimate (3.27), as well as standard algebraic manipulations, one gets Using the stability estimate (2.11) we obtain Step )︀ one can show, with the Cauchy-Schwarz and triangle inequalities, that then, since by definition of the energy (3.22) we have we can simplify (3.31), and we obtain (3.32) Using the stability estimate (2.11) we obtain Step 3. Estimation of Π . The difficulty here is that one can not expect in general to have a uniform bound on ,ℎ 2 ,ℎ ( ) in the norm in . The standard strategy is to use the following equality then, a discrete by part integration in time. The objective is to "exchange space and time derivatives" between the error term and the solution to the semi-discrete problem. The by-part integration in time is done using the following algebraic rule: for all sequences of real numbers { } and { } we have We apply the above equality to the term Π and use property (3.34) as mentioned. We obtain Moreover using the mean value theorem we find that and by the definition of the energy (3.22) one gets Injecting the estimate above as well as estimate (3.37) into (3.36), one obtains after using Cauchy-Schwarz inequality Step 4. Estimation of Λ . A similar strategy than for the estimation of Π can be applied. For that it is essential to observe the following property ,ℎ , that can be proven by diagonalisation of the operators involved on the family of eigenvectors of ,ℎ . Then the same proof as in step 3 can be used. We obtain where is given by (3.25).
Step 5. Final energy estimate and Gronwall's lemma application. Combining inequalities (3.30), (3.33), (3.38) and (3.39) obtained above, we find Then using Young's inequality we write that , for = 1 and = , and, using the above estimation into (3.40) we obtain To conclude let us use the following discrete Gronwall's lemma (see proof in Appendix B): for any real positive sequences { } and any positive scalar numbers and we have, for all ≥ 1, where we use the convention that the sum ∑︀ ( ) Estimate (3.26) shows that it is important to obtain reasonable bounds on the coefficients ℛ and . In particular, if ( ) has some roots then these coefficients may blow up. This is the main difficulty that is addressed in Section 4.3.2 when constructing polynomials for explicit local time discretisation.

Space-time convergence results.
Corollary 3.9. If the assumptions of Theorems 2.3 and 3.8 hold, then, there exists independent of ∆ and ℎ such that, for all ∈ {1, . . . , − 1} Proof. In the proof the notation refers to a positive scalar independent of ∆ and ℎ than can change from one line to another. Our first objective is to estimate by the energy terms. To do so we use the theory developed in [8] that relies on a separation of the unknowns into two orthogonal subspaces spanned by the eigenvectors of the matrices ,ℎ and ,ℎ . It is then shown in [8] that (3.27) To show that a similar identity holds for ( +1 ,ℎ − −1 ,ℎ )/2∆ , we check the hypothesis 2.2 of [8]. More precisely, since ℛ(0) = 1 one can check that there exists > 0, > 0 and > 0 such that then, following the proof of Lemma 2.3 in [8] we obtain Combining the two estimates above we obtain, From this inequality we deduce straightforwardly that Observe now, thanks to the coercivity (2.2) of the bilinear forms (·, ·), that for all ∈ {1, . . . , − 1}, We now study the initial terms 0 ℎ and 1 ℎ . By definition we have where 0 ≤ ♠ ≤ ∆ . Thanks to Assumption 2.2 we have a uniform bound of ℎ in 3 ([0, ], ℎ ) for the supremum norm in time and ‖ · ‖ in space, moreover we have Then, estimation (3.44) can be simplified using Theorem 3.8 and the above estimation of ℰ 1/2 ℎ . We obtain The statement of the corollary is obtained using an adequate decomposition of the difference { ℎ } 1/4 − ( ) and triangle inequalities. More precisely we have where the first term can be estimated by (3.45); the second term is uniformly bounded (with respect to ℎ) by ∆ 2 since ℎ ∈ 2 ([0, ], ℎ ) and (2.11) holds; and the last term can be estimated using Theorem 2.3. Corollary 3.9 means that the "good" quantity that approximates ( ) is̃︀ ℎ = { ℎ } 1/4 . One way to obtain this quantity is by post-processing the obtained solution, a more efficient approach is to compute it directly, indeed, by linearitỹ︀ ℎ can be computed solving (3.1) for ∈ {1, . . . , −2} with source term̃︀ ℎ = { ℎ ( )} 1/4 instead of ℎ ( ) and with initial datã︀ This involves only a small change in the computation of the source terms and the initial data in (3.1) but allows to recover the expected estimate

Derivation of local implicit or explicit time discretisations
In this section we derive three specific local time discretisations that enter the framework described in Section 3. The presented schemes are of increasing complexity and are constructed assuming ,ℎ and ,ℎ known.

Local implicit scheme
Local implicit strategies for wave equations have been developed and analysed by several authors, see for instance [18,21,30,39]. Moreover in [7] a second order and a fourth order local implicit time discretisation adapted to domain decomposition have been constructed. The second order method of [7] fits naturally into the family of discrete problems (3.1) that we have constructed. It is obtained by choosing ( ) = 1 and ( ) = 1.
The complete scheme reads (4.1) Notice that this means that the first equation of (2.8) is discretised with an explicit leap-frog scheme, while the second is discretised with an unconditionally stable implicit -scheme with = 1/4. It has been shown in [7] that at each time iteration, one needs to solve the following problem where ,ℎ is the identity operator in ,ℎ and whereˆℎ andˆℎ are some source terms that depend on previous iterates and of ℎ ( ). The invertibility of the above system is guaranteed if the discrete inf-sup condition (2.10) holds as explained in Section 3.1. Since ( ) and ℛ( ) are independent of then = 0 and ℛ = 1 are obviously independent of that can be arbitrarily high (hence the ratio ,ℎ / ,ℎ can be arbitrarily high). Finally, the application of Corollary 3.9 proves the space-time convergence of (4.1).
Note that, when applied to the wave equation (2.1), solving System (4.2) corresponds to solving the wave equation in Ω with an implicit scheme augmented by some operator acting on boundaries that accounts for the transmission of fluxes between Ω and Ω as well as the equality between and on Σ. This scheme is particularly adapted if a very strong and very local heterogeneity is considered in the propagating medium. In that case (4.2) is not well conditioned but the algebraic system has a small size and can be solved efficiently.

Stabilised explicit scheme
Our objective in this section is to construct a time discretisation that allows to treat situations for which we have that is to say = 2 in Assumption 3.7. Note that we expect ,ℎ ≃ 4 ,ℎ for the scheme to be meaningful and efficient. For instance, in the case of standard P -finite elements on a uniform mesh for the scalar wave equation (2.1), if the mesh size used to discretise Ω is two times smaller than the mesh size used to discretise Ω , we have ,ℎ = 4 ,ℎ . The scheme is constructed by choosing ( ) = 1 − 16 (4.4) and setting

(4.6)
Observe that at each time iteration, computing ℎ requires to solve: Then ℎ is used to compute +1 ,ℎ and +1 ,ℎ explicitly using the first two equations of (4.6). The well-posedness property of (4.7) is a consequence of the discrete inf-sup condition (2.10). To apply Corollary 3.9 one needs to check that Assumption 3.2 holds. Since we have assumed = 2 (i.e. ,ℎ ≤ 4 ,ℎ ), we need to check (3.4), which reads, using (3.5), From the definition of ( ) given by (4.5) (see Fig. 1) it is clear that ( ) ≥ 0 for all positive ≤ 1 (it has a double root at = 8). However from the definition of (4.4) we see that ( ) > 0 only if is strictly less than one, moreover we have and therefore = max This estimate illustrates that the value = 1 is forbidden to apply Theorem 3.8. However we will see that in practice a value really close to 1 gives satisfactory results (to back up this claim, several space-time convergence curves for different values of are presented in Sect. 5.2). To conclude, we have constructed a time discretisation that is stable and convergent if ,ℎ ≤ 4 ,ℎ and ∆ is chosen strictly below the optimal value 2/ √ ,ℎ .

Principle
In the same spirit than Section 4.2, we construct now a method that can be characterised as a conservative local time stepping technique with an implicit treatment of transmission terms. As in Section 4.2, the unknown ,ℎ should be explicitly updated, hence following Remark 3.3, we assume that ( ) satisfies The complete scheme reads (4.10) We do not yet specify the polynomial ( ), but from Assumption 3.1 (consistency assumption) we must have (0) = (0) = 1. Notice that by definition (4.9) we have (0) = 1 ⇒ (0) = 1.
Our objective is then to construct a sequence of functions ,ℓ that satisfy the properties for a monotonically increasing sequence ℓ > 2. To satisfy Assumption 3.2 (stability assumption) it is sufficient to check that The last inequality means that Assumption 3.7 is verified with = ℓ / . Observe that this is an improvement compared to condition (4.3) associated to scheme (4.6). In Section 4.3.2 below we present a procedure to construct the sequence of polynomials that satisfy the property (4.11) for some increasing sequence of ℓ . In Section 4.3.3 we apply the algorithm of Section 4.3.2 below and construct a family of polynomials for which we have 2 ≃ 3, 3 ≃ 4 and 4 ≃ 5. (4.12) Remark 4.1. The stability condition of the scheme (4.10) with ≡ ,ℓ can be rewritten Because of (4.12), this is clearly in improvement compared to the stability condition (3.10).

Construction of a parametrized polynomials sequence
To construct the sequence of polynomials that satisfy property (4.11) for a monotonically increasing sequence ℓ , we start from the polynomials introduced in [24] that correspond to shifted and stretched Chebychev's polynomials. They are given bỹ︀ where ℓ ( ) is the ℓth Chebychev polynomial. The first polynomials are given by It is proven in [24] that the polynomials̃︀ ,ℓ ( ) satisfỹ︀ The polynomials satisfy the good requirements that we have stated in order to construct the local time stepping process, i.e. (4.11), except for the fact that thẽ︀ ,ℓ ( ) do vanish for some ≤ 4(ℓ + 1) 2 . An idea used in [40] and [32] in the context of stabilisation of the Runge-Kutta method is to transform̃︀ ,ℓ ( ) to obtain the required behavior (i.e.̃︀ ,ℓ ( ) > 0). Note that a similar idea is used concurrently in the context of non-linear wave propagation phenomena in [15]. We define the family of functions ,ℓ ( ) parametrized by ( , , ) such that, for positive and sufficiently small, see Figure 2 for an illustration of a representation of the polynomial ,ℓ for various parameters values. Note that if = 1, = 0 and = 0 one recovers ,ℓ =̃︀ ,ℓ .
Let us suppose given 0 < < 4 small enough. We propose a procedure (see again Fig. 2) that computes ≡ and ≡ such that ,ℓ is a well-defined polynomial and consistency as well as stability are ensured. Namely, one should check that equations (4.11) holds.
Step i. From the definition (4.16) one can check that ,ℓ is a polynomial if the a-priori blow-up at = 0 is compensated, this means, that one should have, It can be observed that is a root of polynomial of order ℓ + 1. However since the polynomial̃︀ ,ℓ ( ) behave like the linear function is a neighborhood of = 0 for positive and sufficiently small there exists a real negative solution to (4.17). Hence, we choose as the negative solution to (4.17) with the smallest absolute value. Note that we have → 0 when → 0.  Step ii. To satisfy Assumption 3.1, i.e. the consistency assumption, one must check that ,ℓ (0) = 1. To do so we first differentiate (4.16) with respect to , we obtain Note that for positive and small enough the coefficient is well defined and close to one and we we have → 1 when → 0.
Step iii. The last step concerns the definition of ℓ such that (4.11c) and (4.11d) hold. From the definition (4.16) and the property (4.15) one can see that (4.11) holds for 4 ( ℓ ) 2 = 4(ℓ + 1) 2 − · Notice that ℓ → ℓ + 1 when → 0. Regarding stability, the resulting scheme indeed needs the constants ℛ and to be bounded in order to provide second order accuracy in the ‖·‖ norm. If ,ℓ reaches 0, ℛ and degenerate. For the value of ℓ given above, it can be shown that 4 ( ℓ ) depending on the parity of .

Numerical construction
In this subsection we apply the algorithm explained above for ℓ = 2 to ℓ = 4. Since the proposed algorithm is parametrized by we choose ∈ {0.1, 0.5, 1}. Table 1 gives the values of ( , ) that are computed, along with the corresponding values of ℓ . The script used to compute these values is also provided at the web link [44]. In Figure 3 we have plotted the obtained polynomials.
Remark 4.2. The process we have presented can be used for arbitrarily large ℓ however should be chosen small enough (and presumably smaller and smaller as ℓ increase) and in that case the constants and ℛ will degenerate. Numerical results presented Sections 5.2 and 5.3 confirm the fact that it is necessary to have and ℛ bounded to obtain second order convergence in the norm ‖ · ‖ .

Algorithmic discussion
In practice System (4.10) is solved by in three steps.
(1) Explicit computation of intermediate unknowns ,* ,ℎ and ,* ,ℎ by The evaluation of the term requires ℓ + 1 evaluations of the operator ,ℎ (i.e. matrix-vector product with its algebraic representation). These evaluations can be done in two different ways, either by computing explicitly the coefficient of the polynomials ,ℓ ( ) and using Horner algorithm, or using the second order recurrence relation of the Chebychev polynomials that can be derived from the definitions (4.14) and (4.16). For numerical stability issues, this latter choice should be privileged if ℓ is large, however for small values of ℓ both algorithms perform similarly and Horner algorithm is -in our opinion -easier to implement.
(2) Computation of ℎ by solving Note that in the context of the domain decomposition method, the Lagrange multiplier ℎ corresponds to an unknown discretised on an interface and therefore the size of the system to invert in (4.21) is limited. Note also that, at the algebraic level, the corresponding matrix has a bandwidth that increases with ℓ and is pre-computed before iterating (see Sect. 6.2 for more details on practical algorithmic aspects).  This step involves ℓ evaluations of the operator ,ℎ . This is clearly an over-cost compared to an ideal local time stepping strategy. Note however that the term ,ℎ ℎ may have a small support, this property is used in practice to reduce the amount of computations needed in this step.
To summarize, the main complexity of the algorithm comes from -at each time iteration -the product with the operator ,ℎ , 2ℓ + 1 products with ,ℎ and the solving step of (4.21). The time step being relaxed ℓ ≃ ℓ + 1 times, the proposed algorithm is especially efficient when the operator ,ℎ is more costly to evaluate than the operator ,ℎ . The proposed algorithmic complexity is the same as in the local time stepping introduced in [17]. The strategy presented in [19] required a smaller number of matrix-vector product is required and no system must be solved. Therefore, in the same configuration (same mesh and same time step for which both methods are stable), the algorithm proposed in [19] should be more efficient. However, this method does not rely on a decomposition domain approach and has yet to be extended to transmission problems between non-conformal meshes. For these reasons, definitive comparisons between both methods are difficult. Remark 4.3. Note that although the size of the matrix system corresponding to equation (4.21) may be relatively small (it is reduced to the Lagrange multiplier unknown on the interface) the matrix may become full for large ℓ. In that case the local implicit scheme of Section 4.1 may be preferred.

Numerical illustrations in 1D
The computational code used to obtain the results of this section is available as supplementary material at the web link [44].
In this section we present numerical results in 1D that illustrate the convergence behavior of the schemes we have proposed. We consider the wave equation (2.1), with homogeneous Neumann boundary condition, posed on the domain Ω = (−0.5, 0.5) with Ω = (−0.5, 0), Ω = (0, 0.5) and Σ = {0}. We assume that = 1 and we denote ≡ ≤ 1. Note that is therefore a measure of the contrast between the two subdomains. We consider the propagation problem of a pulse generated by a source term. More precisely, initial data are null and the source term is computed so that the exact solution is given by, for ∈ [0, ] with = 0.5, 2 / 2 ) (a smooth compactly supported pulse), 0 = −0.25, = 0.05 and Qualitatively speaking the sought solution corresponds after time > + 0 to a right propagating pulse that is solution to a transport problem -with no source term -before it reaches the interface Σ.
For the space discretisation we use standard second order Galerkin finite elements with a lumped mass matrix (see for instance [10][11][12]) on a uniform mesh of Ω and Ω and we denote ℎ and ℎ the respective mesh sizes and the refinement rate. We have s ℎ = ℎ, ℎ = ℎ and  The chosen space discretisation provides a second order convergence in space for the 1 -norm (see [23]). We recall that ∆ = 2 / √ ,ℎ for some 0 < ≤ 1. In what follows, we plot space-time convergence curves by setting to some given values, and computing the solution to the discrete problem for some sequence ℎ going to zero (this implies that ∆ goes to zero accordingly). Then the discrete solution ( ,ℎ , ,ℎ ) is compared to the analytic expression (5.1) and we plot where ℐ ,ℎ and ℐ ,ℎ denote the interpolation operators on the nodal finite element spaces.

Local implicit scheme
In order to assess numerically the behavior of local implicit schemes described in Section 4.1, we set ( ) = 1 and ( ) = 1. More specifically, this means that the left-hand side of the domain is discretized with an explicit leap-frog scheme, while the right-hand side of the domain is discretized with an unconditionally stable implicit -scheme with = 1/4. We have chosen the optimal value = 1 and the values for the refinement ratio and the ratio of velocity ( , ) = (10, 1) and ( , ) = (20,2). This configuration is difficult since the constraint on the time step is severe (a pure explicit case would require a time step around 10 to 30 smaller). Nevertheless, the convergence plots represented in Figure 4 show that a second order rate of convergence is achieved. In fact this was to be expected since Corollary 3.8 can be applied with = 0 and ℛ = 1.

Stabilised explicit scheme
In order to assess the behavior of the stabilised explicit scheme described in Section 4.2, we set We first investigate the situation of a homogeneous medium ( = 1) where the subdomain Ω is refined by a factor = 2. We make the value of increase from = 0.9 to 1. As stated in Section 4.2, the value = 1 prevents us from applying Corollary 3.9 since the values of ℛ and blow up when approaches 1. The numerical results displayed in Figure 5a show that values of very close to 1 (up to 0.999) give the expected convergence rate of 2, and that indeed, choosing = 1 does not lead to a second order space/time convergence (the convergence is of order 1 before diverging). As a second example, we consider an inhomogeneous medium with = 0.25, we choose = 4, and we perform the same numerical tests. As observed in Figure 5b, the same conclusions can be drawn.

Local time stepping using the Leap-Frog Chebychev method
In order to assess the behavior of the schemes constructed in Section 4.3, we set, for a given , ( ) = ,ℓ ( ) according to the values given in Table 1 or ( ) =̃︀ ℓ ( ) and we always choose ( ) as in (4.9) meaning that we are using explicit schemes. We consider = 1, and = 3 or = 4. According to equation (5.2) we have ,ℎ = 2 ,ℎ . The polynomial ,ℓ is chosen such that ℓ = − 1 and = 0.1, and the value of ≡ is chosen as where the values of ℓ are given in Table 1. This choices of parameters ensure that the stability condition is fulfilled. In Figure 6, we have displayed the convergence obtained. In all the cases we observe second order space-time convergence with the polynomials ( ) = ,ℓ ( ) while the choice ( ) =̃︀ ℓ ( ) only gives a first order convergence behavior.

Numerical illustrations in 2D
In this section we present numerical results in 2D that illustrate the convergence behavior of the scheme (4.10) as well as the performance gain one could achieved with such method with = ℓ where ℓ ∈ {2, 3} and = 0.1. We systematically use Galerkin high order spectral finite elements on quadrangles (see [10]) for the construction of the spaces ,ℎ and ,ℎ , whereas for the construction of ℎ we use discontinuous nodal finite elements (on Gauss points, see [29]) on the edges defined by the intersection of the boundary Σ and the coarse mesh of Ω . Numerical results are obtained with an in-house solver developed in C++, parallelization is done using OpenMP and Eigen [43] is used as a linear algebra library.

Convergence result
We consider the wave equation (2.1), with homogeneous Neumann boundary condition, posed on the domain Ω = (0, 1) 2 and Ω = (0.6, 0.7) 2 . We assume that = 1 and = 1, so that the domain decomposition is artificial, but allow us to compare the simulations to an accurate reference solution. We consider the propagation of initial data. More precisely, the source term is null and the initial data are computed so that the reference solution is given by, for ∈ The solution can be computed semi-analytically using the expression of the Green's function for the wave equation in R 2 . The meshes of Ω and Ω are discretised uniformly with ℎ /ℎ = 8. The order used in Ω (resp. Ω ) is 8 (resp. 4). Consequently, the meshes are not conform. Order 7 is used to construct the space ℎ . Using the macro-element analysis of [28] one can check that (2.10) holds. Indeed it is sufficient to check an injectivity property on a reference configuration (see Eq. (30) in [28]). In our case, such injectivity property is guaranteed by the high number of test functions of ℎ, on the transmission boundary. The time step is chosen as where = 0.99 and (˜, ℎ ,˜, ℎ ) are approximations of ( ,ℎ , ,ℎ ) computed by power iteration. The approximated spectral radii are close enough to the exact spectral radii such that (6.2) implies the stability condition (4.13). Figure 7 displays a snapshot of the solution at different times. The relative error is computed at the final simulated time. Table 2 shows the relative errors as the discretisation diminishes. A second order convergence rate is observed.

Performance illustration
The chosen configuration is the same as in the previous section expect that the domain Ω differs. Indeed, in this configuration the domain Ω includes a circular hole that is meshed with quadrangles (see Figs 8 and 9). Order 10 finite elements are used in Ω while order 9 are used to construct the space ℎ and second order finite elements are used in Ω . We aim at comparing the computational cost of two schemes: (4.10) with = 2 and with = 3 where = 0.1 (the time step being chosen by formula (6.2)) and (4.10) with equals to the polynomial 1. In this latter case the scheme reads   ,ℎ thanks to the second equation. We give in Table 3 the parameters of the experiments as well as the measured times spent in each steps (in relative values). In the specific chosen configuration the proposed algorithms reaches the final simulation time, approximately 2.27 (resp. 2.63) times faster for ℓ = 2 (resp. ℓ = 3). One can see in Table 3 that the time spent in the computation of ,* ,ℎ is constant for the three methods whereas the time spent to compute ,* ,ℎ , +1 ,ℎ and ℎ decreases by a factor 3 and 4 approximately which is a consequence of the time step increase. Finally note that the time spent in Step 3f (computation of +1 ,ℎ ) increases with ℓ without penalizing too much the performances, so does Step 0 (computation and factorization of the matrix used in Step 2). This latter step representing a marginal effort in term of computations.
To sum-up, this case well illustrates the actual performance gain one can expect with this method. Note that this gain tends to 3 (resp. 4) for ℓ = 2 (resp. 3) as Dim( ,ℎ )/Dim( ,ℎ ) tends to zero.

Comparisons with existing approaches
7.1. The fully explicit Local Time Stepping of [26] The computational code used to obtain the results of this section is available as supplementary material at the web link [44].
In [19,26], an explicit Local Time Stepping Algorithm is proposed and is proved to be second order convergent for the 2 -norm. It is used in the context of solving the following semi-discrete wave equation: From formula (12) of [26] with = 2 one can derive the following scheme where ℎ : ℎ ↦ → ℎ is a restriction operator on a region discretised with a fine grid (with an overlap of one element in our computations). Notice that this algorithm amounts to solving a leap-frog scheme for the kernel of ℎ and to a modified scheme for the complement. Moreover, from Algorithm 1 (page 1000) of [26] we can deduce a variant of (EX-2b), where˜ℎ can be defined in two ways, denoting The choice (7.2a) gives exactly the Algorithm 1 of [26] while the more simple choice (7.2b) gives similar observed convergence behavior and so will be used in what follows. In the following we present a numerical assessment, in a one-dimensional setting, of local time stepping procedures that have the same computational cost: the stabilised explicit scheme (4.6) of Section 4.2, the scheme (EX-2) and its variant (EX-2b).
The considered case is the same as in Section 5 with = 1 (homogeneous medium). More precisely, we solve up to time = 0.5, the equation with homogeneous Neumann boundary condition. The discretisation parameters are = 0.99 and = 2. The purpose of these tests is to quantify the relative 2 and 1 -errors with respect to three chosen continuous analytical solutions associated with adequate source terms or initial data.

Propagating pulse
The first considered case is a propagating pulse as described in Section 5. In Figure 10 are displayed the relative 2 and 1 -errors between the solutions of the three numerical schemes (Scheme (4.6), (EX-2) and (EX-2b)) and the analytical solution, with respect to the mesh size ℎ (note that ∆ goes to zero with ℎ because of Assumption 3.2).

Quasi-static solution
We choose vanishing initial data for the wave equation Note that the solution reaches a static state after time ≥ 0.35. The obtained convergence curves are displayed in Figure 11. One can see that the three schemes behave similarly in terms of convergence in the 2 -norm, however the scheme (EX-2b) is less accurate in the 1 -norm. More precisely, half an order of convergence is lost.

Spatially constant solution
We choose vanishing initial data for the wave equation (7.3) and choose such that the solution is smooth and given by ( , ) = ℎ (︂ − 1 )︂ · with = 0.1 and 1 = 0.8. The analytical solution is therefore constant in space. The obtained convergence curves are displayed in Figure 12. Again, one the one hand, one can see that the three schemes behave similarly in terms of convergence in the 2 -norm. On the other hand, it is this time the scheme (EX-2) which is less accurate in the 1 -norm (half an order of convergence is lost).

The implicit LTS of [17]
In equations (13.74)-(13.76) of [17], we find an algebraic formulation for conservative local time stepping. This formulation is written for the system of elastodynamics written in a first order form in time. However, by elimination of the variable corresponding to the velocity, one can show that the algebraic formulation is This new formulation of system (7.4) shows that the local time stepping proposed in [17] is in fact equivalent to the scheme developed in Section 4.2 (compare the above equations with (4.6)). Therefore the local time stepping proposed in [17] can be seen as a transmission problem between two second order schemes, one of which having a relaxed stability condition by adding stabilising terms. The computational burden of the schemes we propose in Section 4.3 is equivalent to the one of the schemes proposed in [17]. In fact we conjecture that, the local time stepping of [17] can be recast in the formalism of Section 4.3 with ≡̃︀ ,ℓ . It has to be noted that the schemes in [17] are not proven to be second order convergent (in space and time) for the 1 -norm which is in accordance with the convergence results of Section 5.3.

Conclusions
In this work we have presented and analysed a family of second order in time discretisation strategy for linear wave equations. We have shown that they correspond to either locally implicit schemes or to local time stepping. For the analysis we have considered the case of smooth solutions. Then, we have shown that, if a well-defined stability condition holds, then second order space-time convergence property holds in the context of abstract Galerkin approximations of the wave equation. Finally, we have presented 1D numerical convergence results that confirm the obtained theoretical results. As far as local time stepping strategies are concerned, after comparisons with existing methods we have confirmed the interest of the proposed approach since, in terms of accuracy, it yields second order convergence results in the 1 -norm (contrary to [17,19]) and in terms of computational cost, it is equivalent to the method proposed in [17].

Appendix A. Convergence result for the semi-discrete problem
In Section 2.1 we have assumed that the solution to the continuous problem satisfies some extra-regularity in time. One can expect that, in general, existence and uniqueness results of solutions for source term and initial data do not immediately provide the adequate time regularity. Such extra regularity can be obtained in several ways. A simple one is to assume that initial data are null and the source term is smooth and vanishes at the initial time. More precisely, let ∈ N * and assume ∈ ,1 (0, ; ) and ( ) (0) = 0 for ∈ N such that < , then it is expected that the unique solution of Problem 2.4 with null initial data belongs to Such result is obtained by straightforward differentiation in time of the variational formulation. If initial data do not vanish then the question of the regularity of solution in time is more intricate since the initial data should be sufficiently smooth in space (or in a abstract framework, it should belong to the domain of some unbounded operator). Moreover, the same question arises at the semi discrete level, since it is not clear that the continuous result in general transfers to the semi-discrete setting. Existence and uniqueness results for the semi-discrete problem are generally a direct consequence of the results obtained at the continuous level, and the uniform estimates of Assumption 2.2 are obtained by energy techniques.
To obtain smooth semi-discrete solutions one can assume again that initial data are null and the source term is smooth and vanishes at the initial time, then one can define the source term for the semi-discrete problem as follows: the source term ℎ ( ) is chosen at every time as the orthogonal projection of ( ) in ℎ with respect to the scalar product of , i.e., ℎ ( ) ∈ ℎ and ( ℎ ( ), ℎ ) = ( ( ), ℎ ), ℎ ∈ ℎ . (A.2) From this definition we deduce that if ∈ ,1 (0, ; ) then ℎ ∈ ,1 (0, ; ℎ ) and Then, one can prove that the semi-discrete solution is smooth in time by first differentiating in time the semi-discrete variational formulation and second, by using the standard existence, uniqueness and regularity results for wave propagation problems for the successive time derivatives of the solution.
For the sake of completeness we now give the proof of Theorem 2.3 that is rather standard and is inspired from the one given in [1].
The purpose of this section is to provide a proof of the discrete Gronwall's lemma (3.42). Let Since 0 = we get the expected result.