DEVELOPING AND ANALYZING AN EXPLICIT UNCONDITIONALLY STABLE FINITE ELEMENT SCHEME FOR AN EQUIVALENT B´ERENGER’S PML MODEL

. The original B´erenger’s perfectly matched layer (PML) was quite effective in simulating wave propagation problem in unbounded domains. But its stability is very challenging to prove. Later, some equivalent PML models were developed by B´ecache and Joly [ ESAIM: M2AN 36 (2002) 87–119] and their stabilities were established. Hence studying and developing efficicent numerical methods for solving those equivalent PML models are needed and interesting. Here we propose a novel explicit unconditionally stable finite element scheme to solve an equivalent B´erenger’s PML model. Both the stability and convergence analysis are proved for the proposed scheme. Numerical results justifying the theoretical analysis are presented. We also demonstrate the effectiveness of this PML in simulating wave propagation in the free space. To our best knowledge, this is the first explicit unconditionally stable finite element scheme developed for this PML model.


Introduction
In 1994, Bérenger [7] introduced the perfectly matched layer (PML) techqniue to develop efficient numerical absorbing boundary conditions for solving the time-dependent Maxwell's equations in unbounded domains.This PML technique has been shown to be very effective and can absorb all impinging waves over a wide frequency range.Since 1994, many different PML models have been developed and applied to solve Maxwell's equations in both time-domain [1,11,15,21,25,26,32] and frequency domain [3,9,16].More details and references can be found in the classic computational electromagnetic book [30], a review paper on PMLs [31], and our recent book on metamaterials ( [22], Chap.8).Furthermore, the PML technique has been extended to solve other wave propagation problems in different media, such as acoustics, elastodynamics [2,14], elasticity [20], anisotropic dispersive media [6] and metamaterials [5,12,13].Due to many potential applications of metamaterials such as design of invisibility cloaks, recently there has been a growing interest in the study of the Maxwell's equations involving metamaterials (e.g., [8,23,24,33]).
Even though the original splitted Bérenger PML works very well in practical numerical simulations, the mathematical analysis of its energy stability for the general variable damping function case is still an open problem.Until 2002, Bécache and Joly [4] managed to establish a stability result for an equivalent Bérenger PML model.Inspired by their work, recently we [17,18] proposed and analyzed some finite diffence and finite element methods for solving those equivalent PML models.Those schemes we proposed so far are explicit and has a time step constraint.In 2021, we [19] discovered a novel technique in constructing an explicit unconditionally finite element method for solving Maxwell's equations in the free space and in Drude metamaterials.One major contribution of this paper is that we successfully extend that technique to construct and prove an explicit unconditionally stable scheme for solving this complicated equivalent Bérenger's PML model.This new scheme is not only easy in implementation but also very efficient like other explicit leapfrog type schemes.To the best of our knowledge, this is the first explicit unconditionally stable finite element scheme constructed and analyzed for this PML model.
The rest of the paper is organized as follows.In Section 2, we first present the PML model equations and construct a semi-discrete scheme for the PML model.This semi-discrete scheme is a small perturbation of the usual leapfrog scheme for this PML model.This small pertubation plays the magic for the construction of an explicit unconditionally stable scheme.In Section 3, we first develop our fully-discrete finite element scheme.Then we prove the unconditional stability and the optimal error estimate for this scheme.In Section 4, we present some numerical results to confirm our theoretical analysis and further apply our scheme to simulate some practical wave propagation problem to show the long stability of the scheme and the effective wave absorbing property of this equivalent PML model.We conclude the paper in Section 5.
Reducing all 's in (4) by 1, adding the result with (9), and denoting

