Notes on Perfectly Matched Layers (PMLs)

This note is intended as a brief introduction to the theory and practice of perfectly matched layer (PML) absorbing boundaries for wave equations, originally developed for MIT courses 18.369 and 18.336. It focuses on the complex stretched-coordinate viewpoint, and also discusses the limitations of PML.

Authors

• 7 publications
• Complex scaled infinite elements for exterior Helmholtz problems

The technique of complex scaling for time harmonic wave type equations r...
07/23/2019 ∙ by Lothar Nannen, et al. ∙ 0

• Preprint Déjà Vu: an FAQ

I give a brief overview of arXiv history, and describe the current state...
06/13/2017 ∙ by P. Ginsparg, et al. ∙ 0

• Application of a Generalized Secant Method to Nonlinear Equations with Complex Roots

The secant method is a very effective numerical procedure used for solvi...
05/28/2021 ∙ by Avram Sidi, et al. ∙ 0

• Perfectly Matched Layers for nonlocal Helmholtz equations Part II: higher dimensions

Perfectly matched layers (PMLs) are formulated and numerically applied t...
12/03/2020 ∙ by Yu Du, et al. ∙ 0

• Adaptive Fibonacci and Pairing Heaps

This brief note presents two adaptive heap data structures and conjectur...
03/09/2020 ∙ by Andrew Frohmader, et al. ∙ 0

• Generalized Approach to Matched Filtering using Neural Networks

Gravitational wave science is a pioneering field with rapidly evolving d...
04/08/2021 ∙ by Jingkai Yan, et al. ∙ 6

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Whenever one solves a PDE numerically by a volume discretization,111As opposed to a boundary discretization, e.g. in boundary-element methods, where the unknowns are on the interfaces between homogeneous regions, and the homogeneous regions are handled analtyically. In this case, no artificial truncation is required…except in the case of interfaces that extend to infinity, which lead to some interesting unsolved problems in boundary-element methods. one must truncate the computational grid in some way, and the key question is how to perform this truncation without introducing significant artifacts into the computation. Some problems are naturally truncated, e.g. for periodic structures where periodic boundary conditions can be applied. Some problems involve solutions that are rapidly decaying in space, so that the truncation is irrelevant as long as the computational grid is large enough. Other problems, such as Poisson’s equation, involve solutions that vary more and more slowly at greater distances—in this case, one can simply employ a coordinate transformation, such as , to remap from to , and solve the new finite system. However, some of the most difficult problems to truncate involve wave equations, where the solutions are oscillating and typically decay with distance only as in dimensions.222

The square of the solutions are typically related to a rate of energy transport, e.g. the Poynting vector in electromagnetism, and energy conservation requires that this decay be proportional to the surface area

. The slow decay means that simply truncating the grid with hard-wall (Dirichlet or Neumann) or periodic boundary conditions will lead to unacceptable artifacts from boundary reflections. The oscillation means that any real coordinate remapping from an infinite to a finite domain will result in solutions that oscillate infinitely fast as the boundary is approached—such fast oscillations cannot be represented by any finite-resolution grid, and will instead effectively form a reflecting hard wall. Therefore, wave equations require something different: an absorbing boundary that will somehow absorb waves that strike it, without reflecting them, and without requiring infeasible resolution.

The first attempts at such absorbing boundaries for wave equations involved absorbing boundary conditions (ABCs) [1]. Given a solution on a discrete grid, a boundary condition is a rule to set the value at the edge of the grid. For example, a simple Dirichlet boundary condition sets the solution to zero at the edge of the grid (which will reflect waves that hit the edge). An ABC tries to somehow extrapolate from the interior grid points to the edge grid point(s), to fool the solution into “thinking” that it extends forever with no boundary. It turns out that this is possible to do perfectly in one dimension, where waves can only propagate in two directions (). However, the main interest for numerical simulation lies in two and three dimensions, and in these cases the infinite number of possible propagation directions makes the ABC problem much harder. It seems unlikely that there exists any efficient method to exactly absorb radiating waves that strike a boundary at any possible angle. Existing ABCs restrict themselves to absorbing waves exactly only at a few angles, especially at normal incidence: as the size of the computational grid grows, eventually normal-incident waves must become the dominant portion of the radiation striking the boundaries. Another difficulty is that, in many practical circumstances, the wave medium is not homogeneous at the grid boundaries. For example, to calculate the transmission around a dielectric waveguide bend, the waveguide (an inhomogeneous region with a higher index of refraction) should in principle extend to infinity before and after the bend. Many standard ABCs are formulated only for homogeneous materials at the boundaries, and may even become numerically unstable if the grid boundaries are inhomogeneous.

