1. Introduction

Silicon-on-insulator (SOI) waveguides and optical on-chip circuitry rely on the high refractive index contrast between core (silicon, n = 3.5) and cladding (silicon dioxide, n = 1.45) materials to guide light in very compact structures and with small bending radii [1, 2]. SOI photonics is one of the most active areas of ongoing research and large-scale integrated circuits are being designed and fabricated. In many of these proposed circuits, one of the most critical waveguide components is the directional coupler between two parallel waveguides, which is used in microring-based filters [3], Mach-Zehnder interferometers and modulators [4, 5, 6], arrayed waveguide structures [7, 8, 9, 10] etc.

Here, we investigate and quantify the limitations of coupled-mode theory (CMT) in designing high index contrast SOI couplers and arrayed waveguides. CMT is a simple and often reliable approach to the design of such structures, in which the coupling coefficient can be written down in terms of overlap integrals of the individual waveguide modes and the refractive index distribution n(x,y) in the cross-sectional plane [11, 12, 13]. More accurate corrections to CMT for slab waveguides have been investigated by Chiang [14], and Payne [15], among others (see references therein).

In this paper, we compare “exact” results of the coupling coefficients of directional couplers calculated using a fully-vectorial finite difference frequency-domain (FDFD) mode-solving computer program with predictions of CMT. We demonstrate how to solve the inverse problem of reconstructing the coupling matrix from the solutions of the FDFD program. The significance of these results lies, firstly, in the resulting simplification of the design process of high-index contrast SOI waveguide directional couplers and coupled-waveguide structures, providing graphs and rule-of-thumb estimates for the validity and invalidity of simple CMT instead of having to rely on time-consuming ab initio numerical simulations in every case. Secondly, our method of solution of the inverse problem can be easily applied to design couplers from a mode solver rather than a time-consuming propagation simulation.

As a test structure which highlights both the applicability and shortcomings of CMT, we will consider the multi-waveguide coupled-array structure [16] which consists of a number of directional couplers parallel to one another, as is often used in arrayed waveguide gratings and multi-element lasers and amplifiers, and which require an accurate estimate of the coupling coefficients to prevent imaging and phase errors [17]. This structure, also called the “multi-slot waveguide,” has been recently fabricated and demonstrated by us in SOI [18]. We show that, because of the narrow waveguide widths allowed by the high index contrast, multi-waveguide structures can reveal significant next-to-nearest-neighbor coupling and other deviations from the conventional picture of modal coupling, and we find the “critical” waveguide-to-waveguide separation distance at which such terms become significant.

2. Coupled-mode theory (CMT) of the modes of multi-slot waveguides

If many slot waveguides are arranged in a parallel array, as would be encountered in the cross-section of an arrayed waveguide grating or coupler, or coupled-waveguide laser, then their modes can often be adequately described by supermode theory [16], which is one of the fundamental predictions of coupled-mode theory (CMT), and hence can be a test of the applicability of CMT to SOI photonics. The first part of this section will briefly describe an ab initio numerical algorithm we have encoded in MATLAB to calculate the modes accurately. The next part of this section presents the analysis of the modes based on supermode theory. In the following sections, we will compare the predictions of CMT with the FDFD calculations.

2.1. Finite-difference frequency-domain (FDFD) algorithm

In the finite-difference frequency-domain (FDFD) algorithm, the dielectric profile of the waveguide’s cross section is discretized on a rectangular grid. As developed by C.L. Xu et al. [20], the vectorial wave equation,

where E _⊥(x,y) is the transverse electric field vector, is written in matrix form as

and,

In this formulation the displacement vectors n ² E_x and n ² E_y are both continuous across any dielectric discontinuity, and with the graded-index approximation, the central difference equations can be applied directly without any special treatment at the boundaries. Fully discretized versions of these operators can be found in Ref. [21]. We have written a set of routines in MAT-LAB to calculate the modes and their effective indices for arbitrarily-shaped index profiles.

 