2
, we have Reducing all 's in ( 6) and ( 8) by 1, and subtracting the result from (11) and (13), respectively, we have and Then substituting (15) and ( 16) into ( 14), we have Reducing all  ′  of (5) by 1, and adding the result with (10), we obtain Adding ( 6) and (11) together, adding (7) and ( 12) together, and adding (8) and ( 13) together, we have and Using the following average operator and central difference operator in time: we can be rewrite ( 17)-( 21) as follows: Remark 1.It is interesting to remark that the semi-discrete scheme (22a)-(22e) can be seen as a small perturbation of a standard explicit leapfrog scheme for solving (1a)-(1e) with the first equation added by the last ( 2 ) term.With this extra term, we can develop and prove an explicit unconditionally stable scheme out of (22a)-(22e) in the next section.

The fully discrete scheme and its analysis
To solve the problem (1a)-(1e) by a finite element method, we partition the physical domain Ω by a family of regular triangular mesh  ℎ with maximum mesh size ℎ, and adopt the -th ( ≥ 1) order Raviart-Thomas-Nédélec (RTN) mixed finite element spaces  ℎ and  ℎ [27][28][29]: For any  ≥ 1, where p denotes the space of homogeneous polynomials of degree , and   denotes the space of polynomials of degree less than or equal to  in variables , , respectively.To impose the PEC boundary condition (2), we denote the subspace  0 ℎ = { ∈  ℎ :  ×  = 0 on Ω}.Now we construct the following leapfrog type scheme for (1a)-(1e): given initial approximations Note that for theoretical analysis purpose, we multiplied     to both sides of (23d).
The initial approximations  0 ℎ , ̃︀  0 ℎ can be simply obtained by the Nédélec interpolation, while and be done by a Taylor expansion followed by a standard  2 projections.More specifically, we have the following approximations: where Π   0 ∈  ℎ denotes the Nédélec interpolation operator, and  ℎ denotes the standard  2 projection operator into the space  ℎ .It is known that the following interpolation and projection error estimates hold true [27]: Here and in the rest of the paper, we denote ‖•‖ for the  2 norm on Ω.

The unconditional stability analysis
This subsection is devoted to the unconditional stability analysis of our scheme (23a)-(23e).
Theorem 1. Denote the discrete energy at time   as: Then under the time step constraint we have the following stability for the scheme (23a)-(23e): for any  ≥ 0, Remark 2. We want to remark that the stability ( 28) is unconditionally stable, since the time step constraint (27) is independent of mesh size ℎ.Furthermore, the stability (28) has exactly the same form as the following continuous stability established in our previous work ([18], Thm.1): where the energy is denoted as: in (23c), and adding the results, we obtain: Choosing ̃︀  ℎ = ̃︀ in (23e), respectively, we have and Adding ( 29)-( 31), we attain Choosing ̃︀ in (23d), and using it to replace the fourth term in (32), we have Adding ( 33) and ( 34) together, then summing up from  = 0 to , we obtain Similarly, we can obtain where we used the estimate to obtain the first inequality.
It is easy to check that the following identity holds true: Substituting the above estimates (36)-( 38) into (35), we have To get a nice stability result, now we relax the bound for those coefficients on the right hand side of (39) by the following simple estimates: Using the following time step constraints (equivalent to ( 27)): dropping the nonnegative term of the left hand side of (39), then using the discrete Gronwall inequality, we have which completes the proof of (28).