In 1994, however, the problem of absorbing boundaries for wave equations was transformed in a seminal paper by Berenger [2]. Berenger changed the question: instead of finding an absorbing boundary condition, he found an absorbing boundary layer, as depicted in Fig. 1. An absorbing boundary layer is a layer of artificial absorbing material that is placed adjacent to the edges of the grid, completely independent of the boundary condition. When a wave enters the absorbing layer, it is attenuated by the absorption and decays exponentially; even if it reflects off the boundary, the returning wave after one round trip through the absorbing layer is exponentially tiny. The problem with this approach is that, whenever you have a transition from one material to another,333Technically, reflections occur when translational symmetry is broken. In a periodic structure (discrete translational symmetry), there are waves that propagate without scattering, and a uniform medium is just a special case with period . waves generally reflect, and the transition from non-absorbing to absorbing material is no exception—so, instead of having reflections from the grid boundary, you now have reflections from the absorber boundary. However, Berenger showed that a special absorbing medium could be constructed so that waves do not reflect at the interface: a perfectly matched layer, or PML. Although PML was originally derived for electromagnetism (Maxwell’s equations), the same ideas are immediately applicable to other wave equations.

There are several equivalent formulations of PML. Berenger’s original formulation is called the split-field PML, because he artificially split the wave solutions into the sum of two new artificial field components. Nowadays, a more common formulation is the uniaxial PML or UPML, which expresses the PML region as the ordinary wave equation with a combination of artificial anisotropic absorbing materials [3]. Both of these formulations were originally derived by laboriously computing the solution for a wave incident on the absorber interface at an arbitrary angle (and polarization, for vector waves), and then solving for the conditions in which the reflection is always zero. This technique, however, is labor-intensive to extend to other wave equations and other coordinate systems (e.g. cylindrical or spherical rather than Cartesian). It also misses an important fact: PML still works (can still be made theoretically reflectionless) for inhomogeneous media, such as waveguides, as long as the medium is homogeneous in the direction perpendicular to the boundary, even though the wave solutions for such media cannot generally be found analytically. It turns out, however, that both the split-field and UPML formulations can be derived in a much more elegant and general way, by viewing them as the result of a complex coordinate stretching [4, 5, 6].444It is sometimes implied that only the split-field PML can be derived via the stretched-coordinate approach [1], but the UPML media can be derived in this way as well [6]. It is this complex-coordinate approach, which is essentially based on analytic continuation of Maxwell’s equations into complex spatial coordinates where the fields are exponentially decaying, that we review in this note.

In the following, we first briefly remind the reader what a wave equation is, focusing on the simple case of the scalar wave equation but also giving a general definition. We then derive PML as a combination of two steps: analytic continuation into complex coordinates, then a coordinate transformation back to real coordinates. Finally, we discuss some limitations of PML, most notably the fact that it is no longer reflectionless once the wave equation is discretized, and common workarounds for these limitations.

2 Wave equations

There are many formulations of waves and wave equations in the physical sciences. The prototypical example is the (source-free) scalar wave equation:

 ∇⋅(a∇u)=1b∂2u∂t2=¨ub (1)

where is the scalar wave amplitude and is the phase velocity of the wave for some parameters and of the (possibly inhomogeneous) medium. For lossless, propagating waves, and should be real and positive.

Both for computational convenience (in order to use a staggered-grid leap-frog discretization) and for analytical purposes, it is more convenient to split this second-order equation into two coupled first-order equation, by introducing an auxiliary field :

 ∂u∂t = b∇⋅v, (2)
 ∂v∂t=a∇u, (3)

which are easily seen to be equivalent to eq. (1).

