A Novel Boundary Element Method Using Surface Conductive Absorbers for Full-Wave Analysis of 3-D Nanophotonics

A Novel Boundary Element Method Using Surface Conductive Absorbers for Full-Wave Analysis of 3-D Nanophotonics

Lei Zhang, Student Member, IEEE, Jung Hoon Lee, Ardavan Oskooi, Amit Hochman, Member, IEEE,
Jacob K. White, fellow, IEEE and Steven G. Johnson
Manuscript received July 22, 2019.L. Zhang, A. Hochman and J. K. White are with the Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA, 02139, USA. (Email: zhangl@mit.edu; hochman@mit.edu; white@mit.edu)J. H. Lee is with Applied Computer Science and Mathematics, Merck & Co., Inc., Rahway, NJ, 07065, USA (Email: jung_hoon_lee@merck.com)A. Oskooi is with the Center for Materials Science and Engineering, Massachusetts Institute of Technology, Cambridge, MA, 02139, USA. (Email: ardfar@mit.edu)S. G. Johnson is with the Department of Mathematics, Massachusetts Institute of Technology, Cambridge, MA, 02139, USA. (Email: stevenj@math.mit.edu)

Fast surface integral equation (SIE) solvers seem to be ideal approaches for simulating 3-D nanophotonic devices, as these devices generate fields both in an interior channel and in the infinite exterior domain. However, many devices of interest, such as optical couplers, have channels that can not be terminated without generating reflections. Generating absorbers for these channels is a new problem for SIE methods, as the methods were initially developed for problems with finite surfaces. In this paper we show that the obvious approach for eliminating reflections, making the channel mildly conductive outside the domain of interest, is inaccurate. We describe a new method, in which the absorber has a gradually increasing surface conductivity; such an absorber can be easily incorporated in fast integral equation solvers. Numerical experiments from a surface-conductivity modified FFT-accelerated PMCHW-based solver are correlated with analytic results, demonstrating that this new method is orders of magnitude more effective than a volume absorber, and that the smoothness of the surface conductivity function determines the performance of the absorber. In particular, we show that the magnitude of the transition reflection is proportional to , where is the absorber length and is the order of the differentiability of the surface conductivity function.


nanophotonics, surface conductive absorber, boundary element method, surface integral equation, reflections.

I Introduction

In this paper we describe an absorber technique for terminating optical waveguides with the surface integral equation (SIE) method, which otherwise has difficulties with waveguides and surfaces extending to infinity. In order to attenuate waves reflected from truncated waveguides, we append a region with surface absorption to the terminations, as diagrammed in Fig. 1. The transition between the non-absorbing and absorbing regions will generate reflections that can be minimized by making the transition as smooth as possible. We show how this smoothness can be achieved with the surface absorber by smoothly changing integral-equation boundary conditions. Numerical experiments demonstrate that the reflections of our method are orders of magnitude smaller than those of straightforward approaches, for instance, adding a volume absorptivity to waveguide interior. In addition, we show the asymptotic power-law behavior of transition reflections as a function of the length of the surface absorber and demonstrate that the power law is determined by the smoothness of the transition [1].

Fig. 1: Schematic diagram of a photonic device with input and output waveguide channels, which must be truncated in a surface integral equation method.

Unlike the finite-difference or the finite-element volume-discretization methods, SIE methods treat infinite homogeneous regions (and some other cases) analytically via Green’s functions, and therefore often require no artificial truncation of space. Because SIE methods only require surfaces to be discretized, they can be computationally efficient for problems involving piecewise homogeneous media. A popular SIE method is the boundary-element method (BEM) [2, 3, 4, 5, 6], particularly since the development of fast solvers [7, 8, 9, 10]. However, a truncation difficulty arises with unbounded surfaces of infinitely extended channels common in photonics. Fig. 1 is a general photonic device schematic with input and output waveguide channels. In order to accurately simulate and characterize the device, such as calculating its scattering parameters, formally, the transmission channel must be extended to infinity, requiring infinite computational resources. A more realistic option is to truncate the domain with an absorber that does not generate reflections.