The convergence analysis
In this subsection, we will carry out the convergence analysis of our scheme (23a)-(23d).To accomplish that, we need the following lemma.
Theorem 2. Under the following regularity assumptions: then we have: for any 0 ≤  ≤   − 2, where the constant  > 0 is independent of ℎ and  . Proof.
in (49) and  ℎ =  +1 ℎ in (51), then adding the results together, we have where we use the  2 projection property (︁ )︁ = 0, and the identity , respectively, we have and where we dropped the following two zero terms by the  2 projection property: )︁ = 0.
Choosing ̃︀  ℎ =  −2 0  0 ̃︀  +1 ℎ in (52), we have Adding ( 56)-( 59) together, we have Choosing ̃︀  = ̃︀ in (50), we obtain Adding ( 60) and (61) together, then summing up the result from  = 0 to , we obtain Now we need to bound those right hand side terms of (62).First, similar to (36) and (37) in the stability analysis, respectively, we have and Moving the second right hand side term of (62) to combine with some left hand side terms, we have Now we need to estimate the rest right hand side items in (62).Using the arithmetic-geometric mean inequality , Lemma 1 and the interpolation error estimate, we have Let us introduce the wave propagation speed notation   = 1 √ 00 .Using the fact that integration by parts and Lemma 1, we have Using integration by parts, Lemma 1 and the interpolation error estimate, we obtain where in the last step we used the standard inverse estimate ( [22], Chap.3): ‖∇ × ̃︀ By Lemma 1 and the interpolation error estimate, we have Similar to the analysis of Est 3 , we easily obtain Similarly, we can obtain and 3  4 By Lemma 1 and the projection error estimate, we have Similar to Est 10 , we obtain Finally, by Lemma 1 and the interpolation error estimate, we have and Substituting the estimates (63)-( 78) into (62), and combining many like terms, we obtain Now first choosing those   small enough, then choosing  small enough (but independent of mesh size ℎ), , and using the discrete Gronwall inequality and the fact that ( + 1) ≤  , we have where in the last step, we used the initial approximations (24a)-(24d).Then by the triangular inequality, and the interpolation and projection error estimates, we immediately have which completes the proof of (55).