Fig. 1. (a) Refractive index profile in the transverse plane n(x,y) for an N arrayed-waveguide structure. (b) The same refractive index profile may be decomposed, mathematically, into the sum of parts, Δn² _i = n² _i -n² _s , each of which appears in integrals equation for the coupling coefficients. For the structures considered in this paper, h = 500 nm, w = 200 nm, and s varies over the range 50 nm to 1μm. For these waveguide widths and heights (similar to those in Ref. [18]), the polarization direction of the principal transverse component of the electric field is indicated for the (quasi) TE and TM modes.

Download Full Size | PDF

2.2. Coupled mode theory and its predictions

To describe the modes of multislot waveguides, we begin with the wave equation [19],

and consider each polarization in turn. (TE and TM polarization are defined in terms of as the major component of the electric field, which for the structure in question, are polarized vertically and horizontally as shown in Fig. 1.)

2.2.1. TE Polarization

Consider an array of N single mode waveguides, whose refractive index profile is shown schematically in Fig. 1a. As long as the waveguides are not too close to each other (i.e., greater than a separation distance which we will investigate and quantify in a subsequent section), the transverse mode profile of the multislot waveguide structure can be approximated by an expansion of the individual high index waveguide modes

As shown by Fig. 1, the relative dielectric coefficient distribution of the entire N waveguide structure n ² (x,y) can be written as a sum of individual waveguide contributions, so that

where n ² _s (x,y) corresponds to the cladding. Thus, n ² _s (x,y) + Δn ² _l (x,y) would yield the dielectric coefficient profile of the l th waveguide in the absence of the others. Substituting the above two equations into the wave equation, we have,

The modes of the individual waveguides satisfy their respective eigenvalue equations,

and therefore, using Eq. (7), Eq. (6) can be written as,

where,

N equations are formed by multiplying Eq. (8) by 𝓔_j ^* (j = 1,2,…,N), and integrating each of these equations over x and y,

We define the modal overlap integrals as follows:

with the normalization

I_jl is the overlap integral of the modes of two waveguides which are not orthogonal to each other (particularly in the case of small waveguide separation), and κ are the self-coupling and cross-coupling (exchange coupling) coefficients familiar from coupled-mode theory [19, p. 362].

Equation (9) can then be written in matrix form as an eigenvalue problem,

(The matrix on the left-hand side of the above equation will be referred to as M.)

 

Fig. 2. TE Polarization E_y: The modes of an N = 5 coupled waveguide array for λ = 1550 nm, calculated using coupled-mode theory (blue solid lines), and a finite-difference frequency-domain algorithm (black crosses). The coupled-mode theory calculations were done by using the effective index method, calculating the overlap integrals, solving Eq. (11), and reassembling the field. Waveguide height = 500 nm, width = 200 nm, separation = 200 nm, n_core = 3.47, and n_clad =1.46. Under nearest neighbor coupling, the scaling relationship predicted by Eq. (13) adequately predicts the field amplitudes within each waveguide.

Download Full Size | PDF

If we assume that only nearest neighbor coupling is significant, then the integrals in Eq. (10) are nonzero only when l = j -1, j, j +1. M takes the tridiagonal form,

If the waveguides are identical and equally spaced, M can be further simplified by setting β ² ₁ = β ² ₂ … = β ^N ₁ ≡ β ² ₀ and also, I _l,l+1 = I _l-1,l ≡ I,κ _l,l+1 = κ _l-1,l ≡ κ. However, even if the waveguides are identical and equally spaced, κ ₁₁ and κ_NN are not equal to κ ₂₂, κ ₃₃,…,κ _N-1N-1. In fact, Eq. (10b) shows that for those waveguides at the edges (l = 1 and I = N) there are approximately only half as many contributing terms as the other waveguides: there are no waveguides to the left of the l = 1 waveguide, and there are no waveguides to the right of the l = N waveguide, whereas all the other waveguides have contribution terms from both the left and right halves of their modal profiles.

We define κ _self ≡ κ ₂₂, κ ₃₃,…, κ _N-1N-1, κ _self,edge ≡ κ ₁₁ and κ_NN, and δκ _self = κ _self - κ _self,edge.

 

Fig. 3. TM Polarization E_x: The modes of an N = 5 coupled waveguide array, calculated using coupled-mode theory (blue solid lines), and a finite-difference frequency-domain algorithm (black crosses). The coupled-mode theory calculations were done by using the effective index method, calculating the overlap integrals, solving Eq. (11), and reassembling the field. Waveguide height = 500 nm, width = 200 nm, separation = 1μm, n_core = 3.47, and n_clad =1.46. Under nearest neighbor coupling, the scaling relationship predicted by Eq. (13) adequately predicts the field amplitudes within each waveguide.

