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

 

Fig. 1. Schematic of waveguide with wall roughness. The waveguide channel, having nominal width ω and permittivity ∊ ₂, is bounded by media with permittivity ∊ ₁. In the z-direction, the waveguide channel consists of a smooth segment (L _r≤z≤L _t), a roughness segment (0≤z≤L _r), and a semi-infinite substrate region z≤0. The x-dimension of the computational region is bounded by x=±L _x/2 and the shaded regions at the extreme left and right denote the PML absorption layers.

Download Full Size | PDF

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=ξ ²=ξ·ξ. The columns of the N×N matrix S are the eigenvectors of M and the diagonal matrix ξ ² 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=δ ₁, we initiate the recursion with R _ij (δ ₁)=r _ij (δ ₁). If the next z-invariant layer has thickness δ ₂, new r _ij (δ ₂) matrices are computed and the R _ij (δ ₁) are used in Eqs. (26)–(29) to yield R _ij (δ ₁+δ ₂). 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 θ ⁱ. We write the x and y-components of the incident field as

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

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)=(qc/ω)${\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 ∊₁=(2.25, 0.) and ∊₂=(2.50, 0.). In most cases, we assume that the Gaussian beam is normally incident (θ ⁱ=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λ.

 

Fig. 2. Schematic showing modeling of waveguide roughness region. The waveguide channel has nominal width ω. For this example, the roughness region is divided into N _r=5 sublayers of thickness δ and the channel width within each sublayer is varied in accordance with the roughness. By doing this, each sublayer is z-invariant and solutions of the form given in Eqs. (34) and (35) apply throughout the waveguide. In the numerical results of this work, the number of sublayers ranged from 0 to 1600, δ=0.1λ, and 0λ≤L _r≤160λ. The permittivity of the light-shaded and dark-shaded regions are ∊₁=(2.25, 0.) and ∊₂=(2.50, 0.), respectively.

Download Full Size | PDF

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.

 

Fig. 3. Normalized transmitted power versus length of roughness calculated at z=-10000λ. There are four groups of curves where each group has four curves. In each group the four curves, labeled 1–4, correspond to a different roughness realization. The upper and lower two sets of curves are for TE and TM polarization, respectively. The upper two sets of data have been displaced upward 0.4 units for clarity. The horizontal curves are for the normalized reflected power. The downward sloping curves are the normalized transmitted power. The number of discretization points N=299 and the nominal waveguide channel width ω=1λ. All realizations are generated with rms roughness 0.05λ and correlation length 1λ.

Download Full Size | PDF

 

Fig. 4. Normalized transmitted power versus length of roughness calculated at z=-10000λ. These data are analogous to those in Fig. 3 except that all realizations are generated with rms roughness 0.10λ.

Download Full Size | PDF

 

Fig. 5. Normalized transmitted power versus length of roughness calculated at z=-10000λ. These data are analogous to those in Fig. 3 except that all realizations are generated with rms roughness 0.20λ.

Download Full Size | PDF

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.

 

Fig. 6. Convergence of normalized transmitted power for various values of discretization points N. For the same roughness realization 1 as used in Figs. 3–5, the transmitted power is calculated at z=-10000λ versus length of roughness. The upper set of curves, which have for clarity been displaced upward by 0.2, are for TM polarization and the values of N vary from 199 to 449. The lower set of curves are for TE polarization and the values of N vary from 199 to 399. For the corresponding polarization, the curves are essentially converged except for the N=199 curves. The rms roughness is 0.2λ.

Download Full Size | PDF

 

Fig. 7. Convergence of normalized TM transmitted power for various values of discretization points N and length L _x. The ratio of L _x/N=Δx is adjusted so that Δx≈0.055λ is nearly constant. For the same roughness realization 1 that was used in Figs. 3–5, the transmitted power is calculated at z=-10000λ versus length of roughness and the rms roughness is 0.2λ. Note that as L _x increases, there is an increase in transmitted power due to a decrease in the absorption rate for the plane wave solutions. This is not a convergence related issue.

Download Full Size | PDF

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.

 

Fig. 8. Normalized TM transmitted power for various values of PML absorption layers N _PML. For 1≤N _PML≤24, the value A _PML=8. For numerical reasons, the results shown by the curve labeled N _PML=0 was actually calculated with N _PML=1 and A _PML=0.1 (see Eq. (42)). Thus, the N _PML=0 curve only approximates the absence of PML layers. The incident beam is at normal incidence, L _x=16.3λ, N=299, and rms roughness is 0.2λ. These data are for the same roughness realization 1 that was used in Figs. 3–5. The group of horizontal lines are the normalized reflected power. It is clear that if there is sufficient PML absorption, then the solutions converge to a common result.

Download Full Size | PDF

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 ₀10^m1x where y is the normalized transmission and x is L _r/λ. The m ₀ and m ₁ are fitting parameters. In Figs. 3–5, the fitting parameters are m ₀=0.293, 0.286, and 0.277, and m ₁=-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₁₀[y(x)/y(0)] which yields 0.0045L _r/λ, 0.0168L _r/λ, and 0.0542L _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.

 

Fig. 9. TM E-field intensity profile at various depths from the waveguide opening. There is no roughness region in this example and the 8 indicated z values are relative to the waveguide input aperture plane. The width L _x=18.3λ and the angle of incidence is θ ⁱ=30°. For clarity, all curves, except the z=-205λ curve, have been displaced by multiples of 0.01 units. The z=0λ curve is the electric field intensity across the waveguide input aperture plane and the remaining z-value curves show the progression of the intensity curves at indicated depths into the waveguide. It is clear that the refracted part of the incident Gaussian profile is ultimately absorbed by the right-hand set of PML layers without reflection. In addition, the final z=-205λ curve is beginning to take shape as the residual guided wave since much of the transmitted energy has disappeared. The two vertical lines at x=±0.5λ indicate the boundaries of the waveguide channel. The two vertical lines at x=±7.65λ indicate the start of the PML absorbing layers.

Download Full Size | PDF

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.

 

Fig. 10. Analytical curves derived from the data shown in Figs. 3–5. The four curves in each of these figures are each compiled into one analytical curve as a power of 10. This yields an estimate of the dB loss as a function of the length of roughness.

Download Full Size | PDF