Florent Pled
email: florent.pled@univ-eiffel.fr

Christophe Desceliers
email: christophe.desceliers@univ-eiffel.fr

Review and recent developments on the perfectly matched layer (PML) method for the numerical modeling and simulation of elastic wave propagation in unbounded domains

Keywords: Elastic wave propagation, Elastodynamics, Unbounded domain, Perfectly matched layers (PMLs), Absorbing boundary conditions (ABCs), Absorbing layers (ALs), Finite

This review article revisits and outlines the perfectly matched layer (PML) method and its various formulations developed over the past 25 years for the numerical modeling and simulation of wave propagation in unbounded media. Based on the concept of complex coordinate stretching, an efficient mixed displacement-strain unsplit-field PML formulation for second-order (displacementbased) linear elastodynamic equations is then proposed for simulating the propagation and absorption of elastic waves in unbounded (infinite or semi-infinite) domains. Both time-harmonic (frequency-domain) and time-dependent (timedomain) PML formulations are derived for two-and three-dimensional linear elastodynamic problems. Through the introduction of only a few additional variables governed by low-order auxiliary differential equations, the resulting mixed timedomain PML formulation is second-order in time, thereby allowing the use of standard time integration schemes commonly employed in computational structural dynamics and thus facilitating the incorporation into existing displacementbased finite element codes. For computational efficiency, the proposed timedomain PML formulation is implemented using a hybrid approach that couples a mixed (displacement-strain) formulation for the PML region with a classical (displacement-based) formulation for the physical domain of interest, using a standard Galerkin finite element method (FEM) for spatial discretization and a Newmark time scheme coupled with a finite difference (Crank-Nicolson) time scheme for time sampling. Numerical experiments show the performances of the PML method in terms of accuracy, efficiency and stability for two-dimensional linear elastodynamic problems in single-and multi-layer isotropic homogeneous elastic media.

Mathematics Subject Classification

(2010) 35L05 • 35L20 • 35B35 • 65M60 • 65M12 • 74J05 • 74J10 • 74S05 • 74H15 1 Introduction
Large-scale complex wave propagation, radiation or scattering problems are often set on unbounded (or very large) domains. When using classical numerical approximation methods such as finite difference (FD), finite volume (FV), finite element (FE), spectral element (SE) or discontinuous Galerkin (DG) methods, their numerical solutions can be computed in using computational models that are constructed for bounded subdomains (corresponding to the physical domains of interest), thereby requiring the truncation of the unbounded domain with ad hoc modeling artifacts, i.e. appropriate conditions at the artificial truncation boundaries, to properly simulate the unimpeded outward propagation of waves without any spurious reflections coming from the truncation boundaries. Among the numerous methods dealing with the construction of such bounded (or finite) computational domains, such as artificial transparent and absorbing boundary conditions (ABCs), absorbing layers (ALs) methods and alternative techniques (see [137,333,321,[START_REF] Aguilar | A class of analytic perturbations for one-body Schrödinger Hamiltonians[END_REF]149,151,[START_REF] Colonius | Modeling artificial boundary conditions for compressible flow[END_REF]155,140,[START_REF] Antoine | A Review of Transparent and Artificial Boundary Conditions Techniques for Linear and Nonlinear Schrödinger Equations[END_REF]160,127] and the references therein for a comprehensive overview of the existing approaches), the high-order local ABCs (involving only low-order derivatives, i.e. no high-order spatial and temporal derivatives, having recourse to special auxiliary variables on the artificial boundary) [143,144,154,[START_REF] Antoine | Bayliss-Turkel-like Radiation Conditions on Surfaces of Arbitrary Shape[END_REF]148,139,156,141,147,155,152,157,[START_REF] Bécache | High-order absorbing boundary conditions for anisotropic and convective wave equations[END_REF][START_REF] Baffet | Long-time stable highorder absorbing boundary conditions for elastodynamics[END_REF]153,296] and the perfectly matched layers (PMLs) [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF]189,[START_REF] Bérenger | Three-Dimensional Perfectly Matched Layer for the Absorption of Electromagnetic Waves[END_REF]170,163,322,318,[START_REF] Collino | Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF][START_REF] Appelö | Perfectly Matched Layers for Hyperbolic Systems: General Formulation, Well-posedness, and Stability[END_REF][START_REF] Bérenger | Perfectly matched layer (PML) for computational electromagnetics[END_REF]174,[START_REF] Bermúdez | Perfectly Matched Layers for Time-Harmonic Second Order Elliptic Problems[END_REF]200,114] have been extensively studied during the last decades, since they are neither computationally expensive, nor difficult to implement, nor limited to domains with simple geometric shapes (e.g. rectangular, cuboidal, circular, spherical, ellipsoidal, etc.), and allow the outgoing waves to be efficiently absorbed with arbitrarily-high controlled accuracy for a broad class of wave propagation problems.