Download Full Size | PDF

To first order in the perturbation δκ _self, the eigenvectors are

where m is the modal number and l indicates which high-index rib waveguide (or low-index slot) is being described. [The expression for the eigenvalues is written later, Eq. (22).]

For large N, the second term in the above expression, is smaller than the first by (N + 1)^-1 and can be ignored, yielding a simpler expression. The progression of peak-amplitude values (in the high index regions) {A ^(m) _l},l = 1,…,N matches with the numerical calculations shown in Fig. 2. However, we shall see that the agreement is good only at large separation distances between the individual waveguides.

2.2.2. TM Polarization

For the TM polarization (in which the electric field is normal to the waveguide/slot boundary), we again start with the wave equation only now defined in terms of the magnetic field, which we expand in terms of the individual waveguide modes,

Under nearest neighbor coupling, the magnetic field will also obey the scaling relationship of Eq. (13). If it can be assumed that ∣∂H_z/∂_y∣ ≪ ∣∂H_y/∂_z∣, then H_y and E_x are related by

so that, within the high index ribs, we expect that Eq. (13) describes the scaling relationship of the peak electric field amplitudes. Although the peak electric field amplitudes of the entire modal profile are found not in the high-index regions, but in the low-index regions (just inside the core-cladding boundary), they obviously satisfy the same scaling law, as can be seen in Fig. 3.

3. Numerically-assisted CMT: The “Inverse Problem”

CMT offers valuable physical insight into how waveguides couple—in particular, the structure of matrix M in Eq. (11) is revealing—but the quantitative predictions of CMT are in error in high-index-contrast SOI structures at short separation distances. To obtain a numerically-accurate picture of modal coupling, we propose a new extension of CMT, which we call “numerically assisted” CMT, to use the simulation results of the FDFD algorithm to back-calculate the elements of the coupling matrix M. We can thereby check if the assumption of nearest-neighbor coupling is valid at short separation distances, and identify various other interesting coupling phenomena (e.g., non-Hermiticity of M) which have not been pointed out earlier.

To develop NA-CMT, we use the following mathematical procedure, based on a matrix theorem previously developed for coupled-resonator structures [23].

4. Asymptotic accuracy of numerically-assisted CMT

In this paper, the eigenvalues (λ_k) and eigenvectors (A ^(k)) are obtained from a computer simulation. However, they may be obtained from measurements on fabricated structures, in order to test whether the intended coupling matrix was successfully obtained in practice. The experimental procedure to measure eigenvalues and eigenvectors could be similar to that used to image the modes of laser resonators.

We assume that the measurements result in some small, uncorrelated errors in the eigenvalues (Δλ_k) and eigenvectors (Δu_k). The inversion algorithm presented in the previous section can also be used with measured data. In this section we study the accuracy of the nearest-neighbor coupling and next-to-nearest-neighbor coupling coefficients in terms of Δλ_k and Δu_k.

First, we will carry out a simple theoretical estimation. We will assume that the coupling matrix is Hermitian. After some algebra, the error in any element of M can be written to first order as

For simplicity, in this paper, we assume that Δu_k is zero, i.e., the errors are only in the measured eigenvalues, since for identical arrayed waveguide structures, successive eigenvectors look quite different from each other and are easily distinguished [23]. We will also use the simpler form of the eigenvectors, retaining only the first term of Eq. (13), so that the errors in the reconstructed nearest-neighbor coupling and next-to-nearest-neighbor coupling coefficients are

Assuming that Δλ_k are uncorrelated identically-distributed random variables with mean E[Δλ] and variance Var[Δλ], the mean and variance of the nearest-neighbor coupling and next-to-nearest-neighbor coupling coefficients can be calculated. Both ΔM _jj+1 and ΔM _j-1j+1 are zero-mean, since, for example,

because the summation vanishes as a consequence of the orthogonality of the eigenvectors (the sum is equal to u ^(j).u ^(j+1) = 0).

To calculate the variance, we see that

