## Abstract

We consider the propagation of guided waves in a planar waveguide that has smooth walls except for a finite length segment that has random roughness. Maxwell’s equations are solved in the frequency domain for both TE and TM polarization in 2-D by using modal expansion methods. Obtaining numerical solutions is facilitated by using perfectly matched boundary layers and the **R**-matrix propagator.

Varying lengths of roughness segments are considered and numerical results are obtained for guided wave propagation losses due to roughness induced scattering. The roughness on each waveguide boundary is numerically generated from an assumed Gaussian power spectrum. The guided waves are excited by a Gaussian beam incident on the waveguide aperture. Considerable numerical effort is given to determine the stability of the algorithm.

© 2001 Optical Society of America

## 1 Introduction

The work presented here has been motivated by the development of waveguide components for use in fiber optic gyros. The fiber optic gyro requires that two counter rotating guided waves be the same single polarization and mode. These “clean” guided modes are created by a specially designed waveguide that supports and allows only a single mode and polarization to emerge from the output side. The output intensity must be sufficient and this is where loss mechanisms are important. Since the waveguide has wall roughness, radiative scattering losses can occur and the question is how much roughness can be tolerated. Another possibility is that roughness scattering is secondary to scattering due to bulk material inhomogeneity of the waveguide core. In this work scattering from material inhomogeneity is not considered, rather we calculate radiative losses per unit length of wall roughness.

The method used in this work is flexible so that many planar waveguide problems can be considered with either metallic or dielectric walls. This includes step discontinuities, periodic or random roughness, multiple waveguide configurations, multiplexing, beam combining, etc. Using this method, it is straightforward to consider the case of additional scattering due to waveguide core inhomogeneity and this work is in progress.

We numerically solve Maxwell’s equations in the frequency domain for propagation in planar waveguides. Waveguide modes are initiated by a Gaussian beam incident on the waveguide aperture. This is illustrated by the schematic in Fig. 1. Part of the Gaussian incident energy that is not reflected is coupled into one or more guided wave modes that propagate down a smooth portion of the waveguide that extends from *z*=*L*
_{t} to *z*=*L*
_{r}. At the end of the smooth segment begins a segment of roughness that has length extending from *z*=*L*
_{t} to *z*=0. The length of the smooth portion is chosen so that at *z*=*L*
_{r}, the only remaining electromagnetic energy propagating downward is a guided wave mode. While transiting the roughness region, some of the guided wave energy is lost to scattering and after exiting the roughness region at *z*=0, the remaining guided wave mode energy continues to propagate unimpeded in the semi-infinite smooth substrate region (*z*<0) of the waveguide. The 1-D roughness is numerically generated with a Gaussian correlation function.

In Fig. 1, the primary computational domain includes |*x*|<*L*
_{x}/2an d 0≤*z*≤*L*
_{t} and solutions to Maxwell’s equations within this region are obtained in the form of modal expansions that are individually valid over contiguous z-invariant portions within 0≤*z*≤*L*
_{t}. The roughness region 0≤*z*≤*L*
_{r} is approximated by dividing this region into many thin layers each of which are z-invariant and such that the waveguide channel width within each layer is dependent on the roughness and this is discussed in more detail below. All solutions obtained within the z-invariant regions are then related by boundary conditions and the numerically stable **R**-matrix algorithm [1]–[5]. In addition, solutions in the region -∞<*z*≤*L*
_{t} are also obtained as modal expansions. Since the extent of the computational domain is finite, we include absorbing boundary conditions based on the perfectly matched layer (PML) [6, 7] formalism.

Many authors have investigated waveguide scattering losses [8]–[15]. Early theoretical work by Marcuse [8, 9] considered the coupling between guided and radiative wave modes. Based on these works, Payne and Lacey [10, 11] also provided a theoretical analysis of scattering losses from planar optical waveguides. Their analytical results are dependent on numerous waveguide characteristics including the power spectrum of wall roughness that is available to allow coupling into the radiative scattering field. Some experimental results are given in [12]–[14] and these data are also compared with theory similar to that given in [10, 11].