The key challenge is to design an absorber that both has small reflections and is also easily incorporated into an SIE solver. The best-known absorber is a perfectly matched layer (PML) [11, 12, 13, 14, 15], but a PML corresponds to a continuously varying anisotropic absorbing medium, whereas SIE methods are designed for piecewise homogeneous media. A similar problem arises if one were to simply add some absorption within the waveguides—in order to minimize transition reflection, the absorption would need to increase gradually from zero [1], corresponding again to inhomogeneous media. One could also use a volume integral equation (VIE) [16] or a hybrid finite-element method in the inhomogeneous absorbing region, but then one would obtain numerical reflections from the discontinuous change in the discretization scheme from the SIE to the VIE. Moreover, it has been proposed that an integral-equation PML can be obtained by varying the Green’s function instead of the media [17], but a continuously varying Green’s function greatly complicates panel integrations and fast solvers for SIE methods.

With SIE methods, one could, in principle, implement a scalar waveguide absorption in a piecewise homogeneous fashion as a discontinuous increase in absorption, but this will obviously generate large reflections due to the discontinuity of the medium and also the numerical reflection due to the discretization of the interfaces. We demonstrate this with a finite rectangular waveguide, to which an absorber with constant volume electrical conductivity and magnetic conductivity is attached, as shown in Fig. 2. The longitudinal cross-section of this arrangement is shown in Fig. 3(a). To achieve small reflections, the intrinsic impedance of the absorber is matched to that of the waveguide. Hence, the volume electrical conductivity and the volume magnetic conductivity should satisfy , where and are the permittivity and permeability of both the waveguide and absorber media, respectively. We quantify the reflection by use of the standing wave ratio (SWR), the ratio of the maximum field magnitude to the minimum field magnitude in the standing-wave region, evaluated on the waveguide axis. From the SWR, a reflection coefficient is then readily obtained as in a conventional transmission line. We calculated the field reflection coefficients in this way for a range of and , and the smallest reflection coefficient obtained was . This value, which is listed in the first column of Table I, is unacceptable for many design applications. In particular, the design of tapers [18] requires field reflection from terminations in the order of or smaller. Evidently, a more sophisticated way of terminating waveguides is called for.

In this paper, we examine an alternative approach to absorbers, adding electrical conductivity to the waveguide surface rather than to the volume, via a delta-function conductivity on the absorber surface, as shown in Fig. 3(b). The absorber’s interior medium remains the same as the waveguide’s, thus eliminating the need to discretize the waveguide-absorber interface, shown as a solid line in Fig. 3(a) and a dashed line in Fig. 3(b). This surface-conductivity strategy permits an efficient surface-only discretization, but at the same time allows for a smoothly increasing surface conductivity, thereby reducing transition reflections. Specifically, surface conductivity is easily implemented in SIE methods as it corresponds to a jump discontinuity in the field boundary conditions at the absorber surface. Since the SIE explicitly discretizes the surface boundary, continuously varying the field boundary conditions is easily implemented.

As a preview of results to be described below, in Fig. 4, we show the numerically computed complex magnitudes of the electric field along the direction inside a rectangular waveguide with several different surface absorbers. The surface absorber is in the region where in which is the position of the interface and is the absorber length. The surface electrical conductivity in this region is given by , where for constant, linear and quadratic profiles. The constant is chosen so that the total attenuation over the length of the absorbing region matches that of the optimal volume absorber above. The approach for calculating the attenuation along the absorber will be explained in section IV-A. In the complex magnitude plots of Fig. 4, the peak-to-peak magnitudes of ripples are an indication of the amount of reflections. As is easily seen in Fig. 4, there are substantial reflections when using a constant conductivity, smaller reflections when using a linearly increasing conductivity, and almost no reflections for a quadratically increasing conductivity. The magnitudes of the field reflection coefficients are listed in Table I, and show that the reflection coefficient for the quadratically varying surface conductivity is nearly one thousand times smaller than the reflection coefficient for the volume absorber. The results for the constant, linearly and quadratically varying surface absorber verify the results in [1], that the smoothness of any transition in the waveguide largely determines the resulting reflection.

This paper is organized as follows. In section II, the BEM formulation incorporating surface conductivity is derived. In section III, the decay rate due to the surface conductivity is examined using both numerical experiments and calculations using perturbation theory and Poynting’s theorem. In section IV, the asymptotic behavior of the transition reflection with respect to absorber length is presented. In section V, we note the existence of radiation modes originating from the excitation source mismatch with the waveguide mode, and show that the radiation complicates the interpretation of certain numerical results, but does not affect the performance of the surface absorber. Details of the numerical implementation of the SIE solver are given in Appendix, specifically, the detailed matrix construction, acceleration and preconditioning techniques.