Equations (23) can be written more abstractly as:

 (4)

for a linear operator and a 4-component vector (in three dimensions). The key property that makes this a “wave equation” turns out to be that is an anti-Hermitian operator in a proper choice of inner product, which leads to oscillating solutions, conservation of energy, and other “wave-like” phenomena. Every common wave equation, from scalar waves to Maxwell’s equations (electromagnetism) to Schrödinger’s equation (quantum mechanics) to the Lamé-Navier equations for elastic waves in solids, can be written in the abstract form for some wave function and some anti-Hermitian operator .555See e.g. Ref. [7] The same PML ideas apply equally well in all of these cases, although PML is most commonly applied to Maxwell’s equations for computational electromagnetism.

3 Complex coordinate stretching

Let us start with the solution of some wave equation in infinite space, in a situation similar to that in Fig. 1(a): we have some region of interest near the origin , and we want to truncate space outside the region of interest in such a way as to absorb radiating waves. In particular, we will focus on truncating the problem in the direction (the other directions will follow by the same technique). This truncation occurs in three conceptual steps, summarized as follows:

1. In infinite space, analytically continue the solutions and equations to a complex contour, which changes oscillating waves into exponentially decaying waves outside the region of interest without reflections.

2. Still in infinite space, perform a coordinate transformation to express the complex as a function of a real coordinate. In the new coordinates, we have real coordinates and complex materials.

3. Truncate the domain of this new real coordinate inside the complex-material region: since the solution is decaying there, as long as we truncate it after a long enough distance (where the exponential tails are small), it won’t matter what boundary condition we use (hard-wall truncations are fine).

For now, we will make two simplifications:

• We will assume that the space far from the region of interest is homogeneous (deferring the inhomogeneous case until later).

• We will assume that the space far from the region of interest is linear and time-invariant.

Under these assumptions, the radiating solutions in infinite space must take the form of a superposition of planewaves:

 w(x,t)=∑k,ωWk,ωei(k⋅x−ωt), (5)

for some constant amplitudes , where is the (angular) frequency and is the wavevector. (In an isotropic medium, and are related by where is some phase velocity, but we don’t need to assume that here.) In particular, the key fact is that the radiating solutions may be decomposed into functions of the form

 W(y,z)ei(kx−ωt). (6)

The ratio is the phase velocity, which can be different from the group velocity (the velocity of energy transport, in lossless media). For waves propagating in the direction, the group velocity is positive. Except in very unusual cases, the phase velocity has the same sign as the group velocity in a homogeneous medium,666The formulation of PML absorbers when the phase velocity has sign opposite to the group velocity, for example in the “left-handed media” of electromagnetism, is somewhat more tricky [8, 9]; matters are even worse in certain waveguides with both signs of group velocity at the same  [10]. so we will assume that is positive.

3.1 Analytic continuation

The key fact about eq. (6) is that it is an analytic function of . That means that we can freely analytically continue it, evaluating the solution at complex values of . The original wave problem corresponds to along the real axis, as shown in the top panels of Fig. 2, which gives an oscillating solution. However, if instead of evaluating along real axis, consider what happens if we evaluate it along the contour shown in the bottom-left panel of Fig. 2, where for we have added a linearly growing imaginary part. In this case, the solution is exponentially decaying for , because is exponentially decaying (for ) as increases. That is, the solution in this region acts like the solution in an absorbing material.

However, there is one crucial difference here from an ordinary absorbing material: the solution is not changed for , where is no different from before. So, it not only acts like an absorbing material, it acts like a reflectionless absorbing material, a PML.

The thing to remember about this is that the analytically continued solution satisfies the same differential equation. We assumed the differential equation was -invariant in this region, so only appeared in derivatives like , and the derivative of an analytic function is the same along any direction in the complex plane. So, we have succeeded in transforming our original wave equation to one in which the radiating solutions (large |) are exponentially decaying, while the part we care about (small ) is unchanged. The only problem is that solving differential equations along contours in the complex plane is rather unfamiliar and inconvenient. This difficulty is easily fixed.

3.2 Coordinate transformation back to real x

