A Modified PML Acoustic Wave Equation

In this paper, we consider a two-dimensional acoustic wave equation in an unbounded domain and introduce a modified model of the classical un-split perfectly matched layer (PML). We apply a regularization technique to a lower order regularity term employed in the auxiliary variable in the classical PML model. In addition, we propose a staggered finite difference method for discretizing the regularized system. The regularized system and numerical solution are analyzed in terms of the well-posedness and stability with the standard Galerkin method and von Neumann stability analysis, respectively. In particular, the existence and uniqueness of the solution for the regularized system are proved and the Courant-Friedrichs-Lewy (CFL) condition of the staggered finite difference method is determined. To support the theoretical results, we demonstrate a non-reflection property of acoustic waves in the layers.


Introduction
It is quite important to effectively truncate an unbounded domain in wave propagation simulations in open space, where the perfectly matched layer (PML) methods that surround the domain of interest with thin artificial absorbing layers are popularly used in easy and effective ways.After the method was introduced by J. P. Bérenger [1], which involves splitting a field into two nonphysical electromagnetic fields, many studies were conducted regarding the PML method and its modified reformulations in many different wave-type equations.These include Maxwell's equations [2,3], elastodynamics [4,5], linearized Euler equations [6][7][8][9], Helmholtz equations [10], and other types of wave equations [10][11][12].Most PML models by the splitting technique, named a split PML method, yield a hyperbolic system of first order partial differential equations [1,6,[13][14][15].It is known that the split PML models demonstrate excellent overall performance from the viewpoint of applications.However, it was pointed out in [7,16,17] that Bérenger's split, as well as other split models, transform Maxwell's equations from being strongly hyperbolic into weakly hyperbolic.These transforms imply a transition from strong to weak well-posedness in the Cauchy problem and may lead to ill-posedness under certain low-order damping functions in PML layers [18].The authors of [6,19] mention that the use of artificial dissipation is necessary to stabilize the numerical scheme of such formulations for long-time simulations.
The resulting concerns about the well-posedness and stability of the split PML models have prompted the development of other PMLs.Some examples of such developments, without splitting the fields, include un-split PML models using convolution integrals [20,21] and auxiliary variables [17,22,23].In contrast to the split PML models, it is known that the un-split PML wave equations are more effective at time discretization [22] and does not make the use of additional memory for the nonphysical field variables.However, it has also been found that the un-split PML models are susceptible to developing gradual instabilities in long-time simulations [10,19].To overcome this instability issue, various studies are reported: a low-pass filter inside the absorbing layer [6], selective damping coefficients [24], a new layer by regularizing the damping terms [8], a change of variable [25], etc.These issues are the motivation for the mathematical study of the well-posedness and stability for the un-split PML acoustic wave model in various sound speed.A time-domain analysis of PML acoustic wave equation with a constant sound speed is presented with a time-dependent point source in two dimensions using the Cagniard-de-Hoop method [25,26], which includes the time-stability and error estimates.However, it is not easy to extend the analysis to general initial value problems in variable sound speed, because those include not only straight propagating but also evanescent waves [27].There is another approach to demonstrating the well-posedness and stability by investigating the eigenvalues of the Cauchy hyperbolic problems for the PML wave equations [4,7,12,[16][17][18]28].This approach gives a restricted result when the original formulation of the PML wave equation is considered in a bounded domain, in which the solutions should be affected by boundary conditions.
Alternatively, energy techniques are used to analyze the issue of stability for the PML wave equations by presenting the energy behavior for the solution in each model [12,16,29].In general, the restriction of the PML equations to the computational domain coincides with the original problem [12], so that damping terms are required to vanish identically in the computational region.As the constant damping function can be considered as the Heaviside function, the equation [12,16,29] is not valid at the interface between the domain of interest and the layers for the constant damping case from a discontinuity.However, all these approaches only provide its well-posedness, the stability has not been clearly proved in finite PML acoustic wave equations with variable sound speed.
The main contribution of this manuscript is not only to introduce a regularized system of the second order PML acoustic wave equation that exhibits well-posedness without losing the non-reflection property of PMLs, but also to demonstrate its numerical stability.To construct the system, we adopt a regularization technique for the term ∇ • q that has a lower regularity, which is introduced in [8], to regularize the PML model for the Maxwell equation, where q is the auxiliary variable (see (2)).The standard Galerkin approximation and energy estimation of the solution are used to show the well-posedness of the regularized system.A concrete energy estimate yields the boundedness of the solution (see Theorem 1) together with the existence and uniqueness of the solution under the regularity assumption of the damping terms σ x , σ y ∈ L ∞ (Ω) (see Theorem 2).As a numerical scheme for the regularized system, a family of finite difference schemes using half-step staggered grids in space and time is used.All spatial and temporal derivatives are discretized with central finite differences that maintain the second order approximation in both space and time, respectively.A concrete von Neumann stability analysis for the numerical scheme indicates that the scheme is stable under the Courant-Friedrichs-Lewy (CFL) condition between the temporal and spatial grids (see Theorem 3).The novel features of this study include the good performance of the solution that present not only the well-posedness and stability but also the non-reflection property of the wave propagation compared to the classical PML model; even the regularized system does not possess PMLs in the original wave equation.This novelty is numerically illustrated in Section 4.
The remainder of the manuscript is organized as follows.Section 2 describes a regularized system for the un-split PML model of the acoustic wave equation and also contains the well-posedness of its solution based on the energy estimation.In Section 3, we develop a staggered finite difference scheme for the regularized system and determine the CFL condition for the numerical stability.In Section 4, several numerical results are presented to support our theoretical analysis and demonstrate the efficiency of the regularized system.Finally, some discussions are given in Section 5.