The PML method has been introduced in the mid 1990s for simulating the propagation of electromagnetic waves over time in two-dimensional (2D) [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Berenger | Perfectly matched layer for the FDTD solution of wave-structure interaction problems[END_REF] and three-dimensional (3D) [189,[START_REF] Bérenger | Three-Dimensional Perfectly Matched Layer for the Absorption of Electromagnetic Waves[END_REF] unbounded spatial domains. It basically consists in surrounding the bounded physical domain of interest (resulting from the truncation of the unbounded domain) by artificial (non-physical) absorbing nonreflecting layers of finite thickness especially designed to absorb the outgoing waves (leaving the physical domain) without any spurious (unwanted) reflections. In theory, within the continuous framework, a PML model is then constructed such that the outgoing waves are perfectly transmitted at the interface between the physical domain and the absorbing layers whatever their non-zero angles of incidence and their non-zero angular frequencies, and then they decrease exponentially in the direction normal to the interface when propagating inside the absorbing layers. Nevertheless, in practice, within the discrete framework, the numerical approximations introduced in the PML computational model do not allow such a perfect matching (reflectionless) property to be still satisfied due to the numerical dispersion inherent to any discretization method. Furthermore, since the absorbing layers (hereinafter referred to as the PML region) have a limited (finite) thickness, the outgoing waves that are not completely damped while traveling inside the PML region, being reflected at the exterior boundary and again propagating back through the PML region (hence traveling twice through it) re-enter the physical domain and therefore generate spurious waves with a non-zero amplitude, which is relatively small in practical applications since it is damped by an exponential factor that depends on twice the PML thickness. In numerical practice, the layers are then designed to have minimal reflection and strong absorption properties. Besides, due to its simplicity of implementation and versatility (broad applicability), the PML method can be easily implemented in conjunction with any classical numerical approximation method available in standard computer codes and often performs very well compared to other existing ABCs for artificially handling problems defined on unbounded domains [187,75,162,[START_REF] Collino | Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF]193,140,276,354,341,127], especially for arbitrarily heterogeneous materials/media.

Overview of existing PML formulations

The classical PML formulation, henceforth known as split-field PML (sometimes referred to as S-PML), is built by artificially splitting the physical fields solution of the time-domain wave equations (written in Cartesian coordinates) into unphysical components (according to their spatial derivatives) in the PML region and considering specific dissipative anisotropic material properties. These dissipative material properties are introduced through the use of artificial damping terms satisfying a matching impedance condition to achieve a perfect transmission of outgoing propagative waves (without any reflections) at the interface between the physical domain and the absorbing layers, with an exponential decay in the PML region along the direction normal to the interface [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Berenger | Perfectly matched layer for the FDTD solution of wave-structure interaction problems[END_REF][START_REF] Bérenger | Three-Dimensional Perfectly Matched Layer for the Absorption of Electromagnetic Waves[END_REF]. Such a field-splitting PML formulation alters the original formulation of second-order time-domain wave equations and introduces two distinct systems of equations set on the physical domain and the PML region, requiring a special treatment at the interface between the two and making the implementation into existing codes rather cumbersome. In the frequency domain, the PML method can be interpreted as a complex space coordinate stretching approach [69,284,285,[START_REF] Chew | Perfectly matched layers for elastodynamics: A new absorbing boundary condition[END_REF][START_REF] Chew | A 3D perfectly matched medium from modified maxwell's equations with stretched coordinates[END_REF] through an analytic continuation of the real spatial coordinates into the complex space via complex stretching functions, corresponding to a change of metric of the physical coordinate space [308,315,205,207]. As already mentioned in [START_REF] Bermúdez | Perfectly Matched Layers for Time-Harmonic Second Order Elliptic Problems[END_REF], the complex coordinate stretching strategy is closely related to the so-called complex scaling technique (also referred to as analytic dilatation technique) originally introduced in the mid 1970s for the numerical simulation and analysis of the dynamical properties (eigenstates) of atomic systems in quantum mechanics [8,300,168]. The resulting PML formulations based on complex stretched coordinates then involve complex operators, rendering the PML method difficult to implement along with any conventional numerical approximation method. Alternative unsplit-field (or nonsplitting) PML formulations, frequently called uniaxial PMLs (sometimes referred to as U-PML or N-PML), preserve the original (unstretched) formulation of second-order time-harmonic wave equations (involving physical fields and classical real operators), but with specific frequency-dependent complex anisotropic material properties, thereby allowing the PML region to be interpreted as an artificial anisotropic absorbing medium (sometimes called anisotropic material absorber) [291,130,129,356,357,353,363,287,181,310,311,313,315,2]. The unsplit-field PML formulations can then be straightforwardly incorporated into existing codes (with minor modifications) and implemented using any classical numerical method unlike their split-field counterparts. Time-dependent PML formulations are usually derived from their time-harmonic counterparts by taking the inverse Fourier transform in time of time-harmonic PML equations, which can be performed by directly using convolution products that are usually computed by means of a recursive convolution update technique [228] to avoid expensive temporal convolution operations [286,288,329,290,326,[START_REF] Appelö | A new absorbing layer for elastic waves[END_REF]100,192,235,236,234,146,216,240,136,335,270,263,214,215,230] (henceforth known as the convolutional PML or C-PML), or by using a recursive time integration approach [347,[START_REF] Drossaert | A nonsplit complex frequency-shifted PML based on recursive integration for FDTD modeling of elastic waves[END_REF]133,134,135] (henceforth known as the recursive integration PML or RI-PML), or by using digital signal processing (DSP) techniques such as the Z-transform and digital filtering methods [282,279,280,281,210,212,213,328,116,118,120,121] (henceforth known as the matched Z-transform PML or MZT-PML), or by introducing additional fields and auxiliary differential equations [322,171,359,193,278,[START_REF] Basu | Perfectly matched layers for transient elastodynamics of unbounded domains[END_REF]302,325,211,[START_REF] Basu | Explicit finite element perfectly matched layer for transient three-dimensional elastic waves[END_REF]198,275,131,237,350,145,104,105,117,185,[START_REF] Benito | A GFDM with PML for seismic wave equations in heterogeneous media[END_REF]101,107,108,103,292,335,352,[START_REF] Ma | A novel unsplit perfectly matched layer for the second-order acoustic wave equation[END_REF]128,330,331,360,248] (henceforth known as the auxiliary-differential-equation PML or ADE-PML), and thus require additional resources in terms of computational cost and memory storage capacity. Note that ADE-PML formulations provide more flexibility allowing for an easy extension to high-order time integration schemes and are also easier to implement with classical numerical methods as compared to C-PML formulations [131,237,350,292,335,352]. In addition, such ADE-PML formulations can be seen as perturbations of the original (unstretched) formulation of time-domain equations since they reduce to solving the original system of time-domain governing equations over the physical domain by considering all additional fields (auxiliary variables) set to zero and adding the non-zero source terms, boundary and initial conditions imposed over the physical domain of interest.

Some mathematical analysis results on the existence and convergence properties of the PML method can be found in [75,206,207,169,[START_REF] Bao | Convergence Analysis of the Perfectly Matched Layer Problems for Time-Harmonic Maxwell's Equations[END_REF][START_REF] Bramble | Analysis of a finite PML approximation for the three dimensional time-harmonic Maxwell and acoustic scattering problems[END_REF]. The well-posedness and stability properties1 of the PML method have also been studied in numerous research works [255,[START_REF] Abarbanel | A Mathematical Analysis of the PML Method[END_REF]2,[START_REF] Abarbanel | Well-posed Perfectly Matched Layers for Advective Acoustics[END_REF]265,319,[START_REF] Abarbanel | Long Time Behavior of the Perfectly Matched Layer Equations in Computational Electromagnetics[END_REF][START_REF] Bécache | On the analysis of Bérenger's Perfectly Matched Layers for Maxwell's equations[END_REF][START_REF] Diaz | Stabilized Perfectly Matched Layer for Advective Acoustics[END_REF][START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF][START_REF] Bécache | On the long-time behavior of unsplit perfectly matched layers[END_REF][START_REF] Bécache | Perfectly Matched Layers for the Convected Helmholtz Equation[END_REF]277,[START_REF] Appelö | Perfectly Matched Layers for Hyperbolic Systems: General Formulation, Well-posedness, and Stability[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF][START_REF] Appelö | A General Perfectly Matched Layer Model for Hyperbolic-Parabolic Systems[END_REF][START_REF] Diaz | A time domain analysis of PML models in acoustics[END_REF]242,243,244,[START_REF] Dmitriev | Application of M-PML absorbing boundary conditions to the numerical simulation of wave propagation in anisotropic media. Part II: Stability[END_REF]226,294,295,158,[START_REF] Chen | Long-Time Stability and Convergence of the Uniaxial Perfectly Matched Layer Method for Time-Domain Acoustic Scattering Problems[END_REF][START_REF] Bécache | Remarks on the stability of Cartesian PMLs in corners[END_REF]183,104,105,185,195,[START_REF] Demaldent | Perfectly matched transmission problem with absorbing layers: Application to anisotropic acoustics in convex polygonal domains[END_REF]101,103,102,159,[START_REF] Bécache | Perfectly matched layers in negative index metamaterials and plasmas[END_REF][START_REF] Bécache | Stable perfectly matched layers for a cold plasma in a strong background magnetic field[END_REF][START_REF] Bécache | Stable perfectly matched layers for a class of anisotropic dispersive models. Part I: necessary and sufficient conditions of stability[END_REF]184,[START_REF] Benito | Wave propagation in soils problems using the Generalized Finite Difference Method[END_REF]293,126,[START_REF] Bécache | On the analysis of perfectly matched layers for a class of dispersive media and application to negative index metamaterials[END_REF]. On the one hand, for wave propagation problems governed by symmetric hyperbolic systems (e.g. in electromagnetics, acoustics and elastodynamics), classical split-field PML formulations result in weakly hyperbolic systems and lead to only weakly well-posed problems, which may become ill-posed under some low-order (small) perturbations on the split fields [START_REF] Abarbanel | A Mathematical Analysis of the PML Method[END_REF]2,167,[START_REF] Bécache | On the analysis of Bérenger's Perfectly Matched Layers for Maxwell's equations[END_REF]183,195], thus making numerical computations impracticable. On the other hand, unsplit-field PML formulations result in symmetric strongly hyperbolic systems and lead to strongly well-posed problems [2,[START_REF] Abarbanel | Well-posed Perfectly Matched Layers for Advective Acoustics[END_REF]265,277,158,104]. Furthermore, standard PML systems derived from both split-and unsplitfield PML formulations have been proved to be only weakly stable under some suitable conditions (depending on the symbol of the hyperbolic partial differential operator) [START_REF] Abarbanel | Long Time Behavior of the Perfectly Matched Layer Equations in Computational Electromagnetics[END_REF][START_REF] Bécache | On the analysis of Bérenger's Perfectly Matched Layers for Maxwell's equations[END_REF][START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF][START_REF] Bécache | On the long-time behavior of unsplit perfectly matched layers[END_REF]193,[START_REF] Chen | Long-Time Stability and Convergence of the Uniaxial Perfectly Matched Layer Method for Time-Domain Acoustic Scattering Problems[END_REF], which means exponentially growing waves are not supported by the time-domain PML formulation, but linear (or polynomial) growing waves may potentially exist and be triggered for instance by numerical integration (quadrature) or round-off (machine precision) errors [335]. Such growing waves can then lead to spurious instabilities that spread and adversely pollute the solution within the physical domain in long-time simulations. Despite its general success and satisfactory performances in most practical applications, the PML method has some apparent drawbacks and may fail to absorb the outgoing waves in some specific situations [170,307,167,[START_REF] Abarbanel | Well-posed Perfectly Matched Layers for Advective Acoustics[END_REF]171,[START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF][START_REF] Diaz | Stabilized Perfectly Matched Layer for Advective Acoustics[END_REF]254,[START_REF] Correia | 3D-FDTD-PML analysis of left-handed metamaterials[END_REF][START_REF] Cummer | Perfectly matched layer behavior in negative refractive index materials[END_REF][START_REF] Dong | Perfectly matched layer-absorbing boundary condition for left-handed materials[END_REF]298,[START_REF] Diaz | A time domain analysis of PML models in acoustics[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF]303,192,236,242,259,226,[START_REF] Appelö | A high-order super-grid-scale absorbing layer and its application to linear hyperbolic systems[END_REF]216,[START_REF] Dmitriev | Application of M-PML absorbing boundary conditions to the numerical simulation of wave propagation in anisotropic media. Part II: Stability[END_REF]104,195,[START_REF] Dhia | On the use of perfectly matched layers in the presence of long or backward propagating guided elastic waves[END_REF]269,268,270,[START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF][START_REF] Bécache | Perfectly matched layers in negative index metamaterials and plasmas[END_REF]126,214,215,184,[START_REF] Bécache | Stable perfectly matched layers for a cold plasma in a strong background magnetic field[END_REF][START_REF] Bécache | Stable perfectly matched layers for a class of anisotropic dispersive models. Part I: necessary and sufficient conditions of stability[END_REF][START_REF] Bécache | On the analysis of perfectly matched layers for a class of dispersive media and application to negative index metamaterials[END_REF]. In particular, numerical instabilities may occur in time-domain simulations for some wave propagation problems in anisotropic and/or dispersive media (such as negative refractive index metamaterials, also called left-handed or double negative metamaterials, and cold plasmas, also called non-thermal plasmas), e.g. in (visco)elastodynamics [START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF]303,192,303,236,242,[START_REF] Appelö | A high-order super-grid-scale absorbing layer and its application to linear hyperbolic systems[END_REF]216,[START_REF] Dmitriev | Application of M-PML absorbing boundary conditions to the numerical simulation of wave propagation in anisotropic media. Part II: Stability[END_REF]104,[START_REF] Dhia | On the use of perfectly matched layers in the presence of long or backward propagating guided elastic waves[END_REF]269,268,270,[START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF][START_REF] Bécache | Perfectly matched layers in negative index metamaterials and plasmas[END_REF]126,214,215], electromagnetics [START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF][START_REF] Correia | 3D-FDTD-PML analysis of left-handed metamaterials[END_REF][START_REF] Cummer | Perfectly matched layer behavior in negative refractive index materials[END_REF][START_REF] Dong | Perfectly matched layer-absorbing boundary condition for left-handed materials[END_REF]298,259,226,[START_REF] Bécache | Perfectly matched layers in negative index metamaterials and plasmas[END_REF]184,[START_REF] Bécache | Stable perfectly matched layers for a cold plasma in a strong background magnetic field[END_REF][START_REF] Bécache | Stable perfectly matched layers for a class of anisotropic dispersive models. Part I: necessary and sufficient conditions of stability[END_REF][START_REF] Bécache | On the analysis of perfectly matched layers for a class of dispersive media and application to negative index metamaterials[END_REF], aeroacoustics [170,307,167,[START_REF] Abarbanel | Well-posed Perfectly Matched Layers for Advective Acoustics[END_REF]171,[START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF][START_REF] Diaz | Stabilized Perfectly Matched Layer for Advective Acoustics[END_REF][START_REF] Diaz | A time domain analysis of PML models in acoustics[END_REF][START_REF] Demaldent | Perfectly matched transmission problem with absorbing layers: Application to anisotropic acoustics in convex polygonal domains[END_REF], geophysical fluid dynamics [254] and even in anisotropic scalar wave propagation problems [195]. Such instabilities are due to the presence of backward propagating waves (for which the phase and group velocities point in opposite directions with respect to the interface, i.e. they are not oriented in the same way and then have differing signs with respect to the direction normal to the interface) exponentially growing within the PML region and producing strong spurious reflections within the physical domain. Even for isotropic media, the PML method may not only suffer from long-time instabilities but also exhibit instabilities for near-grazing incident propagating waves or evanescent waves [250,[START_REF] Bérenger | An effective PML for the absorption of evanescent waves in waveguides[END_REF][START_REF] Bérenger | Evanescent waves in PML's: origin of the numerical reflection in wavestructure interaction problems[END_REF]75,221,332,123,233,193,124,122,[START_REF] Drossaert | A nonsplit complex frequency-shifted PML based on recursive integration for FDTD modeling of elastic waves[END_REF]100,192,235,236,234,237,216,240,242,243,244,350,352]. Besides, the PML method also performs very poorly in the presence of low-frequency propagating waves or long waves (for which the wavenumber is almost zero) arising near cut-off frequencies and decaying slowly within the PML region [303,[START_REF] Dhia | On the use of perfectly matched layers in the presence of long or backward propagating guided elastic waves[END_REF].

Several methods have been proposed to avoid and remove growing waves and improve the stability and accuracy of PML formulations in the presence of evanescent, near-grazing incident propagating or long waves, e.g. by applying a numerical filtering (low-pass filter) [170,254], by introducing a complex frequency shift (known as the complex-frequency-shifted PML or CFS-PML) [202,129,221,320,286,288,[START_REF] Bérenger | Numerical reflection from FDTD-PMLs: a comparison of the split PML with the unsplit and CFS PMLs[END_REF][START_REF] Bérenger | Application of the CFS PML to the absorption of evanescent waves in waveguides[END_REF][START_REF] Bécache | On the long-time behavior of unsplit perfectly matched layers[END_REF]124,122,[START_REF] Correia | A simple and efficient implementation of CFS-PML in the FDTD analysis of periodic structures[END_REF]210,[START_REF] Appelö | Perfectly Matched Layers for Hyperbolic Systems: General Formulation, Well-posedness, and Stability[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF]326,325,[START_REF] Drossaert | A nonsplit complex frequency-shifted PML based on recursive integration for FDTD modeling of elastic waves[END_REF]100,192,133,235,236,234,275,131,237,350,216,240,136,328,272,273,120,335,104,105,101,108,103,128,336,230], by adding artificial stabilizing parameters (with loss of the perfect matching property) and/or making use of ad hoc space-time coordinate transformations [307,167,322,[START_REF] Abarbanel | Well-posed Perfectly Matched Layers for Advective Acoustics[END_REF]171,[START_REF] Abarbanel | Long Time Behavior of the Perfectly Matched Layer Equations in Computational Electromagnetics[END_REF][START_REF] Abarbanel | Unsplit Variables Perfectly Matched Layers for the Shallow Water Equations with Coriolis Forces[END_REF]150,[START_REF] Diaz | Stabilized Perfectly Matched Layer for Advective Acoustics[END_REF][START_REF] Bécache | Perfectly Matched Layers for the Convected Helmholtz Equation[END_REF]173,[START_REF] Bécache | Perfectly Matched Layers for Time-Harmonic Acoustics in the Presence of a Uniform Flow[END_REF][START_REF] Diaz | A time domain analysis of PML models in acoustics[END_REF]262], or by considering different multi-dimensional damping functions (known as the multiaxial PML or M-PML) [242,243,244,344,[START_REF] Dmitriev | Application of M-PML reflectionless boundary conditions to the numerical simulation of wave propagation in anisotropic media. Part I: Reflectivity[END_REF][START_REF] Dmitriev | Application of M-PML absorbing boundary conditions to the numerical simulation of wave propagation in anisotropic media. Part II: Stability[END_REF]114,269,268,270,126]. To overcome (or at least alleviate) the stability issues, the complex-frequency-shifted PML (CFS-PML) method consists in using a more sophisticated (frequency-shifted) complex coordinate stretching with a frequency-dependent real part (in addition to the standard frequency-dependent imaginary part) by simply shifting the frequency-dependent pole (zero-frequency singularity) of the complex stretching function off the real axis and onto the negative imaginary axis of the complex plane, so that the PML region can act as a low-pass ) filter with a cut-off angular frequency ω c to efficiently improve the absorption of near-grazing incident propagating waves [START_REF] Bérenger | Application of the CFS PML to the absorption of evanescent waves in waveguides[END_REF]124,[START_REF] Drossaert | A nonsplit complex frequency-shifted PML based on recursive integration for FDTD modeling of elastic waves[END_REF]192,235,236,234,237,216,350] and also the attenuation of evanescent waves [221,286,288,[START_REF] Bérenger | Numerical reflection from FDTD-PMLs: a comparison of the split PML with the unsplit and CFS PMLs[END_REF][START_REF] Bérenger | Application of the CFS PML to the absorption of evanescent waves in waveguides[END_REF]122,[START_REF] Correia | Performance of regular PML, CFS-PML, and second-order PML for waveguide problems[END_REF]100]. The PML region then depends on the angular frequency ω since it behaves like the physical medium (losing its absorption properties) at low frequencies (ω ω c ) and like a dissipative layer (losing its stabilizing properties) at high frequencies (ω ω c ), thus switching from a transparent behaviour (without attenuation) to a dissipative behaviour (with an exponential amplitude decay) in a frequency range around the cut-off frequency ω c [124,242,240]. The CFS-PML method can then be interpreted as a low-pass filter allowing to efficiently absorb high-frequency propagating waves (as it reduces asymptotically to the classical PML method at high frequencies) and improve the absorption of near-grazing incident propagating waves as well as evanescent waves, while slightly degrading the absorption of low-frequency propagating waves [288,[START_REF] Bérenger | Numerical reflection from FDTD-PMLs: a comparison of the split PML with the unsplit and CFS PMLs[END_REF][START_REF] Bérenger | Application of the CFS PML to the absorption of evanescent waves in waveguides[END_REF][START_REF] Bécache | On the long-time behavior of unsplit perfectly matched layers[END_REF]124,[START_REF] Correia | Performance of regular PML, CFS-PML, and second-order PML for waveguide problems[END_REF]242,240]. To overcome this potential limitation, a second-order PML method has been proposed in [START_REF] Correia | On the development of a higher-order PML[END_REF][START_REF] Correia | Performance of regular PML, CFS-PML, and second-order PML for waveguide problems[END_REF]227,131,237,116,117,118,119,115,330,331] and extended to a general higher-order PML method in [134,121]. It allows combining the advantages of both classical PML and CFS-PML methods (in terms of absorption performance) by simply considering a more complicated (multipole) stretching function defined as the product of individual CFS stretching functions. Despite a rather complex implementation in the time domain and more computational resource requirements (in terms of computational cost and memory storage capacity), the resulting higher-order PML method is then highly effective in absorbing both low-frequency propagating waves and strong evanescent waves, especially for open-region and periodic problems as well as for waveguide problems defined over elongated domains. More recently, another multiple-pole (or multipole) PML method has been introduced in [135] and is based on an ad hoc multipole complex stretching function simply defined as the sum of individual CFS stretching functions in order to facilitate the implementation and parameter optimization of the multipole PML method when compared to a higher-order PML method. Note that other generalized (multipole) complex stretching functions have also been considered in [START_REF] Bécache | On the analysis of perfectly matched layers for a class of dispersive media and application to negative index metamaterials[END_REF] to derive stable PML formulations for dispersive media (e.g. negative index metamaterials). As another attempt to prevent late-time instabilities occurring in some kinds of anisotropic elastic media, the multiaxial PML (M-PML) method extends the classical PML method by considering different multi-dimensional damping functions (that are proportional to each other in orthogonal directions and in which the ratios of damping functions are additional user-tunable correction parameters for improving stability), leading to a more general complex coordinate stretching involving extra stabilizing correction parameters, so that the outgoing waves are absorbed along all coordinate directions (i.e. not only along the direction normal to the interface) with different amplitude attenuation factors (decay rates). Nevertheless, it has been shown that such a M-PML method results in an improved sponge layer with better absorption and stability properties in long-time simulations for some particular 2D orthotropic (visco)elastic media [242,243,244,236,216,269,268,270,126] or 2D isotropic (visco)elastic media with high Poisson's ratios [344,269,270] or even 2D and 3D general anisotropic elastic media [126], but it is not perfectly matched at the interface with the physical domain even within the continuous framework [START_REF] Dmitriev | Application of M-PML reflectionless boundary conditions to the numerical simulation of wave propagation in anisotropic media. Part I: Reflectivity[END_REF][START_REF] Dmitriev | Application of M-PML absorbing boundary conditions to the numerical simulation of wave propagation in anisotropic media. Part II: Stability[END_REF]324,269,268,270,126]. In addition, let us mention that the CFS-PML and M-PML methods have been recently combined in [352,214,215] (hence termed the hybrid PML or H-PML) to further improve the stability in some strongly anisotropic (2D orthotropic) elastic media without degrading the accuracy (absorption capacity), thus combining the advantages of both CFS-PML and M-PML methods for maximizing both accuracy and stability through the optimization (fine tuning) of the PML parameters involved in the H-PML (i.e. the combined CFS-M-PML) stretching function. Lastly, we refer to [182,[START_REF] Bécache | Perfectly matched layers in negative index metamaterials and plasmas[END_REF] for further discussions about the limitations and failure cases of the PML method as well as some workarounds.

During the last 25 years, the PML method has been widely used for (transient) electromagnetic problems [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF]69,189,253,291,284,285,[START_REF] Berenger | Perfectly matched layer for the FDTD solution of wave-structure interaction problems[END_REF][START_REF] Bérenger | Three-Dimensional Perfectly Matched Layer for the Absorption of Electromagnetic Waves[END_REF]202,130,129,111,356,357,353,287,203,[START_REF] Bérenger | Improved PML for the FDTD solution of wave-structure interaction problems[END_REF][START_REF] Bérenger | An effective PML for the absorption of evanescent waves in waveguides[END_REF]187,363,181,232,338,[START_REF] Chew | A 3D perfectly matched medium from modified maxwell's equations with stretched coordinates[END_REF]309,310,311,312,[START_REF] Abarbanel | A Mathematical Analysis of the PML Method[END_REF]2,337,177,308,315,314,332,317,316,319,318,75,[START_REF] Collino | The Perfectly Matched Layer in Curvilinear Coordinates[END_REF]322,223,165,320,[START_REF] Bérenger | Evanescent waves in PML's: origin of the numerical reflection in wavestructure interaction problems[END_REF]264,265,286,288,110,220,219,[START_REF] Abarbanel | Long Time Behavior of the Perfectly Matched Layer Equations in Computational Electromagnetics[END_REF][START_REF] Bérenger | Numerical reflection from FDTD-PMLs: a comparison of the split PML with the unsplit and CFS PMLs[END_REF][START_REF] Bérenger | Application of the CFS PML to the absorption of evanescent waves in waveguides[END_REF]180,179,278,279,280,281,266,290,[START_REF] Bécache | On the analysis of Bérenger's Perfectly Matched Layers for Maxwell's equations[END_REF][START_REF] Bécache | On the long-time behavior of unsplit perfectly matched layers[END_REF][START_REF] Correia | 3D-FDTD-PML analysis of left-handed metamaterials[END_REF][START_REF] Cummer | Perfectly matched layer behavior in negative refractive index materials[END_REF][START_REF] Dong | Perfectly matched layer-absorbing boundary condition for left-handed materials[END_REF]142,302,306,[START_REF] Correia | A simple and efficient implementation of CFS-PML in the FDTD analysis of periodic structures[END_REF][START_REF] Correia | On the development of a higher-order PML[END_REF][START_REF] Correia | Performance of regular PML, CFS-PML, and second-order PML for waveguide problems[END_REF]227,298,[START_REF] Abarbanel | Non-Linear PML Equations for Time Dependent Electromagnetics in Three Dimensions[END_REF]326,325,260,261,210,211,212,213,[START_REF] Bérenger | Perfectly matched layer (PML) for computational electromagnetics[END_REF][START_REF] Donderici | Mixed Finite-Element Time-Domain Method for Transient Maxwell Equations in Doubly Dispersive Media[END_REF][START_REF] Donderici | Conformal Perfectly Matched Layer for the Mixed Finite Element Time-Domain Method[END_REF]258,133,190,259,226,272,273,[START_REF] Abarbanel | Long-Time Performance of Unsplit PMLs with Explicit Second Order Schemes[END_REF][START_REF] Dosopoulos | Interior Penalty Discontinuous Galerkin Finite Element Method for the Time-Dependent First Order Maxwell's Equations[END_REF]131,134,[START_REF] Alvarez | A Leap-Frog Discontinuous Galerkin Time-Domain Method for HIRF Assessment[END_REF]116,117,118,119,120,121,305,105,107,102,330,331,[START_REF] Bécache | Perfectly matched layers in negative index metamaterials and plasmas[END_REF][START_REF] Bécache | Stable perfectly matched layers for a cold plasma in a strong background magnetic field[END_REF][START_REF] Bécache | Stable perfectly matched layers for a class of anisotropic dispersive models. Part I: necessary and sufficient conditions of stability[END_REF]184,[START_REF] Bécache | On the analysis of perfectly matched layers for a class of dispersive media and application to negative index metamaterials[END_REF]135], then quickly extended and successfully applied to various wave propagation problems [245,222,[START_REF] Bermúdez | Validation of acoustic models for time-harmonic dissipative scattering problems[END_REF][START_REF] Bermúdez | Perfectly Matched Layers for Time-Harmonic Second Order Elliptic Problems[END_REF]239] in diverse scientific and engineering fields including hydrodynamics and aeroacoustics (e.g. arising in radiation and scattering problems) [170,340,225,167,274,307,322,164,[START_REF] Abarbanel | Well-posed Perfectly Matched Layers for Advective Acoustics[END_REF]162,171,188,150,[START_REF] Diaz | Stabilized Perfectly Matched Layer for Advective Acoustics[END_REF]277,[START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF][START_REF] Bécache | Perfectly Matched Layers for the Convected Helmholtz Equation[END_REF]348,172,173,251,252,161,[START_REF] Bécache | Perfectly Matched Layers for Time-Harmonic Acoustics in the Presence of a Uniform Flow[END_REF][START_REF] Diaz | A time domain analysis of PML models in acoustics[END_REF][START_REF] Appelö | Perfectly Matched Layers for Hyperbolic Systems: General Formulation, Well-posedness, and Stability[END_REF]299,343,[START_REF] Bermúdez | An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems[END_REF]174,175,194,262,[START_REF] Appelö | A General Perfectly Matched Layer Model for Hyperbolic-Parabolic Systems[END_REF]361,178,342,283,185,[START_REF] Demaldent | Perfectly matched transmission problem with absorbing layers: Application to anisotropic acoustics in convex polygonal domains[END_REF][START_REF] Ma | A novel unsplit perfectly matched layer for the second-order acoustic wave equation[END_REF]323,128], geophysical fluid dynamics (e.g. arising in atmospheric and oceanic sciences) [START_REF] Darblade | Conditions limites non réfléchissantes pour un modèle de Saint-Venant bidimensionnel barotrope linéarisé[END_REF]254,[START_REF] Abarbanel | Unsplit Variables Perfectly Matched Layers for the Shallow Water Equations with Coriolis Forces[END_REF]208,[START_REF] Barucq | New absorbing layers conditions for short water waves[END_REF]246], elastodynamics (e.g. arising in soil dynamics, soil-structure interaction, earth seismology, geotechnical site characterization and earthquake engineering, geophysical subsurface sensing/probing) [START_REF] Chew | Perfectly matched layers for elastodynamics: A new absorbing boundary condition[END_REF]163,351,221,[START_REF] Collino | Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF][START_REF] Bécache | Fictitious domains, mixed finite elements and perfectly matched layers for 2-D elastic wave propagation[END_REF]359,[START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF]193,123,233,329,224,[START_REF] Basu | Perfectly matched layers for time-harmonic elastodynamics of unbounded domains: theory and finite-element implementation[END_REF]347,[START_REF] Cohen | Mixed Spectral Finite Elements for the Linear Elasticity System in Unbounded Domains[END_REF]124,122,[START_REF] Bindel | Elastic PMLs for resonator anchor loss simulation[END_REF]229,161,[START_REF] Appelö | A new absorbing layer for elastic waves[END_REF]303,[START_REF] Drossaert | A nonsplit complex frequency-shifted PML based on recursive integration for FDTD modeling of elastic waves[END_REF]100,192,236,237,275,125,350,349,[START_REF] Basu | Explicit finite element perfectly matched layer for transient three-dimensional elastic waves[END_REF]216,242,243,244,186,199,200,201,240,344,104,328,191,354,[START_REF] Benito | A GFDM with PML for seismic wave equations in heterogeneous media[END_REF][START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF]335,269,352,114,113,292,108,106,103,[START_REF] Brun | Hybrid Asynchronous Perfectly Matched Layer for seismic wave propagation in unbounded domains[END_REF]263,360,[START_REF] Benito | Wave propagation in soils problems using the Generalized Finite Difference Method[END_REF]293,[START_REF] Benito | Implementations with generalized finite differences of the displacements and velocity-stress formulations of seismic wave propagation problem[END_REF]126,214,215], viscoelastodynamics [START_REF] Basu | Perfectly matched layers for transient elastodynamics of unbounded domains[END_REF][START_REF] Bécache | A mixed finite element approach for viscoelastic wave propagation[END_REF]234,268,270], poroelastodynamics [345,346,304,235,146,166], poroviscoelastodynamics [136], optics and quantum mechanics (e.g. arising in atomic, molecular and laser physics) [START_REF] Collino | Perfectly Matched Absorbing Layers for the Paraxial Equations[END_REF][START_REF] Ahland | Accurate mesh truncation for Schrödinger equations by a perfectly matched layer absorber: Application to the calculation of optical spectra[END_REF]271,209,151,112,[START_REF] Cheng | 3D quantum transport solver based on the perfectly matched layer and spectral element methods for the simulation of semiconductor nanodevices[END_REF]358,[START_REF] Dohnal | Perfectly matched layers in photonics computations: 1D and 2D nonlinear coupled mode equations[END_REF][START_REF] Antoine | A Review of Transparent and Artificial Boundary Conditions Techniques for Linear and Nonlinear Schrödinger Equations[END_REF][START_REF] Dohnal | Perfectly matched layers for coupled nonlinear Schrödinger equations with mixed derivatives[END_REF][START_REF] Antoine | Numerical solution of time-dependent nonlinear Schrödinger equations using domain truncation techniques coupled with relaxation scheme[END_REF]257,241,267,[START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF], fluid-porous medium coupled problems or fluid-solid interaction problems (e.g. arising in marine seismology, in non-destructive testing techniques such as ultrasound evaluation techniques for the characterization of complex biological systems) [START_REF] Bermúdez | Validation of acoustic models for time-harmonic dissipative scattering problems[END_REF][START_REF] Bermúdez | Perfectly Matched Layers for Time-Harmonic Second Order Elliptic Problems[END_REF][START_REF] Calvo | Simulation of acoustic scattering from an aluminum cylinder near a rough interface using the elastodynamic finite integration technique[END_REF]238,239,[START_REF] Assi | A perfectly matched layer formulation for modeling transient wave propagation in an unbounded fluid-solid medium[END_REF]336,297], as well as other general wave-like hyperbolic problems [217,[START_REF] Appelö | Perfectly Matched Layers for Hyperbolic Systems: General Formulation, Well-posedness, and Stability[END_REF][START_REF] Appelö | Application of a perfectly matched layer to the nonlinear wave equation[END_REF][START_REF] Appelö | A General Perfectly Matched Layer Model for Hyperbolic-Parabolic Systems[END_REF]. Initially derived in Cartesian coordinates for simply-shaped domains with straight (planar) artificial boundaries (e.g. squared, rectangular or cuboidal domains), some PML formulations have been extended to other (non-Cartesian) coordinate systems such as polar, cylindrical and spherical coordinates [232,338,[START_REF] Chew | A 3D perfectly matched medium from modified maxwell's equations with stretched coordinates[END_REF]309,310,339,75,[START_REF] Collino | The Perfectly Matched Layer in Curvilinear Coordinates[END_REF]312,165,221,314,355,317,265,359,266,224,[START_REF] Bermúdez | Validation of acoustic models for time-harmonic dissipative scattering problems[END_REF]199,292], and applied to generally-shaped convex domains using local polygonal or (non-)orthogonal curvilinear coordinates [311,308,177,316,319,205,207,220,219,124,364,260,261,343,[START_REF] Donderici | Conformal Perfectly Matched Layer for the Mixed Finite Element Time-Domain Method[END_REF]125,[START_REF] Dosopoulos | Interior Penalty Discontinuous Galerkin Finite Element Method for the Time-Dependent First Order Maxwell's Equations[END_REF]349,352,[START_REF] Alvarez | A Leap-Frog Discontinuous Galerkin Time-Domain Method for HIRF Assessment[END_REF]239,[START_REF] Demaldent | Perfectly matched transmission problem with absorbing layers: Application to anisotropic acoustics in convex polygonal domains[END_REF]323,248]. Such improvements provide more flexibility and choice in the shape and geometry of the underlying physical domain of interest. As already pointed out in [359], recall that the first attempts to extend the classical PML formulation (written in Cartesian coordinates) to non-orthogonal curvilinear coordinates (such as proposed in [253,334,287,203] for instance) were based on an approximate impedance matching condition derived under the assumption that the complex coordinate stretching function is independent of the spatial coordinates, hence resulting in an approximate PML formulation. Furthermore, it should be noticed that a straightforward extension of the classical PML formulation to cylindrical coordinates (such as proposed in [203] for instance) leads a loss of the perfect matching (reflectionless) property for curvilinear (non-planar) interfaces even within the continuous framework, hence resulting in an approximate PML formulation (henceforth known as the quasi-PML) [223,165]. In the present work, for the sake of simplicity, we restrict to rectangular or cuboidal domains with straight (planar) boundaries in Cartesian coordinates, but the proposed PML formulation could be easily extended and applied to more complex domains in other coordinate systems.

Overview of existing PML formulations in linear elastodynamics

Following this general overview of the PML method and its variants addressed in the literature, we now focus on the different PML formulations as well as their implementations in linear (visco)elastodynamics. Classical split-field PML formulations for the linear elastodynamic equations have been initially derived for the first-order mixed (velocity-stress) formulation [START_REF] Chew | Perfectly matched layers for elastodynamics: A new absorbing boundary condition[END_REF]163,351,221,[START_REF] Collino | Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF][START_REF] Bécache | Fictitious domains, mixed finite elements and perfectly matched layers for 2-D elastic wave propagation[END_REF][START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF]123,233,224,[START_REF] Cohen | Mixed Spectral Finite Elements for the Linear Elasticity System in Unbounded Domains[END_REF]229,125,349], then extended to the second-order (displacementbased) formulation (with a velocity field as additional variable to render the resulting mixed formulation second-order in time) [193,354] and to the second-order mixed (displacement-stress) formulation of the viscoelastodynamic equations [START_REF] Bécache | A mixed finite element approach for viscoelastic wave propagation[END_REF]. Classical unsplit-field PML formulations have also been obtained for the linear elastic wave equations formulated as a first-order mixed (velocity-stress) system [329,[START_REF] Benito | Implementations with generalized finite differences of the displacements and velocity-stress formulations of seismic wave propagation problem[END_REF] and as a second-order system in displacement or velocity (with potentially a stress, strain and/or other auxiliary field(s) as additional unknown(s), such as time-integral or memory-like terms, thus resulting in a second-order mixed formulation) [359,[START_REF] Basu | Perfectly matched layers for time-harmonic elastodynamics of unbounded domains: theory and finite-element implementation[END_REF][START_REF] Basu | Perfectly matched layers for transient elastodynamics of unbounded domains[END_REF][START_REF] Basu | Explicit finite element perfectly matched layer for transient three-dimensional elastic waves[END_REF][START_REF] Bindel | Elastic PMLs for resonator anchor loss simulation[END_REF]161,303,186,199,200,201,[START_REF] Benito | A GFDM with PML for seismic wave equations in heterogeneous media[END_REF][START_REF] Benito | Wave propagation in soils problems using the Generalized Finite Difference Method[END_REF]293,[START_REF] Benito | Implementations with generalized finite differences of the displacements and velocity-stress formulations of seismic wave propagation problem[END_REF]114,113,[START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF][START_REF] Brun | Hybrid Asynchronous Perfectly Matched Layer for seismic wave propagation in unbounded domains[END_REF]360]. Later, split-field CFS-PML (resp. M-PML) formulations have been proposed for the linear elastic wave equations written as a first-order system in velocity and stress [124, 122, 236] (resp. [242, 243, 244, 344, 126] and as a secondorder system in displacement [269,268,270]). Lastly, unsplit-field CFS-PML (resp. M-PML) formulations have been developed based on either the first-order mixed (velocity-stress) formulation [START_REF] Appelö | A new absorbing layer for elastic waves[END_REF][START_REF] Drossaert | A nonsplit complex frequency-shifted PML based on recursive integration for FDTD modeling of elastic waves[END_REF]100,192,234,237,275,350,349,103,263] (with an unsplit-field H-PML, i.e. combined CFS-M-PML, formulation proposed in [352,214] and an unsplit-field combined CFS-MZT-PML formulation developed in [328]) or the second-order (displacement-only or mixed displacement-stress) formulation (with potentially additional auxiliary fields) [216, 240, 104, 108, 106, 335] (resp. [114], with an unsplit-field H-PML formulation developed in [215]) of the linear elastodynamic equations. A brief overview of the relations between different split-field and unsplit-field (CFS-)PML formulations and their classification for the linear elastic wave equations formulated as a first-order mixed (velocity-stress) system can be found in [198]. All the aforementioned PML formulations require either the use of convolution products (and related recursive convolution update techniques) or specific (specialized or recursive) time integration schemes with sometimes a large number of additional fields and auxiliary differential equations as well as additional assumptions on the solution fields to derive the PML formulations in the time domain. By and large, the PML formulations differ not only by the use of split or unsplit primary fields (and auxiliary fields if should be the case) and the definition of the complex coordinate stretching function, but also by the choice of numerical approximation methods for space and time discretizations. The interested reader can refer to [199,200,201,114,335] for an overview of the various PML formulations and implementations (based on FD, FE or SE methods) developed for time-domain elastodynamic problems. As already pointed out in [193,192,292], most existing PML formulations based on the linear elastic wave equations written as a first-order system in velocity and stress have been developed for FD methods and cannot be straightforwardly used in classical numerical approximation methods that are based on the linear elastic wave equations written as a second-order system in displacement, such as most FE, SE, DG methods and some FD methods. Conversely, PML formulations based on the second-order (displacement-based) linear elastodynamic equations can be more readily implemented in such numerical schemes [192,216,186,199,200,201,104,335,292,[START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF]360,[START_REF] Benito | Implementations with generalized finite differences of the displacements and velocity-stress formulations of seismic wave propagation problem[END_REF] and are more robust (i.e. have much better discrete stability properties) than the ones based on first-order (velocity-stress) linear elastodynamic equations [104,195,[START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF][START_REF] Benito | Implementations with generalized finite differences of the displacements and velocity-stress formulations of seismic wave propagation problem[END_REF].

Introduction of an efficient PML formulation in linear elastodynamics

The main objective of this paper is to derive both time-harmonic (frequencydomain) and time-dependent (time-domain) PML formulations as well as their numerical implementations for simulating the propagation of elastic waves in both two-and three-dimensional unbounded media. For this purpose, we develop an unsplit-field PML formulation for second-order (displacement-based) linear elastodynamic equations, which is well-suited for numerical implementation using standard FE, SE or DG methods on unstructured meshes that are well adapted to handle computational domains with arbitrarily complex geometries and curved boundaries and allow for a natural application of boundary conditions, as already mentioned in [199,292]. Despite the zero-frequency singularity, we adopt a standard stretching function with careful parameterization, since it leads to a simple and straightforward implementation of the resulting PML formulation and exhibits better performance in the presence of low-frequency propagating waves (and in the absence of strong evanescent waves) when compared to CFS-PML formulations, as clearly explained in [124,240,200,335] and numerically observed in [START_REF] Correia | Performance of regular PML, CFS-PML, and second-order PML for waveguide problems[END_REF] for instance. In order to relax or at least alleviate the temporal complexity arising from the PML-transformed equations and to derive an unsplit-field non-convolutional PML formulation in the time domain for both two-and three-dimensional transient elastodynamic analyses, we first introduce an auxiliary (symmetric) tensor-valued strain field (defined as the difference between the complex stretched tensor-valued strain field and the real classical one) treated as an additional unknown variable (only in the PML region) and governed by a second-order ordinary differential equation in time derived from the strain-displacement kinematic relation. We then propose to construct a space weak formulation of both the equilibrium equations and the strain-displacement kinematic relation, resulting in a mixed (displacement-strain) unsplit-field PML formulation that is second-order in time and can be discretized in space using a standard Galerkin finite element method. We finally employ a classical implicit second-order time integration scheme com-monly used in computational structural dynamics, that is a Newmark time scheme coupled with a finite difference (Crank-Nicolson) time scheme, in order to solve the time-dependent PML equations while preserving second-order accuracy. The resulting mixed PML formulation can be easily implemented using classical numerical tools incorporated in existing (displacement-based) finite element codes. Such codes originally developed for bounded computational domains can then be easily modified to accommodate both time-harmonic (frequency-domain) and time-dependent (time-domain) PML formulations for the propagation of elastic waves in unbounded domains. In addition, it requires no field splitting, no nonlinear solvers (such as Newton-type iterative solvers), no mass lumping and no critical time step size (such as the stability criterion required in explicit time integration schemes) as well as no complicated or costly convolution operations in time and no high-order spatial and temporal derivatives, so that the derivation, resolution and implementation of the proposed PML formulation are made easier compared to most other existing PML formulations. As in [186,199,200,201,114,113,[START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF], the proposed PML formulation preserves the second-order form of the original (unstretched) time-domain linear elastic wave equations while using a very small number of additional fields and auxiliary differential equations as compared to most existing PML formulations based on either first-or second-order equations for time-domain elastic wave propagation in both two-and three-dimensional unbounded domains. For further computational savings, we finally consider a hybrid approach, originally developed for two-dimensional elastodynamic problems in [201] and later extended to three-dimensional elastodynamic problems in [114], that couples a mixed displacement-strain unsplit-field PML formulation in the PML region with a purely (non-mixed) displacement-based formulation (i.e. a standard displacement-only formulation) in the physical domain.

Outline of the paper

The paper is organized as follows. In Section 2, we introduce the linear elastodynamic formulation in both time and frequency domains. In Section 3, we revisit the time-harmonic (frequency-domain) and time-dependent (time-domain) PML formulations based on complex coordinate stretching by introducing an auxiliary tensor-valued strain field and using a space weak formulation of the resulting PML equations. We also discuss the proposed hybrid approach, whereby a mixed displacement-strain unsplit-field PML formulation in the PML region is coupled with a standard (non-mixed) displacement-only formulation in the physical domain of interest, and provide details on its finite element implementation in both the frequency domain and the time domain. In Section 4, we present numerical experiments carried out on two-dimensional linear elastodynamic problems in single-and multi-layer isotropic homogeneous media to numerically simulate the propagation of elastic waves in semi-infinite (unbounded) domains. The capability of the proposed time-dependent PML formulation to efficiently absorb outgoing propagative waves is shown by computing relevant local and global quantities of interest such as the components of the displacement field at some selected receiving points (called receivers) and the kinetic, internal and total energies stored in the physical domain of interest. The performances of the PML formulation in terms of accuracy (absorption capability) are compared to the ones of a classical absorbing layer (CAL) formulation. Finally in Section 5, we conclude with further discussions and give possible perspectives.

2 Problem statement 2.1 Time-domain formulation of linear elastodynamics Within the framework of linear elasticity theory, we consider a two-or threedimensional linear elastodynamic problem. The material is assumed to be linear and elastic, characterized by the fourth-order Hooke's elasticity tensor C and the mass density ρ. The medium is subjected to given external forces represented by a body force field f . In the time domain, the original problem, initially defined on an unbounded domain of the d-dimensional Euclidean physical space R d (d = 2, 3 being the spatial dimension), consists in finding the vector-valued displacement field u(x, t) and its associated tensor-valued Cauchy stress field σ(x, t) satisfying the second-order partial differential equation (derived from linear momentum conservation and hereinafter referred to as the equilibrium equations)

ρ ü -div(σ) = f , (1a)
for spatial position vector x ∈ R d and time t > 0, with given initial conditions at time t = 0 u(x, 0) = u 0 (x) and u(x, 0

) = v 0 (x), (1b)
where div denotes the divergence operator of a tensor-valued field with respect to x, while the dot and the double (or superposed) dots over a variable denote respectively the first-and second-order partial derivatives with respect to time t, i.e. u(x, t) = ∂u ∂t (x, t) and ü(x, t) = ∂ 2 u ∂t 2 (x, t). Body force field f (•, t) is a given function from R d into R d for all t 0. Both fields u and σ are linked through the linear stress-strain constitutive relation (known as the generalized Hooke's law)

σ = C : ε(u), (1c)
in which the colon symbol : denotes the twice contracted tensor product, and where ε(u(x, t)) is the classical linearized tensor-valued strain field associated to displacement field u(x, t), that is

ε(u) = 1 2 ∇ u + (∇ u) T , (1d)
where the superscript T denotes the transpose operator and ∇ (resp. ∇) denotes the gradient operator of a vector-valued (resp. scalar-valued) field with respect to x. Recall that both second-order tensor-valued stress field σ and strain field ε are symmetric, i.e. σ ij = σ ji and ε ij = ε ji . Also, the fourth-order elasticity tensor C satisfies the classical major and minor symmetry properties . Besides, the mass density ρ is a strictly positive-valued field.

C ijkl = C klij = C jikl = C ijlk (
In the case of isotropic homogeneous linear elastic materials, the constitutive relation between σ and ε(u) restricts to

σ = λ tr(ε(u))I + 2µε(u) = λ div(u)I + µ(∇ u + (∇ u) T ), (2)
where tr denotes the trace operator, div denotes the divergence operator of a vector-valued field with respect to x, I is the second-order symmetric identity (or unit) tensor such that I ij = δ ij , with the Kronecker delta δ ij = 1 if i = j and δ ij = 0 if i = j, λ > 0 and µ > 0 are the Lamé's coefficients (that must be strictly positive according to the ellipticity properties of elasticity tensor C [327, 109]) defined by

λ =        Eν (1 + ν)(1 -2ν)
in 2D under plane strain assumption and in 3D

Eν 1 -ν 2
in 2D under plane stress assumption and

µ = E 2(1 + ν)
(also called shear modulus)

with E and ν respectively the Young's modulus and the Poisson's ratio of the isotropic linear elastic material. The fourth-order elasticity tensor C is then constituted of 2 algebraically independent coefficients (λ, µ) (or equivalently, (E, ν))

independent of x such that C ijkl = λδ ij δ kl + µ(δ ik δ jl + δ il δ jk ), leading to σ ij = C ijkl ε kl = λε kk δ ij + 2µε ij (
with summation over indices k and l only). In such an isotropic homogeneous linear elastic medium, the time-dependent elastic wave system (1) of homogeneous equations (under no external body forces) set on an unbounded domain of R d admits bulk (or body) plane wave solutions, namely longitudinal pressure (or compressional) waves (often referred to as P -waves) with characteristic velocity c p = (λ + 2µ)/ρ and whose polarization vector is parallel to the wavefront propagation direction, and transverse shear waves (often referred to as S-waves) with characteristic velocity c s = µ/ρ < √ 2c p and whose polarization vector is perpendicular to the wavefront propagation direction. It also admits surface-or interface-guided wave solutions, such as Rayleigh waves (propagating along a free surface) and Stoneley waves (propagating along the interface between two semi-infinite elastic media), whose amplitudes decay exponentially in the direction normal to and away from the surface (or the interface if should be the case) [START_REF] Basu | Perfectly matched layers for time-harmonic elastodynamics of unbounded domains: theory and finite-element implementation[END_REF]. In the more general case of anisotropic and/or heterogeneous linear elastic media, the fourth-order elasticity tensor C is constituted of 21 algebraically independent coefficients possibly dependent on x. The time-dependent elastic wave system (1) admits further plane wave solutions with different characteristic velocities (depending on the material symmetry properties and on the wavefront propagation direction), such as quasi-longitudinal (or quasi-compressional) waves (often referred to as qP -waves) and quasi-shear (or quasi-transverse) waves (often referred to as qS-waves) for orthotropic elastic media [START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF]242], as well as other wave types in the presence of heterogeneities and interfaces.

For computational purposes, we introduce a bounded domain Ω BD ⊂ R d with sufficiently smooth boundary ∂Ω BD , where the subscript BD refers to "Bounded Domain", which corresponds to the physical domain of interest over which the numerical computations are performed. We assume that the external body force field f and the initial displacement and velocity fields u 0 and v 0 have bounded supports in Ω BD and therefore remain confined within the bounded (or finite) computational domain. The problem of interest then consists in finding the solution fields u and σ satisfying the time-dependent elastic wave system (1) only in the physical bounded domain Ω BD (and not in the entire unbounded domain of R d ).

Frequency-domain formulation of linear elastodynamics

In the frequency domain, the time-harmonic elastic wave system is obtained by taking the Fourier transform in time of ( 1) and reads

-ρω 2 û -div( σ) = f , (3a) σ = C : ε( û), (3b)
ε( û) = 1 2 ∇ û + (∇ û) T , (3c)
where ω > 0 is the angular frequency and the hat (or caret) over a variable denotes the Fourier transform in time, i.e. û(

•, ω) = R u(•, t)e -iωt dt, in which i
denotes the unit imaginary number satisfying i 2 = -1, and where the spatial and frequency dependencies of variables are dropped for the sake of readability and notational brevity. The time-harmonic dependence of the displacement, stress and strain fields is chosen as the multiplicative factor e iωt by convention. Furthermore, all the vector-valued and tensor-valued fields are expressed throughout the paper in a Cartesian coordinate system with canonical basis {e i } 1 i d .

Perfectly matched layers

We now revisit the PML method that aims at constructing artificial layers (of finite thickness) surrounding the physical bounded domain Ω BD and designed in such a way that the incoming waves (entering into the layers) are not reflected at the interface with the physical domain and decay exponentially into them.

In the following, we implicitly assume that all the outgoing waves (leaving the physical domain of interest) correspond to purely forward propagating waves in the direction normal to the interface, i.e. there is no backward propagating waves or evanescent waves.

We consider a bounded Cartesian domain

Ω BD = d × i=1 [-i , i ] ⊂ R d surrounded
by an artificial absorbing layer (the PML region)

Ω PML = d × i=1 [-i -h i , i + h i ] \ Ω BD with a constant thickness h i in each x i -direction. The PML region Ω PML is divided into d overlapping subregions {Ω i } d i=1 such that Ω i = {x ∈ Ω PML ; |x i | > i }
and the outgoing propagative waves are attenuated in the x i -direction within the PML subregion Ω i (see Figure 1). We denote by Γ = ∂Ω BD ∩ ∂Ω PML (resp.

O • Ω BD 1 2 Ω PML Ω 1 Ω 2 h 1 h 2 Γ Γext x 1 x 2
Fig. 1: Physical bounded domain Ω BD surrounded by an artificial PML region

Ω PML = d i=1 Ω i divided into d overlapping subregions Ω i with interface Γ = ∂Ω BD ∩ ∂Ω PML and exterior boundary Γ ext = ∂Ω PML \ Γ Γ i = ∂Ω BD ∩ ∂Ω i )
the interface between the physical domain Ω BD and the PML region Ω PML (resp. the PML subregion Ω i ), and by n = n i e i the unit normal vector to Γ pointing outward Ω BD . Furthermore, we define the entire domain Ω = Ω BD ∪ Ω PML as the disjoint union between Ω BD and Ω PML , and we denote by Γ ext the part of its boundary ∂Ω that coincides with the exterior boundary of Ω PML , that is Γ ext = ∂Ω PML \ (Γ ∪ Γ free ), where Γ free stands for the possible free surfaces in the case where a semi-infinite medium with free-surface boundary conditions is considered.

Complex coordinate stretching

In the frequency domain, the PML method consists in analytically stretching one or more real spatial coordinate(s) of the physical wave equations into the complex plane C, the ones corresponding to the direction(s) along which the outwardly propagating waves have to be attenuated in the absorbing layers. In the following, for any complex number z ∈ C, we denote by Re(z) and Im(z) the real and imaginary parts of z, respectively.

Following the complex coordinate stretching approach [69,284,285,[START_REF] Chew | Perfectly matched layers for elastodynamics: A new absorbing boundary condition[END_REF][START_REF] Chew | A 3D perfectly matched medium from modified maxwell's equations with stretched coordinates[END_REF], we introduce a complex change of spatial variable S :

x ∈ Ω ⊂ R d → S(x) = x ∈ C = S(Ω) ⊂ C d defining a complex stretching (analytic continuation) xi = S i (x i ) ∈ C i of the real spatial coordinate x i along a curve C i in the complex plane C parame- terized by xi = S i (x i ) = x i 0 s i (ξ) dξ, (4)
where the tilde over a variable or an operator denotes its stretched version in the complex space, and s i (x i ) is a complex coordinate stretching function defining the absorption profile along the x i -direction, which is a nowhere zero, continuous, complex-valued function such that

Im (s i (x i )) < 0 in Ω i , s i (x i ) = 1 elsewhere. (5)
Condition [START_REF] Abarbanel | Non-Linear PML Equations for Time Dependent Electromagnetics in Three Dimensions[END_REF]

ensures that xi = S i (x i ) ∈ C \ R in the PML subregion Ω i , whereas xi = S i (x i ) = x i ∈ R
elsewhere and in particular in the physical domain Ω BD , so that the outgoing propagative waves are perfectly transmitted at the interface Γ and attenuated (or damped) in the x i -direction inside the PML subregion Ω i .

As an inverse Fourier transform in time may be required to come back into the time domain by inverting the frequency-domain PML equations, the stretching function s i (x i ) is usually chosen frequency-dependent [69,284,285,[START_REF] Collino | Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF]193,[START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF][START_REF] Bermúdez | An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems[END_REF]303,192,[START_REF] Drossaert | A nonsplit complex frequency-shifted PML based on recursive integration for FDTD modeling of elastic waves[END_REF][START_REF] Bermúdez | Perfectly Matched Layers for Time-Harmonic Second Order Elliptic Problems[END_REF]145,195,247,[START_REF] Bécache | Perfectly matched layers in negative index metamaterials and plasmas[END_REF]360,248] and defined as

s i (x i ) = 1 + d i (x i ) iω , (6)
where d i (x i ) is the so-called damping (or absorption or attenuation) function that is a continuous and monotonically increasing (i.e. non-decreasing) positive-valued function of x i such that

d i (x i ) > 0 in Ω i , d i (x i ) = 0 elsewhere. ( 7
)
The dependence of the stretching function s i (x i ) on the factor iω in ( 6) allows for an easy application of the inverse Fourier transform in time to the time-harmonic (frequency-domain) PML equations, and thus a simple derivation of the timedependent (time-domain) PML equations. It should be noted that the particular form (6) of stretching function s i (x i ) exhibits a zero-frequency singularity (zerofrequency pole). Therefore, the derivation of the PML formulation considered in this work naturally excludes the static case corresponding to ω = 0, that means it does not allow static solutions (displacement and stress fields associated to a given static loading) to be recovered [200]. In theory, any choice of damping function d i (x i ) satisfying [START_REF] Abarbanel | Unsplit Variables Perfectly Matched Layers for the Shallow Water Equations with Coriolis Forces[END_REF] ensure both continuity at the interface Γ i between Ω i and Ω BD and exponential amplitude decay along the x i -direction in the PML subregion Ω i for the outgoing propagative waves. Since the damping function d i (x i ) does not depend on angular frequency ω, it further allows obtaining an exponential amplitude decay of the propagating waves in the PML subregion Ω i that is independent of ω for non-dispersive PML media (see Section 3.2 for a discussion). In practice, typical choices for d i (x i ) in the PML subregion Ω i are low-degree polynomial functions [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Chew | Perfectly matched layers for elastodynamics: A new absorbing boundary condition[END_REF]75,[START_REF] Collino | The Perfectly Matched Layer in Curvilinear Coordinates[END_REF][START_REF] Collino | Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF]359,233,329,193,123,[START_REF] Basu | Perfectly matched layers for time-harmonic elastodynamics of unbounded domains: theory and finite-element implementation[END_REF][START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF]301,306,[START_REF] Cohen | Mixed Spectral Finite Elements for the Linear Elasticity System in Unbounded Domains[END_REF]124,161,[START_REF] Appelö | A new absorbing layer for elastic waves[END_REF]303,[START_REF] Bérenger | Perfectly matched layer (PML) for computational electromagnetics[END_REF]358,192,[START_REF] Drossaert | A nonsplit complex frequency-shifted PML based on recursive integration for FDTD modeling of elastic waves[END_REF]235,236,234,275,350,216,240,344,186,199,200,201,335,269,268,270,[START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF][START_REF] Brun | Hybrid Asynchronous Perfectly Matched Layer for seismic wave propagation in unbounded domains[END_REF]263,360] or singular (non-integrable) hyperbolic or shifted hyperbolic functions [START_REF] Bermúdez | An exact bounded PML for the Helmholtz equation[END_REF][START_REF] Bermúdez | Validation of acoustic models for time-harmonic dissipative scattering problems[END_REF][START_REF] Bermúdez | An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems[END_REF][START_REF] Bermúdez | An Exact Bounded Perfectly Matched Layer for Time-Harmonic Scattering Problems[END_REF][START_REF] Bermúdez | Perfectly Matched Layers for Time-Harmonic Second Order Elliptic Problems[END_REF]185,246,249,247,[START_REF] Cimpeanu | A parameter-free perfectly matched layer formulation for the finite-element-based solution of the Helmholtz equation[END_REF]248] according to the underlying problem. In the numerical experiments shown in Section 4, we adopt the most widely used damping profile (due to its straightforward implementation) by choosing a polynomial function d i (x i ) with a power law dependence on the distance |x i |i to ensure a progressive (gradual)

and smoothly increasing damping of the outgoing propagative waves across the PML subregion Ω i , that is

d i (x i ) = d max i |x i | -i h i p in Ω i , (8)
where p is the polynomial degree (that is a nonnegative integer) and d max i > 0 is the damping coefficient (that is a strictly positive parameter) having the same unit as the angular frequency ω and corresponding to the maximum value of damping function d i (x i ). Both p and d max i are user-chosen scalar parameters allowing the shape (sharpness) and the intensity (strength) of the imposed damping profile of the propagating waves to be controlled within the PML subregion Ω i [200,201]. From a computational viewpoint, choosing a smoothly varying damping profile is of particular interest for the attenuation of outgoing propagative waves within the PML region to be adequately resolved by the spatial discretization in order to minimize the spurious reflections. Furthermore, for any classical numerical method that implicitly assumes continuity of the spatial derivatives of the solution field, choosing a constant or linear function d i (x i ) in the PML subregion Ω i (i.e. taking p = 0 or p = 1 in ( 8)) can generate spurious numerical reflections at the interface Γ i . By contrast, for given values of the damping coefficient d max i and the PML thickness h i , choosing a quadratic or higher-degree polynomial function d i (x i ) in the PML subregion Ω i (i.e. taking p 2 in (8)) avoids those numerical issues and also minimizes the spurious reflections produced at the exterior boundary of the PML subregion Ω i [START_REF] Collino | The Perfectly Matched Layer in Curvilinear Coordinates[END_REF]301,[START_REF] Bindel | Elastic PMLs for resonator anchor loss simulation[END_REF]303]. In the numerical examples presented in Section 4, we choose a quadratic form for d i (x i ) (with a fixed polynomial degree p = 2) so as to enforce a smooth and gradual spatial attenuation of the outgoing propagating waves within the PML subregion Ω i .

In the general case of heterogeneous materials, we further assume that the mass density ρ and the components C ijkl of the fourth-order elasticity tensor C depend analytically on x in the d-dimensional complex manifold C = S(Ω), so that the solution fields u and σ of problem (1) (or equivalently, their Fourier transforms in time û and σ solution of problem ( 3)) are also analytic with respect to x and can be evaluated along the complex path C [259, [START_REF] Bermúdez | Perfectly Matched Layers for Time-Harmonic Second Order Elliptic Problems[END_REF]182].

Reflectionless and absorption properties

The wave propagation properties induced by the complex coordinate stretching in the PML region Ω PML can be merely analyzed by means of a plane wave analysis [START_REF] Collino | Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF][START_REF] Basu | Perfectly matched layers for time-harmonic elastodynamics of unbounded domains: theory and finite-element implementation[END_REF][START_REF] Basu | Perfectly matched layers for transient elastodynamics of unbounded domains[END_REF][START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF]123,124,[START_REF] Appelö | A new absorbing layer for elastic waves[END_REF][START_REF] Bérenger | Perfectly matched layer (PML) for computational electromagnetics[END_REF][START_REF] Bermúdez | Perfectly Matched Layers for Time-Harmonic Second Order Elliptic Problems[END_REF]242,243,244,104,[START_REF] Bécache | Perfectly matched layers in negative index metamaterials and plasmas[END_REF]. We consider harmonic propagative plane waves corresponding to particular solutions of the time-dependent elastic wave system (1) (or equivalently, of its time-harmonic counterpart (3)) of homogeneous equations (under no external body force field, i.e.

with f = 0) of the form u(x, t) = Ae i(ωt-k•x) , ( 9
)
where A is the polarization vector, k = k i e i = kp is the wave vector, with k = k the wavenumber, and p = p i e i is the unit vector denoting the wavefront propagation direction, with p i = k i /k the direction cosinus of k along the x idirection, and k • x = k i x i denotes the usual Euclidean inner product between k and x in R d and k