For convenience, let’s denote the complex contour by , and reserve the letter for the real part. Thus, we have , where is some function indicating how we’ve deformed our contour along the imaginary axis. Since the complex coordinate is inconvenient, we will just change variables to write the equations in terms of , the real part!

Changing variables is easy. Whereever our original equation has (the differential along the deformed contour ), we now have . That’s it! Since our original wave equation was assumed -invariant (at least in the large- regions where ), we have no other substitutions to make. As we shall soon see, it will be convenient to denote , for some function . [For example, in the bottom panel of Fig. 2, we chose to be a step function: zero for and a positive constant for , which gave us exponential decay.] In terms of , the entire process of PML can be conceptually summed up by a single transformation of our original differential equation:

 ∂∂x→11+iσx(x)ω∂∂x. (7)

In the PML regions where , our oscillating solutions turn into exponentially decaying ones. In the regions where , our wave equation is unchanged and the solution is unchanged: there are no reflections because this is only an analytic continuation of the original solution from to , and where the solution cannot change.

Why did we choose , as opposed to just ? The answer comes if we look at what happens to our wave . In the new coordinates, we get:

 eikxe−kω∫xσx(x′)dx′. (8)

Notice the factor , which is equal to , the inverse of the phase velocity in the direction. In a dispersionless material (e.g. vacuum for light), is a constant independent of velocity for a fixed angle, in which case the attenuation rate in the PML is independent of frequency : all wavelengths decay at the same rate! In contrast, if we had left out the then shorter wavelengths would decay faster than longer wavelengths. On the other hand, the attenuation rate is not independent of the angle of the light, a difficulty discussed in Sec. 7.2.

3.3 Truncating the computational region

Once we have performed the PML transformation (7) of our wave equations, the solutions are unchanged in our region of interest (small ) and exponentially decaying in the outer regions (large ). That means that we can truncate the computational region at some sufficiently large , perhaps by putting a hard wall (Dirichlet boundary condition). Because only the tiny exponential tails “see” this hard wall and reflect off it, and even they attenuate on the way back towards the region of interest, the effect on the solutions in our region of interest will be exponentially small.

In fact, in principle we can make the PML region as thin as we want, just by making very large (which makes the exponential decay rate rapid), thanks to the fact that the decay rate is independent of (although the angle dependence can be a problem, as discussed in Sec. 7.2). However, in practice, we will see in Sec. 7.1 that using a very large can cause “numerical reflections” once we discretize the problem onto a grid. Instead, we turn on quadratically or cubically from zero, over a region of length a half-wavelength or so, and in practice the reflections will be tiny.

3.4 PML boundaries in other directions

So far, we’ve seen how to truncate our computational region with a PML layer in the direction. What about other directions? The most important case to consider is the direction. The key is, in the direction we do exactly the same thing: apply the PML transformation (7) with at a sufficiently large negative , and then truncate the computational cell. This works because, for , the radiating waves are propagating in the direction with (negative phase velocity), and this makes our PML solutions (8) decay in the opposite direction (exponentially decaying as ) for the same positive .

Now that we have dealt with , the and directions are easy: just do the same transformation, except to and , respectively, using functions and that are non-zero in the and PML regions. At the corners of the computational cell, we will have regions that are PML along two or three directions simultaneously (i.e. two or three ’s are nonzero), but that doesn’t change anything.

3.5 Coordinate transformations and materials

We will see below that, in the context of the scalar wave equation, the term from the PML coordinate transformation appears as, effectively, an artificial anisotropic absorbing material in the wave equation (effectively changing and

to complex numbers, and a tensor in the case of

). At least in the case of Maxwell’s equations (electromagnetism), this is an instance of a more general theorem: Maxwell’s equations under any coordinate transformation can be expressed as Maxwell’s equations in Cartesian coordinates with transformed materials.777This theorem appears to have been first clearly stated and derived by Ward and Pendry [11], and is summarized in a compact general form by my course notes [12] and in our paper [13]. That is, the coordinate transform is “absorbed” into a change of the permittivity and the permeability (generally into anisotropic tensors). This is the reason why UPML, which constructs reflectionless anisotropic absorbers, is equivalent to a complex coordinate stretching: it is just absorbing the coordinate stretching into the material tensors.

