OSA's Digital Library

Optics Express

Optics Express

  • Editor: C. Martijn de Sterke
  • Vol. 16, Iss. 15 — Jul. 21, 2008
  • pp: 11376–11392

Optics InfoBase > Optics Express > Volume 16 > Issue 15 > Page 11376

« Show journal navigation

The failure of perfectly matched layers, and towards their redemption by adiabatic absorbers

Ardavan F. Oskooi, Lei Zhang, Yehuda Avniel, and Steven G. Johnson  »View Author Affiliations


Optics Express, Vol. 16, Issue 15, pp. 11376-11392 (2008)
http://dx.doi.org/10.1364/OE.16.011376


View Full Text Article

Acrobat PDF (681 KB)


Advanced Search

Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

Although perfectly matched layers (PMLs) have been widely used to truncate numerical simulations of electromagnetism and other wave equations, we point out important cases in which a PML fails to be reflectionless even in the limit of infinite resolution. In particular, the underlying coordinate-stretching idea behind PML breaks down in photonic crystals and in other structures where the material is not an analytic function in the direction perpendicular to the boundary, leading to substantial reflections. The alternative is an adiabatic absorber, in which reflections are made negligible by gradually increasing the material absorption at the boundaries, similar to a common strategy to combat discretization reflections in PMLs. We demonstrate the fundamental connection between such reflections and the smoothness of the absorption profile via coupled-mode theory, and show how to obtain higher-order and even exponential vanishing of the reflection with absorber thickness (although further work remains in optimizing the constant factor).

© 2008 Optical Society of America

1. Introduction

A perfectly matched layer (PML) is an artificial absorbing medium that is commonly used to truncate computational grids for simulating wave equations (e.g. Maxwell’s equations), and is designed to have the property that interfaces between the PML and adjacent media are reflectionless in the exact wave equation [1

J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. 114,185–200 (1994). [CrossRef]

,2

A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech, 2000).

].We describe important cases in which PML fails to be reflectionless, even in the exact (non-discretized) Maxwell equations, most notably in the case of periodic media (photonic crystals [3

J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals: Molding the Flow of Light, 2nd Edition (Princeton Univ. Press, 2008).

])—contrary to previous suggestions of photonic-crystal “PML” absorbers [4–8

M. Koshiba, Y. Tsuji, and S. Sasaki, “High-performance absorbing boundary conditions for photonic crystal waveguide simulations,” IEEE Microwave Wirel. Compon. Lett. 11, 152–154 (2001). [CrossRef]

]. In these cases (similar to PML reflections due to discretization error [2

A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech, 2000).

,9

W. C. Chew and J. M. Jin, “Perfectly matched layers in the discretized space: An analysis and optimization,” Electromagnetics 16, 325–340 (1996). [CrossRef]

]), the remaining approach to reduce reflections is to “turn on” the absorption gradually, asymptotically approaching an “adiabatic” limit of zero reflections [10

S. G. Johnson, P. Bienstman, M. Skorobogatiy, M. Ibanescu, E. Lidorikis, and J. Joannopoulos, “Adiabatic theorem and continuous coupled-mode theory for efficient taper transitions in photonic crystals,”Phys. Rev. E 66, 066608 (2002).

] regardless of whether the medium forms a true PML—here, we provide a deeper understanding of all such adiabiatic absorbers by showing that the reflection’s dependence on the thickness of the absorbing layer is determined by the smoothness of the absorption profile, and can be predicted by coupled-mode theory approximations. For a fixed absorption profile (typically quadratic or cubic in previous work [2

A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech, 2000).

]), the reflection decreases with absorber thickness L proportional to some characteristic power law determined by the smoothness (e.g. 1/L 6 for quadratic absorption). As the absorber becomes thicker, smoother absorptions become favorable, and we show that it is even possible to obtain exponential decrease of the reflection with L by new choices of the absorption profile (although further work remains in optimizing the constant factor). The role of PML (when it works), compared to ordinary absorbing materials, is to improve the constant factor in this reflection convergence, rather than the functional form. For homogeneous materials as in most previous analyses, although some attempts have been made to optimize the PML profile among various polynomial functions [11–14

E. A. Marengo, C. M. Rappaport, and E. L. Miller, “Optimum PML ABC conductivity profile in FDFD,” IEEE Trans. Magn. 35, 1506–1509 (1999). [CrossRef]

], a quadratic or cubic profile works so well [2

A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech, 2000).

] that further attempts at optimization are arguably superfluous. On the other hand, for periodic media—especially when operating in modes with low group velocity—the required absorber thickness can become so large that the choice of absorption profile becomes critical. We also discuss the possibility of other optimizations, such as balancing the “transition” reflection from the absorber interface with the “round-trip” reflection due to the finite absorption, but these optimizations depend more sensitively on the incident-wave medium.

There are several nearly equivalent formulations of PML. Berenger’s original formulation [1

J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. 114,185–200 (1994). [CrossRef]

] split the wave solution into the sum of two new artificial field components. A more common “UPML” (uniaxial-PML) formulation expresses the PML region as the ordinary wave equation with a combination of artificial anisotropic absorbing materials [15

Z. 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 Propag. 43, 1460–1463 (1995). [CrossRef]

]. Both of these formulations were originally derived by laboriously computing the solution for a planewave incident on the absorber interface at an arbitrary angle and polarization, and then solving for the conditions in which the reflection is always zero. Both formulations, however, can also be derived by a complex “stretched-coordinate” approach [16–18

W. C. Chew and W. H. Weedon, “A 3d perfectly matched medium from modified Maxwell’s equations with stretched coordinates,” Microwave Opt. Technol. Lett. 7, 599–604 (1994). [CrossRef]

]—this much simpler and more elegant derivation of PML reveals its underlying meaning and generalizes more easily to inhomogeneous media, other wave equations, and other coordinate systems. In particular, the coordinate-stretching approach derives PML by an analytic continuation of Maxwell’s equations into complex spatial coordinates, where the oscillating fields become exponentially decaying [16–19

W. C. Chew and W. H. Weedon, “A 3d perfectly matched medium from modified Maxwell’s equations with stretched coordinates,” Microwave Opt. Technol. Lett. 7, 599–604 (1994). [CrossRef]

]. (This description can then be converted back into a change of materials via a complex coordinate transformation [18

F. L. Teixeira and W. C. Chew, “General close-form PML constitutive tensors to match arbitrary bianisotropic and dispersive linear media,” IEEE Microwave and Guided Wave Lett. 8, 223–225 (1998). [CrossRef]

, 20

A. J. Ward and J. B. Pendry, “Refraction and geometry in Maxwell’s equations,” J. Mod. Opt. 43, 773–793 (1996). [CrossRef]

].) By viewing PML as an analytic continuation, it can be shown to be reflectionless even for inhomogeneous media such as in Fig. 1(a) [21

W. Huang, C. Xu, W. Lui, and K. Yokoyama, “The perfectly matched layer boundary condition for modal analysis of optical waveguides: leaky mode calculations,” IEEE Photon. Tech. Lett. 8, 652–654 (1996). [CrossRef]

]: for a waveguide entering the PML perpendicularly, complex coordinate stretching is still possible because the material parameters (and hence Maxwell’s equations) are analytic functions (constants) in that direction. The same derivation of PML, however, also immediately points to situations where PML is inapplicable: in any problem where the material parameters are not described by analytic functions in the direction perpendicular to the boundary, a reflectionless absorber cannot be designed by complex coordinate stretching. As discussed in more detail below, this means that “PML” is not reflectionless for photonic crystals as in Fig. 1(b) where the dielectric function varies discontinously in the direction perpendicular to the boundary, or even in cases where a dielectric waveguide hits the PML obliquely [Fig. 1(c)]. (In fact, even for rare cases in which an oscillating dielectric function is analytic in the PML direction, we will explain that the analytic-continuation idea still does not yield a useful PML absorber in the discretized equations.)

Fig. 1. (a) PML is still reflectionless for inhomogeneous media such as waveguides that are homogeneous in the direction perpenendicular to the PML. (b, c) PML is no longer reflectionless when the dielectric function is discontinuous (non-analytic) in the direction perpendicular to the PML, as in a photonic crystal (b) or a waveguide entering the PML at an angle (c).