2 = k•k = k 2 i the associated Euclidean norm of k in R d
, where the Einstein notation convention is used for summation over repeated indices. The complex-valued change of coordinates S i : x i → S i (x i ) = xi introduced in (4) allows a propagative plane wave of the form [START_REF] Ahland | Accurate mesh truncation for Schrödinger equations by a perfectly matched layer absorber: Application to the calculation of optical spectra[END_REF] to be transformed into an attenuated wave

ũ(x, t) = u( x, t) = Ae i(ωt-k•x) = Ae i(ωt-k•x) e -α(x)•p = u(x, t) e -α(x) (10)
exponentially decreasing as |x i | increases in each PML subregion Ω i , where α(x) = α i (x i )e i is the amplitude attenuation vector whose components α i (x i ) are such that

α i (x i ) =    k ω x i 0 d i (ξ) dξ > 0 in Ω i , 0
elsewhere, [START_REF] Antoine | A Review of Transparent and Artificial Boundary Conditions Techniques for Linear and Nonlinear Schrödinger Equations[END_REF] and α(x) = α(x) • p = α i (x i )p i (with summation over index i) is the so-called amplitude attenuation factor (or decay rate). The phase of the outgoing propagative plane waves is maintained as they penetrate the PML subregion Ω i , whereas the amplitude is modulated by an exponential factor e -α(x) that introduces an exponential decay of the outgoing propagative plane waves along the x i -direction once they enter the PML subregion Ω i . As the attenuation factor α(x) depends on the direction of propagation p (characterizing the angle of incidence) and on the phase velocity v p = ω/k but not directly on the angular frequency ω itself, it turns out to be frequency-independent (i.e. independent of ω) for non-dispersive PML media, thus leading to a uniform exponential decay of the outgoing propagative plane waves with respect to the angular frequency ω. Also note that, using ( 11), the stretching function s i (x i ) and damping function d i (x i ) introduced in ( 6) can be respectively expressed as

