A variational mode solver for optical waveguides based on quasi-analytical vectorial slab mode expansion

A variational mode solver for optical waveguides based on quasi-analytical vectorial slab mode expansion

Alyona Ivanova MESA+ Institute for Nanotechnology, University of Twente, Enschede, The Netherlands Remco Stoffer Manfred Hammer MESA+ Institute for Nanotechnology, University of Twente, Enschede, The Netherlands

A flexible and efficient method for fully vectorial modal analysis of 3D dielectric optical waveguides with arbitrary 2D cross-sections is proposed. The technique is based on expansion of each modal component in some a priori defined functions defined on one coordinate axis times some unknown coefficient-functions, defined on the other axis. By applying a variational restriction procedure the unknown coefficient-functions are determined, resulting in an optimum approximation of the true vectorial mode profile. This technique can be related to both Effective Index and Mode Matching methods. A couple of examples illustrate the performance of the method.

Keywords: Integrated optics; Dielectric waveguides; Guided modes; Numerical modelling; Variational Methods
PACS: 42.82.–m; 42.82.Et

1 Introduction

Three-dimensional optical channel waveguides are basic components of integrated optical devices such as directional couplers, wavelength filters, phase shifters, and optical switches. The successful design of these devices requires an accurate estimation of the modal field profiles and propagation constants. Over already some decades several classes of methods for the analysis of dielectric optical waveguides were developed: among these are techniques of more numerical character, like Finite Element and Finite Difference approximations, the Method of Lines, and Integral Equations Methods, but also more analytical approaches like Film Mode Matching (FMM) and the Effective Index Method (EIM). Detailed overviews of these techniques can be found in [1, 2, 3].

In the present paper we propose an extension of the scalar mode solver [4] to vectorial problems. Our method is based on individual expansions of each mode profile component into a set of a priori defined functions of one coordinate axis (vertical), here, field components of some slab waveguide mode. The expansion is global, meaning that the same basis functions are used at any point on the horizontal axis. The unknown expansion coefficients – in our case functions, defined on the horizontal axis – are found by means of variational methods [5, 6].

The present method can be viewed as some bridge between two popular approaches, namely the FMM on the one hand and the EIM on the other. In the standard EIM the 2D problem of finding modes of the waveguide is reduced to consecutive solving two 1D problems: at first, the 1D modes, and their propagation constants, of the constituting slab waveguides are found, and then their propagation constants are used to define effective refractive indices of a reduced 1D problem. In general this is a very quick and easy approach for a rough estimation of mode parameters. However, in case one of the constituting slabs doesn’t support a guided mode (for example, some substrate material with air on top) it is impossible to uniquely define the effective refractive index in that particular region of the reduced problem. Should it be the refractive index of the cladding, refractive index of the air, or something in between? The restriction of the present approach to one-term expansions will answer this question.

The validity of the method was checked on several structures, including waveguides with rectangular and non-rectangular piecewise-constant refractive index distributions, and a diffused waveguide. Comparison shows that the present method is a more consistent and accurate alternative to the standard EIM and also can be pushed to its limits and used for rigorous computations.

The paper is organized as follows. In section 2 the problem of finding vectorial modes of the dielectric waveguides is stated, then some properties of slab modes and the modal field ansatz are described in sections 3 and 4. The equations for the coefficient functions are derived in section 5. Section 6 outlines the numerical solution methods. The relation of the present method to the EIM and FMM is explained in more detail in sections 7 and 8. Then in section 9 numerical results for several waveguide configurations are presented. Finally some concluding remarks are made in section 10.

2 Variational form of the vectorial mode problem

Consider a -invariant dielectric isotropic waveguide defined on its cross-section by a refractive index or relative dielectric permittivity distribution. Figure LABEL:Introductionfig shows two examples.


figure\end@float The propagation of monochromatic light, given by the electric and magnetic components of the optical field, with propagation constant and frequency ,


is governed by the Maxwell equations for the mode profile components and




vacuum permittivity , vacuum permeability , relative permittivity . Here and further in this paper it is assumed that the relative permeability is equal to , as is the case for most materials at optical frequencies.

We will work with a variational formulation of the Maxwell equations. Solutions of the equations (2) correspond to stationary points of the functional [5]