Fig. 2: A discretized dielectric waveguide with an absorber attached.
(a) A waveguide with a volume absorber
(b) A waveguide with a surface absorber
Fig. 3: The 2-D longitudinal section of a waveguide with an absorber. The lengths of the waveguide and absorber are and , respectively, with denoting the wavelength in the waveguide medium. The waveguide cross section size is . The relative permittivities of the waveguide (silicon) and the external medium (air) are and , respectively.
(a) constant surface conductivity
(b) linear surface conductivity along direction
(c) quadratic surface conductivity along direction
Fig. 4: The complex magnitude of the electric field inside a waveguide and a surface absorber. The dashed line indicates the position of the waveguide-absorber interface.
Absorber type Volume Surface
Conductivity profile Constant Constant Linear Quadratic
TABLE I: The Standing wave ratio (SWR) and field reflection versus the conductivity distribution of the absorber

Ii BEM formulations with the Surface Conductive Absorber

In this section, we describe the 3-D BEM formulation for a waveguide truncated with a surface conductive absorber. Fig. 3(b) shows the plane cross-section of an -directed truncated rectangular waveguide, where the surface conductive absorber region is to the right of the dashed line. The permittivity and permeability of the waveguide interior and the exterior media are denoted as , and , , where the subscripts and denote the exterior and the interior, respectively. The electrical surface conductivity is subscripted with as a reminder that only electrical conductivity is being considered, though the generalization of what follows to both electrical and magnetic conductivity could be considered. As is described in section III, using only electrical conductivity can have a saturation phenomenon that can be avoided at the cost of using a longer absorber. The system is excited by a Gaussian beam propagating in direction. The Gaussian beam is generated by a dipole in a complex space [19], where the real part of the dipole position is inside the waveguide, from the left end. In this paper, the convention of the time-harmonic mode is adopted.

In SIE methods, for computing time-harmonic solutions, the unknowns are surface variables. In our case, we use surface electric and magnetic currents on both the interior, and , and exterior, and , of every surface. The currents on surfaces with on the left side of Fig. 3(b), satisfy a simpler set of equations than the currents on surfaces where , the right side of Fig. 3(b). When appropriate, we distinguish between the and currents with superscripts and , respectively.

Invoking the equivalence principle [20] yields a relation between surface currents and fields


where is an exterior-directed normal unit-vector, , , are the electric and magnetic fields of the Gaussian beam in a homogeneous space with material parameters equal to those of the waveguide interior, and , and , are electric and magnetic fields due to the equivalent currents in the exterior and interior, respectively.

On the surface of the waveguide, the continuity of the tangential components of the electric and magnetic fields yields the well-known PMCHW formulation [2, 3],


where and are integral operators described in the Appendix. From (1)-(4) and the tangential field continuity in the surface conductivity free region, the equivalent currents on the two sides of the waveguide surface are of the equal magnitude, but are opposite in direction. Specifically,


Thus, the unknown currents on the waveguide side are reduced to and .

For the surfaces where , a modified surface formulation is needed, one that incorporates the discontinuity due to the surface conductivity. When , the tangential electric field is still continuous across the absorber surface, and therefore


The tangential magnetic field is not continuous as a sheet of surface electric current is induced due to the electrical surface conductivity, thus creating a jump. Therefore,


where is the tangential electric field on the absorber surface, and could choose the field on either side of the absorber surface according to the enforced equality in (9).

As a result of (1), (3) and (9), the interior and exterior magnetic currents can be represented with a single variable


but the discontinuity of tangential magnetic field implies , and (2), (4), and (10) must be combined to generate a local equation


Finally, using the integral operator relation between and and , and substituting into (10), (12), yields


The unknowns, the surface currents , , , and , can be determined by solving equations (5), (6) on the left waveguide surfaces and (9), (13), (14) on the right absorber surfaces. The resulting linear system can be constructed and solved as described in Appendix. It is possible to simplify the formulation by eliminating one extra variable of electric currents on the absorber surface, though there are some accuracy issues [21].

Iii The decay rate of the fields

Fig. 5: The 2-D longitudinal section of a waveguide with uniform surface conductivity. The waveguide length is and cross section size is . The relative permittivity of the waveguide and external medium are 11.9 and 1, respectively.

In this section the exponential decay rate of waves propagating through the surface absorber region is analyzed. We demonstrate the relation between decay rate and surface conductivity using the example of a single dielectric waveguide with uniformly distributed surface conductivity. The longitudinal cross-section is shown in Fig. 5. The behavior of interior fields generated by a Gaussian beam source were computed using a BEM method based on solving (9), (13) and (14).