s i (x i ) = 1 + α i (x i ) ik and d i (x i ) = α i (x i ) ω k , (12)
where the prime over a variable denotes the first derivative with respect to spatial coordinate x i , i.e. α i (x i ) = dα i dx i (x i ). Also, the complex stretched coordinate xi = S i (x i ) defined in ( 4) can be rewritten as [69,285,[START_REF] Collino | Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF]193,[START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF]192] xi = S i (x i ) = x i + 1 iω

x i 0 d i (ξ) dξ = x i + α i (x i ) ik . (13)
Except for near-grazing incident propagative waves (for which k • n ≈ 0), all the outgoing propagative waves (for which k • n > 0) are exponentially damped in the PML region Ω PML since the amplitude attenuation factor α(x) is positive and increases as |x i | increases in each PML subregion Ω i , while all the ingoing ones (for which k • n < 0) are exponentially amplified in Ω PML since the amplitude attenuation factor α(x) is negative and increases as

|x i | decreases in Ω i .
As all outgoing propagative waves are expected to be damped in the PML region Ω PML , an ad hoc (Dirichlet, Neumann, characteristic or mixed) boundary condition can be applied on the exterior boundary Γ ext of the PML region. Let us mention that, provided that the PML model is carefully parameterized, the choice of boundary conditions specified on the exterior boundary of Ω PML does not seem to have a significant impact on the PML absorption performance [START_REF] Collino | The Perfectly Matched Layer in Curvilinear Coordinates[END_REF]264]. In the numerical experiments presented in Section 4 (as in most practical applications), for the sake of convenience and simplicity, we choose to impose homogeneous Dirichlet boundary conditions on the exterior boundary of Ω PML , that is for all t 0, u(•, t) = 0 on Γ ext .

The absorptive and attenuative properties in the PML region Ω PML can be analyzed by studying the reflection of an incident plane wave at the fixed exterior boundary Γ ext of Ω PML . A propagative plane wave traveling outward from the physical domain Ω BD is absorbed into the PML region Ω PML without any reflections at the interface Γ between Ω BD and Ω PML . At the continuous level, the PML region Ω PML is thus said perfectly matched to the physical domain Ω BD . A propagative plane wave entering and traversing the PML subregion Ω i , being reflected at the exterior boundary Γ ext of Ω PML and again propagating back through Ω i is attenuated by an exponential factor [75, [START_REF] Collino | The Perfectly Matched Layer in Curvilinear Coordinates[END_REF][START_REF] Collino | Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF]303,185]

R i = exp(-2α i ( i + h i )p i ) = exp -2 p i v p i +h i i d i (ξ) dξ (15)
before re-entering the physical domain Ω BD , as it travels twice through the PML subregion Ω i (of finite thickness h i ). The theoretical reflection coefficient R i a priori depends on the choice of three user-defined and user-tunable PML parameters, namely the position i of the interface Γ i , the thickness h i of the PML subregion Ω i and the damping function d i (x i ), but it is also influenced by the phase velocity v p = ω/k and the direction of propagation p i of the plane waves along the x i -direction. For the quadratic damping function d i (x i ) previously defined in (8) (with a polynomial degree p = 2), the reflection coefficient R i can be rewritten as [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Bérenger | Three-Dimensional Perfectly Matched Layer for the Absorption of Electromagnetic Waves[END_REF][START_REF] Collino | Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF]303,240]

R i = exp - 2d max i h i p i 3v p . ( 16
)
As a result, the reflected-wave amplitude R i defined in [START_REF] Appelö | A General Perfectly Matched Layer Model for Hyperbolic-Parabolic Systems[END_REF] decreases as the PML thickness h i , the wavefront propagation direction p i (along the x i -direction) or the damping coefficient d max i increases, whereas it increases as the phase velocity v p increases. Note that it does not depend of the size i of the physical domain Ω BD . Also, for non-dispersive media, the reflection coefficient R i is independent of the angular frequency ω since the phase velocity v p = ω/k is equal to the wave velocity [124]. In practical computations, increasing the PML thickness h i allows the outgoing propagative waves to be more attenuated but may induce substantial computational costs. Besides, for a given PML thickness h i , the damping coefficient d max i cannot be chosen arbitrarily large, since large values of d max i would require fine spatial discretizations to well reproduce and approximate the fast attenuation of outgoing waves within the PML region with sufficiently good numerical accuracy, thus leading to high computational costs. Indeed, disregarding any numerical integration or rounding error, the discretization error arising from the numerical approximation method can cause spurious reflections at the interface Γ between Ω BD and Ω PML [START_REF] Bindel | Elastic PMLs for resonator anchor loss simulation[END_REF]. Also, it becomes predominant compared to the truncation error related to the spurious reflections at the exterior boundary of the PML subregion Ω i (of finite thickness h i ) as the damping coefficient d max i increases [START_REF] Bermúdez | An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems[END_REF]303]. Therefore, for fixed values of PML thickness h i and discretization parameters, choosing an optimal damping coefficient d max i is essential for producing minimal reflections from the exterior boundary of the PML region and thus maximizing the attenuation of outgoing waves while maintaining a desired accuracy level. Such an optimal value results from a trade-off between absorption strength (requiring large enough values of d max i ) and numerical accuracy (requiring small enough values of d max i ) [75,[START_REF] Collino | The Perfectly Matched Layer in Curvilinear Coordinates[END_REF][START_REF] Collino | Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF][START_REF] Asvadurov | On Optimal Finite-Difference Approximation of PML[END_REF]161,292]. It a priori depends on the PML thickness h i , the wavefront propagation direction p i , the phase velocity v p and the discretization error introduced by the numerical scheme. In practice, the damping coefficient d max i is usually chosen such that the reflection coefficient R i at normal incidence (i.e. for p i = 1) is equal to an arbitrarily small value R 0 < 1 (typically 10 -3 to 10 -8 ) and defined as [75,[START_REF] Collino | The Perfectly Matched Layer in Curvilinear Coordinates[END_REF][START_REF] Collino | Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF][START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF][START_REF] Bécache | On the long-time behavior of unsplit perfectly matched layers[END_REF]193,124,229,303,192,[START_REF] Drossaert | A nonsplit complex frequency-shifted PML based on recursive integration for FDTD modeling of elastic waves[END_REF]235,234,275,216,186,199,200,201,114,240,344,335,269,268,270,[START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF]263,360]

d max i = 3v p 2h i log 1 R 0 , (17)
or more generally as [123,124,242,243,244]

d max i = A i v p h i , (18)
with A i a given dimensionless constant adapted to the problem being considered. Note that there are only few guidelines reported in the literature and based on either simple heuristics (ad hoc rules of thumb), a priori error estimates or costly general nonlinear optimization algorithms (see [START_REF] Collino | The Perfectly Matched Layer in Curvilinear Coordinates[END_REF][START_REF] Datta | Optimization of uniaxial perfectly matched layer parameters for finite difference time domain simulation and application to coupled microstrip lines with multiple bend discontinuities[END_REF][START_REF] Asvadurov | On Optimal Finite-Difference Approximation of PML[END_REF]301,302,[START_REF] Bindel | Elastic PMLs for resonator anchor loss simulation[END_REF]161,303,245,218,257,292] for instance) to properly select appropriate values for the PML model and discretization parameters. Hence, as already pointed out in [199,200,201], owing to a critical lack of rigorous methodology and simple process for providing proper guidance on the choice of damping functions and on the selection of optimal PML parameters, the derivation of an optimal damping function d i independently of the underlying problem and the discretization method is still an open question.

A possible alternative would consist in using the singular (non-integrable) damping functions proposed in [START_REF] Bermúdez | An exact bounded PML for the Helmholtz equation[END_REF][START_REF] Bermúdez | Validation of acoustic models for time-harmonic dissipative scattering problems[END_REF][START_REF] Bermúdez | An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems[END_REF][START_REF] Bermúdez | An Exact Bounded Perfectly Matched Layer for Time-Harmonic Scattering Problems[END_REF][START_REF] Bermúdez | Perfectly Matched Layers for Time-Harmonic Second Order Elliptic Problems[END_REF]] that exhibit a power-type singularity at the exterior boundary Γ ext of the PML region and allows recovering a reflection coefficient R i = 0 for any angle of incidence and for any angular frequency, since [START_REF] Bermúdez | An exact bounded PML for the Helmholtz equation[END_REF][START_REF] Bermúdez | An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems[END_REF][START_REF] Bermúdez | An Exact Bounded Perfectly Matched Layer for Time-Harmonic Scattering Problems[END_REF]. Numerical computations [START_REF] Bermúdez | An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems[END_REF][START_REF] Bermúdez | An Exact Bounded Perfectly Matched Layer for Time-Harmonic Scattering Problems[END_REF]276,246,249,247,248] have shown that such (shifted) hyperbolic damping functions do not require any tuning of the PML parameters (whatever the numerical method employed) and provide good results that may outperform (or at least be comparable with) those obtained with optimized polynomial damping functions.

d i is a non-integrable function such that i +h i i d i (ξ) dξ = +∞

Frequency-domain PML formulation

The derivation of the time-harmonic (frequency-domain) PML equations classically results from the application of the complex coordinate stretching (4) to the original time-harmonic governing equations (3).

Classical frequency-domain PML formulation based on complex coordinate stretching

From the complex coordinate transformation (4) together with ( 6) and making use of the fundamental theorem of calculus, we have the following relation and derivative rule between the complex stretched coordinate system and the real physical coordinate system

∂ xi ∂x i = s i (x i ) = 1 + d i (x i ) iω and ∂(•) ∂ xi = 1 s i (x i ) ∂(•) ∂x i = iω iω + d i (x i ) ∂(•) ∂x i ,
where the partial derivatives ∂(•) ∂ xi with respect to the complex coordinates xi are interpreted in the usual Cauchy-Riemann sense [289,[START_REF] Bermúdez | Perfectly Matched Layers for Time-Harmonic Second Order Elliptic Problems[END_REF], i.e.

∂(•) ∂ xi = ∂ Re(•) ∂ Re(x i ) - i ∂ Im(•) ∂ Im(x i )
. The jacobian matrix J = ∇ S = ∂ xi ∂x j e i ⊗ e j of the complex coordinate transformation ( 4) is then a complex-valued diagonal matrix that is everywhere non-singular and defined by

J = s i (x i )e i ⊗ e i , (19)
in which the symbol ⊗ denotes the usual tensor product. In the following, the functional (spatial) dependence of stretching function s i and damping function d i will be often omitted for the sake of readability and notational brevity. Therewith, the analytic continuation of the time-harmonic elastic wave system (3) of homogeneous equations (in the absence of body force field, i.e. with f = 0) in the PML region Ω PML is basically obtained by replacing the partial derivatives ∂(•) ∂x i with respect to x i by the partial derivatives

∂(•) ∂ xi with respect to xi in (3)
and reads [START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF][START_REF] Basu | Perfectly matched layers for time-harmonic elastodynamics of unbounded domains: theory and finite-element implementation[END_REF][START_REF] Basu | Perfectly matched layers for transient elastodynamics of unbounded domains[END_REF][START_REF] Bindel | Elastic PMLs for resonator anchor loss simulation[END_REF]200,201,335]

-ρω 2 û -div( σ) = 0, (20a) σ = C : ε( û), (20b) ε( û) = 1 2 ∇ û + ( ∇ û) T , (20c)
which is similar to the original time-harmonic elastic wave system (3) defined over the physical domain Ω BD , except that the external body force field f has been set to zero (since f is non-vanishing only within Ω BD ) while the partial differential operators ∇, ∇, div and ε in the real physical x-coordinate system have been replaced with the frequency-dependent complex stretched operators ∇, ∇ and div and ε in the complex stretched x-coordinate system, respectively defined by

∇û = J -T ∇û = 1 s i ∂ û ∂x i e i , ∇ û = (∇ û)J -1 = 1 s j ∂ ûi ∂x j e i ⊗ e j , div( σ) = ∇ • σ = (J -T ∇) • σ = 1 s j ∂ σij ∂x j e i , ε( û) = 1 2 (∇ û)J -1 + J -T (∇ û) T = 1 2 1 s j ∂ ûi ∂x j + 1 s i ∂ ûj ∂x i e i ⊗ e j .
Note that the complex coordinate stretching (4) allows the symmetry properties of the second-order tensor-valued strain field ε to be preserved, so that the secondorder complex stretched tensor-valued strain field ε is symmetric, i.e. εij = εji . Besides, although by construction the complex coordinate stretching (4) implicitly excludes the static case for which ω = 0 [199,200,201,114,[START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF], we consider continuous extensions of the complex stretched strain tensor ε( û(x, ω)) and stress tensor σ(x, ω) so that in the static case, i.e. for ω = 0, we have for all x ∈ Ω PML , ε( û(x, 0)) = 0 and σ(x, 0) = 0.

Note that condition ( 21) is (often implicitly) stated as an assumption in other PML formulations previously derived in the literature, such as in [START_REF] Basu | Perfectly matched layers for transient elastodynamics of unbounded domains[END_REF][START_REF] Basu | Explicit finite element perfectly matched layer for transient three-dimensional elastic waves[END_REF]186,199,200,201,114,[START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF][START_REF] Brun | Hybrid Asynchronous Perfectly Matched Layer for seismic wave propagation in unbounded domains[END_REF]. Multiplying (20a) by det(J ), where det denotes the determinant operator, and noticing that det(J )/s j = d i=1,i =j s i is independent of x j , the complex stretched (PML-transformed) time-harmonic elastic wave system [START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF] can be rewritten as a system similar to (3) with classical real partial differential operators but for an anisotropic medium made up with a specific frequency-dependent complex anisotropic heterogeneous material. It reads [359,161,216,240,239,335,297]

-ρω 2 û -div( σ) = 0, (22a) σ = C : ∇ û, (22b)
where ρ = ρ det(J ) = ρ d i=1 s i is a frequency-dependent complex scalar-valued field, σ = σ det(J )J -1 is a second-order non-symmetric complex tensor-valued field such that σij = σij det(J )/s j , and C is a fourth-order frequency-dependent complex anisotropic tensor-valued field such that C ijkl = C ijkl det(J )/(s j s l ) with only the major symmetry property preserved, but not the minor ones (due to the non-symmetry of both fields σ and ∇ û)2 . In view of the time-harmonic PML formulation [START_REF] Assimaki | Time-domain forward and inverse modeling of lossy soils with frequency-independent Q for near-surface applications[END_REF], the PML region can be interpreted as an artificial anisotropic heterogeneous absorbing medium with specific complex material properties [291,130,129,356,357,353,363,181,310,311,313,315,2].

Remark 1 Alternatively, left multiplying (or pre-multiplying) (22b) by J -T , the constitutive relation (22b) can be rewritten as [124] σ = C : ∇ û, where σ = J -T σ = det(J )J -T σJ -1 is a second-order symmetric complex tensor-valued field such that σ ij = σij /s i = σij det(J )/(s i s j ), and C is a fourth-order frequency-dependent complex anisotropic tensor-valued field such that C ijkl = C ijkl /s i = C ijkl det(J )/(s i s j s l ) with only the left minor symmetry property preserved (resulting from the symmetry of σ ), but neither the right minor one, nor the major one (due to the non-symmetry of ∇ û).

Remark 2 Another anisotropic medium re-interpretation of the time-harmonic PML formulation [START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF] is proposed in [START_REF] Bindel | Elastic PMLs for resonator anchor loss simulation[END_REF]. It allows for a simple implementation of the time-harmonic (frequency-domain) PML equations within the context of finite element methods. Rewriting the kinematic relation (3c) as ε( û) = L : ∇ û, where L is the fourth-order symmetric identity (or unit) tensor such that L ijpq = (δ ip δ jq + δ iq δ jp )/2, the complex stretched kinematic relation (20c) then becomes ε( û) = L : ∇ û, where L is a fourth-order complex tensor-valued field such that L ijpq = (δ ip (J -1 ) qj + (J -1 ) qi δ jp )/2 = (δ ip δ jq /s j + δ iq δ jp /s i )/2. Note that both tensors L and L satisfy only the left minor symmetry property (due to the symmetry of strain tensors ε( û) and ε( û)), but not the major and right minor ones (due to the non-symmetry of ∇ û). The time-harmonic elastic wave system [START_REF] Assimaki | Time-domain forward and inverse modeling of lossy soils with frequency-independent Q for near-surface applications[END_REF] can then be recast in the same form as the original one (3) and reads [START_REF] Bindel | Elastic PMLs for resonator anchor loss simulation[END_REF] -ρω 2 ûdiv( σ) = 0, (23a)

σ = C : ε( û), ( 23b
)
where C is a fourth-order frequency-dependent complex anisotropic tensor-valued field such that C pqrs = L ijpq C ijkl L klrs det(J ) which inherits the major and minor symmetry properties of C (thanks to the left minor symmetry of L).

Mixed frequency-domain PML formulation with additional auxiliary field

For the stretching function s i chosen as in ( 6), the jacobian matrix J defined in [START_REF] Appelö | Application of a perfectly matched layer to the nonlinear wave equation[END_REF] can be rewritten as

J = I + 1 iω D, (24)
where D = d i (x i )e i ⊗ e i is a diagonal real matrix with positive diagonal terms.

The stress field σ can then be expressed in terms of the original one σ as

σ =        σ I + 1 iω D Σ 2 for d = 2, σ I + 1 iω D Σ 3 - 1 ω 2 D Π 3 for d = 3, (25)
where D Σ 2 , D Σ 3 and D Π 3 are diagonal real matrices expressed in the canonical basis

{e i } 1 i d of R d such that D Σ 2 = d 2 0 0 d 1 =   2 j=1,j =i d j   e i ⊗ e i , D Σ 3 =   d 2 + d 3 0 0 0 d 1 + d 3 0 0 0 d 1 + d 2   =   3 j=1,j =i d j   e i ⊗ e i , D Π 3 =   d 2 d 3 0 0 0 d 1 d 3 0 0 0 d 1 d 2   =   3 j=1,j =i d j   e i ⊗ e i .
Note that, since D is a diagonal matrix, D 20) (or equivalently [START_REF] Assimaki | Time-domain forward and inverse modeling of lossy soils with frequency-independent Q for near-surface applications[END_REF]) could be recast as a single vector equation involving only a single (displacement) field, the resulting system of time-dependent PML equations would involve convolution products or would require specific (specialized or recursive) time integration schemes due to the temporal complexity. To further maintain the secondorder (in time) character of the original linear elastodynamic equations ( 1) while avoiding convolution operations in time, thereby greatly facilitating time integration, we now introduce an auxiliary tensor-valued strain field e(x, t) defined by e(x, t) = ε(u(x, t))ε(u(x, t)) for any time t 0 and spatial position vector x ∈ Ω PML . Its Fourier transform in time ê is then such that