with propagation constant equal to the value of the functional at the stationary point. The inner product used is . The natural interface conditions are the continuity of all tangential field components across the interfaces.

3 Slab modes

In this section we will consider modes of slab waveguides, which we will use in the next section as building blocks to construct approximations of the modes of waveguides with arbitrary 2D cross-sections. Furthermore, we introduce rotations of the slab modes; these rotations will be needed to provide a physical motivation for the particular form of the approximations that we will employ.



A one dimensional TE mode, propagating in the -direction with propagation constant , of the slab waveguide, given by the permittivity distribution (Figure LABEL:SlabModesfig) can be represented as


The principal electric component satisfies the equation


with vacuum wavenumber . The remaining two nonzero components of the mode profile can be derived directly from :


The slab waveguide (Figure LABEL:SlabModesfig) is by definition invariant in the -plane. So if a modal solution of Maxwell equations propagating in the -direction will be rotated in the -plane by an angle (Figure LABEL:RotationAngleFieldTETMfig), it will still remain a modal solution of the Maxwell equations, but now propagating in the direction :



Similarly a 1D TM slab mode, propagating in the -direction with propagation constant


will still be a solution of the Maxwell equations after a rotation around the -axis (Figure LABEL:RotationAngleFieldTETMfig)


The principal magnetic component satisfies the equation


Again the remaining two nonzero components of the mode profile can be derived directly from :


4 Modal field ansatz

We now return to the vectorial modes of the 3D waveguides, as in section 2. Each field component is represented individually as a superposition of a priori known functions , defined on one coordinate axis, times some unknown coefficient-function , defined on the other axis:


For the functions we will take components of slab modes from some reference slice(s). Further in the paper two types of the expansion will be relevant, one which introduces 5 unknown functions per slab mode, and another one, which introduces only 3. These will be called five component approximation (VEIM5) and three component approximation (VEIM3), respectively.

In case of VEIM5, the TE basis mode (5) number  with mode profile components , , contributes to the expansion of components , , , and with the form , , , and . Likewise, the TM basis mode (9) number  with mode profile components , , contributes to the expansion of components , , , and with the form , , , , , such that the complete expansion looks like


This expansion has the drawback that the functions making up some of the components can become linearly dependent; for example, the full set of components from TM modes form a complete set; thus any from a TE mode can be expressed in that complete set of functions. When using a limited number of modes in the expansion, no problems result from this; however, increasing the number of modes will at some point make the problem ill-conditioned. Therefore, we introduce a different expansion, which we call VEIM3, in which we omit contributions of some modal components - making sure that each vector component is only represented by either TE of TM slab mode components. So a TE basis mode (5) number  with mode profile components , , contributes to the expansion of components , and with the form , , . Likewise a TM basis mode (9) number  with mode profile components , , contributes to the expansion of components , and with the form , , , such that the complete expansion looks like


Note that in both expansions each contributing component of a 1D mode is used to represent the field not only in the slab segment where it belongs, as in EIM and FMM methods, but also in the whole waveguide. So even with a single slab mode in both expansions, (14) and (15), it is possible to construct an approximation of the field in the whole structure. In section 7 we will study in detail properties of such one-mode-expansions.

The form of the expansion (14) was inspired by the mode matching techniques that use the physically motivated approach of employing rotated modes (8), (10) to locally expand the total field [8, 9]. In the present approach though, we attribute those parts of the slab mode components that do not depend on to the functions , treating them as unknowns – but the -dependence of the and components is still the same. In the sections 7 and 8 we will study the behavior of these functions .

What concerns the choice of the reference slice(s), it seems that modal components from the slice, where the maximum power is expected to be localized, give the best results. Further in this paper VEIM5 will be used with a few modes only for rough and efficient approximations, while VEIM3 will be used with higher numbers of modes to obtain accurate, converged results. We do not restrict to using modes from one reference slice only; we observe that adding mode(s) from another slice can greatly improve accuracy for lower number of modes in the expansion. However, care must be taken; the problem can become ill-conditioned if the modes become nearly linear dependent.

In the following all the slab mode components , which are used to expand a field component of the complete waveguide, we will denote as (just like in eqn. (13)).

5 Reduced problem