and the same result is obtained for Var[ΔM _j-1j+1].

Summarizing the results,

Numerical calculations, shown in Fig. 4 confirm Eq. (;20). (Numerical calculations show that the same relations are seen to hold in the case of non-identical waveguides, in which case the off-diagonal terms of the coupling matrix are not identical along the sub-diagonals, and also for slightly asymmetric matrices.)

 

Fig. 4. Error versus N: Exact eigenvalues of a tridiagonal symmetric matrix of size N were perturbed by values chosen from a uniform random distribution with variance chosen to be ten percent of the first eigenvalue. The variance and mean of the reconstructed nearest-neighbor coupling and next-to-nearest-neighbor coupling coefficients are plotted, calculated from a distribution of coupling matrices generated by 10⁵ iterations, showing that Eq. (20) is a good predictor of the reconstruction accuracy.

Download Full Size | PDF

These results show that the error in reconstructing the coupling coefficients decreases, rather than increases, as the number of inaccurately-measured eigenvalues increases. This results from (spectral) averaging: each reconstructed coupling coefficient averages over the entire spectrum of eigenvalues, and therefore, benefits from the law of averages. In contrast, directly measuring a coupling coefficient e.g., by a local near-field probe of the field in the coupling region, does not benefit from any ensemble averaging.

5. Discussion

5.1. Next-to-nearest-neighbor coupling

As Table 1 shows (calculated at one specific value of the waveguide separation distance), M contains useful information about non nearest-neighbor coupling. We can read off whichever coupling coefficients are needed: in particular, we calculate the ratio ∣κ ₁₃/κ ₁₂∣, i.e., the ratio of next-to-nearest-neighbor coupling coefficient to the nearest-neighbor coupling coefficient.

Table 1. An example of a reconstructed coupling (M) matrix from FDFD calculations of eigenmodes and eigenvalues. Si/SiO₂, TE polarization, separation s = 350 nm, β ₀ = 2.26128 (2π/λ). Although the nearest-neighbor coupling coefficients dominate, the self-coupling and off-tridiagonal coupling terms are non-zero.

View Table

First we estimate the expected dependency of this ratio of coupling coefficients to the edge-to-edge separation, s. Using Kuznetsov’s solution for the coupling coefficients of two slab waveguides [22], we observe that κ in both the TE and TM cases varies with s as κ ~ e ^-ps where p is the field decay length in the cladding. Therefore, the ratio κ ₁₃/κ ₁₂ for both polarizations has the following expression (to leading order),

i.e., the ratio of next-to-nearest-neighbor coupling coefficient to the nearest-neighbor coupling coefficient should fall off exponentially with increasing separation.

Figure 5 shows the calculations of this ratio using the above algorithm. The exponential fit describes the TE polarization much better than it does the TM polarization, indicating that some of the central assumptions of CMT are starting to fail for the TM polarization at short distances. The next section will describe another symptom of the failure of CMT, obtained by looking at the eigenvalues, i.e., the propagation constants, of the supermodes.

5.2. Eigenvalue fanout: effective index of the supermodes versus separation distance

Another way to evaluate the predictions of CMT is the theory behind the eigenvalues of Eq. (11), which predicts that the effective index of the m-th supermode is given by the equation

where β ₀ is the propagation constant of a single waveguide in isolation. Note that for N = 5, the m = 3 supermode has the special property that the right-hand-side of the above expression , i.e., the index of that supermode does not change with the coupling coefficient κ. Hence, n ⁽³⁾ _eff is only weakly dependent on the separation distance (through the self-coupling coefficients, κ ₁₁,κ ₂₂,…,κ ₅₅).

To verify this prediction, Fig. 6 shows the effective index calculated by FDFD for each of the five supermodes at various separation distances in three different silicon-based material systems. These values of the effective index take into account coupling-induced frequency shifts (CIFS, [25]) because M itself results from a numerical calculation of the supermodes (and their eigenfrequencies), rather than individual waveguide modes and the propagation constants of isolated waveguides.