Previous suggestions to apply PML to photonic crystals by simply overlapping a “PML” anisotropic absorber with the periodic dielectric function [4–8

M. Koshiba, Y. Tsuji, and S. Sasaki, “High-performance absorbing boundary conditions for photonic crystal waveguide simulations,” IEEE Microwave Wirel. Compon. Lett. 11, 152–154 (2001). [CrossRef]

] (including a similar suggestion for integral-equation methods [22

D. Pissoort and F. Olyslager, “Termination of periodic waveguides by PMLs in time-harmonic integral equationlike techniques,” IEEE Antennas and Wireless Propagation Lett. 2, 281–284 (2003). [CrossRef]

]) were therefore not “true” PML media in the sense that the reflection will not go to zero even in the limit of infinite resolution. In this paper, we will refer to such an absorbing layer as a pseudo-PML (pPML). (In the special case of an effectively one-dimensional medium where there is only a single propagating mode, such as a single-mode waveguide surrounded by a complete-bandgap medium, it is possible to arrange an “impedancematched” absorber to approximately cancel that one mode [23

A. Mekis, S. Fan, and J. D. Joannopoulos, “Absorbing boundary conditions for FDTD simulations of photonic crystal waveguides,” IEEE Microwave and Guided Wave Lett. 9, 502–504 (1999). [CrossRef]

], or alternatively to specify analytical boundary conditions of zero reflection for that one mode [24

E. Moreno, D. Erni, and C. Hafner, “Modeling of discontinuities in photonic crystal waveguides with the multiple multipole method,” Phys. Rev. E 66, 036618 (2002).

]. More generally, in a transfer-matrix or scattering-matrix method where one explicitly computes all propagating modes and expands the fields in that basis, it is possible to impose analytically reflectionless boundary conditions [25

Z.-Y. Li and K.-M. Ho, “Light propagation in semi-infinite photonic crystals and related waveguide structures,” Phys. Rev. B 68 (2003).

, 26

D. Pissoort, B. Denecker, P. Bienstman, F. Olyslager, and D. De Zutter , “Comparative study of three methods for the simulation of two-dimensional photonic crystals,” J. Opt. Soc. Am. A 21, 2186–2195 (2004). [CrossRef]

], but such methods become very expensive in three dimensions.) These previous authors were nevertheless able to observe small reflections in a pPML only because they overlapped the pPML with many periods of the crystal and thereby turn on the pPML very gradually. As we explain below, such absorbing layers are more properly understood as adiabatic absorbers rather than PML media, and indeed the “PML” property only improves the constant factor in the long-wavelength limit of an effective homogeneous medium, or in any case where there are large homogeneous-material regions compared to the wavelength. Moreover, as we describe, the reflections worsen rapidly as the group velocity decreases (e.g., as a band edge is approached).

Even in the case of a homogeneous medium (or one uniform in the direction perpendicular to the boundary), where true PML applies, there are well-known numerical reflections due to the finite discretization [2

A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech, 2000).

, 9

W. C. Chew and J. M. Jin, “Perfectly matched layers in the discretized space: An analysis and optimization,” Electromagnetics 16, 325–340 (1996). [CrossRef]

]. It is sometimes claimed that the solutions for a PML converge exponentially to the solution of the open problem as the PML thickness is increased [27

M. Lassas and E. Somersalo, “On the existence and convergence of the solution of PML equations,” Computing 60, 229–241 (1998). [CrossRef]

, 28

L. Zschiedrich, R. Klose, A. Schädle, and F. Schmidt, “A new finite element realization of the perfectly matched layer method for helmholtz scattering problems on polygonal domains in two dimensions,” J. Comput. Appl. Math. 188, 12–32 (2006). [CrossRef]

]. This is true, but only in the limit where the discretization error is negligible. Once the discretization reflections dominate, we show in Sec. 4 that the convergence rate with PML thickness depends on the smoothness of the PML profile in the same way as for any other adiabatic absorber, and the rate is only polynomial for a fixed polynomial profile. That is, there is a universal relationship between smoothness and reflectivity for all adiabatic absorbers, whether discrete or continuous and whether PML or non-PML. Other authors have remarked that the numerical reflection seems to be dominated by the discontinuity in the profile or its derivatives at the PML boundary [2

A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech, 2000).

], but have not presented a precise analysis of the relationship between convergence and smoothness.

The following paper is structured as follows.We begin, in Sec. 2, with a very brief review of the derivation of PML in the simple case of one and two dimensions, and define the key quantities. Then, in Sec. 3, we explain and demonstrate the failure of PML for periodic media, even in the simplest case of one-dimensional structures where only normal-incident, non-evanescent waves are present, and even when the dielectric function varies analytically (sinusoidally). In fact, in this case, pPML may do no better than an ordinary absorbing medium (e.g., a scalar electric conductivity). Next, in Sec. 4, we analyze the relationship of the reflection to the smoothness of the absorption profile, and show via both 1d and 2d numerical calculations that the asymptotic behavior is predicted by coupled-mode theory, as well as the effect of group velocity. In Sec. 5, we describe how the coupled-mode understanding of this transition reflection points the way towards improved absorbing layers—ideally, layers whose reflection decreases exponentially with thickness (not the case even for true PML with a conventional quadratic profile, as mentioned above). Finally, we conclude with some remarks about future directions in Sec. 6.

2. Brief review of PML

Consider Maxwell’s equations in two dimensions (xy) for the TM polarization, in which the electric field (E) is in the z direction and the magnetic field (H) is in the xy plane, for a current source J z and a dielectric function ε(x,y) in natural units (ε0 = μ0 = 1), with time-harmonic fields (time-dependence~e -iωt ) are:

∇×H= ∂Hy ∂x − ∂Hx ∂y= − iωεEz
(1)
∂Ez ∂y= iωHx
(2)
∂Ez ∂x= − iωHy
(3)

One can now derive a PML absorbing boundary in the x direction, assuming for now that ε is a function of y only (e.g., the medium is homogeneous, or a waveguide in the x direction, near the computational cell boundary). In this case, one performs an analytic continuation to complex x coordinates by the transformation:

∂ ∂x→ 1 1+i σ (x)ω ∂ ∂x,
(4)

in terms of a PML profile σ(x), which plays the role of a conductivity or absorption strength. The profile σ(x) can also be a complex function, where the imaginary part corresponds to a real coordinate stretching and enhances the attenuation of purely evanescent waves [2

A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech, 2000).

, 29

J. Fang and Z. Wu, “Generalized perfectly matched layer for the absorption of propagating and evanescent waves in lossless and lossy media,” IEEE Trans. Microwave Theory Tech. 44, 2216–2222 (1998). [CrossRef]

], but in this paper we focus on the case of real σ and the absorption of propagating waves. Maxwell’s equations then become:

∂Hy ∂x− ( 1+ iσω) ∂Hx ∂y= − iωεEz+ σεEz
(5)
∂Ez ∂y= iωHx
(6)
∂Ez ∂x= − iωHy+ σHy
(7)

Note the σεEz and σHy terms, which have the form of electric and magnetic conductivities, respectively. The remaining iσ/ω term becomes an integral or convolution in time-domain and is typically handled by integrating an auxiliary differential Eq. [2

A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech, 2000).

], but is trivial in frequency domain. The extension to PMLs in other directions is straightforward and is not reviewed here.

In a medium independent of x, the wave solutions can be decomposed into normal modes with x dependence exp(ikxx) and kx > 0 for right-going waves in a right-handed [30

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 Wirel. Compon. Lett. 14, 301–303 (2004). [CrossRef]

] medium (e.g. planewaves in a homogeneous medium or waveguide modes in a waveguide). The point of this transformation (4

M. Koshiba, Y. Tsuji, and S. Sasaki, “High-performance absorbing boundary conditions for photonic crystal waveguide simulations,” IEEE Microwave Wirel. Compon. Lett. 11, 152–154 (2001). [CrossRef]

) is that these normal modes are thereby analytically continued to decaying solutions exp⁡ [ ikxx − kxω ∫xσ ( x′)dx′] wherever σ > 0. The 1/ω factor is desirable because, at least in a homogeneous dispersionless medium, the attenuation factor kx /ω is independent of frequency (but not of incidence angle).

Outside the PML regions, where σ = 0, the wave equation and thus the solution are unchanged, and it is only inside the PML (σ > 0) that the oscillating solution becomes exponentially decaying with no reflections (in theory) no matter how fast σ changes, even if σ changes discontinuously. After a short distance L in the PML, the computational cell can then be truncated (e.g. with Dirichlet boundaries), with an exponentially small round-trip reflection

R round−trip~ e −4 kxω ∫0Lσ ( x′) dx′,
(8)

where we have started the PML at x = 0, and the factor of 4 is because the reflection is proportional to the round-trip (2L) field squared.