A theoretical approach different from that given in [8]–[11] is given in [15]. The focus of this work is to obtain a generalized scattering matrix of a cascaded set of abrupt discontinuities. This approach uses the generalized telegraph equation formulation and the modal matching technique. Unlike the present approach where the computational domain is enclosed by perfectly matched absorbing layers, [15] encloses the computational region with infinitely conducting walls. It is not clear how this approach would work if the scattering fields are allowed to reflect from the infinitely conducting walls. This approach has modeling similarities to that used in this work and also appears to be closely related to the **R**-matrix algorithm.

## 2 Theory

#### 2.1 Maxwell’s Equations

We begin with Maxwell’s equations for the electric *E*⃗ and magnetic *B*⃗ fields in the frequency domain.

These equations contain position dependent permittivity *∊*(*x*, *y*, *z*) and permeability *µ*(*x*, *y*, *z*) that will be used to describe the waveguide structure including details of the perfectly matched layers and wall roughness. The angular frequency *ω* is related to the free space wavelength λ by *ω*/c=2π/λ. While not explicit in Eqs. (1) and (2), incorporating the PML formalism actually introduces anisotropy into the permittivity and permeability. However, the anisotropy is only a mathematical construct that is active within the PML absorbing layer region. Outside of these PML regions, the material parameters assume a scalar position-dependent form. In this work we assume non-permeable media and *µ*(*x*, *y*, *z*)=1 everywhere. We consider two-dimensional geometry and assume invariance in the ŷ direction. We then refer to polarization as either TE, where *E*⃗=(0,*E*
_{y}, 0) and *H*⃗=(*H*
_{x}, 0,*H*
_{z}), or TM, where *E*⃗=(*E*
_{x},0,*E*
_{z}) and *H*⃗=(0,*H*
_{y},0).

Reducing Eqs. (1) and (2) to two dimensions and incorporating the PML features [6, 7] yields the following coupled first order equations that can be written in real space as

For reasons that will become clear, Eqs. (3)–(8) have been specifically written in the form as shown and this will be discussed in more detail below. We have defined the quantities

where j=x or z. These position-dependent functions describe the details of the waveguide structure by their dependence on the scalar permittivity *∊*(*x*, *z*), the magnetic conductivity σ_{j}(*x*, *z*)* and the electric conductivity σ_{j}(*x*, *z*). The PML conductivity terms, σ_{j}(*x*, *z*) and σ_{j}(*x*, *z*)^{*}, control absorption for propagation of the electric and magnetic fields in the j-direction, respectively. The electric conductivity σ_{j}(*x*, *z*) is not explicit in Eqs. (9) and (10) because the relationship σ_{j}(*x*, *z*)=∊(*x*, *z*)${\mathrm{\sigma}}_{\mathrm{j}}^{*}$ (*x*, *z*) is required for impedance matching. Ideally, impedance matching eliminates reflection when the PML absorbing layers are encountered. Outside of a PML absorption region, the σ_{j}=${\mathrm{\sigma}}_{\mathrm{j}}^{*}$=0, and then Eqs. (3)–(8) revert to the usual form of Maxwell’s equations.

Equations Eqs. (3)–(5) can be obtained from equations Eqs. (6)–(8) (and vice versa) by making the replacements

For this reason, we will limit the remaining theory discussion to TM polarization only. The case for TE polarization is obtained with Eqs. (11) and (12).

We now convert Eqs. (3)–(5) to a k-space representation and make the transition to a finite sized problem. When these conversions are properly done, equations Eqs. (3)–(5) assume a form that has been shown to have excellent numerical convergence properties [16, 17]. This is in contrast to directly applying the same numerical approach as described here to equations Eqs. (3)–(8) where convergence and numerical stability are very poor. Before making the conversion to k-space, we make two approximations to Eqs. (3)–(5). First, we truncate the x-dimension of the solution space to a finite length *L*
_{x} and discretize the x-coordinate. Second, we assume that the material parameters in Eq. (9) and (10) are z-invariant.