Numerical results
In this section, we present some numerical results to demonstrate the performance of our proposed leapfrog scheme.For simplity, we only implement the lowest order RTN mixed finite element spaces ( = 1) on triangular elements: where   are the barycentric coordinate functions.
To construct an analytical solution, we add extra source terms to the model equations (1a)-(1e).More specifically, we solve the following governing equations: with exact solutions given as follows: )︂]︂ .
In this case, we have the following leapfrog type scheme: given initial approximations  0 ℎ , ̃︀  0 ℎ , We implement this scheme with different time step sizes  = 1 2 ℎ,  = ℎ,  = 2ℎ,  = 4ℎ and varying mesh sizes ℎ from 1  24 to 1 384 .The obtained convergence rates of the  2 errors at final time  = 1 are presented in Table 1, which shows that (ℎ) convergence can be obtained for both the electric field  and the magnetic field  without satisfying the CFL constraint.To test the time convergence rate, we solve this example again by fixing  = √ ℎ and varying mesh sizes ℎ from 1  4 to 1 256 .The convergence rates of  2 errors obtained at  = 1 are presented in Table 2, which justifies the theoretical convergence rate ( 2 + ℎ) for the lowest order RTN spaces.Example 2. This example is used to test the long time stability of our scheme and the wave absorbing capability of this equivalent Bérenger's PML model.For this example, we choose the physical domain Ω = [0, 0.5] m×[0, 0.5] m, which is partitioned by a uniform triangular mesh with mesh size ℎ = 2.5 × 10 −3 m.We surround the physical domain by 20-layer PML cells with thickness  = 20 h.In our simulation, the damping function   is chosen as a fourth-order polynomial function given as: This example is solved by the scheme (23a)-(23e) with zero initial fields and a point source wave located at (0.25, 0.25), the center of the domain.The source wave is imposed on the  field given as  = 0.1 sin(2 ) with frequency  = 3 GHz.
Snapshots of the computed magnetic field  obtained with time step size  = 2.5 × 10 −12 s and up to 10 000 time steps are plotted in Figure 1, which shows a long time stability of the scheme without obvious reflection from the surrounding PML cells.
To see more clearly the performance of the PML model, we solve this example again by stopping the source wave after 200 time steps so that we can see how large the residual wave magnitude can be.Some snapshots of the magnetic fields  are plotted in Figure 2, which shows that the original source wave exits the domain without obvious reflections.The magnitude of the residual wave after 1500 time steps is about 2 × 10 −4 , which is basically the numerical scheme error.Example 3.This example is used to show the wave absorbing performance of our equivalent Bérenger's PML model with a line source wave.For this example, we choose Ω = [0, 2] m×[0, 2] m, which is surrounded by 8-layer PML cells with thickness  = 8 h.The incident source is imposed on the  field given as  = sin(2 ) with frequency  = 1.5 GHz.To make a line source wave, the wave is placed on a line segment located at  = 0.1 m with  ranging from  = 0.5 m to  = 1.5 m.We use ℎ = 0.02 m and  = 10 −12 s for this simulation.Some snapshots of || up to 20 000 time steps are presented in Figure 3, which shows that out scheme enjoys a long stability and the wave propagates in the free space without obvious wave reflection from the truncated PML layers.To measure the PML absorption performance, we solve the same problem in a very larger domain Ω  = [−0.125,0.625] 2 m, which is imposed by the PEC boundary condition and is discretized by a finite element mesh with the same mesh size in Ω  .Then we calculate the errors of magnetic field  at element centers inside Ω  by subtracting the corresponding solutions from those obtained on Ω  .The global error energy is defined by the sum of the squares of those errors at element centers in Ω  .To satisfy the CFL condition, we choose an initial time step  0 = 6.25 × 10 −13 s and denote one time node (shortened as TN in Tab. 3) as 10 steps of  0 , i.e., 1 time node = 10  .We carried out three experiments with PML thicknesses of  = 5 h, 10 h and 20 h to observe the PML absorption capacity.
In Figure 4a, the whole domain plotted is Ω  , while the central square subdomain marked by the red color is Ω  .The Figure 4b compares the global errors versus time for the 5-cell, 10-cell, and 20-cell PMLs.As shown in Figure 4b, the global reflection error is decreasing with the increasing of the PML's thickness, which is consistent with the common performance of PML [34].For each fixed PML thickness, we tested two different mesh sizes: one with ℎ = 2.5 × 10 −3 m plotted by the solid line; and the other one with ℎ/2 = 1.25 × 10 −3 m plotted by the dashed line in Figure 4b.Due to the dominance of the numerical error, we did not achieve the ideal reflection error 10 −7 as chosen in  max .But we did observe the convergence rate ( 2 + ℎ 2 ) in the discrete  2 norm in Table 3.The (ℎ 2 ) is a superconvergence phenomenon often happend for the lowest-order edge elements ( [22], Chap.5).

Conclusion
In this paper, we developed a novel explicit unconditionally stable finite element scheme to solve an equivalent Bérenger's TEz PML model.We rigorously established both the stability and convergence analysis for the proposed scheme.Numerical results are presented to justify the theoretical analysis.We also demonstrated the effectiveness of this PML in simulating wave propagation in the free space.In the future, we will explore the possibility of extending similar idea to develop other explicit unconditionally stable schemes for other PML models [10].
otherwise, where  max = − log() * 5 *  0 *   /(2 * ) with  = 10 −7 and   being the wave propagation speed in vacuum.The damping function   has exactly the same form but in  variable.

Example 4 .
This example is used to illustrate the dependence of the PML absorption capacity on the damping functions, especially on the PML thickness.By adopting a popular numerical strategy (e.g.,[34]), we impose the same source wave given as Example 2 at the center of the physical domain Ω  = [0.125,0.375] 2 m, which is discretized by a 100×100 cells.We surround the domain Ω  by the equivalent Berenger's PML, whose damping functions and related parameters are the same as Example 2.

Figure 4 .
Figure 4. Example 4: The global errors obtained on a 100×100 mesh with a point wave source.(a) Illustration of domains Ω  and Ω  .(b) The global errors obtained with three different PML thicknesses.

Table 3 .
The discrete  2 errors for   of all points in Ω  obtained with different layers PML at several time nodes.