In the exact Maxwell equations, the PML could be made arbitrarily thin by making σ very large, but this is not feasible in practice because, once Maxwell’s equations are discretized (in a finite-difference or finite-element scheme) the reflectionless property disappears. That is, it is not meaningful to analytically continue the discretized equations, and thus in the discretized system there are numerical reflections from the PML boundary that disappear in the limit of high resolution. To reduce these numerical reflections, most authors suggest that the PML be turned on gradually, i.e. that σ(x) be a continuous function starting at zero, typically chosen to grow quadratically or cubically [2

A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech, 2000).

].

More precisely, let us define σ(x) in the PML (x∈[0,L]) by a shape function s(u) ∈ [0,1]:

σ (x)= σ0s ( xL)
(9)

where the argument of s(u) is a rescaled coordinate u = x/L ∈ [0,1] and σ0 is an overall amplitude set to achieve some theoretical maximum round-trip absorption R 0 for normal-incident waves in a medium of index n(kx =ωn). Using Eq. (8) for R 0, we define:

σ0= −lnR0 4Ln ∫01s ( u′) du′.
(10)

For x < 0, outside the PML, σ = 0, i.e. s(u < 0) = 0. As L is made longer and longer for a fixed s(u), the PML profile σ turns on more and more gradually [both because s(u) is stretched out and because σ0 decreases], and the numerical reflections decrease. Several authors have suggested s(u) = u 2 (quadratic) or s(u) = u 3 (cubic) turn-on of the PML, which have discontinuities at u = 0 in the second and third derivatives respectively [2

A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech, 2000).

]. In Sec. 4, we show that there is a simple correspondence between the smoothness of s(u) and the rate of decrease of absorption with L, as a consequence of the adiabatic theorem and coupled-mode theory. Note that the smoothness of s(u) is still relevant in a discretized system—with a fixed resolution and wavelength, as L is increased one samples s(x/L) more and more finely and a discrete version of the adiabatic theorem applies [31

A. Dranov, J. Kellendonk, and R. Seller, “Discrete time adiabatic theorems for quantum mechanical systems,” J. Math. Phys. 39, 1340–1349 (1998). [CrossRef]

].

In fact, we will see that the same adiabatic theorem and the same rate of decrease apply for any absorption, whether or not the absorbing material forms a PML. For example, if we only include σ on the right-hand-side of Eq. (5), and neither on the left-hand-side nor in Eqs. (6) and (7), it corresponds to an ordinary scalar electric conductivity. As we see in Sec. 3, the advantage of PML over this ordinary conductivity is not that the reflection decreases faster with L, but that this decrease is multiplied by a much smaller constant factor (which decreases with increasing resolution) in the case of PML. This advantage mostly disappears for periodic media where analytic continuation fails, but the same relationship between the rate of decrease and the smoothness of s(u) applies.

In general, therefore, we will divide the reflections from PML into two categories: the exponentially small round-trip reflections (above), and transition reflections from the boundary between σ ≠ 0 and σ = 0 (which can arise either from numerical discretization or from other failures of PML as described in the next section). It is possible to obtain exactly zero reflection by balancing the round-trip and transition reflections so that they destructively interfere, but this cancellation can only occur for isolated wavelengths (and incident angles) [12

J. S. Juntunen, N. V. Kantartzis, and T. D. Tsiboukis, “Zero reflection coefficient in discretized PML,” IEEE Microwave Wirel. Compon. Lett. 11, 155–157 (2001). [CrossRef]

] and hence is not useful in general. Instead, we will begin by setting the estimated round-trip reflection R 0 to be negligibly small (10-25) and focus on the transition reflection; we return to the question of balancing round-trip and transition reflections in Sec. 5.

3. Failure of PML

To illustrate the failure of PML in periodic media, we consider a finite-difference frequency-domain simulation (FDFD, with a second-order-accurateYee grid) [32

A. Christ and H. L. Hartnagel, “Three-dimensional finite-difference method for the analysis of microwave-device embedding,” IEEE Trans. Microwave Theory Tech. 35, 688–696 (1987). [CrossRef]

] of the simplest possible case: a periodic dielectric function ε(x) in one dimension [so that we only have the Ez and Hy fields in Eqs. (5) and (7)]. Given a point dipole source at some position (outside the absorber), we then compute the reflection coefficient from a pPML of thickness L as a function of both L and resolution.

Here, pPML (pseudo-PML) is defined by using Eqs. (5) and (7): exactly the same equations as for an ordinary PML, but with an inhomogeneous ε function overlapping the “PML” as in Refs. 4–8. For comparison, we also show a non-PML absorber in which σ is included only in Eq. (5) but not in Eq. (7), i.e. an ordinary electric conductivity only. We consider two dielectric functions: vacuum (ε = 1) for comparison, and a periodic dielectric function ε (x) = 6+5sin(2πx/a) that varies from 1 to 11 with period a. Like all one-dimensional periodic structures, this ε(x) has photonic band gaps that prohibit propagation in certain frequency ranges [3

J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals: Molding the Flow of Light, 2nd Edition (Princeton Univ. Press, 2008).

], but we operate at a vacuum wavelength ≈ a slightly below the first bandgap (at a wavevector kx = 0.9π/a and vacuum wavelength λ = 0.9597a). The reflection is computed as the squared amplitude of the reflected Bloch wave, given by the total field minus the incident Bloch wave (computed by numerically solving for the Bloch waves of the discretized unit cell). Of course, there are two boundaries, at +x and -x, but we make the latter reflection negligible by using an absorber of thickness 5L on the left (and verified that further increasing the left-absorber thickness does not change the result). In this section, we use a quadratic shape function s(u) = u 2 for the absorber profile σ as defined above.

Fig. 2. Reflection coefficient as a function of discretization resolution for both a uniform medium and a periodic medium with PML and non-PML absorbing boundaries (insets). For the periodic medium, PML fails to be reflectionless even in the limit of high resolution, and does no better than a non-PML absorber. Inset: reflection as a function of absorber thickness L for fixed resolution ~ 50pixels/λ : as the absorber becomes thicker and the absorption is turned on more gradually, reflection goes to zero via the adiabatic theorem; PML for the uniform medium only improves the constant factor.

The absorber, here, is a pPML because it is not derived by analytic continuation of the dielectric function, and is instead formed by simply applying the homogeneous-PML equations on top of the inhomogeneous medium, leading to intrinsic reflections. However, in this case the periodic ε(x) function is actually analytic in x, so in principle one could have derived a true PML by using Eq. (5) with the analytically continued dielectric function ε [ x+ iω ∫xσ ( x′)dx′]. Unfortunately, this introduces new problems. Even though the boundary of such a PML would be reflectionless (in the exact Maxwell equations), the field is exponentially growing rather than decaying, and truncating the PML to a finite thickness would be impossible without large reflections from the truncation. The reason is that a propagating Bloch mode generally includes Fourier components with both positive and negative wavevectors, and the latter analytically continue to exponential growth. In any case, this possibility is not applicable in the vast majority of practical periodic structures, which more commonly involve a discontinuous (non-analytic) ε, so we do not consider analytically continuing ε(x) further here and focus only on the pPML case.

Figure 2 shows the results of these one-dimensional FDFD simulations, and the difference between the uniform medium (where PML works) and the periodic medium (where it does not) is stark. In the uniform medium, the reflection from PML rapidly goes to zero as resolution is increased (and in fact, goes to zero quadratically with resolution because FDFD’s center-difference discretization is second-order accurate), whereas the non-PML absorber in the uniform medium goes to a constant nonzero reflection (the Fresnel reflection coefficient from the exact Maxwell equations). For the periodic medium, both the pPML and non-PML absorbers behave roughly the same, going to a constant nonzero-reflection in the high-resolution limit: the pPML is not reflectionless for the exact Maxwell equations.

One way of understanding why pPML is not reflectionless for a periodic medium was described in the previous section: the equations with ’PML’ absorption are no longer derived via analytic continuation of Maxwell’s equations, and so the fundamental justification for PML disappears. This has nothing to do with either evanescent waves or glancing-angle waves, neither of which are present in one dimension, nor is it a numerical reflection from discretization (since it does not vanish as resolution is increased). Another way of understanding this is that the propagatingwaves in a periodic medium are Bloch waves [3

J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals: Molding the Flow of Light, 2nd Edition (Princeton Univ. Press, 2008).

], and consist of a superposition of reflections from all interfaces (all places where ε changes) in the medium—when we absorb waves reflected from interfaces within the “PML,” we have effectively terminated the periodicity and hence see reflections from this termination. (Similar but even stronger reflections are observed if one terminates the periodicity before it enters the absorbing region [23

A. Mekis, S. Fan, and J. D. Joannopoulos, “Absorbing boundary conditions for FDTD simulations of photonic crystal waveguides,” IEEE Microwave and Guided Wave Lett. 9, 502–504 (1999). [CrossRef]

].)