4 PML examples in frequency and time domain

As we have seen, in frequency domain, when we are solving for solutions with time-dependence , PML is almost trivial: we just apply the PML transformation (7) to every derivative in our wave equation. (And similarly for derivatives in other directions, to obtain PML boundaries in different directions.)

In the time domain, however, things are a bit more complicated, because we chose our transformation to be of the form : our complex “stretch” factor is frequency-dependent in order that the attenuation rate be frequency-independent. But how do we express a dependence in the time domain, where we don’t have (i.e. the time-domain wave function may superimpose multiple frequencies at once)? One solution is to punt, of course, and just use a stretch factor for some constant frequency that is typical of our problem; as long as our bandwidth is narrow enough, our attenuation rate (and thus the truncation error) will be fairly constant. However, it turns out that there is a way to implement the ideal dependence directly in the time domain, via the auxiliary differential equation (ADE) approach.

This approach is best illustrated by example, so we will consider PML boundaries in the direction for the scalar wave equation in one and two dimensions. (It turns out that an ADE is not required in 1d, however.)

4.1 An example: PML for 1d scalar waves

Let’s consider the 1d version of the scalar wave equation (23):

 ∂u∂t=b∂v∂x=−iωu
 ∂v∂t=a∂u∂x=−iωv,

where we have substituted an time-dependence. Now, if we perform the PML transformation (7), and multiply both sides by , we obtain:

 b∂v∂x=−iωu+σxu
 a∂u∂x=−iωv+σxv.

The terms have cancelled, and so in this 1d case we can trivially turn the equations back into their time-domain forms:

 ∂u∂t=b∂v∂x−σxu
 ∂v∂t=a∂u∂x−σxv.

Notice that, for , the decay terms have exactly the right sign to make the solutions decay in time if and were constants in space. Similarly, they have the right sign to make it decay in space whereever . But this is a true PML: there are zero reflections from any boundary where we change , even if we change discontinuously (not including the discretization problems mentioned above).

By the way, the above equations reveal why we use the letter for the PML absorption coefficient. If the above equations are interpreted as the equations for electric () and magnetic () fields in 1d electromagnetism, then plays the role of a conductivity, and conductivity is traditionally denoted by . Unlike the usual electrical conductivity, however, in PML we have both an electric and a magnetic conductivity, since we have terms corresponding to currents of electric and magnetic charges. There is no reason we need to be limited to physical materials to construct our PML for a computer simulation!

4.2 An example: PML for 2d scalar waves

Unfortunately, the 1d case above is a little too trivial to give you the full flavor of how PML works. So, let’s go to a 2d scalar wave equation (again for time-dependence):

 ∂u∂t=b∇⋅v=b∂vx∂x+b∂vy∂y=−iωu
 ∂vx∂t=a∂u∂x=−iωvx
 ∂vy∂t=a∂u∂y=−iωvy.

Again performing the PML transformation (7) of in the first two equations, and multiplying both sides by , we obtain:

 b∂vx∂x+b∂vy∂y(1+iσxω)=−iωu+σxu
 a∂u∂x=−iωvx+σxvx.

The the second equation is easy to transform back to time domain, just like for the 1d scalar-wave equation: becomes a time derivative. The first equation, however, poses a problem: we have an extra term with an explicit factor. What do we do with this?

In a Fourier transform,

corresponds to differentiation, so corresponds to integration: our problematic term is the integral of another quantity. In particular, let’s introduce a new auxiliary field variable , satisfying

 −iωψ=bσx∂vy∂y,

in which case

 b∂vx∂x+b∂vy∂y+ψ=−iωu+σxu.

Now, we can Fourier transform everything back to the time-domain, to get a set of four time-domain equations with PML absorbing boundaries in the direction that we can solve by our favorite discretization scheme:

 ∂u∂t=b∇⋅v−σxu+ψ
 ∂vx∂t=a∂u∂x−σxvx
 ∂vy∂t=a∂u∂y
 ∂ψ∂t=bσx∂vy∂y,

where the last equation for is our auxiliary differential equation (with initial condition ). Notice that we have absorption terms in the and equation, but not for : the PML corresponds to an anisotropic absorber, as if were replaced by the complex tensor

 (1a+iσxωa1a)−1.