ε( û) = ε( û) + ê. ( 26
)
Introducing ( 25) into (22a), substituting ( 26) into ( 20b) and ( 20c), left multiplying (or pre-multiplying) (20c) by iωJ T and right multiplying (or postmultiplying) it by iωJ with [START_REF] Baffet | Long-time stable highorder absorbing boundary conditions for elastodynamics[END_REF], rearranging and grouping similar terms, we obtain the following time-harmonic PML formulation that consists in finding the time-harmonic vector-valued displacement field û(x, ω), tensor-valued stress field σ(x, ω) and auxiliary tensor-valued strain field ê(x, ω) satisfying

ρ -ω 2 + iω tr(D) + det(D) û = div σ + 1 iω σD Σ 2 (27a)
in the two-dimensional case (i.e. for d = 2), and

ρ -ω 2 + iω tr(D) + tr(D Π 3 ) + 1 iω det(D) û = div σ + 1 iω σD Σ 3 - 1 ω 2 σD Π 3 ( 27b
-ω 2 ê + iω D T ê + êD + D T êD + iωε D ( û) + D T ε( û)D = 0, (27d)
for any spatial position vector x ∈ Ω PML and angular frequency ω > 0, where ε D ( û) is a second-order symmetric tensor-valued field defined by

ε D ( û) = 1 2 (∇ û)D + D T (∇ û) T = 1 2 d j ∂ ûi ∂x j + d i ∂ ûj ∂x i e i ⊗ e j .
Multiplying (27b) by iω, the PML equilibrium equations (27b) in the threedimensional case (i.e. for d = 3) can be recast as

ρ -iω 3 -ω 2 tr(D) + iω tr(D Π 3 ) + det(D) û = div iω σ + σD Σ 3 + 1 iω σD Π 3 . ( 28
)
Remark 3 Alternatively, multiplying (20c) by det(J ) similarly to [201,114], using ( 26) and ( 24), the time-harmonic kinematic relation (27d) can be replaced with

-ω 2 + iω tr(D) + det(D) ê = -iωε D ( û) -det(D)ε( û) (29a)
in the two-dimensional case (i.e. for d = 2), and

-ω 2 + iω tr(D) + tr(D Π 3 ) + 1 iω det(D) ê = -iωε D ( û)-ε D 3 ( û)- 1 iω det(D)ε( û)
(29b) in the three-dimensional case (i.e. for d = 3), where ε D 3 ( û) is a second-order symmetric tensor-valued field defined by

ε D 3 ( û) = 1 2 (∇ û)DD Σ 3 + (DD Σ 3 ) T (∇ û) T = 1 2   d j   3 k=1,k =j d k   ∂ ûi ∂x j + d i   3 k=1,k =i d k   ∂ ûj ∂x i   e i ⊗ e j .
Since both time-harmonic equilibrium equations and kinematic relation are stretched in the same manner (both multiplied by det(J )), the left-hand sides of the PML equilibrium equations ( 27a)-( 27b) and the PML kinematic relations ( 29) have a similar form, thus leading to a fully-symmetric mixed displacementstrain unsplit-field PML formulation unlike the non-symmetric one ( 27) considered here. Similar mixed unsplit-field PML formulations have already been proposed in the literature, such as the fully-symmetric ones developed in [201,114] or the nonsymmetric ones proposed in [200,114] for instance, but they retain the stress field (or, more precisely, the time-history of stress field) as unknown variable (instead of the auxiliary strain field considered here) in addition to the displacement field.

Time-domain PML formulation

The time-harmonic PML equations can then be easily transformed into the corresponding time-dependent PML equations for transient analysis due to their simple frequency-dependence resulting from the specific choice of stretching function [START_REF] Abarbanel | Long-Time Performance of Unsplit PMLs with Explicit Second Order Schemes[END_REF].

For the subsequent derivations, we introduce the auxiliary fields U (x, t), Σ(x, t) and Π(x, t), corresponding to the time-integrals (or time-histories) of displacement field u, stress field σ and stress time-history field Σ, respectively, defined for any spatial position vector x ∈ Ω PML and time t 0 by

U (x, t) = t 0 u(x, τ ) dτ, Σ(x, t) = t 0 σ(x, τ ) dτ and Π(x, t) = t 0 Σ(x, τ ) dτ = t 0 τ 0 σ(x, θ) dθ dτ.
For any spatial position vector x ∈ Ω PML , we make use of the Fourier transform in time Σ(x, ω) of Σ(x, t) defined by

Σ(x, ω) = 1 iω σ(x, ω) + πδ 0 (ω) σ(x, ω), (30)
where δ 0 (ω) is the Dirac distribution on R at the origin (i.e. at point ω = 0), in order to derive the following inverse Fourier transform in time

F -1 1 iω σ(x, ω) = Σ(x, t) - 1 2 σ(x, 0) = t 0 σ(x, τ ) dτ - 1 2 R + σ(x, τ ) dτ, (31)
where F -1 denotes the inverse Fourier transform operator. Then, using relation [START_REF] Assi | A perfectly matched layer formulation for modeling transient wave propagation in an unbounded fluid-solid medium[END_REF] at angular frequency ω = 0, that is σ(x, 0) = 0, the inverse Fourier transform in time [START_REF] Bécache | Perfectly Matched Layers for the Convected Helmholtz Equation[END_REF] reduces to

F -1 1 iω σ(x, ω) = Σ(x, t). ( 32
)
By applying the inverse Fourier transform in time to the time-harmonic PML equations [START_REF] Basu | Explicit finite element perfectly matched layer for transient three-dimensional elastic waves[END_REF] in which ( 27b) is replaced with [START_REF] Basu | Perfectly matched layers for time-harmonic elastodynamics of unbounded domains: theory and finite-element implementation[END_REF] and with the aid of ( 32), we derive the corresponding time-dependent PML equations. The time-dependent PML formulation then consists in finding the time-dependent vector-valued displacement field u(x, t), tensor-valued stress field σ(x, t) and auxiliary tensor-valued strain field e(x, t) satisfying

ρ ( ü + tr(D) u + det(D)u) = div σ + ΣD Σ 2 (33a)
in the two-dimensional case (i.e. for d = 2),

ρ ... u + tr(D) ü + tr(D Π 3 ) u + det(D)u = div σ + σD Σ 3 + ΣD Π 3 (33b)
in the three-dimensional case (i.e. for d = 3),

σ = C : ε(u) + C : e, (33c) and ë
+ D T ė + ėD + D T eD + ε D ( u) + D T ε(u)D = 0, (33d)
for any spatial position vector x ∈ Ω PML and time t > 0, where the space and time dependencies of variables are implicit and the triple dots over a variable denotes the third-order partial derivative with respect to time t, i.e. ... u(x, t) = ∂ 3 u ∂t 3 (x, t). The set of time-dependent PML equations ( 33) is complemented with the boundary conditions on the exterior boundary Γ ext given by ( 14), the usual continuity conditions on the displacement and normal stress fields at the interface Γ (which are not recalled here for the sake of brevity) and the following homogeneous initial conditions at time t = 0 given for any x ∈ Ω PML by u(x, 0) = u(x, 0) = 0 and e(x, 0) = ė(x, 0) = 0, (33e)

with ü(x, 0) = 0 for d = 3. ( 33f
)
Note that the PML equilibrium equations (33a) derived in the two-dimensional case are second-order in time, while the PML equilibrium equations (33b) derived in the three-dimensional case are third-order in time. As proposed in [114], in order to preserve a fully second-order in time PML formulation in the three-dimensional case, one may analytically integrate (33b) in time over [0 , t] for any time t > 0, using the initial conditions (33e) and (33f). We then obtain the following PML equilibrium equations in the three-dimensional case (i.e. for d = 3)

ρ ü + tr(D) u + tr(D Π 3 )u + det(D)U = div σ + ΣD Σ 3 + ΠD Π 3 . ( 34
)
The time-dependent PML formulation is then second-order in time for both displacement field u and auxiliary strain field e in both two-and three-dimensional cases.

Remark 4 Alternatively, taking the inverse Fourier transform in time of (27b) in the three-dimensional case, as it was done in [START_REF] Basu | Explicit finite element perfectly matched layer for transient three-dimensional elastic waves[END_REF] for instance, requires the use of the Fourier transform in time Û (x, ω) of U (x, t) and the one Π(x, ω) of Π(x, t), respectively defined for any spatial position vector x ∈ Ω PML by

Û (x, ω) = 1 iω û(x, ω) + πδ 0 (ω) û(x, ω), (35a)
Π(x, ω) = - 1 ω 2 σ(x, ω) + iπδ 0 (ω) σ(x, ω), (35b)
where δ 0 (ω) is the first derivative of δ 0 (ω) with respect to ω, in order to derive the following inverse Fourier transforms in time

F -1 1 iω û(x, ω) = U (x, t) - 1 2 û(x, 0) = t 0 u(x, τ ) dτ - 1 2 R + u(x, τ ) dτ, (36a)
F -1 - 1 ω 2 σ(x, ω) = Π(x, t) + i 2 ∂ σ ∂ω (x, 0) + i t σ(x, 0) = t 0 τ 0 σ(x, θ) dθ dτ - 1 2 R + (t -τ )σ(x, τ ) dτ. (36b)
By applying the inverse Fourier transform in time to the time-harmonic PML equilibrium equations (27b) with the aid of ( 32) and ( 36), we may recover [START_REF] Bécache | High-order absorbing boundary conditions for anisotropic and convective wave equations[END_REF] subject to silent initial conditions (33e) at time t = 0 under the additional assumptions û(x, 0) = 0 and ∂ σ ∂ω (x, 0) = 0 for all x ∈ Ω PML .

Space weak formulation

Similarly to [201,114], we adopt a hybrid approach for solving the overall mixed problem defined over the entire domain Ω = Ω BD ∪ Ω PML by considering a mixed displacement-strain formulation in the PML region Ω PML coupled with a purely (non-mixed) displacement-based formulation (i.e. a standard displacement-only formulation) in the physical domain Ω BD for computational efficiency.

Introducing the constitutive relation (27c) into (27a) and (27b), the weak formulation of the system of time-harmonic PML equations ( 27) consists in finding û(x, ω) over the entire domain Ω = Ω BD ∪ Ω PML and ê(x, ω) over the PML region Ω PML such that for all vector-valued test functions δu(x) in an appropriate admissible vector space such that δu(x) = 0 for x ∈ ∂Ω \ Γ free , and for all secondorder symmetric tensor-valued test functions δe(x) in an appropriate admissible vector space,

-ω 2 m( û, δu) + iωc( û, δu) + k( û, δu) + k 2 ( û, δu) + 1 iω g( û, δu) + k1 (ê, δu) + 1 iω g1 (ê, δu) = (δu) for d = 2, ( 37a
) -ω 2 m( û, δu) + iωc( û, δu) + k( û, δu) + k 3 ( û, δu) + 1 iω (g( û, δu) + g 1 ( û, δu)) - 1 ω 2 h( û, δu) + k1 (ê, δu) + 1 iω g1 (ê, δu) - 1 ω 2 h1 (ê, δu) = (δu) for d = 3, ( 37b
) -ω 2 m(ê, δe) + iωc(ê, δe) + k(ê, δe) + iωc 1 ( û, δe) + k 1 ( û, δe) = 0, (37c)
for any angular frequency ω > 0, where m, c, k, k1 , k 2 , k 3 , g, g1 , h, h1 and g 1 = k 2 are sesquilinear forms and is an antilinear (i.e. conjugate-linear) form, respectively defined by

m( û, δu) = Ω ρ û • δu dΩ, c( û, δu) = Ω PML ρ tr(D) û • δu dΩ, k( û, δu) = Ω (C : ε( û)) : ε(δu) dΩ, k1 (ê, δu) = Ω PML (C : ê) : ε(δu) dΩ, k 2 ( û, δu) = Ω PML ρ det(D) û • δu dΩ, k 3 ( û, δu) = Ω PML ρ tr(D Π 3 ) û • δu dΩ, g( û, δu) = Ω PML (C : ε( û)) : ε Σ d (δu) dΩ, g1 (ê, δu) = Ω PML (C : ê) : ε Σ d (δu) dΩ, h( û, δu) = Ω PML (C : ε( û)) : ε Π 3 (δu) dΩ, h1 (ê, δu) = Ω PML (C : ê) : ε Π 3 (δu) dΩ, g 1 ( û, δu) = k 2 ( û, δu), (δu) = Ω BD f • δu dΩ,
where the overline denotes the complex conjugate transpose operator and

ε Σ d (δu) = 1 2 (∇ δu)D Σ d + D Σ d T (∇ δu) T = tr(D)ε(δu) -ε D (δu), ε Π 3 (δu) = 1 2 (∇ δu)D Π 3 + D Π 3 T (∇ δu) T = tr(D Π 3 )ε(δu) -ε D 3 (δu),
and where m, c, k, c 1 and k 1 are sesquilinear forms respectively defined by m(ê, δe) =

Ω PML ê : δe dΩ, c(ê, δe) = Ω PML D T ê + êD : δe dΩ, k(ê, δe) = Ω PML D T êD : δe dΩ, c 1 ( û, δe) = Ω PML ε D ( û) : δe dΩ, k 1 ( û, δe) = Ω PML D T ε( û)D : δe dΩ.
Since the auxiliary strain field ê needs to be introduced only in the PML region Ω PML and the diagonal real matrices D (gathering the damping functions d i ), D Σ 2 , D Σ 3 and D Π 3 are non-zero only in Ω PML , all the space-integrals involved in the sesquilinear forms m, c, c 1 , c, k 1 , k 2 , k 3 , k, k1 , g, g 1 , g1 , h and h1 of (37) restrict to Ω PML . Besides, since the physical domain Ω BD contains the support of the complex vector f associated to external body forces, the space-integral involved in the antilinear form restricts to Ω BD .

In the following, we will consider vectorial representations of second-order symmetric tensor-valued fields e and ê. Also, for the sake of readability and notational convenience, we will use the same notation for both second-order symmetric tensors (or tensor-valued fields) of dimension d(d+1)/2 and their vectorial representations.

Hybrid finite element approximation

For the mixed finite element implementation of space weak formulation [START_REF] Bécache | Fictitious domains, mixed finite elements and perfectly matched layers for 2-D elastic wave propagation[END_REF], both displacement field u and auxiliary strain field e are retained as unknowns and treated as independent variables that can be approximated separately. We then use a classical (displacement-based) finite element method (FEM) [176,362] by introducing finite-dimensional approximation spaces spanned by piecewise polynomial basis functions for the spatial discretization of both displacement field u (resp. its Fourier transform in time û) over Ω = Ω BD ∪ Ω PML and auxiliary strain field e (resp. its Fourier transform in time ê) over Ω PML in the time domain (resp. the frequency domain). For any time t 0, the finite element approximations of displacement field u(•, t) and auxiliary strain field e(•, t) can then be identified with the real vectors U(t) and E(t) gathering the nodal values of u(•, t) in the entire domain Ω = Ω BD ∪ Ω PML and the ones of e(•, t) in the PML region Ω PML , respectively. Similarly, for any angular frequency ω > 0, the finite element approximations of Fourier transforms in time û(•, ω) and ê(•, ω) can be identified with the complex vectors Û(ω) and Ê(ω) gathering the nodal values of û(•, ω) in the entire domain Ω = Ω BD ∪ Ω PML and the ones of ê(•, ω) in the PML region Ω PML , respectively.

Frequency-domain implementation

For the finite element discretization of sesquilinear forms m, c, k, k 2 , k 3 , k1 , g, g 1 , g1 , h, h1 and antilinear form involved in (37a) and (37b), we introduce the finite element (frequency-independent) real matrices M, C, K, K 2 , K 3 , K 1 , G, G 1 , G 1 , H and H 1 , and the finite element (frequency-dependent) complex vector F(ω). Similarly, for the finite element discretization of sesquilinear forms m, c, k, c 1 and k 1 involved in (37c), we introduce the finite element (frequency-independent) real matrices M, C, K, C 1 and K 1 .

In an algebraic setting, the finite element matrix system resulting from the weak formulation of time-harmonic PML equations [START_REF] Bécache | Fictitious domains, mixed finite elements and perfectly matched layers for 2-D elastic wave propagation[END_REF] reads

-ω 2 M Û + iωC Û + (K + K 2 ) Û + K 1 Ê + 1 iω G Û + G 1 Ê = F for d = 2, ( 38a
) -ω 2 M Û + iωC Û + (K + K 3 ) Û + K 1 Ê + 1 iω (G + G 1 ) Û + G 1 Ê - 1 ω 2 H Û + H 1 Ê = F for d = 3, (38b) -ω 2 M Ê + iω C Ê + C 1 Û + K Ê + K 1 Û = 0. (38c)
Then, it comes down to solving the following system of linear algebraic equations for any angular frequency ω > 0:

-ω 2 M Û + iωC Û + K Û + 1 iω G Û = F for d = 2, ( 39a
) -ω 2 M Û + iωC Û + K Û + 1 iω G Û - 1 ω 2 H Û = F for d = 3, ( 39b
)
where the (frequency-independent) real matrices M, C, K, G and H, and the (frequency-dependent) complex vectors F(ω) and Û(ω) are respectively defined by

M = M 0 0 M , C = C 0 C 1 C , K = K + K d K 1 K 1 K , G = G G 1 0 0 for d = 2, G = G + G 1 G 1 0 0 for d = 3, H = H H 1 0 0 only for d = 3, F = F 0 and Û = Û Ê .
Note that all the block-system matrices are frequency-independent square real matrices, but only the mass matrix M is symmetric and has a block-diagonal structure, whereas the other matrices K, C, G and H are non-symmetric. Besides, the block-system (39) has a sparse structure at two levels, coming from the sparsity pattern inherent to each finite element matrix involved in [START_REF] Bécache | On the analysis of perfectly matched layers for a class of dispersive media and application to negative index metamaterials[END_REF], but also from the block-sparsity pattern of matrices

M, C, C 1 , C, K 1 , K 2 , K 3 , K, K 1 , G, G 1 , G 1 ,
H and H 1 that are assembled only from their element-level constituent matrices defined over the PML region Ω PML . For practical implementation in the frequency domain, the FE matrix system (39) can be rewritten as

K dyn (ω) Û = F,
where K dyn (ω) is the complex dynamic (frequency-dependent) stiffness matrix defined by

K dyn (ω) =        K -ω 2 M + i ωC - 1 ω G for d = 2, K -ω 2 M - 1 ω 2 H + i ωC - 1 ω G for d = 3.
Note that K dyn (ω) is a non-symmetric complex matrix due to the non-symmetry of real matrices K, C, G and H. For practical implementation in the time domain, the FE matrix system (39b) in the three-dimensional case can be recast as

-iω 3 M Û -ω 2 C Û + iωK Û + G Û + 1 iω H Û = iω F for d = 3. ( 40
)

Time-domain implementation

The time-harmonic FE matrix system (39a)-( 40) can then be easily inverted back into the time domain for transient analysis due to its simple frequency-dependence. By applying the inverse Fourier transform in time to the time-harmonic FE matrix system (39a)-( 40), we obtain the following second-order and third-order integrodifferential equations (IDEs) in time

M Ü + C U + KU + GI = F for d = 2, (41a) M ... U + C Ü + K U + GU + HI = Ḟ for d = 3, (41b)
for time t > 0, in which the time-dependent real vectors U(t), U(t), Ü(t), ... U(t), I(t), F(t) and Ḟ(t) are respectively defined by

U = U E , U = U Ė , Ü = Ü Ë , ... U = ... U ... E , I = t 0 U dτ, F = F 0 and Ḟ = Ḟ 0 .
The IDEs in time ( 41) are supplemented with the following initial conditions at time t = 0:

U(0) = U 0 = U 0 0 and U(0) = U0 = V 0 0 , (42a) with
Ü(0) = Ü0 = A 0 0 only for d = 3, (42b)
where U 0 and V 0 are the initial displacement and velocity FE vectors corresponding to the spatial discretization of the initial displacement and velocity fields u 0 and v 0 over Ω, and A 0 is the initial acceleration FE vector corresponding to the spatial discretization of the initial acceleration field ü(•, 0) over Ω and such that MA 0 = F 0 -CV 0 -KU 0 with F 0 = F(0) the initial FE vector of nodal forces corresponding to the spatial discretization of the initial body force field f (•, 0) at time t = 0. Note that, since f , u 0 and v 0 have bounded supports in Ω BD (see the last paragraph of Section 2.1), F, U 0 and V 0 have non-zero components only for the degrees of freedom associated to the spatial discretization of f , u 0 and v 0 over the physical domain Ω BD .

In order to recover a second-order IDE in time in the three-dimensional case, one may analytically integrate (41b) in time over [0 , t] for any time t > 0, using the initial conditions [START_REF] Bécache | Perfectly matched layers in negative index metamaterials and plasmas[END_REF]. We then obtain the following second-order IDE in time in the three-dimensional case

M Ü + C U + KU + GI + HJ = F for d = 3, (43)
for time t > 0, in which the time-dependent real vector J(t) is defined by

J = t 0 I dτ = t 0 τ 0 U dθ dτ.
The semi-discrete systems (41a) in the two-dimensional case and ( 43) in the threedimensional case are then both second-order in time, thus allowing the use of standard time integration schemes classically employed in computational structural dynamics. Also, the FE solution U of the time-domain PML formulation obviously satisfies the causality principle.

Remark 5

In the numerical experiments presented in Section 4, the proposed PML formulation is compared to a CAL formulation with ad hoc dissipative material properties introduced in the PML region Ω PML but without a characteristic impedance matching the one of the physical domain Ω BD . In the frequency domain, the time-harmonic CAL formulation consists in solving the following system of linear algebraic equations at any angular frequency ω > 0:

K dyn (ω) Û = F, (44)
where K dyn (ω) is the dynamic (frequency-dependent) stiffness matrix defined by

K dyn (ω) = K -ω 2 M + iωC.
Note that K dyn (ω) is a symmetric complex matrix due to the symmetry of real matrices K, M and C. In the time domain, the CAL formulation leads to the following classical second-order ordinary differential equation (ODE) in time:

M Ü + C U + KU = F (45)
for time t > 0, with the initial conditions U(0) = U 0 and U(0) = V 0 at time t = 0.

Time sampling scheme

We use both a Newmark time scheme [256] and a finite difference time scheme based on the trapezoidal rule for solving the equations in the time domain for a given time sampling of U, U, Ü, I and J. Another possibility would consist in using an extended Newmark time scheme, as it was done in [199,[START_REF] Assimaki | Time-domain forward and inverse modeling of lossy soils with frequency-independent Q for near-surface applications[END_REF]114] for instance.

For the sake of simplicity, we consider a uniform time step ∆t over a given time interval of interest I = [0 , T ], with T = N ∆t the final time and N the number of time steps.

3.7.1 Newmark scheme for the displacement vector U, velocity vector U and acceleration vector Ü

The Newmark scheme may be either explicit or implicit, unstable or (un)conditionally stable (depending on the values of the two parameters γ and β) and it is based on the following recursive formulae expressed in acceleration format as

U n = U n-1 + ∆t Un-1 + ∆t 2 2 (1 -2β) Ün-1 + 2β Ün , (46a) Un = Un-1 + ∆t (1 -γ) Ün-1 + γ Ün , ( 46b
)
where U n , Un and Ün denote respectively the displacement, velocity and acceleration FE vectors evaluated at time t n = n∆t, with ∆t = t nt n-1 the discrete time step. Alternatively, assuming that β = 0, the Newmark time scheme [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF] can also be rewritten in displacement format as

Un = γ 1 U n -U n-1 + (1 + γ 2 ) Un-1 + γ 3 Ün-1 , (47a) Ün = β 1 U n -U n-1 + β 2 Un-1 + (1 + β 3 ) Ün-1 , (47b)
where

β 1 = 1 β∆t 2 , β 2 = - 1 β∆t , β 3 = - 1 2β , γ 1 = γ β∆t , γ 2 = - γ β , γ 3 = 1 - γ 2β ∆t.
The two parameters γ and β allow the degree of implicitness and the stability of the time scheme to be controlled [176]. Recall that the Newmark time scheme [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF] (or ( 47)) is second-order accurate whatever the values of γ and β. Also, when γ < 1/2, it is unstable, whereas when γ 1/2 and 2β γ, it is unconditionally stable for non-active and non-dispersive media [START_REF] Chilton | The discrete origin of FETD-Newmark late time instability, and a correction scheme[END_REF]240], and when γ 1/2 and 2β < γ, it becomes conditionally stable, which means the choice of time step ∆t is limited and bounded by the Courant-Friedrichs-Lewy (CFL) stability condition [START_REF] Courant | Über die partiellen Differenzengleichungen der mathematischen Physik[END_REF][START_REF] Courant | On the Partial Difference Equations of Mathematical Physics[END_REF] for non-dispersive media [240]. Besides, when β = 0, the time integration scheme ( 46) is explicit whatever the value of γ. In particular, when γ = 1/2 and β = 0, one recovers the explicit conditionally stable central-difference scheme. Such an explicit time integration scheme has been recently combined with matrix lumping to avoid the computational burden associated with matrix inversion and applied in (un)split-field (CFS-)PML formulations to perform two-and three-dimensional transient elastodynamic analyses in unbounded (or very large) domains [229,236,[START_REF] Basu | Explicit finite element perfectly matched layer for transient three-dimensional elastic waves[END_REF]240,335,341]. In the numerical experiments shown in Section 4, we use the implicit unconditionally stable constant average acceleration scheme by taking γ = 1/2 and β = 1/4. Such an implicit time integration scheme requires the solution of a linear matrix system at every time step (which may be computationally expensive in very large three-dimensional transient elastodynamic analyses), but it does not suffer from the CFL stability condition in return, which means the choice of time step ∆t is not limited by the smallest element size (spatial resolution), which is of particular interest especially for highly irregular (unstructured) and/or fine spatial FE meshes. It has already been employed in unsplit-field (CFS-)PML formulations to perform transient analyses of second-order elastic wave equations in unbounded domains [START_REF] Basu | Perfectly matched layers for transient elastodynamics of unbounded domains[END_REF]186,199,200,201,240,114,[START_REF] Brun | Hybrid Asynchronous Perfectly Matched Layer for seismic wave propagation in unbounded domains[END_REF]360]. Note that the unconditional stability is guaranteed only for non-active and non-dispersive media [START_REF] Chilton | The discrete origin of FETD-Newmark late time instability, and a correction scheme[END_REF], and therefore it cannot be ensured to hold in the PML region Ω PML that is inherently an anisotropic dissipative (absorbing) medium. Nevertheless, several numerical experiments have shown that the proposed time stepping scheme is actually highly stable [240].