However, the inset of Fig. 2 shows a way in which the reflections can still be made small for the periodic medium: by increasing the thickness L of the absorbing layer. As L is increased, we see that the reflections in all four cases (PML and non-PML, uniform and periodic) go to zero as 1/L 6 asymptotically (although the periodic media take longer to attain this asymptotic power law). The true PML in the uniform medium is only different in that it has a better constant factor (which depends on resolution). The reason for this, as described in the next section, is that all transition reflections can be understood via the same coupled-mode mechanism, and the 1/L 6 rate is a consequence of the second-derivative discontinuity in s(u)=u 2. This reduction of reflection with L is adiabatic absorption, distinct from the PML concept, and it is such adiabatic absorption that one must better understand in order to efficiently truncate periodic media.

4. Smoothness & reflection

In this section, we demonstrate and explain the relationship between the smoothness of the absorber profile’s shape function s(u) and the dependence of reflection on absorber thickness L. The basic principle is that, as L increases, the rate of change of the absorption (PML or otherwise) becomes more and more gradual—as it approaches a perfectly uniform (or perfectly periodic) limit, there is an adiabatic theorem stating that the reflections must go to zero. Such an adiabatic theorem is the well-understood mechanism behind gradual waveguide tapers [33

M. Povinelli, S. Johnson, and J. Joannopoulos, “Slow-light, band-edge waveguides for tunable time delays,” Opt. Express 13, 7145–7159 (2005). [CrossRef] [PubMed]

], and adiabatic theorems also hold in periodic media with slowly varying unit cells [10

S. G. Johnson, P. Bienstman, M. Skorobogatiy, M. Ibanescu, E. Lidorikis, and J. Joannopoulos, “Adiabatic theorem and continuous coupled-mode theory for efficient taper transitions in photonic crystals,”Phys. Rev. E 66, 066608 (2002).

], and there is also an adiabatic theorem for slowly-varying discretized systems [31

A. Dranov, J. Kellendonk, and R. Seller, “Discrete time adiabatic theorems for quantum mechanical systems,” J. Math. Phys. 39, 1340–1349 (1998). [CrossRef]

]. Moreover, as we discuss in the next section, the rate at which the adiabatic (zero-reflection) limit is approached is determined by the smoothness of the transition s(u).

4.1. Numerical results

First, however, let us present the results of numerical experiments using second-order FDFD discretization for four structures: uniform and periodic media in one and two dimensions (with continuous and discontinuous ε, respectively). The reflection versus PML/pPML absorber length L in one dimension is shown for uniform media in Fig. 3 and for a periodic medium (the same medium as for Fig. 2) in Fig. 4, for a variety of shape functions s(u) = u d for exponents d∈{1,2,3,4,5}. In both cases, there is a striking pattern: the reflection asymptotically follows a power law 1/L 2d+2, which we will explain analytically below in terms of the smoothness of s(u).

In two dimensions, we looked at the boundary reflection from a point source at the center of the cell. In this case, defining a single “reflection” coefficient is more difficult because the point source emits waves at multiple angles. Instead, we look at the convergence of the electric field as L is increased, and defined a field convergence factor

∣ Ez ( L+1) ( x,y) − Ez (L) ( x,y)∣2 ∣ Ez (L) ( x,y)∣2
(11)

in terms of the electric field Ez at a point (x,y) (chosen roughly halfway between the point source and the absorbing layer) for two PML/pPML thicknesses L and L+1. This difference should go to zero as L → ∞, assuming that the reflection goes to zero in this limit (and hence the field converges to the solution for open boundaries). Indeed, this adiabatic limit is observed for both the uniform medium (vacuum) in Fig. 5 and for a periodic medium (a square lattice of width-0.7a square air holes in ε = 12) in Fig. 6. Again, there is a simple power-law relationship evident in both plots: when s(u) = ud , the field convergence factor goes as 1/L 2d+4.

Fig. 3. Reflectivity vs. PML thickness L for 1d vacuum (inset) at a resolution of 50pixels/λ for various shape functions s(u) ranging from linear [s(u) = u] to quintic [s(u) = u 5]. For reference, the corresponding asymptotic power laws are shown as dashed lines.

In 1d, we found that the reflection went as 1/L 2d+2 for s(u) =ud , and in 2d we found that the corresponding field convergence factor went as 1/L 2d+4. These two results are mathematically equivalent, for the following reason. Suppose that the reflection coefficient (for waves at any angle) goes asymptotically as 1/L 2α for some exponent α; it follows that the reflected electric field goes as 1/L α , and hence E (L)(x,y) = E (∞)(x,y)+O(1/L α). Substituting this expression into Eq. (11

E. A. Marengo, C. M. Rappaport, and E. L. Miller, “Optimum PML ABC conductivity profile in FDFD,” IEEE Trans. Magn. 35, 1506–1509 (1999). [CrossRef]

) and expanding in powers of 1/L, one finds that the field convergence factor goes as 1/L 2α+2, exactly the difference of 1/L 2 that we observed above. [There is a subtlety in this derivation: it implicitly assumes that the phase of the O(1/L α ) term, i.e. the reflected phase, goes to a constant as L→∞ in order to expand in powers of 1/L. This assumption is confirmed by our numerical results, but it is also predicted analytically by the coupled-mode theory result Eq. (13) in the next section.]

Fig. 4. Reflectivity vs. pPML thickness L for the 1d periodic medium (inset) with period a, as in Fig. 2, at a resolution of 50pixels/a with a wavevector kx = 0.9π/a (vacuum wavelength λ = 0.9597a, just below the first gap) for various shape functions s(u) ranging from linear [s(u) = u] to quintic [s(u) = u 5]. For reference, the corresponding asymptotic power laws are shown as dashed lines.

4.2. Analysis

The natural way to analyze waves propagating along a medium that is slowly varying in the propagation direction (say x) is coupled-mode theory (or coupled-wave theory) [10

S. G. Johnson, P. Bienstman, M. Skorobogatiy, M. Ibanescu, E. Lidorikis, and J. Joannopoulos, “Adiabatic theorem and continuous coupled-mode theory for efficient taper transitions in photonic crystals,”Phys. Rev. E 66, 066608 (2002).

,34

D. Marcuse, Theory of Dielectric Optical Waveguides, 2nd Edition (Academic Press, San Diego, 1991).

]: at each x, one expands the fields in the basis of the eigenmodes (indexed by ℓ) of a uniform structure with that cross-section in terms of expansion coefficients cℓ(x). (The eigenmodes have x-dependence e iβℓx for some propagation constants βℓ.) The expansion coefficients cℓ in this basis are then determined by a set of ordinary differential equations for dc ℓ/dx coupling the different modes, where the coupling coefficient is proportional to the rate of change [here, the derivative s′(x/L)]. In the limit where the structure varies more and more slowly, the solution approaches an “adiabatic” limit in which the c ℓ are nearly constant (i.e. no scattering between modes). Although coupled-mode theory was originally developed for media that are slowly varying in the propagation direction [34

D. Marcuse, Theory of Dielectric Optical Waveguides, 2nd Edition (Academic Press, San Diego, 1991).

], it has been generalized to periodic media with a slowly varying unit cell [10

S. G. Johnson, P. Bienstman, M. Skorobogatiy, M. Ibanescu, E. Lidorikis, and J. Joannopoulos, “Adiabatic theorem and continuous coupled-mode theory for efficient taper transitions in photonic crystals,”Phys. Rev. E 66, 066608 (2002).

], in which very similar coupled-mode equations appear. A similar adiabatic limit has also been derived for slowly varying discrete systems. Using coupled-mode theory, one can derive a universal relationship between the smoothness of the rate of change [s′(u)] and the asymptotic rate of convergence to the adiabatic limit. This relationship, derived below, analytically predicts the convergence rates of the reflection with absorber length observed in the previous section.