In the limit of large separation, the effective indices of all the supermodes tends to that of the single waveguide. As the separation distance is decreased, the coupling coefficients increase, and the effective indexes of the different modes separate [24]. The first three modes remain guided even as s shrinks to zero, since their effective indices are higher than that of the single waveguide. From Fig. 6, one can read off the waveguide separation distance at which conventional CMT is expected to fail, and more accurate design tools, such as FDFD calculations, should be used to accurately predict the coupling coefficients.

 

Fig. 5. Ratio of coupling coefficients for different separation distances extracted from Eq. (11), which was reconstructed using an algorithm described in the text. (a) TE Polarization An exponential fit expected from a simple nearest-neighbor-coupling theory holds throughout this regime. (b) TM Polarization At a separations less than 450 nm, the ratio deviates significantly from the predicted behavior. (c) TE Polarization The ratio of cross coupling coefficients show that the reconstructed coupling matrix M becomes asymmetric as the waveguide separation is reduced. (d) TM Polarization The asymmetry of the coupling matrix begins at a larger separation.

Download Full Size | PDF

An interesting observation obtains from the m = 3 supermode: at a certain (small) waveguide separation, n ⁽³⁾ _eff is no longer independent of s and begins to deviate substantially from a straight line, contrary to the prediction of Eq. (22). This deviation is much more pronounced in the case of the TM polarization.

 

Fig. 6. Left column: TE Polarization, and Right column: TM polarization. Effective index of the five supermodes for different separation distances with n_core = 3.47, and (a,b) n_clad =1.46, (c,d) n_clad = 1 (e,f) n_clad = 2.05. For each case as the separation between the waveguides increases, the effective indexes of the modes converge to that of the single waveguide. These values are (a) n_eff = 2.36 and (b) n_eff = 1.66 for oxide cladding, (c) n_eff = 2.24 and (d) n_eff = 1.07 for air cladding, (e) n_eff = 2.56 and (f) n_eff = 2.26 for nitride cladding. The shaded regions indicate > 5% deviation of n_eff for the m=3 supermode from its theoretical value, which as discussed in the text, is predicted by CMT to be independent of the separation distance.

Download Full Size | PDF

5.3. Field skewing and reshaping

At short separation distances, the reconstructed matrix M can become non-symmetric (non-Hermitian), although the eigenvalues remain stricly real as long as the mode is above cut-off. This can be seen in fact in the matrix written in Table 1 and Fig. 5(c,d): κ ₁₂ ≠ κ ₂₁ and κ ₁₃ ≠ κ ₃₁, etc.

 

Fig. 7. TM Polarization E_x: The field profile of the fifth eigenmode in the first waveguide. When the separation is decreased below 450 nm, the peak of the field in the high-index rib indicated by the dotted red line in (a) is no longer centered, and the mode shape is considerably altered, thereby changing both κ and n _eff. Consequently, CMT can no longer accurately predict the mode coupling.

Download Full Size | PDF

The reason for this asymmetry is that the fields within the individual waveguides are no longer centered between the dielectric boundaries. As shown in Fig. 7. the modal profile starts to deviate in the location of its maxima and minima. For example, the peaks of the field in the outermost ribs are skewed and no longer centered in the middle of the dielectric boundaries, and can even reach the boundaries of the high-index and low-index regions. It is no longer accurate to read off the peak amplitudes of the supermode in order to write the eigenvectors A ^(m) in Eq. (13)—doing so would result in asymmetric M matrices.

At a short separation distance of 80 nm, Fig. 8 shows the coupled-mode theory used to reconstruct the m = 1 and m = 5 supermode (plotted with continuous lines) and the supermode calculation of FDFD (with crosses). Note that in both cases the field is asymmetrically centered within the dielectric boundaries of the outer waveguides. Recall that CMT is based on writing the field as a summation of the scaled individual waveguide modes, Fig. 8(a,d), each of which is centered within its own core-cladding boundaries. At short separation distances, when, for example, there is a significant contribution of the (asymmetric) tail from the field in the second waveguide to the (symmetric) mode of the first waveguide, CMT itself predicts a lateral shift of the peak (of the sum) away from the exact center of the waveguide. The scaling relationships from Eq. (13) will enhance this effect for an multi-waveguide arrayed structure compared to a (twin-waveguide) directional coupler.

 