The next question is how to find corresponding functions , such that the expansion (13) represents the true solution in the best possible way. For this purpose we apply variational restriction [2, 6] of the functional (4). In short it can be outlined as follows. As it was already mentioned the critical points of the functional (4), which satisfy some continuity conditions, are solutions of the Maxwell equations (2) and, vice versa, solutions of the Maxwell equations (2) are critical points of the functional (4).

After insertion of the expansions (14) or (15), variation of the functional (4) with respect to a function , a vector function made up of all functions , results in the following system of first order differential equations for with parameter :


The elements of the matrices include the overlap integrals (here: ) of the functions , their derivatives, and the local permittivity distribution of the waveguide:


Note that the permittivity appears only in , and , hence only these matrices are -dependent.

If the permittivity exhibits discontinuities along the -direction, the functions


are required to be continuous at the respective positions.

It turns out that by algebraic operations the system of first order differential equations (16) can be reduced to a system of second order differential equations for the vector functions and only. Moreover, since the components and , and are approximated by the same functions in the representations (14) and (15), the matrices satisfy the following equalities:


and hence the system (16) reduces to


with and (anti-)block-diagonal matrices of the following form:


Across the vertical interfaces continuity of


is required.

As soon as the function , or in other words and , are known, the functions corresponding to the four other components can be derived as follows:


Note that substituting Equation (23) into the continuity conditions (22) shows that (22) exactly implies the continuity of the relevant electromagnetic components (18).

6 Method of solution

In general the system (20) can be solved by the Finite Element method [10, 11]. It relies on a spatial discretization, i.e. divides the whole computational domain into a number of elements. On each of these elements the unknown function is represented as a superposition of some basis functions. The coefficients of the expansion are found using the weak form of eqn. (20). While this method is very general, it quickly introduces a large number of unknowns.

However, due to common techniques of fabrication many waveguides do not have a completely arbitrary refractive index distribution, but rather one which is piecewise constant along the horizontal axis. The waveguide then can be split in several vertical slices, where the refractive index does not change in the horizontal direction. In each of these layers the general solution of (20) can be written down analytically. Gluing them together across the vertical interfaces will give the desired mode profile.

Both of these methods can be applied to find not only the fundamental, but also higher order modes. In the following we will outline each of these methods in more detail.

6.1 Arbitrary refractive index distribution: Finite Element Method

In case of an arbitrary permittivity distribution (diffused waveguide, waveguide with slanted sidewalls) the matrices depend on , as their elements include overlap integrals with the permittivity . One of the ways to solve the differential equation (20) is by using the Finite Element Method.

By multiplying both sides of (20) from the left by some continuous test vector-function and integrating over one gets the weak form of equation (20):


Then we expand the solution into a finite combination of the basis functions ,


with the dimension of the vector ; the number of consecutive grid points into which the -axis has been divided, and


with, for example, linear basis functions


As eqn. (24) should hold for an arbitrary continuous , we choose it to be one of the basis functions . For and this results in the system of exactly linear equations


where (the subscript here refers to the element of the subvector). Since for any square matrix of dimension


holds, the matrices turn to be of the following form


where the indices and have the same meaning as in the definition of the vector .

The solution of the quadratic eigenvalue problem (28) with as an eigenvalue can be found by introducing an auxiliary vector . (28) can then be transformed into the following linear eigenvalue problem


This is a quite straightforward, but expensive approach, as the dimension of the transformed problem is doubled in comparison to the original one. Other more involved approaches to tackle a quadratic eigenvalue problem can be found e.g. in [12]. We apply standard general eigenvalue solvers as embedded within the LAPACK [13] package. Specialized solvers could be employed, provided that an initial guess for the propagation constant, or a range of possible eigenvalues, are available for the problem at hand. On the other hand, there are situations where all the propagation constants and corresponding functions need to be found together, e.g. if one wants to expand a 3D field in terms of vectorial modes of some channel waveguide, as required for the implementation of transparent boundary conditions [14].

While the entire 2D problem could also be solved directly by means of a Finite Element method, the number of degrees of freedom in such cases would be much higher than when solving the 1D equations (20) using the Finite Element Method; instead of having to use a triangulation of the entire 2D domain, only 1D finite elements are needed; furthermore, the number of degrees of freedom on each node is equal to the number of modes in the expansion, which is typically a small number.