We discretize the x-coordinate into *N* points as *x*
_{m}=*m*Δ*x* where integer *m*=-*N*/2 → *N*/2 and Δ*x*=*L*
_{x}/*N*. Likewise, we discretize the x-component of the wavevector *k*
_{n}=*n*Δ*k* where *n*=-*N*/2 → *N*/2 and Δ*k*=2*π*/*L*
_{x}. Now that we have discretized and truncated the problem, we use matrix-vector notation to represent the Fourier transform and inverse Fourier transform: **F** and **F**
^{-1}. Such transform relations are written as

where the **F** is a *N*×*N* matrix and f is a *N*-element column vector. Note that these are one dimensional transforms applied to the *k*
_{x}-*x* transform pair where *k*
_{x}≡*k*. Applying these transform operators to Eqs. (3)–(5) yields

where **k** is interpreted as a diagonal matrix with elements *k*
_{n} and the **H**
_{y}, **E**
_{x}, and **E**
_{z} are *N*-element column vectors with each element corresponding to *k*
_{n}. We have defined the *N*×*N* square matrices

where the quantities *µ*
_{j}(*x*), *ε*
_{j}(*x*), *β*
_{j}(*x*), and α_{j}(*x*) are interpreted as diagonal matrices with elements *µ*
_{j}(*x*
_{n}), etc.

Using Eqs. (14)–(16), we obtain the second-order differential equation

where **I** is the identity matrix. Since there are *N* discrete *k* values, Eq. (19) represents *N* coupled differential equations. When **M** is independent of z (this is the case when the material parameters are z-invariant), solution of this equation system is straightforward by diagonalization as **S**
^{-1}
**MS**=*ξ*
^{2}=*ξ*·*ξ*. The columns of the *N*×*N* matrix **S** are the eigenvectors of **M** and the diagonal matrix *ξ*
^{2} contains the eigenvalues of **M**. The solution to Eq. (19) is

and the electric field follows from Eq. (14) as

The *C*
_{±} are scalar constants and the *e*
^{±ξz} represent diagonal matrices where the n-th diagonal element is an exponential term as exp(±*ξ*
_{n}
*z*). The *ξ*
_{n} is the n-th diagonal root-eigenvalue element of *ξ*. In Eqs. (20) and (21) we now have solutions to Maxwell’s equations for an inhomogeneous solution space that has an x-dimension of length *L*
_{x} and has material parameters that are z-invariant.

## 2.2 R-matrix Propagator

In order to consider solution spaces that are not z-invariant, we assume that we can build up such structures from numerous z-invariant layers. Thus, for all the z-invariant layers in the computational domain, we need to relate the solutions for each layer as given in Eqs. (20) and (21). A numerically stable relationship is given by the **R**-matrix algorithm which is but one of a class of several types of transfer matrices. First, for solutions within a single z-invariant layer bounded from *z* → *z*+*δ*, a relationship is assumed to be given by

where *δ* is the thickness of the layer. Inserting Eqs. (20) and (21) into Eq. (22) yields the *N*×*N* r-matrices as

These r-matrices relate the electric and magnetic fields at the boundaries across a single z-invariant layer. Generalizing to the case of two or more contiguous layers bounded from *z* → *z*+*z*
_{t}, we use the **R**-matrices which have a similar form as

The parameter *z*
_{t} is the cumulative thickness of all the z-invariant layers as described in Eq. (22). The *N*×*N*
**R**-matrices are calculated recursively by

The recursive algorithm Eqs. (26)–(29) is obtained from Eq. (22) and Eq. (25) and imposing the continuity of **E**
_{x} and **B**
_{y} across the boundaries separating z-invariant layers. It is seen from Eq. (22) and Eq. (25) that for the first z-invariant layer when *z*
_{t}=*δ*
_{1}, we initiate the recursion with R
_{ij}
(*δ*
_{1})=r
_{ij}
(*δ*
_{1}). If the next z-invariant layer has thickness *δ*
_{2}, new r
_{ij}
(*δ*
_{2}) matrices are computed and the **R**
_{ij}
(*δ*
_{1}) are used in Eqs. (26)–(29) to yield **R**
_{ij}
(*δ*
_{1}+*δ*
_{2}). These recursive relations are repeated until the desired thickness is achieved and when this occurs, we have a relation between the electric and magnetic fields at two planar boundaries that encompass the structure of interest.

## 2.3 Superstrate and Substrate

Referring to Fig. 1, we want to relate solutions within 0≤*z*≤*L*
_{t} to solutions in the semi-infinite superstrate *z*≥*L*
_{t} and substrate *z*≤0 and we can use Eq. (25) to provide such a relation. Setting *z*=0 and *z*
_{t}=*L*
_{r} in Eq. (25), we can relate to the superstrate and substrate fields by invoking the continuity of the tangential components. If we have incident and reflected fields in the superstrate and a transmitted field in the substrate, then these continuity conditions are

where the superscripts i, r, and t denote incident, reflected, and transmitted, respectively.

In Fig. 1 we see that the superstrate and substrate regions are z-invariant. In addition, the superstrate region *z*≥*L*
_{t} is homogeneous. In the superstrate region, we let a Gaussian beam be incident on the waveguide aperture and our goal will be to calculate the resulting fields at the *z*=0 and *z*=*L*
_{t} boundaries. This will allow us to further calculate the guided wave transmitted field after it has emerged from the roughness region and propagates into *z*<0. Although we do not do so here, we can also calculate the reflected field for *z*>*L*
_{t}.

## 2.3.1 Transmitted field z≤0

In this section, we discuss the ${\mathbf{E}}_{\mathrm{x}}^{\mathrm{t}}$(*k*, *z*) and ${\mathbf{H}}_{\mathrm{y}}^{\mathrm{t}}$(*k*, *z*) transmitted fields. Solutions of the form given by Eqs. (20) and (21) are still valid for the transmitted region by setting C-=0 and this yields

Since the region *z*≤0 has only downward propagating waves, the eigenvalues associated with Eqs. (32) and (33) must have ℜ*ξ*≥0 and *ℑξ*≤0. The latter equations in Eqs. (32) and (33) give the z-dependence after the guided wave exits the roughness region. These equations can be examined numerically after we have obtained the solution for ${\mathbf{H}}_{\mathrm{y}}^{\mathrm{t}}$(*k*, 0). Equations Eqs. (32) and (33) may be combined to give the relationship between the electric and magnetic transmitted fields in the inhomogeneous region *z*≤0 which is

and this will be used later.

## 2.3.2 Incident and Reflected Fields: *z*≥*L*_{t}

In the homogeneous superstrate region, the Gaussian incident field is assumed to be a superposition of plane waves [18] with half-width determined by *ωσ*/c and the angle of incidence centered about *θ*
^{i}. We write the x and y-components of the incident field as