We omit the derivation of the coupled-mode equations here; their general form is considered in detail elsewhere [10

S. G. Johnson, P. Bienstman, M. Skorobogatiy, M. Ibanescu, E. Lidorikis, and J. Joannopoulos, “Adiabatic theorem and continuous coupled-mode theory for efficient taper transitions in photonic crystals,”Phys. Rev. E 66, 066608 (2002).

, 34

D. Marcuse, Theory of Dielectric Optical Waveguides, 2nd Edition (Academic Press, San Diego, 1991).

].We simply quote the result: in the limit of slow variation (large L), the equations can be solved to lowest order in 1/L in terms of a simple integral. In particular, if the structure is smoothly parameterized by a shape function s(x/L) (e.g. the absorption profile as given here), then the amplitude cr (corresponding to a reflected power |cr |2) of a reflected mode is given to lowest-order (for large L) by [10

S. G. Johnson, P. Bienstman, M. Skorobogatiy, M. Ibanescu, E. Lidorikis, and J. Joannopoulos, “Adiabatic theorem and continuous coupled-mode theory for efficient taper transitions in photonic crystals,”Phys. Rev. E 66, 066608 (2002).

]:

Fig. 5. Field convergence factor [Eq. (11

E. A. Marengo, C. M. Rappaport, and E. L. Miller, “Optimum PML ABC conductivity profile in FDFD,” IEEE Trans. Magn. 35, 1506–1509 (1999). [CrossRef]

)] (~ reflection/L 2) vs. PML thickness L for 2d vacuum (inset) at a resolution of 20pixels/λ for various shape functions s(u) ranging from linear [s(u) = u] to quintic [s(u) = u 5]. For reference, the corresponding asymptotic power laws are shown as dashed lines. Inset: ℜ[Ez ] field pattern for the (point) source at the origin (blue/white/red = positive/zero/negative).
cr (L)= ∫01 s′ (u) M [s (u)] Δβ [s (u)] e iL ∫0uΔβ [s ( u′)] du′du
(12)

Here, M is a coupling coefficient depending on the mode overlap between the incident and reflected field (in the changing part of the structure) and Δβ ≠ 0 is the differenceβ i -β r between the propagation constants of the incident and reflected modes. Both of these are some analytic functions of the shape s(u). In general, there may be more than one reflected mode, and in a periodic structure the coefficient even for a single reflected mode is a sum of contributions of above form from the different Brillouin zones [10

S. G. Johnson, P. Bienstman, M. Skorobogatiy, M. Ibanescu, E. Lidorikis, and J. Joannopoulos, “Adiabatic theorem and continuous coupled-mode theory for efficient taper transitions in photonic crystals,”Phys. Rev. E 66, 066608 (2002).

], but it suffices to analyze the rate of convergence of a single such integral with L. The basic reason for the adiabatic limit is that, as L grows, the phase term oscillates faster and faster and the integral of this oscillating quantity goes to zero.

There are many standard methods to analyze the asymptotic (large L) properties of such an integral. In particular, we apply a technique that is commonly used to analyze the convergence rate of Fourier series: one simply integrates by parts repeatedly until a nonzero boundary term is obtained [35

K. O. Mead and L. M. Delves, “On the convergence rate of generalized fourier expansions,” IMA J. Appl. Math. 12, 247–259 (1973). [CrossRef]

,36

J. P. Boyd, Chebyshev and Fourier Spectral Methods, 2nd Edition (Springer, 1989).

]. Each integration by parts integrates the e iL∫Δβ term, dividing the integrand by iLΔβ(u), and differentiates the s′M/δβ term. After integrating by parts d times, the boundary term at u = 0 is zero if the corresponding derivative s (d)(0+) is zero, whereas the boundary term at u = 1 is always negligible because of the absorption (leading to a complex Δβ and exponential decay), assuming a small round-trip reflection R 0. The dominant asymptotic term is the first (lowest-d) u = 0 boundary term that is nonzero, since all subsequent integrations by parts have an additional factor of 1/L. [Here, we have assumed that s is a smooth function in (0,1) so that there are never delta-function contributions from the interior.] The result is the following asymptotic form for cr (L), independent of the particular details of the geometry or the modes:

Fig. 6. Field convergence factor [Eq. (11)] (~ reflection/L 2) vs. pPML thickness L for the discontinuous 2d periodic medium (left inset: square lattice of square air holes in ε = 12) with period a, at a resolution of 10pixels/a with a vacuum wavelength λ = 0.6667a (not in a band gap) for various shape functions s(u) ranging from linear [s(u) = u] to quintic [s(u) = u 5]. For reference, the corresponding asymptotic power laws are shown as dashed lines. Right inset: ℜ[Ez ] field pattern for the (point) source at the origin (blue/white/red = positive/zero/negative).
cr (L)= s (d) ( 0+) M ( 0+) Δβ ( 0+) [ −iLΔβ ( 0+)] −d +O ( L − ( d+1)),
(13)

where s (d)(0+) is the first nonzero derivative of s(u) at u= 0+, and integrating by parts d times yielded a division by (-iLΔβ) d (flipping sign each time). This result corresponds to what is sometimes called “Darboux’s principle:” the convergence is dominated by the lowest-order singularity [36

J. P. Boyd, Chebyshev and Fourier Spectral Methods, 2nd Edition (Springer, 1989).

], which here is the first discontinuity in the rate of change s′(u) at u = 0. A similar result applies, for example, to the convergence rate of a Fourier series: a function that has a discontinuity in the d-th derivative has a Fourier series whose coefficients cn decrease asymptotically as 1/n (d+1) [35

K. O. Mead and L. M. Delves, “On the convergence rate of generalized fourier expansions,” IMA J. Appl. Math. 12, 247–259 (1973). [CrossRef]

, 36

J. P. Boyd, Chebyshev and Fourier Spectral Methods, 2nd Edition (Springer, 1989).

] (the d +1 instead of d is due to the fact that our integral starts with s′).

Equation (13) would seem to imply that the reflection ~|cr |2 is O(L -2d ), but this is not the case because there is a hidden 1/L factor in the coupling coefficient M, thanks to the 1/L dependence of σ0 in Eq. (10). The coupling coefficient M is a matrix element proportional to the rate of change of the materials [10

S. G. Johnson, P. Bienstman, M. Skorobogatiy, M. Ibanescu, E. Lidorikis, and J. Joannopoulos, “Adiabatic theorem and continuous coupled-mode theory for efficient taper transitions in photonic crystals,”Phys. Rev. E 66, 066608 (2002).

], which in this case is ∂σ ∂u = s′ (u) σ0~ 1L. Therefore, the reflection scales as |M|2/L 2d =O(L -(2d+2)), exactly corresponding to our numerical results above.

Other useful results can be obtained from Eq. (13), and in particular one can show that the reflections due to nonuniformity worsen in a periodic structure as a flat band edge (β0,ω0) is approached [33

M. Povinelli, S. Johnson, and J. Joannopoulos, “Slow-light, band-edge waveguides for tunable time delays,” Opt. Express 13, 7145–7159 (2005). [CrossRef] [PubMed]

]. As a quadratic-shaped band-edgeω -ω 0~(β -β0)2 is approached, the group velocity νg= dωdβ scales proportional to β -β0, while the Δβ between the forward and reflected modes is 2(β -β0) ~ vg . Also, the coupling coefficient M is proportional to 1/vg because of the constant-power normalizations of the incident and reflected modes [10

S. G. Johnson, P. Bienstman, M. Skorobogatiy, M. Ibanescu, E. Lidorikis, and J. Joannopoulos, “Adiabatic theorem and continuous coupled-mode theory for efficient taper transitions in photonic crystals,”Phys. Rev. E 66, 066608 (2002).

, 33

M. Povinelli, S. Johnson, and J. Joannopoulos, “Slow-light, band-edge waveguides for tunable time delays,” Opt. Express 13, 7145–7159 (2005). [CrossRef] [PubMed]

]. Hence, by inspection of Eq. (13), the reflection |cr |2 = O(vg -(2d+4)). For example, the reflection is O(vg -6)for a linear taper s(u) = u [33

M. Povinelli, S. Johnson, and J. Joannopoulos, “Slow-light, band-edge waveguides for tunable time delays,” Opt. Express 13, 7145–7159 (2005). [CrossRef] [PubMed]

]. Because of this unfavorable scaling, an imperfect absorbing layer such as a pPML is most challenging in periodic structures when operating close to a band edge where there are slow-light modes (in the same way that other taper transitions are challenging in this regime [33

M. Povinelli, S. Johnson, and J. Joannopoulos, “Slow-light, band-edge waveguides for tunable time delays,” Opt. Express 13, 7145–7159 (2005). [CrossRef] [PubMed]

]).