6.2 Piecewise constant refractive index distribution

If a waveguide has a piecewise constant rectangular refractive index profile, it can be divided by vertical lines into slices with constant refractive index distribution along the -direction. In each of these slices the matrices do not depend on . Then (20) can be rewritten in a more familiar manner: Inside each of the slices should satisfy a system of second order differential equations with constant coefficients and a parameter


together with the continuity conditions (22). Moreover the matrices and are block-diagonal in such a way that the equations for the functions and decouple inside each of the slices; coupling occurs only across the vertical interfaces.

Inside each slice a particular solution of the system (32) can be readily written as


with some constants and a vector . By substituting (33) into (32) we find a generalized eigenvalue problem with as an eigenvalue:


So inside each of the slices with uniform permittivity along the -axis the function can be represented as


with eigenvalues and corresponding eigenvectors from (34).

By matching the solutions of the each individual slab across the vertical interfaces using (22) and looking only for exponentially decaying solutions for , one can obtain an eigenvalue problem


The vector consists of all unknown coefficients and from the representations of (35) on all individual slices. The matrix depends on in a non-linear, even non-polynomial way. One of the strategies to tackle this is at first to specify a range of admissible values , where solutions are sought. As we are looking only for propagating modes, with decaying field (35) at , should be not smaller than the biggest eigenvalue of (34) in the left-most and the right-most slabs. At the same time we require that there exists at least one oscillating function in at least one vertical slab. So should be smaller than the biggest eigenvalue of (34) of all the constituting slabs, except the left- and the right-most ones. Once this interval is at hand, we scan through it looking for a such that the matrix has at least one zero eigenvalue. Obviously, to find a non-trivial solution with certain accuracy requires some iterations. Moreover a large step size might lead to missing some roots while scanning the interval.

Once a nontrivial solution , of (36) is at hand, can be reconstructed using (35). And then all field components can be obtained according to expressions (23) together with (14) or (15).

7 Relation with the Effective Index Method

In the following section we are going to show what happens if only a single, TE or TM, slab mode is taken into account in VEIM5 (14). Using the variational reasoning we will rigorously derive an analog to the Effective Index method.

7.1 TE polarization

Let us take only one TE slab mode with propagation constant from a reference slice r with permittivity distribution , and use it to represent the vectorial field profile of the complete waveguide as in eqn. (14). Due to the fact that , according to (20) the unknown function satisfies the eqn.


After some manipulations, using the relations between the modal components , and of the slab mode the above relation can be rewritten as follows




This looks exactly as a TM mode equation, similar to the standard Effective Index Method. In the reference slice one has , and the effective permittivity is equal to the squared effective index of the mode of the reference slice . In other slices this squared effective index is modified by the difference between the local permittivity and that of the reference slice, weighted by the local intensity of the fundamental component of the reference mode profile. Hence, on the contrary to the EIM, even in slices where no guided mode exist, the effective permittivity can still be rigorously defined.

Now it is instructive to see how the mode profile adjusts both in the reference slabs and elsewhere. Inside a slice with constant permittivity , eqn. (38) permits solutions of the form


for arbitrary constants and and with


With the abbreviation from (23) it follows that


By introducing an angle such that , one can write


If we use the principal square roots of and for and , and the principal inverse cosine for eq. 43 can be interpreted as follows. In the slice where the reference slab mode lives , and we find that functions act as a rotation of the slab mode, such that the projection of the propagation constant of this mode onto the -axis will match the global propagation constant . In other slices, in addition to the rotation of the and components of the slab mode, the component is scaled by .

7.2 TM polarization

Analogously, the eqn. (20) can be rewritten for a single TM mode, with a field template as in (14). We now have in eqn. (20) and using the properties of the TM slab mode, the original equation for the unknown function ,


can be rewritten as




This appears to be neither a standard TE nor a TM mode equation, but something in between, with the local refractive index distribution appearing both in the terms with and without derivative. In the reference slice with , the effective permittivity is equal to the squared effective index of the mode of the reference slice and . Contrary to the EIM, even in slices where no guided mode exists quantities that act like effective indices can still be rigorously defined.

What concerns the mode profile