$$\times \phantom{\rule{.9em}{0ex}}\mathrm{exp}\left[-{\left(\frac{\omega \sigma}{2c}\right)}^{2}{\left(\theta -{\theta}^{i}\right)}^{2}\right]\mathrm{exp}\left[i(\omega /c)\left(x\mathrm{sin}\theta -z\mathrm{cos}\theta \right)\right]$$

where the superstrate medium is assumed to be vacuum. As in Eqs. (30) and (31), we want the Fourier transforms of these fields and Fourier inversion of these components yields the column vectors ${\mathbf{E}}_{\mathrm{x}}^{\mathrm{i}}$(*k*, *z*) and ${\mathbf{H}}_{\mathrm{y}}^{\mathrm{i}}$(*k*, *z*) where the n-th element is given by

when *k*
_{n}<(*ω*/c) and zero when *k*
_{n}>(*ω*/c). The wavevector *k*
_{n}=2πn/*L*
_{x} and sin *θ*
_{n}=*k*
_{n}/(*θ*/c).

The reflected fields are written as a superposition of plane waves as

$$=\frac{1}{2\pi}\int \text{d}k\left(\begin{array}{c}{E}_{x}^{r}(k,z)\\ {H}_{y}^{r}(k,z)\end{array}\right)\mathrm{exp}\left(ikx\right)$$