There is one thing missing from the above analysis, and that is the discretized-space adiabatic case. In a slowly varying discrete system [i.e, sampling some slow change sn = s(nΔx/L) as L grows larger], there is still a proof of the adiabatic theorem (cr → 0), but the only published proof is currently for the lossless case (unitary evolution) [31

A. Dranov, J. Kellendonk, and R. Seller, “Discrete time adiabatic theorems for quantum mechanical systems,” J. Math. Phys. 39, 1340–1349 (1998). [CrossRef]

]. Also, an analogous integral form of the lowest-order reflection has not been presented, nor has the rate of convergence to the adiabatic limit been analyzed in the discrete case. So, our prediction of the asympototic convergence rate is rigorously proven only for the case of the continuous-space wave propagation. However, our numerical results demonstrate that a slowly-changing discretized system exhibits exactly the same scaling (e.g. in the PML case for uniform media, where the only reflections are due to discretization). (This seems analogous to the fact that the discretization error of a discrete Fourier transform converges at the same rate as the decay of the coefficients of the continuous-space Fourier series [36

J. P. Boyd, Chebyshev and Fourier Spectral Methods, 2nd Edition (Springer, 1989).

].) In future work, we hope to further validate our numerical result for the convergence rate in discretized space with a proper generalization of the coupled-mode analysis.

5. Towards better absorbers

From the previous section, there is a close relationship between the smoothness of the absorption profile and the asymptotic convergence rate of the reflections R(L) as a function of absorber thickness L: if the profile s(u) has a discontinuity in the d-th derivative (e.g. for s = ud ), then the reflection coefficient goes as 1/L 2d+2 for a fixed round-trip reflection. This result raises several interesting questions. Can one do better than polynomial convergence? What is the optimal shape s(u)? And what if the round-trip reflection is not fixed?

The above result relating smoothness and convergence has a natural corollary: if s(u) is C ∞, i.e. all of its derivatives are continuous, then the reflection goes to zero faster than any polynomial in 1/L. This is similar to a well-known result for the convergence of Fourier series of C ∞ functions [36

J. P. Boyd, Chebyshev and Fourier Spectral Methods, 2nd Edition (Springer, 1989).

]; the exact rate of faster-than-polynomial convergence again depends on the strongest singularity in s(u). For example, for s(u) = (tanh(u)+1)/2, which goes exponentially to zero as u→-∞ and to one as u→+∞, the reflection should decrease exponentially with L, as determined by contour integration from the residue of the pole at u= ±iπ/2 that is closest to the real axis (similar to the analysis for the convergence of a Fourier series for an analytic function [36

J. P. Boyd, Chebyshev and Fourier Spectral Methods, 2nd Edition (Springer, 1989).

, 37

D. Elliott, “The evaluation and estimation of the coefficients in the Chebyshev series expansion of a function,” Mathematics of Computation 18, 274–284 (1964). [CrossRef]

]). However, such an absorption taper would require an infinitely thick absorber in order to avoid discontinuously truncating the exponential tail of tanh(u). To have a C ∞ function with a finite absorber, with s(u)=0 for u≤0, the s(u) function must be non-analytic; a standard example of such a function is s(u) = e 1-1/u for u> 0 (all of whose derivatives go to zero as u→0+, where there is an essential singularity). Because s(u)=e 1-1/u is C ∞, its reflection R(L) must decrease faster than any polynomial. Exactly how much faster than polynomial is determined by asymptotically evaluating the integral of Eq. (12) by a saddle-point method [38

J. P. Boyd, “The optimization of convergence for Chebyshev polynomial methods in an unbounded domain,” J. Comput. Phys. 45, 43–79 (1982). [CrossRef]

, 39

H. Cheng, Advanced Analytic Methods in Applied Mathematics, Science and Engineering (Luban Press, Boston, 2006).

]: the result is that R(L) decays asymptotically as e -α√L for some constant α > 0 [38

J. P. Boyd, “The optimization of convergence for Chebyshev polynomial methods in an unbounded domain,” J. Comput. Phys. 45, 43–79 (1982). [CrossRef]

]. This is confirmed by Fig. 7, which plots the PML/pPML reflection for the 1d uniform and periodic cases on a semilog scale versus √L, and results clearly approach a straight line as expected.

Fig. 7. Reflectivity vs. PML thickness L for 1d vacuum (blue circles) at a resolution of 50pixels/λ, and for pPML thickness L in the 1d periodic medium of Fig. 4 (red squares) with period a at a resolution of 50pixels/a with a wavevector kx = 0.9π/a (vacuum wavelength λ = 0.9597a. In both cases, a C ∞ (infinitely differentiable) shape function s(u) = e 1-1/u for u > 0 is used, leading to asymptotic convergence as e -α√L for some constants α.

Although s(u) = e 1-1/u yields an exponential convergence of the absorption in Fig. 7, the constant factor and the exponential rate are almost certainly suboptimal for this arbitrary choice of C ∞ function. If we compare Fig. 7 to Fig. 3 for the uniform case and Fig. 4 for the periodic case, we see that this C ∞ s(u) is superior to the polynomial s(u) for the periodic case where PML is not perfect, but inferior for the uniform case until the reflection becomes inconsequential (~ 10-20). This is still a useful result in the sense that one mainly needs to improve pPML for the periodic case, whereas PML is already good enough for uniform media. However, one would ideally prefer a shape function that is consistently better than the polynomial s(u), regardless of the dielectric function, so further exploration of the space of possible absorption profiles seems warranted.

Finally, in the above analysis we fixed the round-trip reflection R 0, via the estimate of Eq.(8), to approximately 10-25 in order for our calculations to isolate the effect of the transition reflection. Obviously, in a real application, one is unlikely to require such low reflections and one will set R 0 to a larger value, corresponding to a larger σ0 ~ lnR 0 in Eq. (10). This will also reduce the transition reflection [as seen from Eq. (13)], but only by a logarithmic constant factor. The best choice to minimize reflection for a given absorber length, in principle, is to set R 0 to be roughly equal to the transition reflection for that length. (Another reason to make them equal is the possibility of destructive interference between the round-trip and transition reflection [12

J. S. Juntunen, N. V. Kantartzis, and T. D. Tsiboukis, “Zero reflection coefficient in discretized PML,” IEEE Microwave Wirel. Compon. Lett. 11, 155–157 (2001). [CrossRef]

], but such destructive interference is inherently restricted to narrow bandwidths and ranges of incident angles and so we do not concern ourselves with this possibility.) In order to make them roughly equal, one needs an estimate of the transition reflection; for example, one could simply numerically fit the power law of Eq. (13). The result of such matching is shown in Fig. 8 for a quadratic profile s(u) = u 2 in 1d uniform media, and the overall reflection is reduced by a factor of 3-400 compared to a fixed R 0 = 10-16. This is a significant reduction, but is not overwhelming (especially for smaller L) and changes the asymptotic convergence rate Eq. (13) only by a factor of lnR 0~lnL. The drawback of this optimization is that it is difficult to determine the transition reflection analytically for inhomogeneous media, and so one is generally forced to make a conservative estimate of R 0, which reduces the advantage gained.

Fig. 8. Reflectivity vs. PML thickness L for 1d vacuum (inset) at a resolution of 50pixels/λ for s(u) =u 2, with the round-trip reflection either set to R 0 = 10-16 (upper blue line) or set to match the estimated transition reflection from Fig. 3 (lower red line). By matching the round-trip reflection R 0 to the estimated transition reflection, one can obtain a substantial reduction in the constant factor of the total reflection, although the asymptotic power law is only changed by a lnL factor

6. Conclusion

Perfectly matched layers are an extraordinarily powerful technique to absorb waves incident on the boundaries of wave-equation simulation, but they are not a panacea. In particular, for cases such as photonic crystals where the medium is not analytic in the direction perpendicular to the boundary, the fundamental coordinate-stretching idea behind PML breaks down, and the interface has intrinsic reflections (even for simple 1d cases with only normal-incident nonevanescent waves). However, one can still obtain small reflections by gradually ramping up the “pseudo-PML” (pPML) absorption, similar to the idea behind the quadratic PML profiles commonly used to circumvent discretization-based reflections in uniform media, forming an adiabatic absorber. In fact, for both cases (pPML in periodic media and PML in discretized uniform media), we show that the basic mechanism behind the reflection is determined in the same way by the smoothness of the absorption profile, which can be predicted analytically by coupled-mode theory. More generally, an adiabatic absorber is applicable in any situation where a true PML is inconvenient or impossible to implement.

The same theory then predicts that an exponential absorber, one whose reflections decrease exponentially with some power of the absorber thickness L, is possible, for example by using an infinitely differentiable absorption profile. (In contrast, ordinary PML in a uniform medium with a quadratic profile is not an exponential absorber: its numerical-discretization reflections decrease as 1/L 6.) We gave a simple C ∞ example profile that led to such exponential absorption, but much future work remains to be done in identifying profiles with both exponential absorption and good constant factors. In particular, one possibility that we will examine in a subsequent manuscript is an absorption profile whose smoothness increases with L, so that it matches simple quadratic profiles for small L but becomes exponentially smoother with large L. (Such L-varying profiles require a more careful convergence analysis, however, in order to ensure that they approach the adiabatic zero-reflection limit. A closely related mathematical idea is explored in Ref. 38.)

Acknowledgments

This work was supported in part by the U.S. Army Research Office under contract W911NF- 07-D-0004.We are also grateful to Profs. J. White and H. Cheng at MIT for helpful discussions.

References and links

1.

J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. 114,185–200 (1994). [CrossRef]

2.

A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech, 2000).

