{'fraction_non_alphanumeric': 0.03378064445598404, 'fraction_numerical': 0.013774631719915822, 'mean_word_length': 4.3307109557109555, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On energy-stable and high order finite element methods for the wave equation in heterogeneous media with perfectly matched layers

This paper presents a stable finite element approximation for the acoustic wave equation on second-order form, with perfectly matched layers (PML) at the boundaries. Energy estimates are derived for varying PML damping for both the discrete and the continuous case. Moreover, a priori error estimates are derived for constant PML damping. Most of the analysis is performed in Laplace space. Numerical experiments in physical space validate the theoretical results.


Introduction
Accurate and effective domain truncation schemes are necessary for efficient numerical simulation of wave propagation problems in many applications. This is because waves typically propagate over large domains, while the domain of interest usually is much smaller. Solving the problem over the whole domain may be computationally costly or even impossible. It is necessary to truncate the large domain so that the problem only has to be solved over the smaller domain of interest. When truncating the domain, an artificial boundary is introduced, and careful treatment of this artificial boundary is necessary to get reliable results from numerical simulations.
There are now two powerful and equally appealing approaches for realizing an effective domain truncation scheme. The first is to use a non-reflecting boundary condition (NRBC), [11,25,23], which is a boundary condition that is enforced at the artificial boundary such that spurious reflections from the artificial boundary are sufficiently minimized. The second one is adding an absorbing layer where the underlying equations are transformed such that waves decay rapidly as they propagate in the layer. The perfectly matched layer (PML), [7,16,13], is an absorbing layer with the desirable property that all waves enter the layer and are absorbed without any reflections, regardless of frequency and angle of incidence.
The PML technology was initially introduced for electromagnetic waves by Bérenger in [7]. The approach has since then been extended to other wave propagation problems such as elastic wave propagation, see for example [3,21,14], and acoustic wave propagation, for example, [24,1]. Please see [20] for a recent review. In this paper, we consider the PML formulation for the second-order scalar wave equation presented in [24].
The PML approach has grown popular over the last decades due to its versatility, efficiency, and ease of implementation using standard numerical methods such as finite difference, finite volume, finite element, spectral, and discontinuous Galerkin methods. One caveat is that in practice, once the PML-formulation is discretized, the PML approach is no longer completely non-reflecting. This is due to the discretization, which disturbs the perfect matching, and the finite PML width, which means that the residual outgoing waves will reflect from the outer boundary and travel back through the layer.
However, stability for the chosen continuous PML and the corresponding discretization is imperative for the approach to be efficient and robust and for results from simulations to be reliable. There is a large body of work treating the stability and well-posedness of PML formulations at the continuous PDE level. For instance, in [6,12] the Cagniard-De Hoop method was used to derive analytical solutions of the continuous PML for the 2D and 3D acoustic wave equation in a homogeneous medium. We note that the stability analysis of the PML at the discrete level has attracted less attention in the literature. This is attributable to the complex and asymmetrical forms of the PML equations, which make applying standard analytical tools such as the energy methods technically difficult.
The stability of a PML model has been studied from a mode perspective. In [5] a necessary geometric stability condition was formulated for the initial value problems of general firstorder hyperbolic systems. When violated, exponentially growing modes in time are present, making the layer useless. However, even if the Cauchy problem does not have growing modes, there can still be stability issues related to boundaries and interfaces. Mode analysis has to some degree been extended to include such issues, see [14,13,18,19,15]. However, the mode-based stability results are also cumbersome to extend to discrete settings.
A standard approach to the stability of a numerical scheme is to first derive an energy estimate for the corresponding continuous problem and then mimic the energy analysis in the discrete setting. In physical space, there are severe technical difficulties in deriving energy estimates for problems with PMLs, especially when including boundaries, heterogeneous medium parameters, and allowing for spatially varying damping parameters. In [2] firstorder constant coefficient Cauchy problems are considered. A temporally decaying energy in physical space could be found by requiring several algebraic inequalities, which must be verified in Fourier-Laplace space. However, the energy involves combinations of higher-order derivatives of the unknowns and is therefore difficult to work within a discrete setting. An important step was recently taken in [4], where energy estimates, valid in both continuous and discrete settings, are derived in physical space for the PML formulation for the secondorder wave equation presented in [24], but only for constant media and damping functions.
This paper aims to formulate a stable and accurate finite element method (FEM) for the wave equation on a domain surrounded by a PML. We take the approach of deriving an energy estimate in Laplace space, as in [17], where stability results are derived for a discontinuous Galerkin spectral element method (DGSEM) for the wave equation as a firstorder system. The analysis in [17] allows constant damping parameters, while the stability analysis in this paper extends the analysis to also include continuously variable damping coefficients. The key step taken in [17] was to construct physics-based numerical fluxes with appropriate penalty weights, which allowed the derivation of energy estimates in the Laplace space. For continuous Galerkin FEM, we show that the choice of the weak form and solution space are crucial for proving numerical stability and convergence of high order FEM for the PML for the wave equation in general domains. The chosen weak form allows us to prove energy estimates for heterogeneous medium parameters with piecewise continuous wave speed and continuously varying damping functions, both in continuous and discrete settings. In addition, we prove a convergence result. These are the main results of this paper. Also worth noting is that the paper demonstrates the strength of the strategy of doing energy analysis in Laplace space. We believe that this strategy can facilitate stability results for many other PML problems.
The outline of this paper is as follows. In Section 2 we present the underlying wave problem in three space dimensions. In Section 3 we introduce the PML formulation. In Section 4 we derive continuous energy estimates. In Section 5 the weak form is presented. Numerical stability and a priori error estimates are derived in Section 6. The numerical experiments are presented in Section 7. Finally, in Section 8 we summarise the results and draw conclusions.

Equations
We consider the time dependent acoustic wave equation on second order form in a source free heterogeneous medium where x = (x, y, z) ∈ Ω ⊂ R 3 , u : Ω × (0, T ) → R is the solution to the wave equation (1), and T > 0 is the final time. The medium is defined by the wave speed c = κ/ρ, where κ > 0 is the bulk-modulus and ρ > 0 is the density. The wave equation (1) is subjected to the initial conditions where we assume that the initial data is compactly supported. The domain Ω is bounded by the boundary ∂Ω. On ∂Ω we augment the wave equation with the well-posed linear boundary conditions where n is the outward unit normal on the boundary ∂Ω, r with |r| ≤ 1 is the reflection coefficient and Z = ρc > 0 is the wave impedance. Note that r = −1 corresponds to a Dirichlet boundary condition r = 1 yields a Neumann boundary condition and r = 0 corresponds to the classical first order absorbing boundary condition Remark 1 We have tacitly considered homogeneous boundary data under the usual assumptions of compatibility between boundary and initial data. The analysis follows through for the case of non-homogeneous boundary data, but this would complicate the algebra.
We will now derive the standard energy estimate, on which a proof of well-posedness can be based. To begin, for real functions u, v defined over Ω we introduce the inner products and for the boundary together with the corresponding norm If b ≡ 1 we will omit the subscript b. We define the energy as To begin with assume that the material parameters are continuous and there is a sufficiently smooth solution. We will now use the standard technique to derive an energy estimate. Multiply equation (1) byu and integrate by parts in space to get The boundary term vanishes for the cases when |r| = 1, which corresponds to either the Dirichlet boundary condition (4) or the Neumann boundary condition (6).
where c > 0 is the wave speed. It follows that 1 2 where In particular we note that BT ≥ 0 since |r| ≤ 1. Thus which after integration in time gives us the energy estimate If |r| = 1 then we have energy conservation, that is E(t) = E(0) for all t ≥ 0. If |r| < 1 then energy is lost through the boundaries. We can include piece-wise continuous material properties in the above. At interfaces Γ where the wave speed c > 0 is discontinuous we prescribe the interface condition where n is a unit normal at the interface. In the derivation of the energy estimate interface terms appear, but vanish due to this condition, and the same energy estimate holds also for this case. In the following analysis we will include the possibility of material interfaces.

PML
In this section we re-derive the PML problem used in [24] using the well-known complex coordinate stretching technique [9].
We begin with introducing the Laplace transform in time. Let u = L (u) be the Laplace transform of u defined by where s ∈ C with Re(s) = a > 0. Laplace transforming the wave equation (1) in time, with homogeneous initial conditions, yields We introduce the PML coordinate transformation where d η (η) are the damping functions, and s is dual Laplace time variable. Note that In the Laplace space we apply the coordinate transformation to the wave equation (16) yielding where To invert the Laplace transforms we multiply (19) with S x S y S z , we have Introducing the diagonal matrices and inverting the Laplace transforms results in the PML equation in the time-domain We will use boundary and interface conditions, which in Laplace space have the form and Note that (23) can be written as and since S x S y S z = 0 we must also have Similarly the second part of (24) can be written as Corresponding to (23) and (24) we have relations in physical space in terms of the solution and the auxilary variables. The corresponding boundary condition is and the corresponding interface conditions are where The conditions (28) and (29) are non-trivial modifications of (3) and (14) respectively. The modifications are chosen to ensure the energy estimate, and to allow for energy dissipation at the boundary. Note that r = −1 corresponds to the Dirichlet boundary condition and r = 1 corresponds to the Neumann boundary condition

Energy analysis for the PML problem
In this section we will derive the energy estimate for the PML problem in Laplace space. Consider the inhomogeneaus system with vanishing initial data and subject to the boundary and interface conditions (28) and (29). Problems with homogeneous initial and boundary data can be reformulated as (30) by a suitable change of variables. Next we Laplace transform (30) in time and get Using (31b)-(31c) and eliminating the auxiliary variables φ, ψ in (31a) yields with the boundary condition and the interface condition Here and Remark 2 In realistic settings we typically have f φ = f ψ = 0. However, for generality, we choose to include them in the analysis in this section.
In [17], an acoustic first-order system with a surrounding PML is considered. After eliminating the velocity field the pressure is shown to satisfy (32)-(34). This means that we can take advantage of the analysis there. Following [17], we note that for s = a + ib = 0, a > 0, and d η ≥ 0, with η = x, y, z, we have Next define an energy for the solution in Laplace space and an energy for the data The following theorem gives an energy estimate in Laplace space. It is similar to Theorem 1 in [17], but more general since it allows continuously varying damping.
Theorem 1 Consider the Laplace transformed PML equation (32) in a piece-wise continuous heterogeneous media with continuously varying damping coefficients d x , d y , d z ≥ 0, subject to boundary condition (33), |r| ≤ 1, and interface conditions (34) at discontinuities in material parameters. Then for any where η=x,y,z a(a+dη)+b 2 (a+dη) 2 +b 2 n 2 η dS ≥ 0, r = −1, E u and E f are defined in (36) and (37), n η are the components of the normal, and b = Im(s).

Proof 1
We multiply (32) with (s u) * /(S x S y S z ) from the left and integrate over our domain Ω, we have Integrating by part gives us where n η is the η component of the normal. Note that the last term in the first row of the above vanishes if the damping is constant in space. At discontinuities in ρ the interface terms vanishes due to the interface conditions (24), when we take into account that the damping is continous. At the boundaries we enforce the boundary conditions (23). Note that if r = −1 we have (s u) = 0 and the boundary integral vanishes. And if r = −1 we enforce the boundary conditions and we have Note also that the boundary integral vanishes when r = 1. Adding the conjugate transpose yields Here we have used give us the desired estimate, which completes the proof.

Remark 3
In particular, constant damping coefficients gives us β = 0 and the proof holds for any a > 0.
The problem is asymptotically stable in the sense that no exponentially growing solutions are supported.
Transforming our estimate in Theorem 1 from Laplace domain into time domain is hard to do while keeping the sharpness of the bound. However we can straight-forwardly achieve a bound that increases exponentially with time. This does not prove asymptotic stability but it is enough to assure well-posedness of the PML problem. Let L −1 denote the inverse Laplace transform and introduce the following norms and Theorem 2 Consider the Laplace transformed PML equation (19) with constant damping d η ≥ 0 and subject to homogeneous initial data. Let the energy norms E 2 for any a > β and T > 0.

Proof 2
The proof can be found in [17].
From the results in this section we can conclude that our PML problem is well posed. In particular, for constant damping (41) is valid for any a > 0, and we can choose a > 0 in relation to the time interval. A bound in physical space is obtained, which grows only algebraically in time. With variable damping the bound grows exponentially. In section 6 we will see that we can follow the same approach to prove stability in the corresponding sense also for the discrete problem.

The finite element approximation
In this section we formulate the PML problem in weak form, and state the corresponding FEM. For the weak form we will consider u, ψ ∈ V , where and φ ∈ [L 2 (Ω)] 3 . Note that u ∈ V is the standard choice for the wave equation without damping. In the PML u ∈ V in (31) implies φ ∈ [L 2 (Ω)] 3 and ψ ∈ V . Here H 1 0 (Ω) and H 1 (Ω) are the standard Sobolev spaces. We begin by multiplying the PML system in (31) by test functions v, q ∈ V and p ∈ [L 2 (Ω)] 3 and integrating over the domain. We obtain When r = −1 we have a Dirichlet boundary condition, which is imposed strongly in the function space for u and v. When r = −1 we need to consider the boundary condition multiplied by the test function v and integrated over the boundary, Similarly, we will use the interface condition multiplied by the test function v, and integrated over the interface, Next we perform integration by parts on the right hand side of equation (42a). If we take (44) into account the result is ψ , q = (u + f ψ , q) .
When r = −1 we use (43) to express the boundary term inu and θ, and thus the boundary term can be written as When r = −1 the term vanishes due to the strong Dirichlet condition imposed on the function space.
Remark 4 Note that both the first and second term of the right hand side of (42a) are integrated by parts. This weak form will allow the discrete energy analysis in the next section to follow the steps in Section 4. We also note that with this formulation the interface condition is imposed weakly.
Remark 5 From experience we know that naively only integrating the first term of the right hand side of (42a) by parts, yields an unstable numerical method.
Note that for |r| = 1 the boundary term vanishes identically. Since we derived this weak form from the strong form consistency follows. It remains to formulate the Galerkin approximation. For this purpose we discretise our domain Ω by introducing a tesselation K Ω , with no hanging nodes, consisting of elements, K, that completely cover Ω. More precisely Next we introduce our discrete spaces and where N p (K) denotes an element with polynomial order p ≥ 1, over K.
The Galerkin approximation for the PML problem can now be stated. For all v, q ∈ V h and p ∈ W h , find u h , ψ h ∈ V h and φ h ∈ W h such that The initial data (2) is L 2 -projected onto the finite element space.
The numerical experiments will be performed in 2D. There the approximation reads as in 3D, but the number of auxiliary variables is reduced to two, corresponding to the first two components of φ in the 3D case.

Discrete Analysis
In this section we derive energy and convergence results for particular cases. The energy result is quite general and holds for continuously varying damping coefficient, while the convergence proof is more restrictive and relies on constant damping coefficients. For simplicity we will consider Dirichlet boundary condition, but the analysis can easily be extended to the more general boundary conditions.

Discrete energy analysis
Here, we derive an energy estimate for the discrete problem with a particular choice of quadrature, which corresponds to the polynomial order. In the final part of the derivation we will assume that in each element the quadrature points can be used to interpolate any function to the local polynomial space. For later use we introduce For example, with rectangular elements we can use a tensor product Lagrange basis with piece-wise polynomials of degree p in each spatial direction based on Lagrange-Gauss-Lobatto nodes, and a corresponding Gauss-Legandre quadrature with p+1 quadrature points in each direction, to satisfy the Interpolation Condition. Such a quadrature rule can integrate polynomials of order 2p + 1 exactly. If a higher order quadrature involving more points is used, the Interpolation Condition will not be satisfied. Comments on how to extend the result to the other quadratures or exact integration will follow in the end of the section.
We will analyze the discretization in Laplace space, and therefore we consider the Laplace transform of the discrete weak form, The effects of initial data is incorporated in the forcing functionF u h , but here we have for brevity assumed forcing only in the main equation.
Here v ∈ V h and w ∈ W h . Further, the discrete inner product used in (53a)-(53c), and the corresponding norm are defined by quadrature instead of integrals. The next step is to derive a discrete reduced weak form involving onlyû h and v. We would like to usep(x) = (s * I + d(x)) −1 ∇v as test function in (53b). Here I is the identity matrix and d is a diagonal matrix defined in section 3. However, for variable dampingp is not in W h even if v ∈ V h . If the Interpolation Condition (see definition 1) is satisfied, the corresponding Quadrature Point Interpolant can be used, and we have Theorem 3 Letû h ,φ h , andψ h satisfy the discrete weak form in Laplace space, (53a)-(53c). For constant damping functions the solutionû h satisfies the reduced weak form If the Interpolation Condition (see definition 1) is fulfilled, then the solutionû h satisfies (54) also in the case of spatially varying damping functions.
Proof 3 Consider the case of spatially varying damping functions. The constant case follows in the same way. In (53c) we note thatψ h andû h are in the same space, which is spanned by the test functions q. Thus we haveψ We can thus eliminateψ h from (53b) and (53a). Next we let for any v ∈ V h . By the Interpolation Condition p h ∈ W h . With p h as test function in (53b) we get φ h , ∇v By the definitions of the diagonal matrices Γ and Υ, this relation can be rewritten as and we can also eliminateφ h from (53a), arriving at the problem in reduced weak-form, (54).
If the Interpolation Condition is not satisfied we need to find a suitable test functioñ p h ∈ W h that will allow eliminatingφ h from (53a). We will use a projection ofp(x) = (s * I + d(x)) −1 ∇v defined by Now introduce p h = Π p (s * I + d) −1 ∇v as test function in (53b), and get By the definitions of the diagonal matrices Γ and Υ, this relation can be rewritten similarly as above, and we can again eliminateφ h from (45), arriving at the problem in reduced weak-form. This proves the more general theorem Theorem 4 Letû h ,φ h , andψ h satisfy the discrete weak form in Laplace space (53a)-(53c). Thenû h also satisfies the reduced weak form (57) Here Π p is defined in Definition 2.
For simplicity we will in the remainder of this section assume that the Interpolation Condition is satisfied. It follows that the reduced weak form is (54).
We will use a similar energy for the approximate solution in Laplace space,û h , as for the exact solution and an energy for the data The discrete norms are defined using the same quadrature points as above. We can now formulate a stability result.
Theorem 5 Consider the discrete weak form (54) in a piece-wise continuous heterogeneous media and with continuous damping d η (η) ≥ 0, with homogeneous initial data. Assume that the Interpolation Condition (see Definition 1) is satisfied. There exists a β h = β + O(h p−1 ), with p being the polynomial order of the discrete space, such that for any Re(s) = a > β h the Here Eû h , E f h are defined in (58)-(59) and β is defined in Theorem 1.
i , and that for variable damping parameters the interpolation is needed to ensure that v is in the test space. It follows that Thus adding the conjugate transpose of (61) and using Cauchy-Schwarz inequality with gives us the desired estimate, which completes the proof.

Remark 6
There is a corresponding stability estimate for the discrete solution in physical space. As for the continuous problem, the estimate excludes exponential growth when the damping is constant, while with variable damping it allows bounded, exponential growth.

Remark 7
If some other quadrature is used the more general reduced form, (57), is valid. Note that the difference between integral values when using different quadratures is proportional to h p for some p > 0. Therefore we expect that when using other quadratures, or exact integration, the discrete solutionû h will satisfy a similar, but slightly perturbed energy estimate.

Convergence
Next, we turn our attention to the question of convergence. We will prove a Theorem which is valid for constant damping coefficients, d η ≥ 0. A discussion on how to extend the convergence proof to include continuously varying damping is given at the end of this subsection. However, the lemmas in this subsection are valid also for continuously varying damping functions. In this subsection we will use the Sobolev norm · H k = · W k 2 as well as the corresponding semi-norm |·| H k = |·| W k 2 . Let I h : H 1 (Ω) → V h , where V h is the space of piece-wise polynomials of degree k, be the standard Clément interpolation operator, see for example [22]. We will use the following variant of the standard interpolation result where C I is a constant independent of u, h and s, and 1 ≤ l ≤ k. Note that the first inequality in (62) holds since | 1 Sη | ≤ 1. The second inequality holds as long as u ∈ H l+1 (Ω). Next we introduce a bi-linear form that will be useful in our convergence proof. Let u and v be complex valued functions and define Lemma 1 Consider Re(s) = a > 0 and a given function f ∈ L 2 . The bilinear form (63) can be used to formulate two problems. Find u, w ∈ V s.t.
respectively. The following inequality holds as well as the following elliptic regularity results Here C R is a constant independent of u, h and s.
Thanks to (35) the following Lemma can straightforwardly be proven.
We can define finite element approximations corresponding to the two elliptic problems defined in Lemma 1. In particular we have Lemma 3 Let u h ∈ V h ∈ V be the standard finite element approximation corresponding to the first problem defined in (64), with A being the bi-linear form (63). Then u h satisfies the orthogonality Proof 6 Subtracting the variational formulation from its corresponding finite element formulation gives us the result.
We can then define our Ritz-projection R h : V → V h as: The Ritz projection satisfies the following best approximation result For all Re(s) = a > 0 the Ritz projection, R h u, defined in 69, satisfies the following best approximation result where we used the orthogonality result, the Ritz projection and the inequalities (65) and (67). This concludes the proof.
In the next Lemma we state the convergence of the Ritz projection.
Lemma 5 Assume Re(s) = a > 0. For all u ∈ H k+1 (Ω) with k ≥ 1, R h u in (69) satisfies the following two estimates and Here C 1 and C 2 are constants independent of u, h and s.

Proof 8
To begin with we note that the error estimate (71) follows from Lemma 4 and the interpolation result (62).
Next we prove the estimate (72). Let e = R h u − u and consider the elliptic problem defined in Lemma 1. Find w ∈ H 2 (Ω) such that Note that the forcing is chosen to be the error e. Choosing v = e and applying the definition of our Ritz projection yields where we used the interpolation result (62) together with result (71) and elliptic regularity (66). The constants C * I and C 2 = C * I C R are constants independent of u, s and h. C * I is obtained from applying the interpolation result (62) two times and C 2 = C * I C R , where C R comes from applying (66). This concludes the proof.
Finally we will show convergence of the reduced weak form (54). Our convergence result is for the case of constant damping. Thus the reduced weak form (54) is equivalent to Clearly the solutionû of the continuous problem satisfies the same reduced weak form for ∀v ∈ V . For this weak form we have Theorem 6 For any s ∈ C, with Re(s) = a > 0 the solutionû h to (74) with constant damping satisfies the a priori estimate Hereû is assumed to be a sufficiently smooth solution of the corresponding continuous problem and C is a constant independent of h,û and s.

Proof 9
Let e =û h −û and introduce ξ =û h − R hû and γ = R hû −û so that e = ξ + γ. For γ we have the error estimate in Lemma 5. What remains is to find a bound for ξ as well.
We have We let v = sξ and get The last term disappear by the Ritz projection (69). Adding conjugate to both sides yields a( sξ 2 + ∇ d ξ 2 ) ≤ s 2 γ sξ which further gives us Note that sξ = |s| ξ and s 4 u H k+1 = |s| s 3û H k+1 . We now have our estimate of Combining (77) and Lemma 5 together with the triangle inequality we get our desired result.
Remark 8 Theorem 6 was proven for constant damping coefficients d η . To show the corresponding convergence result for varying damping coefficients we need a discrete stability result, which is based on exact integration. As discussed above such a stability result is likely to be a perturbation of the stability estimate for the particular discretization discussed in Section 6.1. It is therefore not surprising that good convergence properties are seen in Section 7 for our discretizations of the variable damping case.

Numerical Experiments
In this section, we present some numerical experiments to validate the theoretical results presented in Section 6, as well as to quantify the numerical PML errors. We use the spatial discretisation presented in Section 5. In the experiments, we use standard continuous square quadrilateral Lagrange elements, Q p , for the solution u, where p = 1, 2, 3 denotes the order of the Lagrange polynomials. To match the theory we use discontinuous square quadrilateral Lagrange elements for the auxiliary variables. Gauss-Legendre quadrature nodes are used for integration over the elements. As boundary condition, we use the Dirichlet boundary condition defined in equation (4), which is imposed strongly. The initial condition is set to zero and we have a forcing function in the form of a Gaussian pulse centered at (0, 0). For all experiments, we let the damping profile be defined as where d 0 ≥ 0 is the damping strength coefficient and it is chosen as in [16] to balance the modeling error and the discretisation error. It is given by where ∆ pml is the PML-width, c is the wave-speed and tol is chosen as Here p is the order of the polynomial approximation, h is the element size, and C 0 = 2. The PML error tolerance, tol, is controlled by the PML with ∆ pml , the order of the polynomial approximation p, and the element size h. For a PML of finite thickness the PML error will converge to zero at the rate of the accuracy of the underlying numerical method [16,17,15]. For time integration we will use the fourth-order Runge-Kutta method.

Homogeneous problem
In this subsection we consider the homogeneous case i.e. κ = ρ = 1, and first investigate the numerical PML-error in maximum norm. In order to measure the numerical PML-error, we simulate both the PML problem as described in this section, as well as a reference problem. The reference problem is a problem defined on the larger domain Ω r = [−12, 12] × [−12, 12], on which we solve the wave equation without any PML. The size of Ω r is chosen so that the reflections from the outer boundary will not reach Ω in during the simulation. We use the same initial condition, the same source, the same time step, and the same spatial resolution. The numerical PML-error is defined as the difference between the PML solution and the reference solution, and is calculated at every node in Ω in , and at every time step. The simulation is performed for t ∈ [0, 10]. At the end of this time interval a reflected wave in the reference problem will just have reached the outer boundary of Ω in . The computations are done with Q1, Q2, and Q3 elements, as well as for different element sizes. We use a time step ∆t = 0.01 for all the elements and all spatial resolutions. This time step is chosen so that the dominant error for all resolutions and all orders will come from the space discretisation.
In Figures 2a, 2b, and 2c we have plotted the maximum PML-error in Ω in as a function of time for different resolutions and for different polynomial orders. In Figure 2d the convergence of the PML-error in maximum norm can be seen. The error is measured at the end time t = 10. We can see that the solution converges with at least the expected rates as we refine the element size.
To verify the stability results from Section 6 we look at the long-time behavior of our simulated PML-problem. For this purpose solve the same problem as before but for times t ∈ [0, 150], so that we can see what happens when the simulation is allowed to run for a longer time. We then measure the maximum amplitude of our solution in Ω in . The result can be seen in Figure 1a. Note that no long-time growth of the solution can be observed.
In Figure 3 we can see snapshots of the solution at times t = 2, t = 4, t = 6 and t = 12. We can follow how the Gaussian pulse propagates in Ω in and vanishes as it enters the PML.

Heterogeneous problem
Next we consider the heterogeneous problem. The setup is the same as for the homogeneous problem except that we let the wave speed vary in the domain and we have chosen our end time to be t = 14 to make up for the lower wave speed. The wave speed is chosen as  Note that the element boundaries are chosen so that they aline with the material interfaces.
We perform the same experiments for the heterogeneous case as for the homogeneous case. In Figures 4a, 4b and 4c we can see the maximum PML error as a function of time for different resolutions and for different elements Q1, Q2 and Q3.
In Figure 4d we can see the convergence of the PML-error in maximum norm at the end time t = 14. We can note that the convergence rate is slightly lower than expected for the Q3 elements.
In Figure 1b we can see the long-time behavior of the solution to the heterogeneous problem. To capture the long-time behavior we run the simulation until t = 150, just as in the homogeneous case. No asymptotic growth is observed, which again agrees with the stability analysis in Section 6. Finally, we can see snapshots of the solution at different times in Figure 5. We can follow how the Gaussian pulse propagates and vanishes as it enters the PML. One can also see that the wave spreads faster when y < −2.4 and slower when y > 2.4.

Summary and Discussion
We have presented a stable FEM for the acoustic wave equation on second-order form with a PML. To begin with, we prove an energy estimate in Laplace space, which extends the result in [17] from constant damping to continuously varying damping. We also formulate the problem weakly in terms of a bi-linear form. In this formulation, we use the standard H 1 space for the physical variable and a combination of H 1 and L 2 spaces for the auxiliary variables. Our FEM is then obtained by restricting to the corresponding continuous and discontinuous piece-wise polynomial finite-dimensional subspaces. The chosen weak form allows us to mimic the energy estimate in a discrete setting.
An a priori error estimate for the case of constant damping follows by applying standard   Figures 2a, 2b, 2c for elements Q1, Q2 and Q3 respectively. In Figure 2d we can see the PML error convergence for Q1, Q2 and Q3 elements at end time t = 10 (c) Snapshot at t = 6 (d) Snapshot at t = 12   Figures 4a, 4b, 4c for elements Q1, Q2 and Q3 respectively. In Figure 2d we can see the convergence of the PML error for Q1, Q2 and Q3 elements at end time t = 14 techniques in Laplace space. It indicates convergence of order p + 1, where p is the local polynomial order. For variable damping, a technical difficulty is related to the fact that the discrete stability relies on quadrature while exact integration is used in the continuous case. This prevents the straightforward extension of the convergence result to spatially varying damping.
Numerical experiments for the case of varying damping verify the stability and convergence results. We compute in 2 space dimensions and use square, Cartesian elements. Note that the theoretical results presented in Section 6 are more general and cover the 3D problem with no hanging nodes as the only restriction on element shape. Therefore we are confident that when the method is applied to a 3D problem, it will be stable and have good accuracy properties.
Note that the method uses a discontinuous function space for some of the auxiliary variables. Experiments using continuous function spaces for all variables were also performed, but with less accurate results. The degradation of accuracy was pronounced for piece-wise linears, in analogy with results in [10].
A standard approach to stability is to derive energy estimates in physical space. The analysis in this paper is done in Laplace space, both for the continuous and the discrete problems. It demonstrates how an energy estimate in Laplace space for the continuous problem can be transformed to stability and accuracy results for a discrete method. By performing analysis in Laplace space, the energy approach is applicable to many more problems, and it is our expectation that the techniques presented in this paper will be useful for many different problems.