Regularized System
The aim of this section is to introduce a modified PML system using a regularization technique in a classical PML model for the acoustic wave equation.For the sake of argument, we let H 1 (Ω) := {ϕ : ϕ, ∂ x ϕ, ∂ y ϕ ∈ L 2 (Ω)} and H −1 (Ω) be the Sobolev space and dual space of H 1 0 (Ω), respectively.
The target problem we consider with here is a general second order acoustic wave equation with a variable sound speed c(x) > 0 described by with initial conditions u(•, 0) = f and u t (•, 0) = 0, where supp( f ) ⊂ Ω 0 with a domain Here, T > 0 and the sound speed c(x) is bounded by Let the domain surrounded by PML layers, where a, b, L x , L y > 0. Using a complex coordinate stretch, we consider the following system of the PML wave equation which is introduced in [28 with the initial conditions and the zero Dirichlet boundary condition u(x, •)| ∂Ω = 0, where Here, the damping terms σ x := σ x (x) and σ y := σ y (y) are assumed to be nonnegative C 0 functions which vanish in the computational domain in the sense of the analytical continuation of the PML.
Please note that a weak solution (u, q) of ( 2) is in H 1 0 (Ω) × L 2 (Ω), i.e., ∇ • q ∈ H −1 (Ω), which regularity is not enough to show the existence.In order to provide regularity on the term by an operator, we introduce a mollifier ρ .
|x| and satisfies . Then, we consider the operator where δ * ε : , as ε → 0 in the strong operator topology (see, for detail, Theorem 3 on page 7 in [30]).Then, we obtain Please note that δ ε is a linear and bounded operator from H −1 (Ω) to L 2 (Ω).Now, following [8,31], we introduce a regularized system of the classical PML model (2) by using with initial and boundary conditions The remainder of this section details the analysis of the well-posedness of the solution to the regularized system (4) based on the energy estimation under the assumption that the dampings σ x and σ y are in L ∞ (Ω).

Energy Estimate of Weak Solution
We assume that the damping functions σ x , σ y satisfy σ x , σ y ∈ L ∞ (Ω), which implies that under the condition of c(x) = 1 in the layers of the PML model ( 2), where • ∞ denotes the L ∞ -norm.Under these assumptions, the aim of this subsection is to provide an energy estimation of the weak solution of (4) in the sense that with , and almost everywhere 0 ≤ t ≤ T and the initial data satisfy Here, < •, • > denotes the duality pairing between H −1 (Ω) and H 1 0 (Ω), and (•, •) is the inner product in L 2 (Ω).In addition, the time derivatives are understood in a distributional sense.
To investigate the weak solution of (4) that satisfies ( 7) and ( 8), we use the standard Galerkin approximation and estimate the energy of the solution, which will be used to show the well-posedness of the regularized system (4) in the subsequent subsection.Let {w j |j ∈ N} be an c −2 -weighted orthonormal basis in L 2 (Ω), i.e., (c −2 w j , w k ) = δ jk , where the Kronecker delta is given by δ jk = of the eigenfunctions of the eigenvalue problem Let U k be the subspace generated by the orthonormal system {w 1 , Then, one can see that U k also becomes the c −2 -weighted orthogonal basis of H 1 0 (Ω) in the sense that Let us also denote Q k , which is the space generated by the smooth functions whose coefficients are satisfied for all theory of ordinary differential equations guarantees the existence of the approximation u k (t), q k (t) satisfying ( 9) and (10).
The following theorem gives a uniform bound of energy of the approximate solutions (9), which allows us to send k → ∞.Theorem 1.There exists a constant C T > 0 that depends only on σ x , σ y , Ω, and T such that for k ≥ where the energy E k (t) is defined by Proof.Please note that u k t ∈ U k and q k ∈ Q k .Hence, we apply u k t and q k in the first and second equations of (10), respectively, to obtain for almost everywhere 0 ≤ t ≤ T. Combining the two equations with the equality 1 where Based on the linear bounded operator ϕ −→ δ ε (ϕ), Hölder's inequality, assumptions for σ x , σ y , and Poincaré inequality, it can be noted that E k (t) satisfies the inequality Furthermore, by applying Gronwall's inequality, Poincaré inequality, and (1) in the above equation, one can obtain for some where From ( 9) and ( 10), we have Thus, we have .
Consequently, we obtain The proof is carried out by combining (11) and (12).

Existence and Uniqueness
In this subsection, we will discuss the well-posedness of the regularized system by demonstrating the existence and uniqueness of the solution (6) based on the result of Theorem 1.

Remark 3.
The most important concern in the proof is the estimation of the term δ ∇ • q in the regularized system, which has roles of a convolution, improving the stability of the system from the regularization of the term from H −1 (Ω) to L 2 (Ω).

Numerical Scheme
The aim of this section is to introduce a staggered finite difference method for discretizing the regularized system and to find a stability condition for the numerical scheme.For the staggered finite difference method, we use a family of finite difference schemes [33] with half-step staggered grids in space and time.All spatial derivatives are discretized with the centered finite differences over two or three cells, which guarantees a second order approximation in space.For the time discretization, we also use the centered finite differences for the first and second order time derivatives on a uniform mesh, which is also of the second order approximation in time.Based on the standard von Neumann stability analysis technique, we analyze the stability of the numerical scheme and obtain its CFL condition.

Staggered Finite Differences
Let t > 0 denote the time step size and x > 0 and y > 0 denote the spatial mesh sizes in the x and y directions, respectively.In addition, we also introduce the time step t n = n t and the spatial nodes x i = i x and y j = j y for n ∈ N ∪ {0} and i, j ∈ Z.We also define staggered nodes in the time direction and the x and y directions, respectively, as x, and y j± 1 2 = y j ± 1 2 y for n, i, j ∈ N. To simplify the notation, we denote u n i,j := u(t n , x i , y j ) and , y j+ 1 2 ) for q = (q x , q y ), α = x, y.For the discretization of the regularization defined in Remark 1 for the regularized system, the smooth function ρ (x, y) chosen in the following examples is constant on a rectangle centered at zero, where For a given two-dimensional finite difference grid with spatial sizes x and y, a possible choice of k is 1 = n x x and 2 = n y y with n x , n y ∈ N.For instance, with n x = n y = 1 and the usual integration formula (see Chapter 3 in [34]), we discretize the regularized term δ ε (v) i,j := (ρ * v) i,j , using the 9-point central difference formula, as follows: Let us now introduce new notations and for k = i, j, α = x, y, ).
Based on these notations, the staggered finite difference scheme for discretizing the regularized system is defined in the following steps.
Step 1. Compute q 2 where the cell averages ∂x u n i+ The definition of the cell averages allows us to compute the regularized term in (3) and ∂ y q n y i,j = 1 2 ∂y q n+ 1 2 , where the cell averages of the derivatives of the function (q .

Stability Analysis
To obtain the stability condition of the staggered finite difference scheme defined above, we restrict our concern to the constant damping case with σ x = σ y = σ 0 ≥ 0 for simplicity in our analysis.The stability condition for the scheme in the computational domain is as follows.
Remark 4. The CFL condition of scheme (13)-( 14) in the computational area (i.e., for x = y = h from the standard von Neumann stability analysis technique.
Generally the stability condition for the staggered finite difference scheme developed in Section 3.1 can be obtained as follows.
To prove Theorem 3 and use the technique of the standard von Neumann stability analysis, we recall the definition of the simple von Neumann polynomial and some of its properties as follows.
Definition 1.A polynomial is a simple von Neumann polynomial if all its roots, r, lie on the unit disk (|B(0, r)| < 1) and its roots on the unit circle are simple roots.
The following theorem demonstrates that a simple von Neumann polynomial can be a sufficient stability condition.Theorem 4. A sufficient stability condition is that φ be a simple von Neumann polynomial, where φ is the characteristic polynomial (see [35] for the proof).
With Theorem 4, the stability condition for a polynomial is presented in the following.Theorem 5. Let φ be a polynomial of degree p written as and the conjugate polynomial φ is defined as where c is the complex conjugate of c.The main ingredient in the proof of the theorem is Rouché's theorem; the proof is detailed in [36].Now, we can computationally verify the stability condition (15) in Theorem 3 using Theorems 4 and 5.
Proof of Theorem 3 .Assume that σ x = σ y = σ 0 in scheme ( 13)-( 14) and we rewrite the scheme as the second order central difference scheme of the variables u and q.
By von Neumann analysis, we can assume a spatial dependence of the following form in the field quantities: i,j = ûn+1 (k x , k y )e ik x x i +ik y y j , u n i,j = ûn (k x , k y )e ik x x i +ik y y j , where k x , k y , is the component of the wave vector k, i.e., k = (k x , k y ) T , and the wave number is k = Then, we have the system ûn+1 , ûn , qn+ , where the amplification matrix G of scheme ( 16), ( 17) is given by , where C qx and C qy satisfy . Then, it is noted that the characteristic function of G is given by Please note that |η| < This inequality gives the CFL condition (15), which completes the proof.
Remark 5. From the proof of Theorem 3, we notice that the characteristic function φ of the amplification matrix G does not depend on any quantity related to the regularized term.That is, the staggered finite difference scheme corresponding to the classical PML model (2) with a constant damping in the layers is stable under the CFL condition (15).

Numerical Result
The aim of this section is to provide numerical evidence of the well-posedness of the regularized system and the non-reflection properties of the acoustic wave in the layers of the classical PML model.For the discussion of the non-reflection properties, we demonstrate the behavior of the maximum error at t n defined as the maximum of the differences between the numerical solution and a reference solution in the computational domain Ω 0 := [0, 1] × [0, 1].Here, the reference solution is taken in the same computational domain instead of the layers with an additional large domain, for example, 15 times wider in the x and y directions in our experiment, causing the wave in the computational domain to be unaffected by the wave propagating from outside in the chosen long-time step.Furthermore, we use the energy method introduced in [37] and numerically examine the well-posedness or stability of the model ( 4) by observing the long-time behavior of the acoustic wave energy defined by For the numerical simulation, we use the same initial condition defined by (4) and, in the absorbing layer, the damping function of the form given by where β = 0, 1, 2, σ 0 is a given constant and L denotes the thickness of the layers.
For the comparisons of non-reflection property, we first demonstrate the maximum error for both Formulas ( 2) and ( 4) with two sets of thickness and damping as (L, σ 0 ) = (0.25, 30) and (L, σ 0 ) = (0.1875, 30).The numerical results are displayed in Figure 1.The classical PML has slightly smaller errors than the modified one in both cases, as shown in Figure 1, but it can be observed that these errors of the modified one can be reduced by simply increasing small amounts of thickness or damping such as L = 0.27 or σ 0 = 35.To see the influence of absorbing property by incidence angle, we demonstrate both formulas with different positions of source function.The resulted differences between reference and computed values of the solution during simulation at one point within the computational domain are plotted in Figure 2. The errors of the classical PML have relatively smaller than the modified one and both formulas have slightly better absorbing property when the angle of incidence to the interface between the computational domain and the layers is bigger.Next, to investigate the energy E (t) behavior, we choose a time step size t of x/3, which satisfies the CFL condition (15) to guarantee the stability of the staggered finite difference scheme (see Remark 4).Here, the first order backward and second order central finite differences in time and space, respectively, are used to discretize the energy E (t n ) of (18) at each time step t n .We investigate the behavior of the energy for a long-time simulation at time t n =10,000 according to the thickness of the layers and magnitude of the damping.The numerical results are displayed in Figure 3: (a) the energy with various dampings σ 0 = 40, 50, 50, 60, 70 for a fixed thickness L = 0.0625 and (b) the energy with various thicknesses L = 0.0625, 0.1, 0.125, 0.15 for a fixed damping σ 0 = 50.The results indicate that the numerical stability of the modified formula is consistently stable in the long-time simulation regardless of the magnitudes of damping and thickness of the layer.This provides proof of the well-posedness of the developed system and numerical stability for the finite difference method.Lastly, in order to illustrate this visual investigation, we consider the damping β = 2 and display the snap shots of the wave propagation at times t n = 1, 30, 60, 100, 130, 150, 200, 300, 500 with σ 0 = 35, L = 0.25 in Figure 4.One can see that the regularized system displays a good property of non-reflection in the layers, which is the purpose of building the layers .It is remarkable that from a mathematical point of view, the analytical well-posedness without losing the non-reflection property in the layers of that the classical PML model.

Discussion
We have introduced a new and efficient formulation related to the acoustic wave equation based on the regularization of the un-split PML wave equation.By regularizing the lower order regularity term in the original equation and the standard von Neumann stability analysis, we have achieved well-posedness as well as numerical stability of the solution in the new formulation.We summarize the main novelty and results of this study as follows: (1) We have proved the analytical well-posedness of our formulation without any restriction of damping terms; (2) a staggered finite difference scheme for the formulation is introduced and numerical stability is also analyzed; (3) several numerical tests are exhibited to show the numerical stability and a non-reflection property.
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Figure 2 .
Figure 2. Comparison of the difference at a point from different positions of source function with σ 0 = 35 and L = 0.1