The plots in Fig. 6 show the complex magnitudes of the electric fields along the axis inside the waveguide for two cases, S and S. As expected, the complex magnitude decays exponentially with distance from the source with a surface-conductivity dependent rate. Also, as can be seen, waves reflect back from the right end and presumably these reflections decay as they travel to the left.

Fig. 6: The complex magnitude of the electric field along inside the waveguide in Fig. 5 with uniform surface conductivity.

An approximation to the rate of exponential decay can be determined by fitting the field plots. The fitted decay rates for a range of are shown in Fig. 7 and denoted with a dashed star curve. The decay rate does not monotonically increase with the surface conductivity. The curve shape can be explained as follows. When is small, the propagating wave is able to penetrate the lossy surface and is absorbed, with the absorption increasing with as expected. However, for large , the surface conductor itself becomes reflecting, forming essentially an enclosed metallic waveguide; as the tangential electric field vanishes at the surface and therefore there is no absorption. The practical implications of this upper bound on effective values for are limited, and are described in section IV.

The following sections introduce two alternative approaches to calculate the decay rate from surface conductivities, to confirm and further illustrate the numerical observations above.

Fig. 7: A comparison of three methods for computing the rate of field exponential decay along the propagation direction versus surface conductivity.

Iii-a Decay rate calculation by perturbation theory

In this section, a first-order closed-form decay rate formula, valid for low surface conductivity, is derived using perturbation theory.

Assume the electric field of the fundamental mode of a lossless rectangular waveguide is given, and the superscript denotes the unperturbed quantity. When electrical surface conductivity is put on the waveguide surface, it is equivalent to a perturbation of permittivity, denoted as , where is the Dirac delta function across the waveguide surface. According to [20], [22], the first-order variance of angular frequency due to the perturbation of permittivity is


where is the whole volume domain and the superscript denotes a first-order approximation. After applying the triple product rule to partial derivatives of the three interdependent variables , , , we obtain a first-order change in propagation constant due to the frequency change in (15), denoted as , where is the group velocity, as in which the subscript indicates is held fixed. Combining (15) and equations of and above, the integral in the numerator of (15) is reduced to a surface integral of the tangential components of the electric field, therefore


where is the surface of the waveguide. As expected, the perturbation in the propagation constant is imaginary, which in turn represents the decay rate . With a uniform cross section, the volume and surface integrals in (16) can be further reduced to surface and line integrals on the cross section, respectively. The electric field before the perturbation , along with , can be obtained numerically, for example, with a planewave method [23].

The decay rate calculated using the perturbation is plotted in Fig. 7 with a dashed diamond curve. Note that the curve overlaps with decay rates computed with other methods when the surface conductivity is small, and deviates for larger conductivity as should be expected given the first-order approximation.

Iii-B Decay rate calculation using Poynting’s theorem

Fig. 8: Illustration of the approach using Poynting’s theorem to calculate the decay rate of a waveguide with surface conductivity. The plot of the surface conductivity distribution along the longitudinal direction is aligned with the waveguide.

The perturbation method above predicts the decay rate when the surface conductivity is small. An alternative approach, based on Poynting’s theorem, can be used to calculate the decay rate for the entire range.

Figure 8 shows a waveguide with surface conductivity and also plots the conductivity function with -axis aligned with the -axis of the waveguide. Since this approach requires integrating the fields of source-excited propagating modes in the exterior region, some inevitably excited modes, such as radiation modes that will be discussed in section V, must be suppressed. For this reason, the surface conductivity starts with a large constant value ( S). The large surface conductivity leads to the saturation as seen in Fig. 7, and therefore, the Gaussian-beam source excites metallic waveguide modes at the beginning, propagating in direction in the closed interior region. The surface conductivity is then smoothly reduced to a smaller value, with which the decay rate is to be calculated. In this way, the metallic waveguide modes can smoothly change to the desired decaying dielectric waveguide modes with minimal radiation modes excited.

The decaying propagating mode in the domain of interest, shown in Fig. 8, is assumed in the form


where is the real propagation constant, and is the unknown decay rate. The Poynting vector in frequency domain is given by , and together with the assumed forms of in (17) and in (18), the derivative with respect to is given by


As illustrated in Fig. 8, apply Poynting’s theorem in the closed volume


where is the surface of the volume , is an exterior-directed normal unit-vector, and is the waveguide surface within . In the limit as , the closed integration surface becomes a surface , and one component of the integrand of the left side of (20) becomes . Combining (19) and this limit of (20) yields a closed-form representation of the decay rate


where denotes the boundary of the surface , and denotes the integral line on the waveguide surface within the surface .