This is an example of the general theorem alluded to in Sec. 3.5 above.

5 PML in inhomogeneous media

The derivation above didn’t really depend at all on the assumption that the medium was homogenous in for the PML layer. We only assumed that the medium (and hence the wave equation) was invariant in the direction for sufficiently large . For example, instead of empty space we could have a waveguide oriented in the direction (i.e. some -invariant cross-section). Regardless of the dependence, translational invariance implies that radiating solutions can be decomposed into a sum of functions of the form of eq. (6), . These solutions are no longer plane waves. Instead, they are the normal modes of the -invariant structure, and is the propagation constant. These normal modes are the subject of waveguide theory in electromagnetism, a subject extensively treated elsewhere [14, 15]. The bottom line is: since the solution/equation is still analytic in , the PML is still reflectionless.888There is a subtlety here because, in unusual cases, uniform waveguides can support “backward-wave” modes where the phase and group velocities are opposite, i.e. for a right-traveling wave [16, 17, 18, 19]. It has problems even worse than those reported for left-handed media [8, 9], because the same frequency has both “right-handed” and “left-handed” modes; a deeper analysis of this interesting case can be found in our paper [10].

6 PML for evanescent waves

In the discussion above, we considered waves of the form and showed that they became exponentially decaying if we replace by for . However, this discussion assumed that was real (and positive). This is not necessarily the case! In two or more dimensions, the wave equation can have evanescent solutions where is complex, most commonly where is purely imaginary. For example consider a planewave in a homogeneous two-dimensional medium with phase velocity , i.e. . In this case,

 k=kx=√ω2c2−k2y.

For sufficiently large (i.e. high-frequency Fourier components in the direction), is purely imaginary. As we go to large , the boundary condition at implies that we must have so that is exponentially decaying.

What happens to such a decaying, imaginary- evanescent wave in the PML medium? Let . Then, in the PML:

 e−κx→e−κx−iσxωx. (9)

That is, the PML added an oscillation to the evanescent wave, but did not increase its decay rate. The PML is still reflectionless, but it didn’t help.

Of course, you might object that an evanescent wave is decaying anyway, so we hardly need a PML—we just need to make the computational region large enough and it will go away on its own. This is true, but it would be nice to accelerate the process: in some cases may be relatively small and we would need a large grid for it to decay sufficiently. This is no problem, however, because nothing in our analysis required to be real. We can just as easily make complex, where corresponds to a real coordinate stretching. That is, the imaginary part of will accelerate the decay of evanescent waves in eq. (9) above, without creating any reflections.

Adding an imaginary part to does come at a price, however. What it does to the propagating (real ) waves is to make them oscillate faster, which exacerbates the numerical reflections described in Sec. 7.1. In short, everything in moderation.

7 Limitations of PML

PML, while it has revolutionized absorbing boundaries for wave equations, especially (but not limited to) electromagnetism, is not a panacea. Some of the limitations and failure cases of PML are discussed in this section, along with workarounds.

7.1 Discretization and numerical reflections

First, and most famously, PML is only reflectionless if you are solving the exact wave equations. As soon as you discretize the problem (whether for finite difference or finite elements), you are only solving an approximate wave equation and the analytical perfection of PML is no longer valid.

What is left, once you discretize? PML is still an absorbing material: waves that propagate within it are still attenuated, even discrete waves. The boundary between the PML and the regular medium is no longer reflectionless, but the reflections are small because the discretization is (presumably) a good approximation for the exact wave equation. How can we make the reflections smaller, as small as we want?

The key fact is that, even without a PML, reflections can be made arbitrarily small as long as the medium is slowly varying. That is, in the limit as you “turn on” absorption more and more slowly, reflections go to zero due to an adiabatic theorem [20]. With a non-PML absorber, you might need to go very slowly (i.e. a very thick absorbing layer) to get acceptable reflections [21]. With PML, however, the constant factor is very good to start with, so experience shows that a simple quadratic or cubic turn-on of the PML absorption usually produces negligible reflections for a PML layer of only half a wavelength or thinner [1, 21]. (Increasing the resolution also increases the effectiveness of the PML, because it approaches the exact wave equation.)