where the z-component of the wavevector is $q=\sqrt{{(\omega /c)}^{2}-{k}^{2}}$. In the latter equation, the z-dependence has been absorbed into the Fourier coefficient. It is straightforward to show that ${E}_{\mathrm{x}}^{\mathrm{r}}$(*k*, *z*)=(*q*c/*ω*)${\mathrm{H}}_{\mathrm{y}}^{\mathrm{r}}$(*k*, *z*). With this and consistent with the notation of Eq. (31), we relate the column vectors associated with the reflected field as

where the n-th element of the diagonal matrix **Z**(k) is

## 2.3.3 Surface Fields: *z*=0 and *z*=*L*_{t}

We now have the relationships to calculate the fields at the upper and lower boundaries of the computational domain as shown in Fig. 1. Using Eqs. (30), (31), (34), (39), and (25) yields

and this linear equation system may be solved for the surface fields ${\mathbf{H}}_{\mathrm{y}}^{\mathrm{t}}$(*k*, 0) and ${\mathbf{H}}_{\mathrm{y}}^{\mathrm{r}}$(*k*,*L*
_{t}). Along with the eigenvectors **S** and eigenvalues *ξ* for the substrate region, the first of these two solutions may be then used in Eqs. (32) and (33) or Eq. (34) to obtain the transmitted fields for *z*<0.

## 3 Numerical Results

The theory discussed in Section 2i s applied to the problem of guided wave propagation in rough waveguides. Numerical results of transmission losses due to the wall roughness are given in Figs. 3–5. In Figs. 6–8 we present some results showing numerical convergence characteristics of the algorithm used in this work. In Fig. 9, we illustrate how effective the PML layers can be in preventing reflection from the hard limits of the computational domain. Finally, in Fig. 10, the transmission curves are fit to analytical model to estimate power loss per unit length of wall roughness.

## 3.1 General Comments

There are some parameters common to the numerical data and reference to Fig. 1 and Fig. 2 may be helpful in parts of the following discussion. All realizations for the waveguide wall roughness were numerically generated for a Gaussian autocorrelation function with a correlation length of 1.0λ and various values of rms roughness: 0.05λ, 0.1λ and 0.2λ. Four realizations have been used and numerical results for each realization are displayed separately rather than an attempt at ensemble averaging. Except in Fig. 10, the initial smooth portion of the waveguide has length *L*
_{t}-*L*
_{r}=3000λ and nominal waveguide channel width *ω*=1λ. This length *L*
_{t}-*L*
_{r} is sufficient to allow the remaining transmitted energy that reaches the *z*=*L*
_{r} plane to only be one or more guided wave modes. With this, only a guided wave mode is incident on the roughness region. The length *L*
_{r} of the roughness region is a parameter of study in this work that varies over many values. In the numerical algorithm, the length *L*
_{r}=*N*
_{r}
*δ* is the cumulative thickness of *N*
_{r} z-invariant layers each having thickness *δ*. This is illustrated in Fig. 2. Within each layer of thickness *δ*, the waveguide channel width is dependent on the roughness at each waveguide boundary. We have let *δ*=0.1λ for all calculations presented here. The waveguide permittivity values are ∊_{1}=(2.25, 0.) and ∊_{2}=(2.50, 0.). In most cases, we assume that the Gaussian beam is normally incident (*θ*
^{i}=0) and the ratio *σ*/λ=5 (see Eq. (35)). In Figs. 3–5, the extent of the computational region in the x-dimension is *L*
_{x}=16.3λ and this length is divided into *N*=299 segments of length Δ*x*=*L*
_{x}/*N*=0.055λ.

In the theory described in Section 2, some of the PML absorption parameters have a z-subscript that identifies them as pertaining to absorption for propagation in the z-direction. However, in the problem described here, absorption for propagation in the z-direction is not applicable and this means that we set ${\mathrm{\sigma}}_{\mathrm{z}}^{*}$=0 everywhere. However, in order to prevent reflection and simulate infinite extent in the x-direction, absorption of fields having propagation components in the x-direction is applicable. As shown in Fig. 1, there are two PML absorption regions when |*x*| in the vicinity of *L*
_{x}/2eac h region has a total thickness denoted by *δ*
_{PML}. If *η* is the magnitude of the penetration depth into the PML region, where *η* ranges over 0 → *δ*
_{PML}, we assume that the PML absorption increases quadratically with penetration depth as

which will be used in Eq. (9). The PML regions are divided into *N*
_{PML} layers each with thickness Δ*x* and the total thickness may then be written as *δ*
_{PML}=*N*
_{PML}Δ*x*. In all cases except the data shown in Fig. 9, we have set *N*
_{PML}=24 and the value of *A*
_{PML}=8. Since the x coordinate is discretized, we calculate the *N*
_{PML }values of Eq. (42) with *η* set to the middle value of the appropriate layer in a PML region.

## 3.2 Transmission Loss

The normalized transmitted power is calculated as the ratio of *P*
_{t}/*P*
_{i} where

Shown in Figs. 3–5 is transmitted power lost as *L*
_{r} increases. These three plots consider both TE and TM polarization and each plot includes four roughness realizations. The difference between the three plots is the rms roughness used to generate the realizations where the rms roughness values are 0.05λ, 0.10λ, and 0.20λ, respectively. The transmitted power is only plotted for every 10*δ*=1λ increase in *L*
_{r}. Thus, each curve in Figs. 3–5 contain 160 data points. However, since the numerical algorithm actually calculates in increments of *δ*=0.1λ, the total number of sublayers over the roughness region ranges from 0 to 1600 over the range of *L*
_{r}. To include the smooth region *L*
_{r}≤*z*≤*L*
_{t}, we use four sublayers with *δ*=750λ for a total thickness *L*
_{t}-*L*
_{r}=3000λ. As would be expected, it is seen in Figs. 3–5 that the rate of scattered power lost as *L*
_{r} increases becomes greater as the rms roughness increases. These curves are summarized later in Section 4.

## 3.3 Convergence

We now look at the stability of the numerical solution as certain parameters are varied. In Fig. 6, the number of digitization points parameter *N* is varied with all other parameters held constant. Figure 6 is similar to Fig. 3 through Fig. 5 where the residual transmitted power versus *L*
_{r} is shown, but only for realization 1. The main purpose of Fig. 6 is to show convergence as the number of digitization points *N* increases for constant *L*
_{x}=16.3λ. The data in Figs. 3–5 were generated for *N*=299 whereas in Fig. 6, the seven TM curves cover *N*=199 → 449 and the five TE cover *N*=199 → 399. These data indicate numerical convergence as *N* increases, or equivalently, as Δ*x* decreases.

In Fig. 7 convergence is now considered for the case where the increment Δ*x* is held fixed. We choose *L*
_{x} and *N* such that *L*
_{x}/*N*=Δ*x*≈0.055λ. Values of *L*
_{r} range from 14.3λ to 24.3λ and it is seen that as *L*
_{r} increases, the transmission at *z*=-10000λ is noticeably greater for the larger values of *N* and *L*
_{x}. The short answer for this discrepancy is that for larger *L*
_{x}, the plane wave solutions are absorbed by the PML layers at a lower rate as *z* → -∞. As an example, when *L*
_{x}=16.3λ one of the 299 calculated complex eigenvalues is *ξ*=(*ω*/*c*)(0.9340×10^{-4}-1.5i). When *L*
_{x}=2 4.3λ there are 441 calculated eigenvalues and the corresponding eigenvalue is now *ξ*=(*ω*/*c*)(0.1217×10^{-4},-1.5i). In the transmission region, the solutions are proportional to exp(*ξz*) (see Eqs. (32) and (33)) and it is clear that these particular eigenvalues pertain to the downward propagating plane wave solution in the *∊*=2.25 bounding medium of the waveguide structure. The absorptive real part of ξ is generated by the presence of the PML layers. Note that at *z*=-10000λ and the two *ξ* values above, the products ℜ(*ξ*)*z* are -0.9340(*ω*/*c*) and -0.1217(*ω*/*c*). For these values it is clear that exp(*ξz*) is still contributing a non-negligible solution value to the transmitted field. In hindsight, this discrepancy could have been avoided simply by choosing *z*<<-10000λ such that ℜ(*ξ*)*z*<<0 and all plane wave solutions would have been sufficiently decayed to be a negligible contribution to the transmitted field. We conclude that the discrepancy noted in Fig. 7 is not a convergence problem.

Shown in Fig. 8 are convergence characteristics as the number of PML layers *N*
_{PML} and the absorption coefficient *A*
_{PML} are varied. For curves labeled *N*
_{PML}=1 through 24, we set *A*
_{PML}=8 and vary *N*
_{PML} as indicated. For the curve labeled *N*
_{PML}=0, we have for numerical reasons actually set *N*
_{PML}=1 with *A*
_{PML}=0.1. This latter case does not completely eliminate PML absorption, but the effect is relatively minimal. We see that the *N*
_{PML}=0 curve shows little absorption at *z*=-10000λ and examination of the plane wave eigenvalues as was done in Fig. 7 shows a very low rate of absorption coefficient. As the PML absorption is increased, the curves quickly converge to a common result.

The dataset shown in Fig. 9 is the TM E-field intensity profile across the width *L*
_{x} of the computational domain. Unlike the previous figures, there is no initial smooth region followed by a roughness region. Here, all z-values are measured relative to the input aperture to the waveguide. It is seen that the PML layers perform quite well. The onset of PML layers is shown by the vertical lines near the x-limits of the computational region. In this example the total thickness of a PML region is 1.47λ. The intensity curve labeled *z*=-205λ will eventually become a single mode guided wave (as is the case in Figs. 3–5) after the plane wave and evanescent modes have vanished. This guided wave mode has eigenvalue *ξ*=(*ω*/*c*)(0,-1.552i).