3.7.2 Finite difference scheme for the time-integral vectors I and J

The finite difference scheme may be either explicit or implicit, (un)conditionally stable (depending on the value of the parameter α) and it is based on the following recursive formulae

I n = I n-1 + ∆t (1 -α)U n-1 + αU n , (48a)
J n = J n-1 + ∆t (1 -α)I n-1 + αI n , ( 48b
)
where I n and J n denote respectively the time-integral I of displacement FE vector U and the time-integral J of FE vector I evaluated at time t n = n∆t. Note that the time scheme ( 48b) is introduced only in the three-dimensional case (i.e. for d = 3), since J is involved in the left-hand side of ( 43) in the three-dimensional case (i.e. for d = 3), whereas it does not appear in the left-hand side of (41a) in the two-dimensional case (i.e. for d = 2). Without loss of generality and for the sake of simplicity, we consider here the same parameter α for both time schemes (48a) and ( 48b). The parameter α allows the degree of implicitness and the stability of the time scheme to be controlled. Recall that when α 1/2, the finite difference scheme ( 48) is unconditionally stable, whereas when α < 1/2, it becomes conditionally stable. Besides, when α = 0 (resp. α = 1), one recovers the explicit conditionally stable forward Euler time scheme (resp. the implicit unconditionally stable backward Euler time scheme), which is first-order accurate, whereas when α = 1/2, it leads to the implicit unconditionally stable Crank-Nicolson time scheme, which is second-order accurate. In the numerical experiments shown in Section 4, we use the Crank-Nicolson time scheme by taking α = 1/2 to preserve the intrinsic second-order accuracy and the unconditional stability of the Newmark time scheme [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF] (or [START_REF] Berenger | Perfectly matched layer for the FDTD solution of wave-structure interaction problems[END_REF]) with γ = 1/2 and β = 1/4 (for non-active and non-dispersive media).

Time sampling scheme for the PML formulation

Introducing the discretized version of the second-order IDEs in time (41a) and ( 43) at time t n and using the above formulae ( 47) and [START_REF] Bérenger | Three-Dimensional Perfectly Matched Layer for the Absorption of Electromagnetic Waves[END_REF], we derive the time stepping scheme for the FE discretized wave equations (41a) and ( 43) expressed in displacement format that reads as follows for n 1,

K eff U n = F n -M Ün β -C Un γ -GI n α for d = 2, ( 49a
)
K eff U n = F n -M Ün β -C Un γ -(G + α∆tH) I n α -HJ n α for d = 3, ( 49b
)
where Un

γ = -γ 1 U n-1 + (1 + γ 2 ) Un-1 + γ 3 Ün-1 , Ün β = -β 1 U n-1 + β 2 Un-1 + (1 + β 3 ) Ün-1 , I n α = I n-1 + (1 -α)∆tU n-1 and J n α = J n-1 + (1 -α)∆tI n-1
are the respective predictors of unknown FE vectors Un , Ün , I n and J n at time t n such that Un = Un

γ + γ 1 U n , Ün = Ün β + β 1 U n , I n = I n α + α∆tU n and J n = J n α + α∆tI n ,
F n is the right-hand side FE vector F evaluated at time t n , and K eff is the effective stiffness matrix defined by

K eff = β 1 M + γ 1 C + K + α∆tG for d = 2, β 1 M + γ 1 C + K + α∆t (G + α∆tH) for d = 3.
Note that K eff is a non-symmetric real matrix (due to the non-symmetry of real matrices K, C, G and H) which depends on the two parameters γ and β involved in the Newmark time scheme [START_REF] Berenger | Perfectly matched layer for the FDTD solution of wave-structure interaction problems[END_REF] though γ 1 and β 1 and also on the parameter α involved the finite difference time scheme [START_REF] Bérenger | Three-Dimensional Perfectly Matched Layer for the Absorption of Electromagnetic Waves[END_REF]. For practical purposes, introducing the predictors Un γ , Ün β , I n α and J n α into ( 49) and using again the discretized version of the second-order IDEs in time (41a) and ( 43) at previous time t n-1 , the linear matrix system [START_REF] Bérenger | Improved PML for the FDTD solution of wave-structure interaction problems[END_REF] to be solved at current time t n for n 1 can be recast as

K eff ∆U n = F n -F n-1 -∆tGU n-1 -(β 2 M + γ 2 C) Un-1 -(β 3 M + γ 3 C) Ün-1 for d = 2, ( 50a
)
K eff ∆U n = F n -F n-1 -∆t (G + α∆tH) U n-1 -(β 2 M + γ 2 C) Un-1 -(β 3 M + γ 3 C) Ün-1 -∆tHI n-1 for d = 3, (50b) in which ∆U n = U n -U n-1
is the displacement increment FE vector. The time stepping scheme boils down to computing the displacement FE vector U n = U n-1 + ∆U n solution of the linear matrix system ( 50), then updating the corresponding velocity and acceleration FE vectors Un and Ün as well as I n (only for d = 3) through ( 47a), ( 47b) and (48a) at each time step t n . Note that there is no need to compute and store the FE vectors I n in the two-dimensional case (i.e. for d = 2) and J n in the three-dimensional case (i.e. for d = 3). Even though the computation and storage of FE vector I n at each time t n are required in the three-dimensional case (i.e. for d = 3), it should be noted that I n has non-zero components only for the degrees of freedom associated to the spatial discretization of the time-histories

U (•, t n ) = t n 0 u dτ and E(•, t n ) =
t n 0 e dτ (evaluated at time t n ) of displacement field u and auxiliary strain field e, respectively, over the PML region Ω PML .

Remark 6 Alternatively, considering the acceleration-based formula [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF] of the Newmark time scheme, we derive an equivalent reformulation of the time stepping scheme [START_REF] Bérenger | Improved PML for the FDTD solution of wave-structure interaction problems[END_REF] expressed in acceleration format that reads as follows for n 1,

M eff Ün = F n -C Un γ -(K + α∆tG) U n β -GI n α for d = 2, ( 51a
)
M eff Ün = F n -C Un γ -(K + α∆t (G + α∆tH)) U n β -(G + α∆tH) I n α -HJ n α for d = 3, ( 51b
)
where

U n β = U n-1 + ∆t Un-1 + (1 -2β) ∆t 2 2
Ün-1 and Un γ = Un-1 + (1γ)∆t Ün-1 are the respective predictors of unknown FE vectors U n and Un at time t n such that U n = U n β + β∆t 2 Ün and Un = Un γ + γ∆t Ün , and M eff is the effective mass matrix defined by

M eff = β∆t 2 K eff = M + γ∆tC + β∆t 2 (K + α∆tG) for d = 2, M + γ∆tC + β∆t 2 (K + α∆t (G + α∆tH)) for d = 3.
In the same way as K eff , M eff is a non-symmetric real matrix that depends on both parameters γ and β and also on parameter α. For practical purposes, introducing the predictors U n β , Un γ , I n α and J n α into ( 51) and using the discretized version of the second-order IDEs in time (41a) and ( 43) at previous time t n-1 , the linear matrix system [START_REF] Bérenger | Evanescent waves in PML's: origin of the numerical reflection in wavestructure interaction problems[END_REF] to be solved at current time t n for n 1 can be recast as

M eff ∆ Ün = F n -F n-1 -∆tGU n-1 -∆t (K + α∆tG) Un-1 -∆t C + ∆t 2 (K + α∆tG) Ün-1 (52a)
in the two-dimensional case (i.e. for d = 2), and

M eff ∆ Ün = F n -F n-1 -∆t (G + α∆tH) U n-1 -∆t (K + α∆t (G + α∆tH)) Un-1 -∆t C + ∆t 2 (K + α∆t (G + α∆tH)) Ün-1 -∆tHI n-1 (52b)
in the three-dimensional case (i.e. for d = 3), in which ∆ Ün = Ün -Ün-1 is the acceleration increment FE vector. The time stepping scheme boils down to computing the acceleration FE vector Ün = Ün-1 + ∆ Ün solution of the linear matrix system ( 52), then updating the displacement and velocity FE vectors U n and Un as well as I n (only for d = 3) through ( 46a), ( 46b) and (48a) at each time step t n . Note that once again there is no need to compute and store the FE vectors I n in the two-dimensional case (i.e. for d = 2) and J n in the three-dimensional case (i.e. for d = 3), whereas it is required to compute and store I n at each time t n in the three-dimensional case (i.e. for d = 3).

Remark 7 Following Remark 5, recall that the time-domain PML formulation is compared to a CAL formulation governed by the second-order ODE in time [START_REF] Benito | Wave propagation in soils problems using the Generalized Finite Difference Method[END_REF] in the numerical experiments shown in Section 4. In the time domain, such a CAL formulation comes down to solving the following linear matrix system at each time step t n for n 1:

K eff ∆U n = F n -F n-1 -(β 2 M + γ 2 C) Un-1 -(β 3 M + γ 3 C) Ün-1 , (53)
with the initial conditions given by U 0 and V 0 , and where

∆U n = U n -U n-1
is the displacement increment FE vector and K eff is the effective stiffness matrix defined by

K eff = β 1 M eff = β 1 M + γ 1 C + K,
with M eff = M + γ∆tC + β∆t 2 K the corresponding effective mass matrix.

Numerical results

In order to assess the performances and the robustness of the proposed PML formulation, we present different numerical experiments for simulating the propagation and absorption of elastic waves in two-dimensional semi-infinite (unbounded) single-or multi-layer isotropic homogeneous elastic media (without any viscous damping) subjected to a directional force point-source. We study numerically the accuracy and stability of the proposed time-domain PML formulation governed by the second-order IDE in time (41a) and compare it to the time-domain CAL formulation governed by the second-order ODE in time [START_REF] Benito | Wave propagation in soils problems using the Generalized Finite Difference Method[END_REF]. In all the numerical experiments, an explosive source term f (x, t), acting as a downward vertical force, is applied to a source point S located in the physical domain Ω BD at the position x S and defined as

f (x, t) = r(t)g(x)e 2 ,
where the source spatial function g(x) = δ(xx S ) is the Dirac delta function δ (also known as the impulse function) on R d at the origin (i.e. at point 0) shifted by the source position x S , and the source time function r(t) is a second-order Ricker wavelet 3 (sometimes called a Mexican hat wavelet) given by [START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF]242,243,244,240,[START_REF] Brun | Hybrid Asynchronous Perfectly Matched Layer for seismic wave propagation in unbounded domains[END_REF]341,215,166]

r(t) = 2(πf d ) 2 (t -t d ) 2 -1 exp(-(πf d ) 2 (t -t d ) 2 ), ( 54
)
with a central (or dominant) frequency f d = 1/3 Hz, and a source time delay t d = 1/f d = 3 s. The time-Fourier transform r(ω) of the wavelet function r(t) is written as

r(ω) = -4 √ π ω 2 ω 3 d exp - ω 2 ω 2 d exp(-iωt d )
with a central (or dominant) angular frequency ω d = 2πf d = 2π/3 ≈ 2.0944 rad/s. Figure 2 shows the evolutions of the Ricker wavelet r(t) with respect to time t and of the amplitude |r(ω)| of its time-Fourier transform r(ω) with respect to angular frequency ω = 2πf . For the time integration scheme used for solving the time-domain equations, a constant time step ∆t = 0.025 s is chosen to correctly reproduce the source time variation and adequately captured the numerical solutions for all the numerical tests. The final sampled time of the numerical simulations is T = 20 s, thus requiring N = 800 time steps. Besides, the elastic medium is assumed to be at rest at initial time t = 0, so that zero initial conditions u(x, 0) = u 0 (x) = 0 and u(x, 0) = v 0 (x) = 0 hold for any x in R d . As regards the spatial discretization, all the finite element meshes considered in the numerical experiments are generated with Gmsh [132]. Finally, all material parameters and physical quantities are expressed in terms of the International System of Units (universally abbreviated SI units), such as times in [s], frequencies in [Hz], angular frequencies in [rad.s -1 ], lengths and displacements in [m], velocities in [m.s -1 ], mass densities in [kg.m -3 ], elastic moduli in [Pa], energies in [J], etc., and we further drop units for the sake of readability. 3 The Ricker wavelet r(t) defined in [START_REF] Bérenger | Perfectly matched layer (PML) for computational electromagnetics[END_REF]

is such that r(t) = d 2 p dt 2 (t), where p(t) = 1 2(πf d ) 2 exp(-(πf d ) 2 (t -t d ) 2
) corresponds to a modified second-order derivative of the Gaussian probability density function t → 1 We first consider the well-known Lamb's problem [204], which consists in modeling the propagation of elastic waves in an isotropic homogeneous half-space with freesurface boundary conditions and a directional force point-source applied onto the top free surface of the semi-infinite (unbounded) domain.

σ √ 2π exp - (t -t d ) 2
We aims at numerically simulating the propagation of elastic waves in the halfplane ]-∞ , +∞[×]-∞ , 0] ⊂ R 2 (with a top free surface Γ free located at x 2 = 0) generated by a vertical force point-source S located at the center of the top free surface Γ free of the semi-infinite domain, i.e. at the position x S = (0, 0) of the frame origin (see Figure 3a). The material properties of the isotropic homogeneous elastic medium are characterized by a mass density ρ = 1, a Young's modulus E = 2.5 and a Poisson's ratio ν = 0.25 under the 2D plane strain assumption, so that the Lamé's coefficients are λ = µ = 1, the longitudinal P -wave velocity (resp. wavelength) is c p ≈ 1.7321 (resp. λ p = c p /f d ≈ 5.1962) and the transverse S-wave velocity (resp. wavelength) is c s = 1 (resp. λ s = c s /f d = 3). This academic test case has already been studied in [START_REF] Brun | Hybrid Asynchronous Perfectly Matched Layer for seismic wave propagation in unbounded domains[END_REF]341]. The exact analytical solution involves bulk P -and S-waves as well as Rayleigh waves generated by the point-source.

We introduce a rectangular (two-dimensional) bounded domain Ω BD = [-, ] × [-, 0] with length 2 = 8 and width = 4, surrounded by a PML region Ω PML = {x ∈ R 2 ; < |x 1 | + h and -h x 2 < -} corresponding to three overlapping absorbing layers (on the right, left and bottom sides of the physical domain of interest Ω BD ) with the same finite thickness h = /2 = 2 for each absorbing layer (see Figure 3a). The entire computational domain

Ω = Ω BD ∪ Ω PML = [--h , + h] × [--h , 0]
is then a rectangular bounded domain with length 2( + h) = 12 and width + h = 6 over which the calculations are performed. As already mentioned in Section 3.2, due to the expected absorption properties in the PML region Ω PML , we apply homogeneous Dirichlet boundary conditions u = 0 on the exterior boundary Γ ext of Ω PML , that is at |x 1 | = + h = 6 and at x 2 = -h = -6, and we consider stress-free boundary given by ( 18) with a phase velocity v p = c p (being equal to the longitudinal P -wave speed) and an amplitude A i = 6 corresponding to a small theoretical reflection coefficient R i = 10 -4 in [START_REF] Appelö | Perfectly Matched Layers for Hyperbolic Systems: General Formulation, Well-posedness, and Stability[END_REF]. Note that the PML thickness h = 2 is slightly lower than a half (1/2) of the longitudinal P -wave wavelength λ p ≈ 5.1962 and equal to two thirds (2/3) of the shear S-wave wavelength λ s = 3.

