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,2].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])—contrary to previous suggestions of photonic-crystal “PML” absorbers [4–8]. In these cases (similar to PML reflections due to discretization error [2,9]), the remaining approach to reduce reflections is to “turn on” the absorption gradually, asymptotically approaching an “adiabatic” limit of zero reflections [10] 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]), the reflection decreases with absorber thickness L proportional to some characteristic power law determined by the smoothness (e.g. 1/L ⁶ 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], a quadratic or cubic profile works so well [2] 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] 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]. 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]—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]. (This description can then be converted back into a change of materials via a complex coordinate transformation [18, 20].) 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]: 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).

Download Full Size | PDF

Previous suggestions to apply PML to photonic crystals by simply overlapping a “PML” anisotropic absorber with the periodic dielectric function [4–8] (including a similar suggestion for integral-equation methods [22]) 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], or alternatively to specify analytical boundary conditions of zero reflection for that one mode [24]. 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, 26], 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, 9]. 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, 28]. 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], 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 (ε₀ = μ₀ = 1), with time-harmonic fields (time-dependence~e ^-iωt) are:

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:

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, 29], but in this paper we focus on the case of real σ and the absorption of propagating waves. Maxwell’s equations then become:

Note the σεE_z and σH_y 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], 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(ik_xx) and k_x > 0 for right-going waves in a right-handed [30] medium (e.g. planewaves in a homogeneous medium or waveguide modes in a waveguide). The point of this transformation (4) is that these normal modes are thereby analytically continued to decaying solutions $\exp [{\mathit{ik}}_{x}x-\frac{k_x}{\omega }{\int }^x\sigma \left(x\prime \right)\mathit{dx}\prime ]$ wherever σ > 0. The 1/ω factor is desirable because, at least in a homogeneous dispersionless medium, the attenuation factor k_x/ω 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

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].

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

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

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 σ₀ decreases], and the numerical reflections decrease. Several authors have suggested s(u) = u ² (quadratic) or s(u) = u ³ (cubic) turn-on of the PML, which have discontinuities at u = 0 in the second and third derivatives respectively [2]. 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].

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] and hence is not useful in general. Instead, we will begin by setting the estimated round-trip reflection R ₀ 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] of the simplest possible case: a periodic dielectric function ε(x) in one dimension [so that we only have the E_z and H_y 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], but we operate at a vacuum wavelength ≈ a slightly below the first bandgap (at a wavevector k_x = 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 ² 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.

Download Full Size | PDF

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 $\epsilon [x+\frac{i}{\omega }{\int }^x\sigma \left(x\prime \right)\mathit{dx}\prime ]$. 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], 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].)

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 ⁶ 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 ⁶ rate is a consequence of the second-derivative discontinuity in s(u)=u ². 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], and adiabatic theorems also hold in periodic media with slowly varying unit cells [10], and there is also an adiabatic theorem for slowly-varying discretized systems [31]. 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

in terms of the electric field E_z 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) = u^d, 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 ⁵]. For reference, the corresponding asymptotic power laws are shown as dashed lines.

Download Full Size | PDF

In 1d, we found that the reflection went as 1/L ^2d+2 for s(u) =u^d, 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) 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 ² 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 k_x = 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 ⁵]. For reference, the corresponding asymptotic power laws are shown as dashed lines.

Download Full Size | PDF

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,34]: 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], it has been generalized to periodic media with a slowly varying unit cell [10], 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, 34].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 c_r (corresponding to a reflected power |c_r|²) of a reflected mode is given to lowest-order (for large L) by [10]:

 

Fig. 5. Field convergence factor [Eq. (11)] (~ reflection/L ²) 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 ⁵]. For reference, the corresponding asymptotic power laws are shown as dashed lines. Inset: ℜ[E_z] field pattern for the (point) source at the origin (blue/white/red = positive/zero/negative).

Download Full Size | PDF

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], 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,36]. 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 ₀. 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 c_r(L), independent of the particular details of the geometry or the modes:

 

Fig. 6. Field convergence factor [Eq. (11)] (~ reflection/L ²) 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 ⁵]. For reference, the corresponding asymptotic power laws are shown as dashed lines. Right inset: ℜ[E_z] field pattern for the (point) source at the origin (blue/white/red = positive/zero/negative).

Download Full Size | PDF

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], 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 c_n decrease asymptotically as 1/n ^(d+1) [35, 36] (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 ~|c_r|² 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 σ₀ in Eq. (10). The coupling coefficient M is a matrix element proportional to the rate of change of the materials [10], which in this case is $\frac{\partial \sigma }{\partial u}=s^{\prime }\left(u\right){\sigma }_0~\frac{1}{L}$. Therefore, the reflection scales as |M|²/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 (β₀,ω₀) is approached [33]. As a quadratic-shaped band-edgeω -ω ₀~(β -β₀)² is approached, the group velocity ${\nu }_g=\frac{\mathrm{d\omega }}{\mathrm{d\beta }}$ scales proportional to β -β₀, while the Δβ between the forward and reflected modes is 2(β -β₀) ~ v_g. Also, the coupling coefficient M is proportional to 1/v_g because of the constant-power normalizations of the incident and reflected modes [10, 33]. Hence, by inspection of Eq. (13), the reflection |c_r|² = O(v_g ^-(2d+4)). For example, the reflection is O(v_g ^-6)for a linear taper s(u) = u [33]. 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]).

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 s_n = s(nΔx/L) as L grows larger], there is still a proof of the adiabatic theorem (c_r → 0), but the only published proof is currently for the lossless case (unitary evolution) [31]. 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].) 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 = u^d), 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]; 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, 37]). 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, 39]: the result is that R(L) decays asymptotically as e ^-α√L for some constant α > 0 [38]. 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 k_x = 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 α.

Download Full Size | PDF

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 ₀, 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 ₀ to a larger value, corresponding to a larger σ₀ ~ lnR ₀ 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 ₀ 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], 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 ² in 1d uniform media, and the overall reflection is reduced by a factor of 3-400 compared to a fixed R ₀ = 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 ₀~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 ₀, 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 ², with the round-trip reflection either set to R ₀ = 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 ₀ 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

Download Full Size | PDF

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 ⁶.) 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.)