7.2 Angle-dependent absorption

Another problem is that the PML absorption depends on angle. In particular, consider eq. (8) for the exponential attenuation of waves in the PML, and notice that the attenuation rate is proportional to the ratio . But , here, is really just , the component of the wavevector in the direction (for a planewave in a homogeneous medium). Thus, the attenuation rate is proportional to , where is the angle the radiating wave makes with the axis. As the radiation approaches glancing incidence, the attenuation rate goes to zero! This means that, for any fixed PML thickness, waves sufficiently close to glancing incidence will have substantial “round-trip” reflections through the PML.

In practice, this is not as much of a problem as it may sound like at first. In most cases, all of the radiation originates in a localized region of interest near the origin, as in Fig. 1. In this situation, all of the radiation striking the PML will be at an angle in the limit as the boundaries move farther and farther away (assuming a cubic computational region). So, if the boundaries are far enough away, we can guarantee a maximum angle and hence make the PML thick enough to sufficiently absorb all waves within this cone of angles.

7.3 Inhomogeneous media where PML fails

Finally, PML fails completely in the case where the medium is not -invariant (for an boundary) [21]. You might ask: why should we care about such cases, as if the medium is varying in the direction then we will surely get reflections (from the variation) anyway, PML or no PML? Not necessarily.

There are several important cases of -varying media that, in the infinite system, have reflectionless propagating waves. Perhaps the simplest is a waveguide that hits the boundary of the computational cell at an angle (not normal to the boundary)—one can usually arrange for all waveguides to leave the computational region at right angles, but not always (e.g. what if you want the transmission through a bend?). Another, more complicated and perhaps more challenging case is that of a photonic crystal: for a periodic medium, there are wave solutions (Bloch waves) that propagate without scattering, and can have very interesting properties that are unattainable in a physical uniform medium [22].

For any such case, PML seems to be irrevocably spoiled. The central idea behind PML was that the wave equations, and solutions, were analytic functions in the direction perpendicular to the boundary, and so they could be analytically continued into the complex coordinate plane. If the medium is varying in the direction, it is most likely varying discontinously, and hence the whole idea of analytic continuation goes out the window.

What can we do in such a case? Conventional ABCs don’t work either (they are typically designed for homogeneous media). The only fallback is the adiabatic theorem alluded to above: even a non-PML absorber, if turned on gradually enough and smoothly enough, will approach a reflectionless limit. The difficulty becomes how gradual is gradual enough, and in finding a way to make the non-PML absorber a tractable thickness [21].

There is also another interesting case where PML fails. The basic analytic-continuation idea is valid in any -invariant medium, regardless of inhomogeneities in the plane. However, certain inhomogeneous dielectric patterns in the plane give rise to unusual solutions: “left-handed” solutions where the phase velocity is opposite to the group velocity in the direction, while at the same the medium also has “right-handed” solutions where the phase velocity and group velocity have the same sign. Most famously, this occurs for certain “backward-wave” coaxial waveguides [16, 17, 18, 19]. In this case, whatever sign one picks for the PML conductivity , either the left- or right-handed modes will be exponentially growing and the PML fails in a spectacular instability [10]. There is a subtle relationship of this failure to the orientation of the fields for a left-handed mode and the anisotropy of the PML [10]. In this case, one must once again abandon PML absorbers and use a different technique, such as a scalar conductivity that is turned on sufficiently gradually to adiabatically absorb outgoing waves.

References