Ω PML h = 2 h S • R • Ω BD = 4 Γ Γext Γ free
For the spatial discretization, the finite element mesh is an unstructured triangulation of the entire domain Ω that is composed of 6-nodes quadratic triangular elements with a characteristic element size ∆x = 0.1 such that ∆x = λ s /30 ≈ λ p /50 (see Figure 3b). It thus contains 38 209 nodes (i.e. 76 418 degrees of freedom) and 18 924 elements. We use quadratic finite elements for the spatial discretization of both displacement field u and auxiliary strain field e, since they usually outperform linear finite elements in terms of accuracy [START_REF] Bindel | Elastic PMLs for resonator anchor loss simulation[END_REF]. Also, standard isoparametric quadratic finite elements have already been used and found to be numerically stable in other PML formulations [186,199,200,[START_REF] Assimaki | Time-domain forward and inverse modeling of lossy soils with frequency-independent Q for near-surface applications[END_REF]201,114,113,[START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF]. We validate the proposed PML approach by comparing the numerical PML solution u to a reference solution u ref computed by using a standard purely displacement-based formulation on the same time interval I = [0 , T ] with the same time step ∆t = 0.025 s, but on an extended (or enlarged) computational domain

Ω ref = [-L , L]×[-L , 0]
(containing the physical domain of interest Ω BD ) with a top free boundary and left, right and bottom fixed exterior boundaries, and whose side length L = 6 = 24 is equal to six times the size of the physical domain Ω BD , so that any outgoing waves (exiting Ω BD ) remain within Ω ref \ Ω BD and do not travel back and do not interfere with the wave solution in Ω BD during the time interval of interest I = [0 , T ] (see the first column in Figure 5). This extended domain is discretized using the same characteristic element size ∆x = 0.1 as the one considered for the spatial discretization of domain Ω. The reference finite element mesh then comprises 610 113 nodes (i.e. 1 220 226 degrees of freedom) and 304 336 quadratic triangular elements. The extended computational domain Ω ref is naturally made of the same isotropic homogeneous linear elastic material as the physical domain Ω BD . Note that the reference solution u ref computed within Ω ref is completely independent from the PML solution u computed within Ω. Both reference and PML solutions are compared only within the physical domain Ω BD .

A receiver R is located on the top free surface Γ free , at a distance of 2 from the source S in the x 1 -direction, i.e. at the position x R = (2, 0) (see Figure 3a), and records the two components of the displacement vector u = (u 1 , u 2 ) over the time period I = [0 , 20] s. Figure 4 represents the time evolutions of the horizontal and vertical displacements u 1 (x R , t) and u 2 (x R , t) at receiver R for both PML and CAL solutions compared to the reference solution. For a more quantitative analysis, the corresponding absolute errors with respect to the reference solution are plotted in a logarithmic scale. The PML transient dynamical response is very close to the reference one, whereas the CAL transient dynamical response is clearly not acceptable from approximately time t = 7 s due to multiple reflections occurring at the interface Γ between the elastic physical domain Ω BD and the absorbing unphysical layers constituting Ω PML . 2 taken at different times are displayed in Figure 5 for the reference solution computed over the extended computational domain Ω ref and truncated over the computational domain Ω = Ω BD ∪ Ω PML as well as for the PML and CAL solutions over Ω. We observe that the vertical force point-source has first generated two weak bulk waves, a longitudinal P -wave followed by a transverse S-wave (the former being faster with a lower amplitude than the latter) and then produced two strong Rayleigh waves propagating at the free surface Γ free with higher amplitudes than the bulk waves. The PML solution is in very good agreement with the reference solution inside the physical bounded domain Ω BD without any visible reflections, unlike the CAL solution that suffers from spurious reflections at the interface Γ between Ω BD and Ω PML . For comparison purposes, Figure 6 shows snapshots of the absolute error on the displacement field magnitude x → | u(x, t)u ref (x, t) | at different times for the PML and CAL solutions compared to the reference solution u ref over both the entire computational domain Ω = Ω BD ∪ Ω PML and the physical domain Ω BD . The errors obtained for the PML solution remain very small, between one and three orders of magnitude lower than that obtained for the CAL solution. The PML method is clearly more efficient than the CAL method for absorbing both bulk and Rayleigh waves.

Snapshots of the displacement field magnitude

x → u(x, t) = u(x, t) T u(x, t) = u 1 (x, t) 2 + u 2 (x, t)
To further illustrate the accuracy and efficiency of the proposed PML method, we compute the total mechanical energy E(t) stored in the physical bounded domain Ω BD as the sum of kinetic and internal energies at each time t where E k (t) is the kinetic energy and E i (t) is the internal energy, respectively defined by

E(t) = E k (t) + E i (t),
E k (t) = 1 2 Ω BD ρ u(x, t) 2 dΩ = 1 2 UBD (t) T M BD UBD (t), E i (t) = 1 2 Ω BD σ(x, t) : ε(u(x, t)) dΩ = 1 2 U BD (t) T K BD U BD (t), where u(x, t) = u(x, t) T u(x, t) = u1 (x, t) 2 + u2 (x, t) 2
is the amplitude of the vector-valued velocity field u(x, t) at time t 0 and position x ∈ Ω BD , U BD (t) and UBD (t) are the finite element real vectors associated to displacement and velocity fields u(•, t) and u(•, t) over the physical domain Ω BD , and K BD and M BD are the corresponding standard finite element stiffness and mass matrices over Ω BD . Figure 7 shows the evolutions of the kinetic energy E k (t), internal energy E i (t) and total energy E(t) (contained in the physical domain Ω BD ) with respect [START_REF] Abarbanel | Long Time Behavior of the Perfectly Matched Layer Equations in Computational Electromagnetics[END_REF][START_REF] Abarbanel | Long-Time Performance of Unsplit PMLs with Explicit Second Order Schemes[END_REF]8,[START_REF] Alvarez | A Leap-Frog Discontinuous Galerkin Time-Domain Method for HIRF Assessment[END_REF][START_REF] Antoine | Bayliss-Turkel-like Radiation Conditions on Surfaces of Arbitrary Shape[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF][START_REF] Appelö | A General Perfectly Matched Layer Model for Hyperbolic-Parabolic Systems[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF][START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF] [START_REF] Abarbanel | Long Time Behavior of the Perfectly Matched Layer Equations in Computational Electromagnetics[END_REF][START_REF] Abarbanel | Long-Time Performance of Unsplit PMLs with Explicit Second Order Schemes[END_REF]8,[START_REF] Alvarez | A Leap-Frog Discontinuous Galerkin Time-Domain Method for HIRF Assessment[END_REF][START_REF] Antoine | Bayliss-Turkel-like Radiation Conditions on Surfaces of Arbitrary Shape[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF][START_REF] Appelö | A General Perfectly Matched Layer Model for Hyperbolic-Parabolic Systems[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF][START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF] 4, 6, 8, 10, 12, 14, 16, 18, 20 s (from top to bottom) for the PML solution over the entire computational domain Ω = Ω BD ∪ Ω PML (first column) and its restriction to the physical bounded domain Ω BD (second column) and the CAL solution over Ω and Ω BD (third and fourth columns) to time t, as well as the ones of the respective absolute errors

|E k (t) -E ref k (t)|, |E i (t) -E ref i (t)| and |E(t) -E ref (t)
| plotted in a logarithmic scale, where the superscript ref over a variable denotes its reference counterpart computed using the reference solution u ref over Ω BD . We observe that the vertical force pointsource injects energy in the physical domain Ω BD from initial time t = 0 s to approximately time t = 4.5 s. The energy is then carried by bulk (P -and Swaves) and Rayleigh waves until they progressively leave the physical domain Ω BD between approximately time t = 4.5 s for the first P -wave and time t = 9 s for the last Rayleigh wave (see Figure 5). The errors on kinetic, internal and total energies computed using the CAL solution reach much higher values than the ones calculated using the PML solution due to the amount of energy carried by the spurious waves reflected at the interface Γ and coming back into the physical domain of interest Ω BD . The energy decay of the PML solution within the physical domain Ω BD is gradual and closely follows that of the reference solution without any distinguishable discrepancy, thus indicating the efficiency of the proposed PML method.

4.2 Directional force point-source buried within a two-dimensional semi-infinite isotropic homogeneous elastic medium We again study the propagation of elastic waves in an isotropic homogeneous halfspace with free-surface boundary conditions but the directional force point-source is now buried within the semi-infinite (unbounded) domain at a given depth.

We perform the numerical simulation of elastic waves propagating in the halfplane ]-∞ , +∞[×]-∞ , ] ⊂ R 2 (with a top free surface Γ free located at x 2 = = 4) submitted to a vertical force point-source buried within the semi-infinite domain at a depth , i.e. at the position x S = (0, 0) of the frame origin (see Figure 8a). We keep the same isotropic homogeneous elastic material properties as the ones considered in the previous numerical example. A similar test case has already been investigated in [240] with different geometric features and material properties. As for the previous test case presented in Section 4.1, the exact analytical solution brings bulk P -and S-waves into play as well as Rayleigh waves induced by the vertical force point-source.

We consider a squared (two-dimensional) physical domain Ω BD = [-, ] × [-, ] with side length 2 = 8, bounded by a PML region Ω PML = {x ∈ R 2 ; < |x 1 | + h and -h x 2 < -} gathering three overlapping absorbing layers (on the right, left and bottom sides of the physical domain of interest Ω BD ) with the same finite thickness h = 2 as in the first numerical experiment (see Figure 8a). The entire computational domain

Ω = Ω BD ∪ Ω PML = [--h , + h]×[--h , ]
is then a rectangular domain with length 2( + h) = 12 and width 2 + h = 10 over which the numerical simulations are run. Recall that the point-source S is located at the center of the physical domain Ω BD . Therefore, bulk and Rayleigh waves are not expected to be formed instantly as in the first numerical example, as already mentioned in [240]. We also expect the reflection and conversion of waves at the top free surface. For the same reasons as in the first numerical experiment shown in Section 4.1, we apply homogeneous Dirichlet boundary conditions u = 0 on the exterior boundary Γ ext of Ω PML , that is at |x 1 | = + h = 6 and at x 2 = -h = -6. We also impose stress-free boundary conditions σ • n = 0 at the top surface Γ free of the entire domain Ω = Ω BD ∪ Ω PML , that is at x 2 = . We use again a quadratic damping function d i (x i ) in the PML region Ω PML but with a higher amplitude A i = 12 corresponding to a very small theoretical reflection coefficient R i = 10 -8 according to [START_REF] Appelö | Perfectly Matched Layers for Hyperbolic Systems: General Formulation, Well-posedness, and Stability[END_REF].

At the spatial level, the entire computational domain Ω is discretized using an unstructured finite element mesh made of 6-nodes quadratic triangular elements 8a). The monitored time evolutions of the two components of the displacement vector for both PML and CAL solutions and their respective absolute errors with respect to the reference solution are displayed in Figures 9 and10. The PML solution closely matches the reference solution for each of the three receivers R 1 , R 2 and R 3 placed at the interface Γ between the physical domain Ω BD and the PML region Ω PML , while the CAL solution is highly distorted and exhibits large spurious oscillations from approximately time t = 6 s for receivers R 1 and R 3 and time t = 8 s for receiver R 2 due to significant spurious reflections at the interface Γ .

Ω PML h = 2 h S • R 1 • R 2 • R 3 • Ω BD = 4 Γ Γext Γ free
Snapshots of the displacement field magnitude x → u(x, t) taken at different times from time t = 4 s to final time t = T = 20 s are shown in Figure 11 for the three performed simulations, namely for the reference solution computed It can be observed that two bulk waves have been formed by the driving buried source, a primary faster longitudinal P -wave followed by a second larger shear Swave. Both bulk waves propagate towards the boundaries of the physical domain Ω BD and enter into the PML region Ω PML at normal incidence. When reaching and impinging upon the top free surface Γ free , the direct incident bulk waves give rise to reflected bulk and Rayleigh waves that interfere together. The large Rayleigh waves can be easily distinguished based upon their elliptical polarizations and large amplitudes in comparison with the bulk waves characterized by circular polarizations and smaller amplitudes. The large-amplitude Rayleigh waves can be observed on the transient dynamical responses of the vertical displacement u 2 recorded at receivers R 2 and R 3 around time t = 6 s and t = 8 s, respectively (see Figures 10c and10e). We observe that no reflected waves are visible for the PML solution (neither any spurious reflections from the interface Γ , nor any residual reflections from the fixed exterior boundary Γ ext ), whereas multiple reflections coming from the interface Γ arise for the CAL solution. For further analysis, Figure 12 shows some snapshots of the absolute errors on the displacement field magnitude x → | u(x, t)u ref (x, t) | obtained by comparing the PML and CAL solutions to the reference one. It confirms that the PML method is clearly superior to the CAL method for the absorption of both bulk and Rayleigh waves in this numerical experiment.

Figure 13 displays the time evolutions of the kinetic energy E k (t), internal energy E i (t) and total energy E(t) stored in the physical domain Ω BD , as well as that of the corresponding absolute errors

|E k (t) -E ref k (t)|, |E i (t) -E ref i (t)| and |E(t)-E ref (t)|.
We can see that some energy is injected by the point-source into the physical domain Ω BD between initial time t = 0 s and approximately time t = 4 s. Then, the total energy transported by bulk (P -and S-waves) and Rayleigh waves decreases gradually as the outgoing waves leave the physical domain Ω BD and are progressively absorbed inside the PML region Ω PML starting from around time t = 4.5 s for the first P -wave to time t = 16 s for the last reflected bulk wave (see Figure 11). The discrepancies between PML and CAL solutions in terms of energy decay are caused by the spurious waves that send energy back into the physical domain Ω BD . The numerical results highlight the poor performances of the CAL solution compared to the PML solution that is in almost perfect agreement with the reference solution.

4.3 Directional force point-source applied onto the free surface of a two-dimensional semi-infinite multi-layer isotropic homogeneous elastic medium We finally investigate the propagation of elastic waves in a two-layer isotropic homogeneous half-space excited by a directional force point-source placed at the top free surface in order to assess the accuracy and efficiency of the PML method in the presence of an interface separating two homogeneous media with different elastic material properties.

The semi-infinite (unbounded) domain corresponding to the half-plane ]-∞ , +∞[ × ]-∞ , 0] ⊂ R 2 (with a top free surface Γ free located at x 2 = 0 as in the first test case) is composed of two elastic layers (with different mechanical [START_REF] Abarbanel | Long Time Behavior of the Perfectly Matched Layer Equations in Computational Electromagnetics[END_REF][START_REF] Abarbanel | Long-Time Performance of Unsplit PMLs with Explicit Second Order Schemes[END_REF]8,[START_REF] Alvarez | A Leap-Frog Discontinuous Galerkin Time-Domain Method for HIRF Assessment[END_REF][START_REF] Antoine | Bayliss-Turkel-like Radiation Conditions on Surfaces of Arbitrary Shape[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF][START_REF] Appelö | A General Perfectly Matched Layer Model for Hyperbolic-Parabolic Systems[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF][START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF] 4, 6, 8, 10, 12, 14, 16, 18, 20 s (from top to bottom) for the PML solution over the entire computational domain Ω = Ω BD ∪ Ω PML (first column) and its restriction to the physical bounded domain Ω BD (second column) and the CAL solution over Ω and Ω BD (third and fourth columns) properties) separated by an horizontal interface located at a depth of 2 below the top free surface, i.e. at x 2 = -2 (see Figure 14a). The linear elastic material is assumed to be isotropic and homogeneous in each of the two horizontal layers. The top layer (denoted Layer 1 in Figure 14a) has the same material properties as in the two previous test cases, namely a mass density ρ 1 = 1, a Young's modulus E 1 = 2.5 and a Poisson's ratio ν 1 = 0.25 in a state of 2D plane strain, leading to the following Lamé's coefficients λ 1 = µ 1 = 1, P -and S-wave velocities c p,1 ≈ 1.7321 and c s,1 = 1, and corresponding wavelengths λ p,1 = c p,1 /f d ≈ 5.1962 and λ s,1 = c s,1 /f d = 3. The point-source S is located at the center of the free surface Γ free of the top layer (Layer 1), i.e. at the position x S = (0, 0) of the frame origin. A receiver R is also placed at Γ free away from the source S, the distance between the source S and the receiver R being equal to 2. The bottom layer (denoted Layer 2 in Figure 14a) is characterized by the following material properties: mass density ρ 2 = ρ 1 = 1, Young's modulus E 2 = 5E 1 = 12.5, Poisson's ratio ν 2 = ν 1 = 0.25, Lamé's coefficients λ 2 = 5λ 1 = 5 and µ 2 = 5µ 1 = 5, P -and S-wave velocities c p,2 = √ 5c p,1 ≈ 3.8730 and c s,2 = √ 5c s,1 = √ 5 ≈ 2.2361, and corresponding wavelengths λ p,2 = √ 5λ p,1 ≈ 11.619 and λ s,2 = √ 5λ s,1 = 3 √ 5 ≈ 6.7082. Thus, the Young's and shear moduli (resp. the wave velocities and associated wavelengths) are five times (resp. more than two times) higher in the bottom layer (Layer 2) than in the top layer (Layer 1). As a consequence, the two-layer elastic system is a dispersive medium and the propagative waves are expected to undergo multiple reflections and refractions at the material interface between the two elastic layers. As already pointed out in [START_REF] Brun | Hybrid Asynchronous Perfectly Matched Layer for seismic wave propagation in unbounded domains[END_REF], this wave propagation problem is similar to an open elastic waveguide system in which the elastic waves are guided and tightly confined within the top layer (Layer 1) with possible radiations in the bottom layer (Layer 2), due to the contrast in the material properties between the two elastic layers. In addition to bulk waves (such as compressional P -waves and shear S-waves propagating within the two elastic media) and Rayleigh waves (propagating along the free surface), Stoneley waves (propagating along the material interface between the two elastic media) can also occur.

We consider the same physical domain Ω BD = [-, ]×[-, 0] (with = 4) as in the first test case, bounded at the top by a free surface Γ free and on the other three remaining sides by a three-overlapping-layer PML region Ω PML = {x ∈ R 2 ; < |x 1 | + h and -h x 2 < -} with a thickness h = 3 /2 = 6 three times higher than in the two previous test cases (see Figure 14a). As a consequence, the PML thickness h = 6 is slightly lower than a half (1/2) of the largest longitudinal P -wave wavelength λ p,2 ≈ 11.619 and close to the shear Swave wavelength λ s,2 ≈ 6.7082 in the underneath layer (Layer 2) corresponding to a medium with higher velocity. Note that both the physical domain Ω BD and the nearby PML region Ω PML feature different material properties on either side of the horizontal interface between the two elastic layers in order to fulfill the perfect matching property at the interface Γ between Ω BD and Ω PML . The computations are performed over the entire rectangular computational domain Ω = Ω BD ∪ Ω PML = [-h , + h]×[-h , 0] with length 2( + h) = 20 and width + h = 10. We apply the same boundary conditions as in the two previous test cases, namely clamped conditions u = 0 on the right, left and bottom exterior boundaries Γ ext of Ω (i.e. at |x 1 | = + h = 10 and at x 2 = -h = -10) and stress-free conditions σ • n = 0 at the top surface Γ free (i.e. at x 2 = 0). The damping functions d i (x i ) allowing the propagative waves to be absorbed in the PML region Ω PML are identical to those used in the previous (second) numerical experiment with a quadratic polynomial degree p = 2 and a rather high amplitude A i = 12.

We use again a standard (displacement-based) finite element method based on piecewise quadratic basis functions on a triangulation of the entire domain Ω with the same characteristic element size ∆x = 0.1 as in the two previous examples. The unstructured finite element mesh of the entire domain In Figure 15, we plot the time evolutions of the horizontal and vertical components u 1 (x R , t) and u 2 (x R , t) of the displacement vector u(x R , t) at the position x R of the receiver R. Despite the numerous reflections and refractions arising from the material interface between the two elastic layers, the PML solution fits the reference solution very well at receiver R during the whole time interval [0 , 20] s, whereas the CAL solution exhibits some discrepancies. According to the time evolutions of the corresponding errors shown in Figures 15b and15d, the PML solution presents higher match with the reference solution than the CAL solution at receiver R.

Ω = Ω BD ∪ Ω PML Ω PML h = 6 h S • R • Ω BD = 4 /2 /2 Γ Γext Γ free x 1 x 2 Layer 1 (ρ 1 , E 1 , ν 1 ) Layer 2 (ρ 2 , E 2 , ν 2 ) (a)
Figure 16 displays snapshots of the displacement field magnitude x → u(x, t) at several times for all the three numerical simulations performed up to T = 20 s over the extended computational domain Ω ext (for the reference solution) and the computational domain Ω (for the PML and CAL solutions). In both PML and CAL numerical simulations, the propagative waves leave the physical domain then damp out with an exponential decay inside the absorbing layers so as to mimic the propagation inside a semi-infinite two-layered medium. Both direct bulk and Rayleigh waves can be observed as well as multiple reflected and transmitted waves at the material interface between the two elastic media. Note that a fast bulk Pwave with a low amplitude has already entered the (high velocity) bottom layer (Layer 2) and left the physical domain Ω BD after only t = 4 s. As expected, most of the waves are guided in the top layer (Layer 1) with minor radiations in the bottom layer (Layer 2). Figure 17 depicts some snapshots of the absolute errors on the displacement field magnitude obtained by comparing the PML and CAL solutions to the reference one. Despite no spurious reflections can be clearly seen from the snapshots in Figure 16, we observe that the errors computed on the PML solution inside the physical domain Ω BD are at least around one order of magnitude lower than that computed on the CAL solution. All the outgoing waves are then slightly better absorbed in the PML region Ω PML using the PML approach than using the CAL approach.

In order to provide further insight into the performances of the PML approach, Figure 18 shows the evolutions of the kinetic, internal and total energies, namely E k (t), E i (t) and E(t), accumulated in the physical domain Ω BD as functions of time t for the PML, CAL and reference solutions. Although all energy curves are almost perfectly superimposed due to the relatively large thickness h = 6 of the PML region Ω PML , the PML solution is more accurate than the CAL solution anyway. Indeed, the errors committed on the kinetic energy E k (t), internal energy E i (t) and total energy E(t) are obviously lower for the PML solution than the CAL solution, especially before all the outgoing waves have left the physical domain Ω BD , that is up to approximately time t = 14 s. Overall, the performances of the proposed PML method in terms of absorptive capacity are very good, with no discernible spurious reflections or numerical instabilities, even in the presence of material heterogeneities.

Conclusion