3.

J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals: Molding the Flow of Light, 2nd Edition (Princeton Univ. Press, 2008).

4.

M. Koshiba, Y. Tsuji, and S. Sasaki, “High-performance absorbing boundary conditions for photonic crystal waveguide simulations,” IEEE Microwave Wirel. Compon. Lett. 11, 152–154 (2001). [CrossRef]

5.

Y. Tsuji and M. Koshiba, “Finite element method using port truncation by perfectly matched layer boundary conditions for optical waveguide discontinuity problems,” J. Lightwave Technol. 20, 463–468 (2002). [CrossRef]

6.

E. P. Kosmidou, T. I. Kosmani, and T. D. Tsiboukis, “A comparative FDTD study of various PML configurations for the termination of nonlinear photonic bandgap waveguide structures,” IEEE Trans. Magn. 39, 1191–1194 (2003). [CrossRef]

7.

A. R. Weily, L. Horvath, K. P. Esselle, and B. C. Sanders, “Performance of PML absorbing boundary conditions in 3d photonic crystal waveguides,” Microwave Opt. Technol.Lett. 40, 1–3 (2004). [CrossRef]

8.

N. Kono and M. Koshiba, “General finite-element modeling of 2-D magnetophotonic crystal waveguides,” IEEE Photon. Tech. Lett. 17, 1432–1434 (2005). [CrossRef]

9.

W. C. Chew and J. M. Jin, “Perfectly matched layers in the discretized space: An analysis and optimization,” Electromagnetics 16, 325–340 (1996). [CrossRef]

10.

S. G. Johnson, P. Bienstman, M. Skorobogatiy, M. Ibanescu, E. Lidorikis, and J. Joannopoulos, “Adiabatic theorem and continuous coupled-mode theory for efficient taper transitions in photonic crystals,”Phys. Rev. E 66, 066608 (2002).

11.

E. A. Marengo, C. M. Rappaport, and E. L. Miller, “Optimum PML ABC conductivity profile in FDFD,” IEEE Trans. Magn. 35, 1506–1509 (1999). [CrossRef]

12.

J. S. Juntunen, N. V. Kantartzis, and T. D. Tsiboukis, “Zero reflection coefficient in discretized PML,” IEEE Microwave Wirel. Compon. Lett. 11, 155–157 (2001). [CrossRef]

13.

Y. S. Rickard and N. K. Nikolova, “Enhancing the PML absorbing boundary conditions for the wave equation,” IEEE Trans. Antennas Propag. 53, 1242–1246 (2005). [CrossRef]

14.

Z. Chen and X. Liu, “An adaptive perfectly matched layer technique for time-harmonic scattering problems,” SIAM J. Num. Anal. 43, 645–671 (2005). [CrossRef]

15.

Z. 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 Propag. 43, 1460–1463 (1995). [CrossRef]

16.

W. C. Chew and W. H. Weedon, “A 3d perfectly matched medium from modified Maxwell’s equations with stretched coordinates,” Microwave Opt. Technol. Lett. 7, 599–604 (1994). [CrossRef]

17.

C. M. Rappaport, “Perfectly matched absorbing boundary conditions based on anisotropic lossy mapping of space,” IEEE Microwave and Guided Wave Lett. 5, 90–92 (1995). [CrossRef]

18.

F. L. Teixeira and W. C. Chew, “General close-form PML constitutive tensors to match arbitrary bianisotropic and dispersive linear media,” IEEE Microwave and Guided Wave Lett. 8, 223–225 (1998). [CrossRef]

19.

F. Collino and P. B. Monk, “Optimizing the perfectly matched layer,” Comput. Methods Appl. Mech. Engrg. 164, 157–171 (1998). [CrossRef]

20.

A. J. Ward and J. B. Pendry, “Refraction and geometry in Maxwell’s equations,” J. Mod. Opt. 43, 773–793 (1996). [CrossRef]

21.

W. Huang, C. Xu, W. Lui, and K. Yokoyama, “The perfectly matched layer boundary condition for modal analysis of optical waveguides: leaky mode calculations,” IEEE Photon. Tech. Lett. 8, 652–654 (1996). [CrossRef]

22.

D. Pissoort and F. Olyslager, “Termination of periodic waveguides by PMLs in time-harmonic integral equationlike techniques,” IEEE Antennas and Wireless Propagation Lett. 2, 281–284 (2003). [CrossRef]

23.

A. Mekis, S. Fan, and J. D. Joannopoulos, “Absorbing boundary conditions for FDTD simulations of photonic crystal waveguides,” IEEE Microwave and Guided Wave Lett. 9, 502–504 (1999). [CrossRef]

24.

E. Moreno, D. Erni, and C. Hafner, “Modeling of discontinuities in photonic crystal waveguides with the multiple multipole method,” Phys. Rev. E 66, 036618 (2002).

25.

Z.-Y. Li and K.-M. Ho, “Light propagation in semi-infinite photonic crystals and related waveguide structures,” Phys. Rev. B 68 (2003).

26.

D. Pissoort, B. Denecker, P. Bienstman, F. Olyslager, and D. De Zutter , “Comparative study of three methods for the simulation of two-dimensional photonic crystals,” J. Opt. Soc. Am. A 21, 2186–2195 (2004). [CrossRef]

27.

M. Lassas and E. Somersalo, “On the existence and convergence of the solution of PML equations,” Computing 60, 229–241 (1998). [CrossRef]

28.

L. Zschiedrich, R. Klose, A. Schädle, and F. Schmidt, “A new finite element realization of the perfectly matched layer method for helmholtz scattering problems on polygonal domains in two dimensions,” J. Comput. Appl. Math. 188, 12–32 (2006). [CrossRef]

29.

J. Fang and Z. Wu, “Generalized perfectly matched layer for the absorption of propagating and evanescent waves in lossless and lossy media,” IEEE Trans. Microwave Theory Tech. 44, 2216–2222 (1998). [CrossRef]

30.

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 Wirel. Compon. Lett. 14, 301–303 (2004). [CrossRef]

31.

A. Dranov, J. Kellendonk, and R. Seller, “Discrete time adiabatic theorems for quantum mechanical systems,” J. Math. Phys. 39, 1340–1349 (1998). [CrossRef]

32.

A. Christ and H. L. Hartnagel, “Three-dimensional finite-difference method for the analysis of microwave-device embedding,” IEEE Trans. Microwave Theory Tech. 35, 688–696 (1987). [CrossRef]

33.

M. Povinelli, S. Johnson, and J. Joannopoulos, “Slow-light, band-edge waveguides for tunable time delays,” Opt. Express 13, 7145–7159 (2005). [CrossRef] [PubMed]

34.

D. Marcuse, Theory of Dielectric Optical Waveguides, 2nd Edition (Academic Press, San Diego, 1991).

35.

K. O. Mead and L. M. Delves, “On the convergence rate of generalized fourier expansions,” IMA J. Appl. Math. 12, 247–259 (1973). [CrossRef]

36.

J. P. Boyd, Chebyshev and Fourier Spectral Methods, 2nd Edition (Springer, 1989).

37.

D. Elliott, “The evaluation and estimation of the coefficients in the Chebyshev series expansion of a function,” Mathematics of Computation 18, 274–284 (1964). [CrossRef]

38.

J. P. Boyd, “The optimization of convergence for Chebyshev polynomial methods in an unbounded domain,” J. Comput. Phys. 45, 43–79 (1982). [CrossRef]

39.