The decay rate calculated using (21) is plotted in Fig. 7 with a dashed circle curve. It shows good agreement with the decay rate computed using fitting for the entire range, verifying that the surface conductivity is handled correctly by the BEM in accordance with Maxwell’s equations.

Iv Reflections and Asymptotic Convergence

In section I, we showed that a smoothly varying surface conductive absorber easily implemented in the SIE method is orders of magnitude more effective at eliminating reflections than a volume absorber of comparable length. And, since computational cost increases with the length of the absorber, it is worth examining the relationship between reflections and the absorber length. In order to do this, we first require a more careful classification of reflections and a more sensitive numerical measure of reflections than the standing wave ratio method. In the following subsections, we will define round-trip and transition reflections and present a measure of the power-law asymptotic convergence of transition reflections.

Iv-a Round-trip and transition reflections

Reflections in a domain of interest can be divided into a round-trip reflection, , and a transition reflection, . The round-trip reflection is generated by waves entering into the absorber, propagating to the end without being completely absorbed, reflected off the end of the absorber, and eventually propagating back into the waveguide. The round-trip reflection coefficient is proportional to


where the coefficient is determined by the reflection at the end of the absorber and is the absorber length. A factor of in the exponent of (22) represents the effect of the round trip, and another factor of indicates that the power is considered.

The transition reflection is the reflection generated by the change in material properties at the waveguide-absorber interface. Therefore, smoother material change at the interface produces smaller transition reflection as predicted by coupled-mode theory [1].

If the round-trip reflection in (22) is held fixed, the decay rate is inversely related to the absorber length . In another words, for a longer absorber, a smaller decay rate yields the same round-trip reflection. Since is proportional to the surface conductivity for small , using a longer absorber implies a smaller transition reflection. It was further shown in [1] that, given a fixed round-trip reflection, the transition reflection decreased as a power law with increasing absorber length . The power-law exponent is determined by the order of the differentiability of the medium (conductivity) function. Suppose a conductivity function along is , where the interface is at , and is a monomial function with order in the region


With the surface conductivity function and a fixed round-trip reflection, the asymptotic behavior of the transition reflection in terms of the length of the absorber is


The power-law behavior in (24) indicates that, with a higher-order conductivity function, the transition reflection decreases faster with increasing the absorber length. It does not follow that should be made arbitrarily large, however, there is a tradeoff in which increasing eventually delays the onset of the asymptotic regime in which (24) is valid [1].

Iv-B Asymptotic Convergence with

We present numerical results to verify the asymptotic power-law convergence of the transition reflection of the surface absorber. Because it is hard to explicitly measure the transition reflection in the integral equation method, instead, we measure alternative electric field expressions as in [1].

First, we define as the electric field at a fixed position in the waveguide when the length of the absorber is , with unit . Thus includes the incident field, the round-trip reflection and the transition reflection. With a small fixed round-trip reflection, the difference of and is the difference of the transition reflections, which in the limit of large approaches to zero, so . Therefore, , and similarly , can be a measure of transition reflection, and specifically, derived from (24), they are subject to the following asymptotic behavior


In the example of a rectangular waveguide attached with a surface absorber, the asymptotic convergence of (25) is shown in Fig. 9(a) in a log-log scale for constant, linear, quadratic, and cubic surface conductivity profiles (). The curve for slowly approaches to the curve. The other three curves align with the expected curves respectively. The curves would eventually converge to a small quantity, which is the difference of the small round-trip reflections due to a phase difference. Fig. 9(b) shows the asymptotic convergence of (26) with the same waveguide example for constant, linear, quadratic, and cubic profile surface conductivities. Again, the curve for the constant conductivity profile converges slowly, and is expected to be aligned with the curve. The other three power-law convergence are the same as the predicted . The agreement of the figures verifies the power-law behavior of the transition reflection of the surface absorber and further illustrates that a higher-order conductivity function leads to smaller transition reflections.

(a) The asymptotic convergence of the transition reflection measured using .
(b) The asymptotic convergence of the transition reflection measured using .
Fig. 9: Asymptotic power-law convergence of the transition reflection with the length of the surface absorber. The length of the waveguide is , with denoting the wavelength in the waveguide medium. The waveguide cross section size is . The relative permittivities of the waveguide (silicon) and the external medium (air) are and , respectively.

V Radiation in the surface absorber

Fig. 10: The complex magnitude of the electric field inside a waveguide and a long surface absorber excited by a dipole source and a Gaussian beam, respectively. The lengths of the waveguide and the absorber are and , respectively. The surface conductivity on the absorber increases quadratically.