Fig. 8. TE Polarization E_y: Using the exact solution from a FDFD simulation of a single waveguide, the horizontal cross section is extracted and five copies are shifted from one another so that their separation corresponds to a waveguide separation of 80 nm. (a) These individual waveguide modes are scaled in accordance with Eq. (13) for the fundamental mode (m=1). (b) The sumation of the individual waveguide modes; superimposed is the FDFD solution of the entire five waveguide structure. (c) Zoomed in to just the first waveguide. CMT and FDFD show a shift of the mode towards the center of the waveguide structure. (d-e) The fifth mode, both CMT and FDFD show a shift towards the edge of the waveguide structure however FDFD shows a shift of greater magnitude.

Download Full Size | PDF

For the fundamental mode, Fig. 8(a-c), the summation of the fields associated with the first (blue) and second (green) waveguides results in the peak shifting towards the center of the five waveguide structure, which qualitatively agrees with the FDFD simulation. But the FDFD result for the fifth mode shows a shift of greater magnitude, now towards the outer edge of the waveguide structure, indicating that CMT no longer accurately predicts the modal profile of the supermode. TM polarized modes start to shift at a larger separation, due to the field discontinuities at the waveguide boundaries and electric field enhancement in the cladding regions.

Note also, as shown in Fig. 8(e), that FDFD predicts a different exponential decay constant of the field wings, compared to CMT. This is a fundamental failure of CMT in the sense that the eigenmode of the composite structure can no longer be written as the sum of modes of individual waveguides (in isolation from each other). Until a satisfactory theory of mode evolution in this strong-coupling regime can be obtained, we recommend that designers rely on direct numerical simulations to obtain field profiles.

5.4. Polarization hybridization

Alongside this phenomenon, we also observe that the polarization becomes strongly hybridized as the separation distance is reduced. Figure 9 shows the major and minor E field components at large and small separation distances, along with the cross section of the Poynting vector, which indicates power flow. Notice that at small separation distances, the polarization component that was previously negligible has become, in fact, the dominant one. Furthermore, the power is actually carried above and below the waveguide structure at the outer edges, rather than within the inner slots, contrary to the original intention of slot waveguides.

These calculations suggest that at small separation distances, CMT of high-index contrast waveguides should consider both polarization components of the electric field. Mathematically, the dimensionality of the basis set (of eigenvectors) depends on the separation distance. Practically, one should pay careful attention to where the power flow is concentrated in waveguide couplers at short separation distances, since both polarization components of the electric field (with different coupling lengths due to different modal effective indices) play a role in the exchange of power between the constituent waveguides.

 

Fig. 9. TM Polarization: finite element simulations (COMSOL) of the transverse field amplitudes, E_x and E_y, and Power flow, P_z for a separation of 1 μm (first Row) and 100 nm (second row) for the fifth mode. For small separations the polarization becomes strongly hybridized and the power is carried at the outer edges.

Download Full Size | PDF

6. Conclusion

In this paper, we have studied the validity of coupled-mode theory CMT for high-index contrast (e.g., silicon-based) optical waveguiding structures, in particular directional couplers and multi-slot waveguides. We have used a fully vectorial finite-difference frequency-domain (FDFD) algorithm to obtain modal profiles and effective indices of the supermodes in a non-perturbative way. When the modal profiles can be “discretized” to read off peak amplitudes within each of the waveguide cores, a theorem from matrix algebra can be employed to solve the “inverse problem”—to reconstruct accurately the matrix of coupling coefficients M (a procedure we have called numerically-assisted coupled-mode theory, NA-CMT). The NA-CMT framework can be used to find out when the nearest-neighbor-coupling approximation breaks down.

Aside from the inversion procedure, the results of FDFD calculations also directly address this problem. We have presented the eigenvalue fanout graphs for three different silicon-based material platforms and indicated the separation distances at which the CMT predictions break down. We have also shown that there are significant deviations in the modal shape as the waveguide separation distance decreases and the state of polarization of the field may change between quasi-TE or quasi-TM to a strongly-hybridized polarization.

We suggest that the terminology “strong coupling” in the context of waveguide directional couplers be used to describe this coupling regime, for which a simple theory is not yet available, rather than simply large values of the length-integrated coupling coefficient, which can approach unity even with small coupling coefficients per unit length for long couplers, as this latter regime has already been well studied.