• [1] A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method. Norwood, MA: Artech, 2000.
• [2] J.-P. Bérenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys., vol. 114, no. 1, pp. 185–200, 1994.
• [3] Z. S. Sacks, D. M. Kingsland, R. Lee, and J. F. Lee, “A perfectly matched anisotropic absorber for use as an absorbing boundary condition,” IEEE Trans. Antennas and Propagation, vol. 43, no. 12, pp. 1460–1463, 1995.
• [4] W. C. Chew and W. H. Weedon, “A 3d perfectly matched medium from modified Maxwell’s equations with stretched coordinates,” Microwave and Optical Tech. Lett., vol. 7, no. 13, pp. 599–604, 1994.
• [5] C. M. Rappaport, “Perfectly matched absorbing boundary conditions based on anisotropic lossy mapping of space,” IEEE Microwave and Guided Wave Lett., vol. 5, no. 3, pp. 90–92, 1995.
• [6] F. L. Teixeira and W. C. Chew, “General closed-form PML constitutive tensors to match arbitrary bianisotropic and dispersive linear media,” IEEE Microwave and Guided Wave Lett., vol. 8, no. 6, pp. 223–225, 1998.
• [7] S. G. Johnson, “Notes on the algebraic structure of wave equations.” Online at http://math.mit.edu/~stevenj/18.369/wave-equations.pdf.
• [8] S. A. Cummer, “Perfectly matched layer behavior in negative refractive index materials,” IEEE Antennas and Wireless Propagation Lett., vol. 3, pp. 172–175, 2004.
• [9] X. T. Dong, X. S. Rao, Y. B. Gan, B. Guo, and W. Y. Yin, “Perfectly matched layer-absorbing boundary condition for left-handed materials,” IEEE Microwave and Wireless Components Lett., vol. 14, no. 6, pp. 301–303, 2004.
• [10] P.-R. Loh, A. F. Oskooi, M. Ibanescu, M. Skorobogatiy, and S. G. Johnson, “Fundamental relation between phase and group velocity, and application to the failure of perfectly matched layers in backward-wave structures,” Physical Review E, vol. 79, p. 065601(R), 2009.
• [11] A. J. Ward and J. B. Pendry, “Refraction and geometry in Maxwell’s equations,” J. Mod. Opt., vol. 43, no. 4, pp. 773–793, 1996.
• [12] S. G. Johnson, “Coordinate transformation and invariance in electromagnetism: Notes for the course 18.369 at mit.” Online at http://math.mit.edu/~stevenj/18.369/coordinate-transform.pdf.
• [13] C. Kottke, A. Farjadpour, and S. G. Johnson, “Perturbation theory for anisotropic dielectric interfaces, and application to sub-pixel smoothing of discretized numerical methods,” Physical Review E, vol. 77, p. 036611, 2008.
• [14] A. W. Snyder, “Radiation losses due to variations of radius on dielectric or optical fibers,” IEEE Trans. Microwave Theory Tech., vol. MTT-18, no. 9, pp. 608–615, 1970.
• [15] D. Marcuse, Theory of Dielectric Optical Waveguides. San Diego: Academic Press, second ed., 1991.
• [16] P. J. B. Clarricoats and R. A. Waldron, “Non-periodic slow-wave and backward-wave structures,” J. Electron. Contr., vol. 8, pp. 455–458, 1960.
• [17] R. A. Waldron, “Theory and potential applications of backward waves in nonperiodic inhomogeneous waveguides,” Proc. IEE, vol. 111, pp. 1659–1667, 1964.
• [18] A. S. Omar and K. F. Schunemann, “Complex and backward-wave modes in inhomogeneously and anisotropically filled waveguides,” IEEE Trans. Microwave Theory Tech., vol. MTT-35, no. 3, pp. 268–275, 1987.
• [19] M. Ibanescu, S. G. Johnson, D. Roundy, C. Luo, Y. Fink, and J. D. Joannopoulos, “Anomalous dispersion relations by symmetry breaking in axially uniform waveguides,” Phys. Rev. Lett., vol. 92, p. 063903, 2004.
• [20] S. G. Johnson, P. Bienstman, M. Skorobogatiy, M. Ibanescu, E. Lidorikis, and J. D. Joannopoulos, “Adiabatic theorem and continuous coupled-mode theory for efficient taper transitions in photonic crystals,” Phys. Rev. E, vol. 66, p. 066608, 2002.
• [21] A. F. Oskooi, L. Zhang, Y. Avniel, and S. G. Johnson, “The failure of perfectly matched layers, and towards their redemption by adiabatic absorbers,” Optics Express, vol. 16, pp. 11376–11392, July 2008.
• [22] J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals: Molding the Flow of Light. Princeton University Press, second ed., 2008.