A review of the current state-of-the-art of the perfectly matched layer (PML) method and its variants (leading to different PML formulations) developed over the past 25 years has been first presented with a focus on the elastic wave propagation in unbounded media. An efficient PML formulation has then been proposed for the numerical simulation and modeling of second-order linear elastodynamic equations in two-and three-dimensional unbounded domains. Both time-harmonic (frequency-domain) and time-dependent (time-domain) PML formulations have been addressed. The frequency-domain PML formulation is based on a classical complex coordinate stretching and involves a specific auxiliary strain field. The time-domain PML formulation is then obtained by applying an inverse Fourier transform in time. The resulting unsplit-field PML weak formulation is discretized in space using a standard Galerkin finite element method and integrated in time through classical second-order accurate implicit time schemes, namely a Newmark time scheme combined with a Crack-Nicolson (finite difference) time scheme. Such a PML formulation turns out to be intrusive but it can be easily implemented into existing finite element codes with minor modifications. Furthermore, it does not require any field splitting, any nonlinear solver, any mass lumping, any limitation in the choice of time step size, and it does not involve any complicated (or costly) convolution operations in time or any high-order derivatives with respect [START_REF] Abarbanel | Long Time Behavior of the Perfectly Matched Layer Equations in Computational Electromagnetics[END_REF][START_REF] Abarbanel | Long-Time Performance of Unsplit PMLs with Explicit Second Order Schemes[END_REF]8,[START_REF] Alvarez | A Leap-Frog Discontinuous Galerkin Time-Domain Method for HIRF Assessment[END_REF][START_REF] Antoine | Bayliss-Turkel-like Radiation Conditions on Surfaces of Arbitrary Shape[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF][START_REF] Appelö | A General Perfectly Matched Layer Model for Hyperbolic-Parabolic Systems[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF][START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF] [START_REF] Abarbanel | Long Time Behavior of the Perfectly Matched Layer Equations in Computational Electromagnetics[END_REF][START_REF] Abarbanel | Long-Time Performance of Unsplit PMLs with Explicit Second Order Schemes[END_REF]8,[START_REF] Alvarez | A Leap-Frog Discontinuous Galerkin Time-Domain Method for HIRF Assessment[END_REF][START_REF] Antoine | Bayliss-Turkel-like Radiation Conditions on Surfaces of Arbitrary Shape[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF][START_REF] Appelö | A General Perfectly Matched Layer Model for Hyperbolic-Parabolic Systems[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF][START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF] 4, 6, 8, 10, 12, 14, 16, 18, 20 s (from top to bottom) for the PML solution over the entire computational domain Ω = Ω BD ∪ Ω PML (first column) and its restriction to the physical bounded domain Ω BD (second column) and the CAL solution over Ω and Ω BD (third and fourth columns) The performances of the proposed PML method have been illustrated on twodimensional linear elastic wave propagation problems defined in semi-infinite (unbounded) domains and stated in Cartesian coordinates. Numerical examples include (horizontally-layered) isotropic homogeneous elastic half-spaces subjected to a directional force point-source applied onto the surface (Lamb's problem) or buried within the volume. The accuracy of the PML approach is compared to a classical absorbing layer (CAL) approach with respect to a reference solution computed on an extended computational domain. Numerical results highlight that the PML acts as a highly efficient absorbing boundary condition (compared to the CAL) used to attenuate the amplitude of the waves traveling outside the physical domain of interest. Eventually, the proposed PML method could be applied to three-dimensional wave propagation problems and readily extended to other coordinate systems, such as polar, cylindrical or spherical coordinates. Also, it would be interested to extend the PML formulation to fluid-solid (or other multiphysics) interaction problems. Another ambitious perspective would consist in investigating the intrinsic stability issues arising in certain kinds of anisotropic (visco)elastic materials and exhibited by other PML formulations [START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF]192,236,242,216,[START_REF] Dmitriev | Application of M-PML absorbing boundary conditions to the numerical simulation of wave propagation in anisotropic media. Part II: Stability[END_REF]104,195,352,269,270,[START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF][START_REF] Bécache | Perfectly matched layers in negative index metamaterials and plasmas[END_REF]126,214,215]. To date, the construction of a stable PML formulation for the linear elastodynamic equations in arbitrary anisotropic elastic media (whatever the level of material anisotropy, i.e. for all the material symmetry classes ranging from isotropy to pure anisotropy) remains an open problem.

Σ 2 =

 2 tr(D)I -D = adj(D) in the two-dimensional case (i.e. for d = 2), while D Σ 3 = tr(D)I -D and D Π 3 = adj(D) in the three-dimensional case (i.e. for d = 3), where adj(D) denotes the adjugate of square matrix D (that is the transpose of its cofactor matrix). Accordingly, D Σ 2 (resp. D Π 3 ) is such that DD Σ 2 = D Σ 2 D = det(D)I = d 1 d 2 I (resp. DD Π 3 = D Π 3 D = det(D)I = d 1 d 2 d 3 I) in the two-dimensional (resp. three-dimensional) case. Besides, D Σ 3 and D Π 3 are linked through the relation DD Σ 3 = tr(D Π 3 )I-D Π 3 . Even though the system of time-harmonic PML equations (

  ) in the three-dimensional case (i.e. for d = 3), σ = C : ε( û) + C : ê (27c)

Fig. 2 :

 2 Fig. 2: Evolutions of (a) the Ricker wavelet r(t) with respect to time t over the time interval [0 , 12] s and (b) the amplitude |r(ω)| of its Fourier transform in time r(ω) with respect to angular frequency ω over the angular frequency band [0 , 4π] rad/s (corresponding to the frequency band [0 , 2] Hz)

  Bounded domain Ω BD and PML region Ω PML separated by interface Γ (b) Finite element mesh of domain Ω = Ω BD ∪ Ω PML with interface Γ = ∂Ω BD ∩ ∂Ω PML represented by white solid lines

Fig. 3 :

 3 Fig. 3: Lamb's problem: (a) two-dimensional semi-infinite domain (half-space) modeled by a physical bounded domain Ω BD surrounded by an artificial PML region Ω PML , subjected to a downward vertical force point-source S, with a receiver R, both located at the top free surface Γ free of Ω BD , and (b) associated finite element mesh

Fig. 4 :

 4 Fig. 4: Lamb's problem: evolutions of (a) the horizontal displacement u 1 , (b) the absolute error |u 1u ref 1 | on the horizontal displacement, (c) the vertical displacement u 2 and (d) the absolute error |u 2u ref 2 | on the vertical displacement at receiver R with respect to time t for the PML and CAL solutions compared to the reference solution

Fig. 5 :

 5 Fig. 5: Lamb's problem: snapshots of the displacement field magnitude x 7 ! ku(x, t)k taken at di↵erent times t =[START_REF] Abarbanel | Long Time Behavior of the Perfectly Matched Layer Equations in Computational Electromagnetics[END_REF][START_REF] Abarbanel | Long-Time Performance of Unsplit PMLs with Explicit Second Order Schemes[END_REF] 8,[START_REF] Alvarez | A Leap-Frog Discontinuous Galerkin Time-Domain Method for HIRF Assessment[END_REF][START_REF] Antoine | Bayliss-Turkel-like Radiation Conditions on Surfaces of Arbitrary Shape[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF][START_REF] Appelö | A General Perfectly Matched Layer Model for Hyperbolic-Parabolic Systems[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF][START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF] s (from top to bottom) for the reference solution over an extended computational domain ⌦ ref ⌦ = ⌦ BD [ ⌦ PML (first column) with interface = @⌦ BD \ @⌦ PML indicated by white solid lines and the exterior boundary ext of ⌦ indicated by red solid lines, the reference solution over the computational domain ⌦ = ⌦ BD [⌦ PML (second column), the PML solution over ⌦ (third column) and the CAL solution over ⌦ (fourth column)

Fig. 5 :Fig. 6 :

 56 Fig. 5: Lamb's problem: snapshots of the displacement field magnitude x 7 ! ku(x, t)k taken at di↵erent times t =[START_REF] Abarbanel | Long Time Behavior of the Perfectly Matched Layer Equations in Computational Electromagnetics[END_REF][START_REF] Abarbanel | Long-Time Performance of Unsplit PMLs with Explicit Second Order Schemes[END_REF] 8,[START_REF] Alvarez | A Leap-Frog Discontinuous Galerkin Time-Domain Method for HIRF Assessment[END_REF][START_REF] Antoine | Bayliss-Turkel-like Radiation Conditions on Surfaces of Arbitrary Shape[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF][START_REF] Appelö | A General Perfectly Matched Layer Model for Hyperbolic-Parabolic Systems[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF][START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF] s (from top to bottom) for the reference solution over an extended computational domain⌦ ref ⌦ = ⌦ BD [ ⌦ PML (first column) with interface = @⌦ BD \ @⌦ PML indicatedby white solid lines and the exterior boundary ext of ⌦ indicated by red solid lines, the reference solution over the computational domain ⌦ = ⌦ BD [⌦ PML (second column), the PML solution over ⌦ (third column) and the CAL solution over ⌦ (fourth column) Fig. 5: Lamb's problem: snapshots of the displacement field magnitude x → u(x, t) taken at different times t =[START_REF] Abarbanel | Long Time Behavior of the Perfectly Matched Layer Equations in Computational Electromagnetics[END_REF][START_REF] Abarbanel | Long-Time Performance of Unsplit PMLs with Explicit Second Order Schemes[END_REF] 8,[START_REF] Alvarez | A Leap-Frog Discontinuous Galerkin Time-Domain Method for HIRF Assessment[END_REF][START_REF] Antoine | Bayliss-Turkel-like Radiation Conditions on Surfaces of Arbitrary Shape[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF][START_REF] Appelö | A General Perfectly Matched Layer Model for Hyperbolic-Parabolic Systems[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF][START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF] s (from top to bottom) for the reference solution over an extended computational domainΩ ref ⊃ Ω = Ω BD ∪ Ω PML (first column) with interface Γ = ∂Ω BD ∩ ∂Ω PML indicated by white solid lines and the exterior boundary Γ ext of Ω indicated by red solid lines, the reference solution over the computational domain Ω = Ω BD ∪Ω PML (second column), the PML solution over Ω (third column) and the CAL solution over Ω (fourth column)

  Fig. 6: Lamb's problem: snapshots of the absolute error on the displacement field magnitude x 7 ! |ku(x, t)k ku ref (x, t)k| taken at di↵erent times t =[START_REF] Abarbanel | Long Time Behavior of the Perfectly Matched Layer Equations in Computational Electromagnetics[END_REF][START_REF] Abarbanel | Long-Time Performance of Unsplit PMLs with Explicit Second Order Schemes[END_REF] 8,[START_REF] Alvarez | A Leap-Frog Discontinuous Galerkin Time-Domain Method for HIRF Assessment[END_REF][START_REF] Antoine | Bayliss-Turkel-like Radiation Conditions on Surfaces of Arbitrary Shape[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF][START_REF] Appelö | A General Perfectly Matched Layer Model for Hyperbolic-Parabolic Systems[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF][START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF] s (from top to bottom) for the PML solution over the entire computational domain ⌦ = ⌦ BD [ ⌦ PML (first column) and its restriction to the physical bounded domain ⌦ BD (second column) and the CAL solution over ⌦ and ⌦ BD (third and fourth columns) Fig. 6: Lamb's problem: snapshots of the absolute error on the displacement field magnitude x → | u(x, t)u ref (x, t) | taken at different times t =[START_REF] Abarbanel | Long Time Behavior of the Perfectly Matched Layer Equations in Computational Electromagnetics[END_REF][START_REF] Abarbanel | Long-Time Performance of Unsplit PMLs with Explicit Second Order Schemes[END_REF] 8,[START_REF] Alvarez | A Leap-Frog Discontinuous Galerkin Time-Domain Method for HIRF Assessment[END_REF][START_REF] Antoine | Bayliss-Turkel-like Radiation Conditions on Surfaces of Arbitrary Shape[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF][START_REF] Appelö | A General Perfectly Matched Layer Model for Hyperbolic-Parabolic Systems[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF][START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF] s (from top to bottom) for the PML solution over the entire computational domain Ω = Ω BD ∪ Ω PML (first column) and its restriction to the physical bounded domain Ω BD (second column) and the CAL solution over Ω and Ω BD (third and fourth columns)

Fig. 7 :

 7 Fig. 7: Lamb's problem: evolutions of (a) the kinetic energy E k (t), (b) the absolute error |E k (t) -E ref k (t)| on the kinetic energy, (c) the internal energy E i (t), (d) the absolute error |E i (t) -E ref i (t)| on the internal energy, (e) the total energy E(t) and (f) the absolute error |E(t) -E ref (t)| on the total energy stored in the physical domain Ω BD with respect to time t for the PML and CAL solutions compared to the reference solution

  Fig. 8: Buried source problem: (a) two-dimensional semi-infinite domain (halfspace) modeled by a physical bounded domain Ω BD surrounded by an artificial PML region Ω PML , subjected to a downward vertical force point-source S buried in Ω BD , with three receivers R 1 , R 2 and R 3 located along the interface Γ between Ω BD and Ω PML , and (b) associated finite element mesh

3 Fig. 9 : 3 Fig. 10 :

 39310 Fig. 9: Buried source problem: evolutions of (a-c-e) the horizontal displacement u 1 and (b-d-f) the corresponding absolute error |u 1u ref 1 | at three receivers R 1 (first row), R 2 (second row) and R 3 (third row) with respect to time t for the PML and CAL solutions compared to the reference solution

Fig. 11 :Fig. 11 :Fig. 12 :

 111112 Fig. 11: Buried source problem: snapshots of the displacement field magnitude x 7 ! ku(x, t)k taken at di↵erent times t =[START_REF] Abarbanel | Long Time Behavior of the Perfectly Matched Layer Equations in Computational Electromagnetics[END_REF][START_REF] Abarbanel | Long-Time Performance of Unsplit PMLs with Explicit Second Order Schemes[END_REF] 8,[START_REF] Alvarez | A Leap-Frog Discontinuous Galerkin Time-Domain Method for HIRF Assessment[END_REF][START_REF] Antoine | Bayliss-Turkel-like Radiation Conditions on Surfaces of Arbitrary Shape[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF][START_REF] Appelö | A General Perfectly Matched Layer Model for Hyperbolic-Parabolic Systems[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF][START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF] s (from top to bottom) for the reference solution over an extended computational domain ⌦ ref ⌦ = ⌦ BD [ ⌦ PML (first column) with interface = @⌦ BD \ @⌦ PML indicated by white solid lines and the exterior boundary ext of ⌦ indicated by red solid lines, the reference solution over the computational domain ⌦ = ⌦ BD [⌦ PML (second column), the PML solution over ⌦ (third column) and the CAL solution over ⌦ (fourth column) Fig. 11: Buried source problem: snapshots of the displacement field magnitude x → u(x, t) taken at different times t =[START_REF] Abarbanel | Long Time Behavior of the Perfectly Matched Layer Equations in Computational Electromagnetics[END_REF][START_REF] Abarbanel | Long-Time Performance of Unsplit PMLs with Explicit Second Order Schemes[END_REF] 8,[START_REF] Alvarez | A Leap-Frog Discontinuous Galerkin Time-Domain Method for HIRF Assessment[END_REF][START_REF] Antoine | Bayliss-Turkel-like Radiation Conditions on Surfaces of Arbitrary Shape[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF][START_REF] Appelö | A General Perfectly Matched Layer Model for Hyperbolic-Parabolic Systems[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF][START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF] s (from top to bottom) for the reference solution over an extended computational domain Ω ref ⊃ Ω = Ω BD ∪ Ω PML (first column) with interface Γ = ∂Ω BD ∩ ∂Ω PML indicated by white solid lines and the exterior boundary Γ ext of Ω indicated by red solid lines, the reference solution over the computational domain Ω = Ω BD ∪Ω PML (second column), the PML solution over Ω (third column) and the CAL solution over Ω (fourth column)

  Fig. 12: Buried source problem: snapshots of the error on the displacement field magnitude x 7 ! |ku(x, t)k ku ref (x, t)k| taken at di↵erent times t =[START_REF] Abarbanel | Long Time Behavior of the Perfectly Matched Layer Equations in Computational Electromagnetics[END_REF][START_REF] Abarbanel | Long-Time Performance of Unsplit PMLs with Explicit Second Order Schemes[END_REF] 8,[START_REF] Alvarez | A Leap-Frog Discontinuous Galerkin Time-Domain Method for HIRF Assessment[END_REF][START_REF] Antoine | Bayliss-Turkel-like Radiation Conditions on Surfaces of Arbitrary Shape[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF][START_REF] Appelö | A General Perfectly Matched Layer Model for Hyperbolic-Parabolic Systems[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF][START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF] s (from top to bottom) for the PML solution over the entire computational domain ⌦ = ⌦ BD [ ⌦ PML (first column) and its restriction to the physical bounded domain ⌦ BD (second column) and the CAL solution over ⌦ and ⌦ BD (third and fourth columns) Fig. 12: Buried source problem: snapshots of the error on the displacement field magnitude x → | u(x, t)u ref (x, t) | taken at different times t =[START_REF] Abarbanel | Long Time Behavior of the Perfectly Matched Layer Equations in Computational Electromagnetics[END_REF][START_REF] Abarbanel | Long-Time Performance of Unsplit PMLs with Explicit Second Order Schemes[END_REF] 8,[START_REF] Alvarez | A Leap-Frog Discontinuous Galerkin Time-Domain Method for HIRF Assessment[END_REF][START_REF] Antoine | Bayliss-Turkel-like Radiation Conditions on Surfaces of Arbitrary Shape[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF][START_REF] Appelö | A General Perfectly Matched Layer Model for Hyperbolic-Parabolic Systems[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF][START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF] s (from top to bottom) for the PML solution over the entire computational domain Ω = Ω BD ∪ Ω PML (first column) and its restriction to the physical bounded domain Ω BD (second column) and the CAL solution over Ω and Ω BD (third and fourth columns)

Fig. 13 :

 13 Fig. 13: Buried source problem: evolutions of (a) the kinetic energy E k (t), (b) the absolute error |E k (t) -E ref k (t)| on the kinetic energy, (c) the internal energy E i (t), (d) the absolute error |E i (t) -E ref i (t)| on the internal energy, (e) the total energy E(t) and (f) the absolute error |E(t) -E ref (t)| on the total energy stored in the physical domain Ω BD with respect to time t for the PML and CAL solutions compared to the reference solution

  Fig. 14: Multi-layer problem: (a) two-dimensional (two-)layered semi-infinite domain (half-space) modeled by a physical bounded domain Ω BD surrounded by an artificial PML region Ω PML , subjected to a downward vertical force point-source S with a receiver R, both located at the top free surface Γ free of Ω BD , and (b) associated finite element mesh

Fig. 15 :

 15 Fig. 15: Multi-layer problem: evolutions of (a) the horizontal displacement u 1 , (b) the absolute error |u 1u ref 1 | on the horizontal displacement, (c) the vertical displacement u 2 and (d) the absolute error |u 2u ref 2 | on the vertical displacement at receiver R with respect to time t for the PML and CAL solutions compared to the reference solution

Fig. 16 :

 16 Fig. 16: Multi-layer problem: snapshots of the displacement field magnitude x 7 ! ku(x, t)k taken at di↵erent times t =[START_REF] Abarbanel | Long Time Behavior of the Perfectly Matched Layer Equations in Computational Electromagnetics[END_REF][START_REF] Abarbanel | Long-Time Performance of Unsplit PMLs with Explicit Second Order Schemes[END_REF] 8,[START_REF] Alvarez | A Leap-Frog Discontinuous Galerkin Time-Domain Method for HIRF Assessment[END_REF][START_REF] Antoine | Bayliss-Turkel-like Radiation Conditions on Surfaces of Arbitrary Shape[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF][START_REF] Appelö | A General Perfectly Matched Layer Model for Hyperbolic-Parabolic Systems[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF][START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF] s (from top to bottom) for the reference solution over an extended computational domain ⌦ ref ⌦ = ⌦ BD [ ⌦ PML (first column) with interface = @⌦ BD \ @⌦ PML indicated by white solid lines and the exterior boundary ext of ⌦ indicated by red solid lines, the reference solution over the computational domain ⌦ = ⌦ BD [ ⌦ PML (second column), the PML solution over ⌦ (third column) and the CAL solution over ⌦ (fourth column)

Fig. 17 :

 17 Fig.16: Multi-layer problem: snapshots of the displacement field magnitude x 7 ! ku(x, t)k taken at di↵erent times t =[START_REF] Abarbanel | Long Time Behavior of the Perfectly Matched Layer Equations in Computational Electromagnetics[END_REF][START_REF] Abarbanel | Long-Time Performance of Unsplit PMLs with Explicit Second Order Schemes[END_REF] 8,[START_REF] Alvarez | A Leap-Frog Discontinuous Galerkin Time-Domain Method for HIRF Assessment[END_REF][START_REF] Antoine | Bayliss-Turkel-like Radiation Conditions on Surfaces of Arbitrary Shape[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF][START_REF] Appelö | A General Perfectly Matched Layer Model for Hyperbolic-Parabolic Systems[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF][START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF] s (from top to bottom) for the reference solution over an extended computational domain ⌦ ref ⌦ = ⌦ BD [ ⌦ PML (first column) with interface = @⌦ BD \ @⌦ PML indicated by white solid lines and the exterior boundary ext of ⌦ indicated by red solid lines, the reference solution over the computational domain ⌦ = ⌦ BD [ ⌦ PML (second column), the PML solution over ⌦ (third column) and the CAL solution over ⌦ (fourth column) Fig. 16: Multi-layer problem: snapshots of the displacement field magnitude x → u(x, t) taken at different times t = 4, 6, 8, 10, 12, 14, 16, 18, 20 s (from top to bottom) for the reference solution over an extended computational domain Ω ref ⊃ Ω = Ω BD ∪ Ω PML (first column) with interface Γ = ∂Ω BD ∩ ∂Ω PML indicated by white solid lines and the exterior boundary Γ ext of Ω indicated by red solid lines, the reference solution over the computational domain Ω = Ω BD ∪ Ω PML (second column), the PML solution over Ω (third column) and the CAL solution over Ω (fourth column)

  Fig. 17: Multi-layer problem: snapshots of the error on the displacement field magnitude x 7 ! |ku(x, t)k ku ref (x, t)k| taken at di↵erent times t =[START_REF] Abarbanel | Long Time Behavior of the Perfectly Matched Layer Equations in Computational Electromagnetics[END_REF][START_REF] Abarbanel | Long-Time Performance of Unsplit PMLs with Explicit Second Order Schemes[END_REF] 8,[START_REF] Alvarez | A Leap-Frog Discontinuous Galerkin Time-Domain Method for HIRF Assessment[END_REF][START_REF] Antoine | Bayliss-Turkel-like Radiation Conditions on Surfaces of Arbitrary Shape[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF][START_REF] Appelö | A General Perfectly Matched Layer Model for Hyperbolic-Parabolic Systems[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF][START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF] s (from top to bottom) for the PML solution over the entire computational domain ⌦ = ⌦ BD [ ⌦ PML (first column) and its restriction to the physical bounded domain ⌦ BD (second column) and the CAL solution over ⌦ and ⌦ BD (third and fourth columns) Fig. 17: Multi-layer problem: snapshots of the error on the displacement field magnitude x → | u(x, t)u ref (x, t) | taken at different times t =[START_REF] Abarbanel | Long Time Behavior of the Perfectly Matched Layer Equations in Computational Electromagnetics[END_REF][START_REF] Abarbanel | Long-Time Performance of Unsplit PMLs with Explicit Second Order Schemes[END_REF] 8,[START_REF] Alvarez | A Leap-Frog Discontinuous Galerkin Time-Domain Method for HIRF Assessment[END_REF][START_REF] Antoine | Bayliss-Turkel-like Radiation Conditions on Surfaces of Arbitrary Shape[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF][START_REF] Appelö | A General Perfectly Matched Layer Model for Hyperbolic-Parabolic Systems[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF][START_REF] Assi | Compact second-order time-domain perfectly matched layer formulation for elastic wave propagation in two dimensions[END_REF] s (from top to bottom) for the PML solution over the entire computational domain Ω = Ω BD ∪ Ω PML (first column) and its restriction to the physical bounded domain Ω BD (second column) and the CAL solution over Ω and Ω BD (third and fourth columns)

Fig. 18 :

 18 Fig. 18: Multi-layer problem: evolutions of (a) the kinetic energy E k (t), (b) the absolute error |E k (t) -E ref k (t)| on the kinetic energy, (c) the internal energy E i (t), (d) the absolute error |E i (t) -E ref i (t)| on the internal energy, (e) the total energy E(t) and (f) the absolute error |E(t) -E ref (t)| on the total energy stored in the physical domain Ω BD with respect to time t for the PML and CAL solutions compared to the reference solution

  C : ζ = C ijkl ζ ij ζ kl and ζ 2 = ζ 2 ij (with summation over all indices)

	d(d + 1)/2, where ζ :
	without summation) and the following ellipticity (hence,
	positive-definiteness) and boundedness properties on R d (see e.g. [176]) : there
	exist constants 0 < α β < +∞ such that 0 < α ζ 2 ζ : C : ζ β ζ 2 < +∞
	for all non-zero second-order symmetric real tensors ζ = (ζ ij ) 1 i,j d of dimension

The interested reader can refer to[197,[START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF] 196,[START_REF] Appelö | Perfectly Matched Layers for Hyperbolic Systems: General Formulation, Well-posedness, and Stability[END_REF][START_REF] Appelö | A new absorbing layer for elastic waves[END_REF][START_REF] Bécache | Perfectly matched layers in negative index metamaterials and plasmas[END_REF] for the definition of well-posedness, hyperbolicity and stability of initial boundary value problems for hyperbolic systems.

As already mentioned in [161,[START_REF] Bermúdez | Perfectly Matched Layers for Time-Harmonic Second Order Elliptic Problems[END_REF] 239] and unlike stated in[216, 240], the fourth-order tensor-valued field C does have the major symmetry property, i.e. C ijkl = C klij , but both minor ones are lost i.e. C ijkl = C jikl and C ijkl = C ijlk .

Acknowledgements The authors gratefully acknowledge and thank Christian Soize, Professor at Université Gustave Eiffel, Laboratoire MSME, for very helpful discussions, constructive remarks and valuable comments.

Conflict of interest

The authors declare that they have no conflict of interest regarding the publication of this article.