In Fig. 10, the TM curves in Figs. 3–5 are averaged and modeled by fitting a curve of the form *y*(*x*)=*m*
_{0}10^{m1x} where *y* is the normalized transmission and *x* is *L*
_{r}/λ. The *m*
_{0} and *m*
_{1} are fitting parameters. In Figs. 3–5, the fitting parameters are *m*
_{0}=0.293, 0.286, and 0.277, and *m*
_{1}=-0.00045, -0.00168, and -0.00542, respectively. These analytical curves are plotted in Fig. 10 and the loss can be expressed in dB as -10Log_{10}[*y*(*x*)/*y*(0)] which yields 0.0045*L*
_{r}/λ, 0.0168*L*
_{r}/λ, and 0.0542*L*
_{r}/λ, respectively.

## 4 Concluding Remarks

We have presented a frequency domain calculation method and numerical results that predict scattering losses in planar waveguides due to wall roughness. The wall roughness is numerically generated assuming a Gaussian autocorrelation function with root mean square values of 0.05λ, 0.10λ, and 0.20λ and a correlation length of 1.0λ. We have also considered several aspects regarding convergence of the analytical/numerical method where various numerical parameters are varied.

The numerical trend of the transmission energy loss results obtained here are intuitive and entirely predictable: the rougher the random surface and the longer its length, the greater the guided wave propagation losses. However, the actual numerical amount of loss per unit length of roughness is not predictable and this can be an important question that depends on many parameters. For example, suppose the intrinsic material absorption is negligible and the measured scattering losses are unacceptably large. It may not be obvious as to what the dominant source or sources of scattering losses are. Should the focus be on reducing the wall roughness or decreasing the inhomogeneity of the waveguide core? These kinds of questions can apply to many electro-optic devices and theoretical/analytical tools can, among other things, help interpret experimental results and aid in development/research decisions. Regarding scattering due to waveguide core inhomogeneity, the work presented here can, with some minor modifications, be used to predict losses resulting from this source.

The objective of this paper is part methodological as well. The method is sufficiently general that many other interesting problems can be examined, such as photonic crystal structure, for example. The main criteria is that the structure of interest be described by spatially varying linear media. With the use of PML absorbing layers, the computational domain is finite and not dependent on periodic boundary conditions. Reflection from the “hard limits” at the truncation boundary of the computational domain is prevented by the PML algorithm which is commonly used in finite difference time domain algorithms. Of course, the method is easily adaptable to be used for infinitely periodic structure as well. The **R**-matrix algorithm allows very large computational domains to be considered with numerical stability.

## References and links

**1. **J. M. Elson and P. Tran, “R-matrix propagator with perfectly matched layers for the study of integrated optical components,” J. Opt. Soc. Am. A **16**2983–2989 (1999). [CrossRef]

**2. **J. M. Elson and P. Tran, “Coupled-mode calculation with the R-matrix propagator for the dispersion of surface waves on a truncated photonic crystal,” Phys. Rev. B **54**1711–1715 (1996). [CrossRef]

**3. **J. M. Elson and P. Tran, “Band structure and transmission of photonic media: a real-space finite-difference calculation with the R-matrix propagator,” *NATO ASI Series E: Applied Sciences on Photonic Band Gap Materials*Vol. 315 341–354 Crete, Greece June 15–29 (1995).

**4. **L. Li, “Multilayer modal method for diffraction gratings of arbitrary profile, depth, and permittivity,” J. Opt. Soc. Am. A **11**2816–2828 (1993). [CrossRef]

**5. **L. Li, “Formulation and comparison of two recursive matrix algorithms for modeling layered diffraction gratings,” J. Opt. Soc. Am. **13**1024–1035 (1996). [CrossRef]

**6. **J. P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comp. Phys. **114**185–200 (1995). [CrossRef]

**7. **J. P. Berenger, “Three-dimensional perfectly matched layer for the absorption of electromagnetic waves,” J. Comp. Phys. **127**363–379 (1996). [CrossRef]

**8. **D. Marcuse, “Mode conversion caused by surface imperfections in a dielectric slab waveguide,” Bell Sys. Tech. J. **48**, 3187–3215 (1969).

**9. **D. Marcuse, “Radiation losses of dielectric waveguides in terms of the power spectrum of the wall distortion function,” Bell Sys. Tech. J. **48**, 3233–3242 (1969).

**10. **F. P. Payne and J. P. R. Lacey, “A theoretical analysis of scattering loss from planar optical waveguides,” Opt. and Quantum Elec. **26**977–986 (1994). [CrossRef]

**11. **J. P. R. Lacey and F. P. Payne, “Radiation loss from planar waveguides with random wall imperfections,” IEE Proc. J. **137**282–288 (1990).

**12. **K. K. Lee, D. R. Lim, H. Luan, A. Agarwal, J. Foresi, and L. Kimerling, “Effect of size and roughness on light transmission on a Si/SiO_{2} waveguide: Experiment and model,” Appl. Phys. Lett. **77**1617–1619 (2000). [CrossRef]

**13. **F. Ladouceur, J. D. Love, and T. J. Senden, “Effect of side wall roughness in buried channel waveguides,” IEE Proc.-Optoelectron. **141**242–248 (1994). [CrossRef]

**14. **F. Ladouceur, J. D. Love, and T. J. Senden, “Measurement of surface roughness in buried channel waveguides,” Electron. Lett. **28**1321–1322 (1992). [CrossRef]

**15. **J. Rodríguez, R. D. Crespo, S. Fernández, J. Pandavenes, J. Olivares, S. Carrasco, I. Ibáñez, and J. Virgós, “Radiation losses on discontinuities in integrated optical waveguides,” Opt. Engr. **38**1896–1906 (1999). [CrossRef]

**16. **P. Lalanne and G. Morris, “Highly improved convergence of the coupled-wave method for TM polarization,” J. Opt. Soc. Am. **13**779–784 (1996). [CrossRef]

**17. **L. Li, “Use of Fourier series in the analysis of discontinuous periodic structures,” J. Opt. Soc. Am. **13**1870–1876 (1996). [CrossRef]

**18. **A. A. Maradudin, T. Michel, A. R. McGurn, and E. Mendez, “Enhanced backscattering of light from a random grating,” Annals of Physics **203**225–307 (1990). [CrossRef]