A surface absorber with a conductivity that rises smoothly with distance from the waveguide-absorber interface would be expected to have interior fields whose magnitude decays with distance at an accelerating exponential rate. Instead, the field magnitudes decay exponentially near the interface, but then switch to an decay, as shown in Fig. 10. The cause of this switch in decay rate is due to coupled radiation. The radiation is generated by the inevitable mismatch between the excitation source and waveguide modes. The situation is similar to a source radiating in a lossy half-space, in which the dominant field contribution is due to a lateral wave that decays algebraically [24, 25]. Fig. 10 shows the complex magnitude of the interior electric field in a waveguide attached with a long quadratic-profile surface absorber. The waveguide system is excited by a dipole source and a Gaussian beam, respectively, located inside the waveguide. In the semilog plot, the dominant waveguide mode decays at an accelerating exponential rate, and then because the guided mode decays faster than the radiation, the radiation dominates at a certain distance from the waveguide-absorber interface. The coupled radiation results in an floor, as is clearly shown in the inset with a log-log scale. Note that using a Gaussian beam source results in a lower floor than using a dipole source, because the Gaussian beam is more directional, and generates less radiation.

The coupled radiation does not affect the performance of the surface absorber, because the coupled radiation itself is several orders of magnitude smaller than the propagating modes in the waveguide, and little will be reflected. The asymptotic convergence of the transition reflections in section IV and the reflection coefficient in section I consistently showed the excellent performance of surface absorbers, obviously unhampered by the effects of radiation.

Vi Conclusion

We presented a novel numerical technique, a surface conductive absorber that is easily combined with the boundary element method, to eliminate the reflections due to the truncation of infinite channels. We illustrated this technique with a dielectric optical waveguide and described the modified BEM formulation needed to allow varying the surface conductivity. We further discussed the non-monotonically increasing decay rate with the surface conductivity and presented methods to calculate the decay rate using perturbation theory and Poynting’s theorem. We demonstrated that the surface conductive absorber is orders of magnitude more effective than the volume conductive absorber, and showed an asymptotic power-law convergence of the transition reflection with respect to the length of the absorber to verify that the smoothness of conductivity function determines the transition reflection, also noted in [1].

The major advantages of the surface conductive absorber are: (1) the varying surface conductivity is easily implemented in BEM and can significantly reduce the transition reflections; and (2) the volume properties of the absorber are the same as the waveguide, so there is no interior cross section to discretize, eliminating a potential source of numerical reflections.

In this paper, the surface conductive absorber is demonstrated with a simple example, a rectangular waveguide. Demonstrating that the method is as effective on more general structures is the subject of future work.

Vii Acknowledgement

This work was supported by the Singapore-MIT Alliance (SMA), the MARCO Focus Center on Interconnect, in part by Dr. Dennis Healy of DARPA MTO, under award N00014-05-1-0700 administered by the Office of Naval Research, and the Army Research Office through the Institute for Soldier Nanotechnologies under contract DAAD-19-02-D0002.

Viii Appendix

Viii-a Construction of a linear system

From section II, the five equivalent currents, , on the waveguide surfaces, and , , on the absorber surfaces, are to be determined by solving the equations (5), (6), (9), (13) and (14). The currents are approximated with the RWG basis function [4] on triangular-meshed surfaces,


where is the RWG function on the th triangle pair, and , are the corresponding coefficients for the electric and magnetic currents, respectively. These unknown coefficients of the five equivalent currents assemble a vector to be solved for, specifically,


Electric and magnetic fields are represented using the mixed-potential integral equation (MPIE)[5] for a low-order singularity, with integral operators and as in [6]


where is the intrinsic impedance, the subscript denotes the exterior or interior region, and the superscript denotes the waveguide or absorber surfaces. The integral operators on the th RWG function are given by


where is the surface of the th triangle pair, is the propagation constant in region , and is the free-space Green’s function in region


where and are target and source positions, respectively.

We employ Galerkin method [26] on the waveguide by using the RWG function as the testing function on target triangle pairs. The tested , operators on the th target triangle pair due to the th source triangle pair become


where is the surface of the th target triangle pair. Substituting the tested field operators into equations (5), (6) yields a matrix due to the currents , , and a matrix due to the currents , , .

On the absorber surfaces, the term in (13) and (14) requires another testing procedure in order to incorporate the surface conductivity


Similarly, substituting the above four tested integral operators into (9), (13) and (14) generate a matrix due to the currents , , and a matrix due to the currents , , .

Assembling the four matrices according to (5), (6), (9), (13) and (14) yields a dense linear system




The right-hand-side vector contains tested incident electric and magnetic fields


in which the th entry of each subvector is given by


Viii-B Acceleration and Preconditioning with FFT

The linear system (39) can be solved with iterative algorithms, for instance, GMRES for this non-symmetrically dense system. In each iteration, the matrix-vector product takes time, where is the number of unknowns. Moreover, to explicitly store the matrix is expensive, requiring memory. In fact, there have been many well-developed fast algorithms to reduce the costs of the integral equation solvers [7, 8, 9, 10]. In this paper, we use a straightforward and easily-implemented FFT-based fast algorithm to accelerate the SIE method on periodic guided structures.

As shown in Fig. 11, the waveguide is discretized into a periodically repeating set of the RWG triangle pairs. Due to the mesh periodicity and the space invariance of the operators (35), (36), the matrices and are block Toeplitz, requiring explicit calculation and storage of only a block row and a block column, reducing memory to approximately . A Toeplitz matrix can be embedded in a circulant matrix, and the circulant matrix-vector product can be computed with the FFT [27, 28]. In this way, the computational costs are reduced approximately to .

Because the surface conductivity must vary with distance from the waveguide-absorber interface, the tested potential operators (37), (38) are not space invariant. Therefore, accelerating the parameterized matrices and is not as straightforward as and . Typically, the integral of (37), (38) is calculated numerically using Gauss quadrature [29], summing up the tested potentials at quadrature (target) points with Gauss weights. The space invariance of the potential operators (32), (33) and periodicity of the mesh allows assembling a matrix of the potentials at target points by explicitly calculating only a block row and block column. Then the potentials at target points are summed after testing and multiplications with Gauss weights and surface conductivity, and eventually stamped into matrices and .

Fig. 11: A discretized waveguide with a periodic unit.

Another great advantage of working with a Toeplitz or a block Toeplitz matrix is the existence of a highly efficient preconditioner [30, 31, 32]. A circulant matrix is approximated from the Toeplitz matrix, and then can be easily inverted with the FFT. We use this method to calculate a preconditioner for , and use the block-diagonal preconditioner [10] for .


  • [1] A. F. Oskooi, L. Zhang, Y. Avniel, and S. G. Johnson, “The failure of perfectly matched layers, and towards their redemption by adiabatic absorbers,” Optics Express, vol. 16, pp. 11 376–11 392, 2008.
  • [2] R. F. Harrington, “Boundary integral formulations for homogeneous material bodies,” Journal of Electromagnetic Waves and Applications, vol. 3, no. 1, pp. 1–15, 1989.
  • [3] K. Umashankar, A. Taflove, and S. M. Rao, “Electromagnetic scattering by arbitrary shaped three-dimensional homogeneneous lossy dielectric objects,” IEEE Trans. Antennas Propagat., vol. 34, no. 6, pp. 758–765, June 1986.
  • [4] S. M. Rao, D. R.Wilton, and A. Glisson, “Electromagnetic scattering by surfaces of arbitrary shape,” IEEE Trans. Antennas Propagat., vol. 30, pp. 409–418, May 1982.
  • [5] C. F. Wang, F. Ling, and J. M. Jin, “A fast full-wave analysis of scattering and radiation from large finite arrays of microstrip antennas,” IEEE Trans. Antennas Propagat., vol. 46, pp. 1467–1474, Oct. 1998.
  • [6] X. Q. Sheng, J.-M. Jin, J. Song, W. C. Chew, and C.-C. Lu, “Solution of combined-field integral equation using multilevel fast multipole algorithm for scattering by homogeneous bodies,” IEEE Trans. Antennas Propagat., vol. 46, no. 11, pp. 1718–1726, 1998.
  • [7] J. R. Phillips and J. K. White, “A precorrected-FFT method for electrostatic analysis of complicated 3-D structures,” IEEE Trans. Computer-Aided Design, vol. 16, pp. 1059–1072, Oct. 1997.
  • [8] Z. Zhu, B. Song, and J. White, “Algorithms in Fastimp: a fast and wideband impedance extraction program for complicated 3-D geometries,” in Proceedings of DAC, 2003, pp. 712–717.
  • [9] L. Zhang, N. Yuan, M. Zhang, L. W. Li, and Y. B. Gan, “RCS computation for a large array of waveguide slots with finite wall thickness using the mom accelerated by p-fft algorithm,” IEEE Trans. Antennas Propagat., vol. 53, no. 9, pp. 3101–3105, 2005.
  • [10] J. Song, C.-C. Lu, and W. C. Chew, “Multilevel fast multipole algorithm for electromagnetic scattering by large complex objects,” IEEE Trans. Antennas Propagat., vol. 45, no. 10, pp. 1488–1493, 1997.
  • [11] J. P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comp. Phy., vol. 114, pp. 185–200, 1994.
  • [12] W. C. Chew and W. H. Weedon, “A 3D perfectly matched medium from modified Maxwell’s equations with stretched coordinates,” Microwave Opt. Technol. Lett., vol. 7, no. 13, pp. 599–604, 1994.
  • [13] Z. S. Sacks, D. M. Kingsland, R. Lee, and J. F. Lee, “A perfectly matched anisotropic absorber for use as an absorbing boundary condition,” IEEE Trans. Antennas Propagat., vol. 43, no. 12, pp. 1460–1463, 1995.
  • [14] S. D. Gedney, “An anisotropic perfectly matched layer-absorbing medium for the truncation of FDTD lattices,” IEEE Trans. Antennas Propagat., vol. 44, no. 12, pp. 1630–1639, 1996.
  • [15] A. Taflove, Computational Electrodynamics: the Finite-difference Time-Domain Method.     Artech House, 1995.
  • [16] D. H. Schaubert, D. R. Wilton, and A. W. Glisson, “A tetrahedral modeling method for electromagnetic scattering by arbitrarily shaped inhomogeneous dielectric bodies,” IEEE Trans. Antennas Propagat., vol. 32, no. 1, pp. 77–85, 1984.
  • [17] D. Pissoort and F. Olyslager, “Termination of periodic waveguides by PMLs in time-harmonic integral equation-like techniques,” IEEE Antennas and Wireless Propagat. Letters, vol. 2, pp. 281–284, 2003.
  • [18] M. L. Povinelli, S. G. Johnson, and J. D. Joannopoulos, “Slow-light, band-edge waveguides for tunable time delays,” Optics Express, vol. 13, no. 18, pp. 7145–7159, 2005.
  • [19] E. Erez and Y. Leviatan, “Electromagnetic scattering analysis using a model of dipoles located in complex space,” IEEE Trans. Antennas Propagat., vol. 42, pp. 1620–1624, 1994.
  • [20] R. F. Harrington, Time-Harmonic Electromagnetic Fields.     McGraw-Hill Book Company, 1961.
  • [21] L. Zhang, S. G. Johnson, and J. K. White, “Calculation of Mie scattering from a sphere with surface conductivity using the boundary element method,” In preparation, 2010.
  • [22] J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals: Molding the Flow of Light, 2nd ed.     Princeton University Press, 2008.
  • [23] S. G. Johnson and J. D. Joannopoulos, “Block-iterative frequency-domain methods for Maxwell’s equations in a planewave basis,” Optics Express, vol. 8, no. 3, pp. 173–190, 2001.
  • [24] W. C. Chew, Waves and Fields in Inhomogeneous Media.     New York: Van Nostrand, 1990.
  • [25] A. Ishimaru, Electromagnetic Wave Propagation, Radiation, and Scattering.     Englewood Cliffs, New Jersey: Prentice-Hall, 1991.
  • [26] R. F. Harrington, Field computation by moment methods.     Macmillan, 1968.
  • [27] C. V. Loan, Computational Frameworks for the Fast Fourier Transform.     Society for Industrial and Applied Mathematics, 1992.
  • [28] Y. Zhuang, K. L. Wu, C. Wu, and J. Litva, “A combined full wave CGFFT method for rigorous analysis of large microstrip antenna arrays,” IEEE Trans. Antennas Propagat., vol. 44, pp. 102–109, Jan. 1996.
  • [29] L. N. Trefethen and D. Bau, Numerical Linear Algebra.     SIAM, 1997.
  • [30] G. Strang, “A proposal for Toeplitz matrix calculations,” Studies in Applied Mathematics, vol. 74, pp. 171–176, 1986.
  • [31] T. F. Chan, “An optimal circulant preconditioner for Toeplitz systems,” SIAM journal on scientific and statistical computing, vol. 9, no. 4, pp. 766–771, 1988.
  • [32] T. F. Chan and J. A. Olkin, “Circulant preconditioners for Toeplitz-block matrices,” in SIAM Conference on Linear Algebra in Signals, Systems, and Control, 1990, pp. 89–101.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description