H. Cheng, Advanced Analytic Methods in Applied Mathematics, Science and Engineering (Luban Press, Boston, 2006).

OCIS Codes
(000.3860) General : Mathematical methods in physics
(000.4430) General : Numerical approximation and analysis
(050.1755) Diffraction and gratings : Computational electromagnetic methods
(160.5298) Materials : Photonic crystals

ToC Category:
Diffraction and Gratings

History
Original Manuscript: May 6, 2008
Revised Manuscript: May 29, 2008
Manuscript Accepted: May 29, 2008
Published: July 14, 2008

Citation
Ardavan F. Oskooi, Lei Zhang, Yehuda Avniel, and Steven G. Johnson, "The failure of perfectly matched layers, and towards their redemption by adiabatic absorbers," Opt. Express 16, 11376-11392 (2008)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-16-15-11376


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. J.-P. Berenger, "A perfectly matched layer for the absorption of electromagnetic waves," J. Comput. Phys. 114, 185-200 (1994). [CrossRef]
  2. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech, 2000).
  3. J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals: Molding the Flow of Light, 2nd Edition (Princeton Univ. Press, 2008).
  4. M. Koshiba, Y. Tsuji, and S. Sasaki, "High-performance absorbing boundary conditions for photonic crystal waveguide simulations," IEEE Microwave Wirel. Compon. Lett. 11, 152-154 (2001). [CrossRef]
  5. Y. Tsuji and M. Koshiba, "Finite element method using port truncation by perfectly matched layer boundary conditions for optical waveguide discontinuity problems," J. Lightwave Technol. 20, 463-468 (2002). [CrossRef]
  6. E. P. Kosmidou, T. I. Kosmani, and T. D. Tsiboukis, "A comparative FDTD study of various PML configurations for the termination of nonlinear photonic bandgap waveguide structures," IEEE Trans. Magn. 39, 1191-1194 (2003). [CrossRef]
  7. A. R. Weily, L. Horvath, K. P. Esselle, and B. C. Sanders, "Performance of PML absorbing boundary conditions in 3d photonic crystal waveguides," Microwave Opt. Technol. Lett. 40, 1-3 (2004). [CrossRef]
  8. N. Kono and M. Koshiba, "General finite-element modeling of 2-D magnetophotonic crystal waveguides," IEEE Photon. Tech. Lett. 17, 1432-1434 (2005). [CrossRef]
  9. W. C. Chew and J. M. Jin, "Perfectly matched layers in the discretized space: An analysis and optimization," Electromagnetics 16, 325-340 (1996). [CrossRef]
  10. S. G. Johnson, P. Bienstman, M. Skorobogatiy, M. Ibanescu, E. Lidorikis, and J. Joannopoulos, "Adiabatic theorem and continuous coupled-mode theory for efficient taper transitions in photonic crystals," Phys. Rev. E 66, 066608 (2002).
  11. E. A. Marengo, C. M. Rappaport, and E. L. Miller, "Optimum PML ABC conductivity profile in FDFD," IEEE Trans. Magn. 35, 1506-1509 (1999). [CrossRef]
  12. J. S. Juntunen, N. V. Kantartzis, and T. D. Tsiboukis, "Zero reflection coefficient in discretized PML," IEEE Microwave Wirel. Compon. Lett. 11, 155-157 (2001). [CrossRef]
  13. Y. S. Rickard and N. K. Nikolova, "Enhancing the PML absorbing boundary conditions for the wave equation," IEEE Trans. Antennas Propag. 53, 1242-1246 (2005). [CrossRef]
  14. Z. Chen and X. Liu, "An adaptive perfectly matched layer technique for time-harmonic scattering problems," SIAM J. Num. Anal. 43, 645-671 (2005). [CrossRef]
  15. Z. 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 Propag. 43, 1460-1463 (1995). [CrossRef]
  16. W. C. Chew and W. H. Weedon, "A 3d perfectly matched medium from modified Maxwell???s equations with stretched coordinates," Microwave Opt. Technol. Lett. 7, 599-604 (1994). [CrossRef]
  17. C. M. Rappaport, "Perfectly matched absorbing boundary conditions based on anisotropic lossy mapping of space," IEEE Microwave and Guided Wave Lett. 5, 90-92 (1995). [CrossRef]
  18. F. L. Teixeira and W. C. Chew, "General close-form PML constitutive tensors to match arbitrary bianisotropic and dispersive linear media," IEEE Microwave and Guided Wave Lett. 8, 223-225 (1998). [CrossRef]
  19. F. Collino and P. B. Monk, "Optimizing the perfectly matched layer," Comput. Methods Appl. Mech. Engrg. 164, 157-171 (1998). [CrossRef]
  20. A. J. Ward and J. B. Pendry, "Refraction and geometry in Maxwell???s equations," J. Mod. Opt. 43, 773-793 (1996). [CrossRef]
  21. W. Huang, C. Xu,W. Lui, and K. Yokoyama, "The perfectly matched layer boundary condition for modal analysis of optical waveguides: leaky mode calculations," IEEE Photon. Tech. Lett. 8, 652-654 (1996). [CrossRef]
  22. D. Pissoort and F. Olyslager, "Termination of periodic waveguides by PMLs in time-harmonic integral equation-like techniques," IEEE Antennas and Wireless Propagation Lett. 2, 281-284 (2003). [CrossRef]
  23. A. Mekis, S. Fan, and J. D. Joannopoulos, "Absorbing boundary conditions for FDTD simulations of photonic crystal waveguides," IEEE Microwave and Guided Wave Lett. 9, 502-504 (1999). [CrossRef]
  24. E. Moreno, D. Erni, and C. Hafner, "Modeling of discontinuities in photonic crystal waveguides with the multiple multipole method," Phys. Rev. E 66, 036618 (2002).
  25. Z.-Y. Li and K.-M. Ho, "Light propagation in semi-infinite photonic crystals and related waveguide structures," Phys. Rev. B 68 (2003).
  26. D. Pissoort, B. Denecker, P. Bienstman, F. Olyslager, and D. De Zutter, "Comparative study of three methods for the simulation of two-dimensional photonic crystals," J. Opt. Soc. Am. A 21, 2186-2195 (2004). [CrossRef]
  27. M. Lassas and E. Somersalo, "On the existence and convergence of the solution of PML equations," Computing 60, 229-241 (1998). [CrossRef]
  28. L. Zschiedrich, R. Klose, A. Sch¨adle, and F. Schmidt, "A new finite element realization of the perfectly matched layer method for helmholtz scattering problems on polygonal domains in two dimensions," J. Comput. Appl. Math. 188, 12-32 (2006). [CrossRef]
  29. J. Fang and Z. Wu, "Generalized perfectly matched layer for the absorption of propagating and evanescent waves in lossless and lossy media," IEEE Trans. Microwave Theory Tech. 44, 2216-2222 (1998). [CrossRef]
  30. 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 Wirel. Compon. Lett. 14, 301-303 (2004). [CrossRef]
  31. A. Dranov, J. Kellendonk, and R. Seller, "Discrete time adiabatic theorems for quantum mechanical systems," J. Math. Phys. 39, 1340-1349 (1998). [CrossRef]
  32. A. Christ and H. L. Hartnagel, "Three-dimensional finite-difference method for the analysis of microwave-device embedding," IEEE Trans. Microwave Theory Tech. 35, 688-696 (1987). [CrossRef]
  33. M. Povinelli, S. Johnson, and J. Joannopoulos, "Slow-light, band-edge waveguides for tunable time delays," Opt. Express 13, 7145-7159 (2005). [CrossRef] [PubMed]
  34. D. Marcuse, Theory of Dielectric Optical Waveguides, 2nd Edition (Academic Press, San Diego, 1991).
  35. K. O. Mead and L. M. Delves, "On the convergence rate of generalized fourier expansions," IMA J. Appl. Math. 12, 247-259 (1973). [CrossRef]
  36. J. P. Boyd, Chebyshev and Fourier Spectral Methods, 2nd Edition (Springer, 1989).
  37. D. Elliott, "The evaluation and estimation of the coefficients in the Chebyshev series expansion of a function," Mathematics of Computation 18, 274-284 (1964). [CrossRef]
  38. J. P. Boyd, "The optimization of convergence for Chebyshev polynomial methods in an unbounded domain," J. Comput. Phys. 45, 43-79 (1982). [CrossRef]
  39. H. Cheng, Advanced Analytic Methods in Applied Mathematics, Science and Engineering (Luban Press, Boston, 2006).

Cited By

 Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited