Accurate vectorial finite element mode solver for magneto-optic and anisotropic waveguides

Paolo Pintus

doi:10.1364/OE.22.015737

1. Introduction

Dielectric waveguides are fundamental components of optoelectronic and microwave devices. They have been extensively studied in the past and different kinds of materials have been considered. To perform the waveguide modal analysis, several numerical methods are commonly used to compute the electromagnetic modes such as the finite element method (FEM) [1], the finite difference method (FDM) [2], the method of lines (MoL) [3], and the film mode matching (FMM) [4]. Due to the possibility of using adaptive meshing, FEM shows several advantages. It usually provides a better approximation and requires less memory to store the stiffness matrix with respect to FDM, while it is more appropriate than MoL and FMM for modal analysis of graded index waveguides, such as ion-implanted waveguides [5], and in general for the analysis of waveguides with complex cross-section geometry and refractive index profiles [6]. In addition, it is the most suited to solve deformation and stress problems in solids, like for the case of stress-induced effects in optical waveguides (e.g., silicon waveguides on silicon-on-insulator platform [7]).

In the finite element analysis, the solution can be numerically computed as a linear combination of basis functions which, in electromagnetism, are usually of two kinds: node elements (also called Lagrangian elements) [1, 8–14] and edge elements (also called Nédélec elements) [1,8,15–19]. The use of the edge elements shows mainly three important benefits: i) the spurious solutions can be effectively removed in several electromagnetic problem formulations, ii) the boundary conditions at material interface and conducting surface can be easily imposed, iii) there are no difficulties in treating conducting and dielectric edges and corners related to the field singularities [1,20]. On the other hand, the node elements are more efficient as regards to the storage requirements and the number of floating point operations (FLOPs). In addition, the solutions computed using node elements provide higher accuracy when extremely flat or elongated elements are used in the mesh [20]. It is worth noting that in order to completely get rid of any spurious solutions introduced by the node elements within the range of interest of the guided modes, the penalty function can be added to the functional [8, 21].

In this work, the node elements have been used, since we do not consider waveguides with field singularities (e.g., by using magnetic-field formulation) and assume a zero-field condition on the border. However, the method can be implemented also with edge elements. The mode solver has been developed for straight waveguides, as the one shown in Fig. 1 where the cross-section is in the xy-plane and the light propagates only along the z-direction.

Fig. 1 Cross section of a uniform waveguide loaded with inhomogeneity.

Download Full Size | PDF

Using Cartesian coordinate system, we have considered the following relative permittivity tensor

ε_{r} = (\begin{matrix} ε_{x x} & ε_{x y} & - ε_{x z} \\ ε_{x y} & ε_{y y} & ε_{y z} \\ ε_{x z} & - ε_{y z} & ε_{z z} \end{matrix}),

where the diagonal blocks are symmetric matrices, while the off-diagonal blocks are antisymmetric ones. All the matrix entries are complex. This tensor generalizes the one presented by Konrad [9], permitting to consider also lossy material. In addition, it includes the case presented by Lu and Fernandez [22] and allows us to investigate optical-isolators based on magneto-optical materials and anisotropic waveguides.

For the sake of clarity, let us consider the case of a magneto-optic garnet. Its permittivity tensor can be written as

ε_{r} = (\begin{matrix} ε_{x x} & 0 & 0 \\ 0 & ε_{y y} & 0 \\ 0 & 0 & ε_{z z} \end{matrix}) + K (\begin{matrix} 0 & M_{z} & - M_{y} \\ - M_{z} & 0 & M_{x} \\ M_{y} & - M_{x} & 0 \end{matrix}),

where M_x, M_y and M_z are the magnetization of the garnet and K is a complex parameter which depends on the material [23, 24]. Note that when KM_x or KM_y are not zero, the forward and backward mode propagation constants have different values. More specifically, when KM_x or KM_y are purely imaginary numbers, the forward and backward modes differ only in the phase constant (i.e., the imaginary part of the propagation constant) and a nonreciprocal phase shift effect (NRPS) arises. By exploiting this effect, optical isolators have been designed using Mach-Zehnder Interferometer (MZI) configuration in order to generate constructive interference for forward light and destructive interference for backward light [25,26]. On the other hand, when the real part of KM_x or KM_y is not zero, also the propagation loss with respect to the two directions is different. Optical isolating function can be achieved combining the nonreciprocal loss (NRL) with the gain of an optical amplifier, which compensates for the loss of the forwarding wave, whereas keeps a large loss for the reverse direction [27].

While M_x and M_y are responsible for NRPS and NRL of the transverse magnetic (TM) and transverse electric (TE) modes, respectively, M_z gives rise to TE-TM mode coupling. In integrated optics M_z is usually negligible for several applications and the magnetization vector is all in the xy-plane. As a result, the permittivity tensor of Eq. (2) has the form described by Eq. (1).

Considering the previous hypothesis, the augmented Rayleight-Ritz functional with the penalty function has been derived. By using the node elements, we have computed the eigenvalue problem as a function of the angular frequency ω and propagation constant γ. The γ-formulation, where the frequency is provided as the input parameter and the propagation constant is the output eigenvalue, has been explicitly derived. This formulation avoids the time-consuming iterations necessary to evaluate the propagation constant at a given frequency [8].

By exploring two examples of optical isolators based on magneto-optical waveguides for silicon compatible platform, we show how the shift between forward and backward propagation constants can be directly computed instead of using the pertubation theory [28].

2. Mathematical model

2.1. Differential electromagnetic problem

In a linear, instantaneous and time invariant medium, the magnetic and electric fields can be written as a linear combination of harmonic waves

H (r, t) = H (r) e^{j ω t}, E (r, t) = E (r) e^{j ω t},

where t is the time and r is the position. By replacing the new expressions for the fields in Maxwell’s equations, and combining them with the constitutive relations, we have

\nabla \cdot μ_{0} μ_{r} (r) H (r) = 0,

\nabla \cdot ε_{0} ε_{r} (r) E (r) = ρ (r),

\nabla \times E (r) = - j ω μ_{0} μ_{r} (r) H (r),

\nabla \times H (r) = j ω ε_{0} ε_{r} (r) E (r) + J (r),

where J and ρ are the electric current and charge density, respectively, the constants ε₀ and μ₀ are the electric permittivity and the magnetic permeability in the vacuum, while the tensors ε_r and μ_r are the relative electric permittivity and the relative magnetic permeability of the materials. In addition, to determine the solution of the electromagnetic problem, a set of boundary conditions associated with the domain must be fixed. Here, we refer to them with ℬ(E, H).

By computing E(r) from (4d) and replacing it into (4c), the result is a differential problem for H(r)

{\begin{array}{l} \nabla \times [ε_{r}^{- 1} \nabla \times H] - {(\frac{ω}{c})}^{2} μ_{r} H = \nabla \times ε_{r}^{- 1} J, & in V, \\ ℬ (H) = 0, & on S, \end{array}

where V is the domain and S is the boundary. In Eq. (5), we used the relation

c = 1 / \sqrt{ε_{0} μ_{0}}

for the speed of light in the vacuum and we omit the explicit dependence on the space variable r for the sake of simplicity. In a similar way, it is possible to derive an equivalent problem for E.

2.2. Variational formulation: non-self-adjoint problem

As it is known from the variational methods, it is possible to construct a functional linked to the differential problem such as the stationary point of that functional is the exact solution of the differential problem [1]. In our case, the relative permittivity tensor defined in eq. (1) is neither symmetric (i.e., $ε_{r} = ε_{r}^{T}$ ) nor hermitian (i.e., $ε_{r} = ε_{r}^{†}$ ). As a result, the differential operator $\nabla \times [ε_{r}^{- 1} (r) \nabla \times \cdot]$ is not symmetric or self-adjoint [29]. To compute the related functional an auxiliary problem is introduced, which is also called adjoint problem [30]. Considering the original problem (5), its adjoint problem is

{\begin{array}{l} \nabla \times [{(ε_{r}^{a})}^{- 1} \nabla \times H^{a}] - {(\frac{ω}{c})}^{2} μ_{r}^{a} H = \nabla \times {(ε_{r}^{a})}^{- 1} J^{a}, & in V, \\ ℬ^{a} (H^{a}) = 0, & on S, \end{array}

where H^a is the adjoint magnetic fields, J^a is the adjoint source distribution, and ℬ^a(H^a) = 0 are the adjoint boundary conditions. Note that J^a and ℬ^a(H^a) = 0 are chosen according to the direct problem, while the relative tensors

ε_{r}^{a}

and

μ_{r}^{a}

are equal to

ε_{r}^{T}

and

μ_{r}^{T}

when we consider the real inner product, or

ε_{r}^{†}

and

μ_{r}^{†}

in case we use the complex inner product [30].

By considering the perfect electric conductor (PEC) or perfect magnetic conductor (PMC) boundary conditions

(\nabla \times H) \times i_{n} = 0, PEC,

H \times i_{n} = 0, PMC,

where i_n is the normal unit vector with respect to the surface. It can be proved that the solutions of the differential equation (5) is the stationary point of the functional

\begin{array}{r} F (H, H^{a}) = \frac{1}{ε_{0}} ∭_{V} [\nabla \times H^{a} \cdot ε_{r}^{- 1} \nabla \times H] d V - ω^{2} μ_{0} ∭_{V} μ_{r} H \cdot H^{a} d V + \\ - \frac{1}{ε_{0}} ∭_{V} \nabla \times ε_{r}^{- 1} J \cdot H^{a} d V - \frac{1}{ε_{0}} ∭_{V} H \cdot \nabla \times {(ε_{r}^{T})}^{- 1} J^{a} d V, \end{array}

where the real inner product has been considered. If the complex inner product is used, the adjoint field (H^a) and current (J^a) must be replaced by their complex conjugate values (H^a* and J^a*), while

ε_{r}^{†}

and

μ_{r}^{†}

take the place of

ε_{r}^{T}

and

μ_{r}^{T}

in the previous formula [30]. A similar formulation can be derived for E and the electric adjoint field E^a.

2.3. Waveguide mode problem

In our work, we consider that the materials are non-magnetic (μ_r = 1). As a result, the normal and tangent component of H are continuous across any boundary separating two different media. For this reason, the formulation in term of the magnetic field is preferred. Moreover, because a mode is a transverse electromagnetic wave which propagates in the waveguide without sources, we set J = 0 and J^a = 0.

By considering an infinite straight waveguide, like the one shown in Fig. 1, the geometry is invariant along the z-axis and the electric and magnetic fields can be written as

E = E (x, y) e^{j ω t - γ z}, H = H (x, y) e^{j ω t - γ z},

where γ = jβ + α is the propagation constant, β the phase constant and α the attenuation. Similarly, the adjoint fields are

E^{a} = E^{a} (x, y) e^{j ω t - γ^{a} z}, H^{a} = H^{a} (x, y) e^{j ω t - γ^{a} z},

where γ^a = −γ in the real inner product case, while γ^a = −γ^* when the complex inner product is used. Indeed, in order to have a useful form for the functional and to get an expression independent of z, the exponential (z-dependent) terms of the original and adjoint vector fields must cancel out [30–32]. To compute the functional in Eq. (8), we need to find the relationship between the direct and the adjoint fields. Equation (4c) and (4d) for the direct field are

(\begin{matrix} 0 & γ & \frac{\partial}{\partial y} \\ - γ & 0 & - \frac{\partial}{\partial x} \\ - \frac{\partial}{\partial y} & \frac{\partial}{\partial x} & 0 \end{matrix}) (\begin{matrix} E_{x} \\ E_{y} \\ E_{z} \end{matrix}) = - j ω μ_{0} (\begin{matrix} H_{x} \\ H_{y} \\ H_{z} \end{matrix}),

(\begin{matrix} 0 & γ & \frac{\partial}{\partial y} \\ - γ & 0 & - \frac{\partial}{\partial x} \\ - \frac{\partial}{\partial y} & \frac{\partial}{\partial x} & 0 \end{matrix}) (\begin{matrix} H_{x} \\ H_{y} \\ H_{z} \end{matrix}) = j ω ε_{0} (\begin{matrix} ε_{x x} & ε_{x y} & - ε_{x z} \\ ε_{x y} & ε_{y y} & ε_{y z} \\ ε_{x z} & - ε_{y z} & ε_{z z} \end{matrix}) (\begin{matrix} E_{x} \\ E_{y} \\ E_{z} \end{matrix}),

where for the curl operator we use the matrix expression and the partial derivative with respect to z has been replaced by −γ. Adopting the real inner product (i.e., γ^a = −γ), the adjoint problem is

(\begin{matrix} 0 & - γ & \frac{\partial}{\partial y} \\ γ & 0 & - \frac{\partial}{\partial x} \\ - \frac{\partial}{\partial y} & \frac{\partial}{\partial x} & 0 \end{matrix}) (\begin{matrix} E_{x}^{a} \\ E_{y}^{a} \\ E_{z}^{a} \end{matrix}) = - j ω μ_{0} (\begin{matrix} H_{x}^{a} \\ H_{y}^{a} \\ H_{z}^{a} \end{matrix}),

(\begin{matrix} 0 & - γ & \frac{\partial}{\partial y} \\ γ & 0 & - \frac{\partial}{\partial x} \\ - \frac{\partial}{\partial y} & \frac{\partial}{\partial x} & 0 \end{matrix}) (\begin{matrix} H_{x}^{a} \\ H_{y}^{a} \\ H_{z}^{a} \end{matrix}) = j ω ε_{0} (\begin{matrix} ε_{x x} & ε_{x y} & ε_{x z} \\ ε_{x y} & ε_{y y} & - ε_{y z} \\ - ε_{x z} & ε_{y z} & ε_{z z} \end{matrix}) (\begin{matrix} E_{x}^{a} \\ E_{y}^{a} \\ E_{z}^{a} \end{matrix}) .

The matrix blocks which have been highlighted in the adjoint problem (12) have opposite sign with respect to the ones in the original problem (11). Using the following relations

(\begin{matrix} E_{x}^{a} \\ E_{y}^{a} \\ E_{z}^{a} \end{matrix}) = (\begin{matrix} - E_{x} \\ - E_{y} \\ E_{z} \end{matrix}), (\begin{matrix} H_{x}^{a} \\ H_{y}^{a} \\ H_{z}^{a} \end{matrix}) = (\begin{matrix} H_{x} \\ H_{y} \\ - H_{z} \end{matrix}),

in the equation system (12) we obtain the system (11). The previous relations allow us to link the components of the original fields with those of the adjoint ones. At this point, it is easy to compute the Rayleigh-Ritz functional (8). The relationship (13) is also verified when both ε_r and μ_r have the form of Eq. (1), which is the case of magnetized ferrites [33]. In that case, the variational formulation with respect to E is preferred.

Let us note that from Eq. (8), the solution of the problem H and H^a belongs to the Sobolev space H(curl).

2.4. Interpolating function

Because the electromagnetic problem for a straight waveguide is independent of the direction of propagation (z-axis), only the transverse cross-section is discretized. Considering a generic mesh element e, the interpolating function on it can be written as

H_{n}^{e} (x, y) = \sum_{i = 1}^{m} [h_{x i}^{e} N_{i}^{e} (x, y) i_{x} + h_{y i}^{e} N_{i}^{e} (x, y) i_{y} + h_{z i}^{e} N_{i}^{e} (x, y) i_{z}] e^{j ω t - γ z},

where m is the number of nodes in the single element,

N_{i}^{e} (x, y)

is the basis function for the node i = 1,...,m, and (

h_{x i}^{e}, h_{y i}^{e}, h_{z i}^{e}

) are the values of

H_{n}^{e}

at the interpolation node i. The approximate field H_n in the whole cross-section is then

H (x, y) ≃ H_{n} (x, y) = \sum_{e = 1}^{N} H_{n}^{e} (x, y),

where N is the total number of mesh elements. In our work, we have considered a triangular mesh and quadratic interpolating function (i.e., m = 6).

As we already mentioned, one of the drawbacks associated with the node element approach is the appearance of spurious or nonphysical solutions. In fact, when the magnetic field H solves the curl-curl equation (5), it also solves the constraint

\nabla \cdot μ_{0} H (r) = 0 .

When we are looking for an approximate solution, this is not always true and some non-feasible solutions are introduced. To force the numerical solutions to satisfy all equations (5), the constraint should be taken into account using the augmented functional F̃(H, H^a)

\tilde{F} (H, H^{a}) = F (H, H^{a}) + \frac{α_{p}}{ε_{0}} ∭_{V} \nabla \cdot H^{a} \nabla \cdot H d V .

The added term is called the penalty function and has been introduced by Rahman and Davies [21]. Equation (17) implies that H, H^a belong to the Sobolev space H(curl) ∩ H(div). The constant α_p is a free parameter and it is usually chosen equal to 1.

2.5. Discretized functional

Using the approximation of the previous paragraph and introducing the tensor $p = ε_{0}^{- 1} ε_{r}^{- 1}$ , we have calculated the functional (17). Note that the integral along the z-axis is neglected because the argument is independent of this variable. To make the computation clearer, let us consider the integrals on a generic element e with surface S^e. In addition, we assume p is constant in each element. At this point, let us introduce a more compact matrix notation by defining the vectors

h_{x}^{e} = (\begin{matrix} h_{x 1}^{e} \\ h_{x 2}^{e} \\ ⋮ \\ h_{x m}^{e} \end{matrix}), h_{y}^{e} = (\begin{matrix} h_{y 1}^{e} \\ h_{y 2}^{e} \\ ⋮ \\ h_{y m}^{e} \end{matrix}), h_{z}^{e} = (\begin{matrix} h_{z 1}^{e} \\ h_{z 2}^{e} \\ ⋮ \\ h_{z m}^{e} \end{matrix}),

and the matrices R^e, J^e, N^e, D^e, E^e, and Z^e with elements

\begin{array}{l} r_{i j}^{e} = \iint_{S^{e}} N_{i}^{e} N_{j}^{e} d x d y, & d_{i j}^{e} = \iint_{S^{e}} \frac{\partial N_{i}^{e}}{\partial y} \frac{\partial N_{j}^{e}}{\partial y} d x d y, \\ j_{i j}^{e} = \iint_{S^{e}} N_{i}^{e} \frac{\partial N_{j}^{e}}{\partial y} d x d y, & e_{i j}^{e} = \iint_{S^{e}} \frac{\partial N_{i}^{e}}{\partial x} \frac{\partial N_{j}^{e}}{\partial x} d x d y, \\ n_{i j}^{e} = \iint_{S^{e}} N_{i}^{e} \frac{\partial N_{j}^{e}}{\partial x} d x d y, & z_{i j}^{e} = \iint_{S^{e}} \frac{\partial N_{i}^{e}}{\partial y} \frac{\partial N_{j}^{e}}{\partial x} d x d y, \end{array}

where the capital letter is used for the matrix and its corresponding lower case letter indicates its generic element [9]. The matrix entry with index (i, j) depends only on the basis functions

N_{i}^{e} (x, y)

and

N_{j}^{e} (x, y)

on the element e. Note that the matrices R^e, D^e, and E^e are symmetric.

From Eq. (17), we can see that the functional F̃(H, H^a) is the sum of three integrals. By computing the first one, we obtain

\begin{array}{l} \iint_{S^{e}} [\nabla \times H^{a} \cdot p \nabla \times H] d x d y = \\ {(\begin{matrix} h_{x}^{e} \\ h_{y}^{e} \\ h_{z}^{e} \end{matrix})}^{T} [γ^{2} (\begin{matrix} - p_{y y} R^{e} & p_{y x} R^{e} & 0 \\ p_{x y} R^{e} & - p_{x x} R^{e} & 0 \\ 0 & 0 & 0 \end{matrix}) + \\ + γ (\begin{matrix} p_{z y} {J^{e}}^{T} - p_{y z} J^{e} & p_{y z} N^{e} - p_{z x} {J^{e}}^{T} & p_{y x} J^{e} - p_{y y} N^{e} \\ p_{x z} J^{e} - p_{z y} {N^{e}}^{T} & p_{z x} {N^{e}}^{T} - p_{x z} N^{e} & p_{x y} N^{e} - p_{x x} J^{e} \\ p_{x y} {J^{e}}^{T} - p_{y y} {N^{e}}^{T} & p_{y x} {N^{e}}^{T} - p_{x x} {J^{e}}^{T} & 0 \end{matrix}) \\ + (\begin{matrix} p_{z z} D^{e} & - p_{z z} Z^{e} & p_{z y} Z^{e} - p_{z x} D^{e} \\ - p_{z z} {Z^{e}}^{T} & p_{z z} E^{e} & p_{x z} {Z^{e}}^{T} - p_{z y} E^{e} \\ p_{x z} D^{e} - p_{y z} {Z^{e}}^{T} & p_{y z} E^{e} - p_{x z} Z^{e} & L_{z z}^{e} \end{matrix})] (\begin{matrix} h_{x}^{e} \\ h_{y}^{e} \\ h_{z}^{e} \end{matrix}), \end{array}

where L_zz = p_xyZ^e + p_yxZ^eT − p_xxD^e − p_yyE^e.

Similarly, we computed the second term

{\iint_{S^{e}} H^{a} \cdot H d x d y = (\begin{matrix} h_{x}^{e} \\ h_{y}^{e} \\ h_{z}^{e} \end{matrix})}^{T} (\begin{matrix} R^{e} & 0 & 0 \\ 0 & R^{e} & 0 \\ 0 & 0 & - R^{e} \end{matrix}) (\begin{matrix} h_{x}^{e} \\ h_{y}^{e} \\ h_{z}^{e} \end{matrix}),

while the last integral becomes

\begin{array}{l} \iint_{S^{e}} \nabla \cdot H^{a} \nabla \cdot H d x d y = \\ {(\begin{matrix} h_{x}^{e} \\ h_{y}^{e} \\ h_{z}^{e} \end{matrix})}^{T} [γ^{2} (\begin{matrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & R^{e} \end{matrix}) - γ (\begin{matrix} 0 & 0 & {N^{e}}^{T} \\ 0 & 0 & {J^{e}}^{T} \\ N^{e} & J^{e} & 0 \end{matrix}) + (\begin{matrix} E^{e} & {Z^{e}}^{T} & 0 \\ Z^{e} & D^{e} & 0 \\ 0 & 0 & 0 \end{matrix})] (\begin{matrix} h_{x}^{e} \\ h_{y}^{e} \\ h_{z}^{e} \end{matrix}) . \end{array}

Considering the previous results, the augmented functional over the element e is approximated by the quadratic form

{\tilde{F} (H, H^{a}) |}_{e} \approx {(\begin{matrix} h_{x}^{e} \\ h_{y}^{e} \\ h_{z}^{e} \end{matrix})}^{T} [γ^{2} M^{e} + γ C^{e} + K^{e} - ω^{2} μ_{0} T^{e}] (\begin{matrix} h_{x}^{e} \\ h_{y}^{e} \\ h_{z}^{e} \end{matrix}) .

where we have introduced the 18 × 18 matrices M^e, C^e, K^e and T^e. Note that M^e, C^e and T^e are the sum of the matrices multiplied by γ², the sum of the matrices multiplied by γ, and the sum of the matrices multiplied by ω², respectively, while the sum of the remaining ones is K^e.

At this point we move from the local computation performed on the generic element e to the global one that defines the matrices associated to the discretized problem. The total functional is the sum of the contributions of all of the elements

\tilde{F} (H, H^{a}) = \sum_{e = 1}^{N} {\tilde{F} (H, H^{a}) |}_{e} .

For this purpose, let us extend the previous notation defining three vectors with unknown coefficients

h_{x} = (\begin{matrix} h_{x 1} \\ h_{x 2} \\ ⋮ \\ h_{x \hat{N}} \end{matrix}), h_{y} = (\begin{matrix} h_{y 1} \\ h_{y 2} \\ ⋮ \\ h_{y \hat{N}} \end{matrix}), h_{z} = (\begin{matrix} h_{z 1} \\ h_{z 2} \\ ⋮ \\ h_{z \hat{N}} \end{matrix}),

where N̂ is the number of the mesh nodes. Therefore, the discrete functional computed in the cross-section is

\tilde{F} (H, H^{a}) \approx {(\begin{matrix} h_{x} \\ h_{y} \\ h_{z} \end{matrix})}^{T} [γ^{2} M + γ C + K - ω^{2} μ_{0} T] (\begin{matrix} h_{x} \\ h_{y} \\ h_{z} \end{matrix}),

where M, C, K and T are 3N̂ × 3N̂ matrices. Those matrices can be easily computed “assembling” all of the matrices M^e, C^e, K^e, and T^e for e = 1,...,N. It is worth noting that the matrix entry with index (i, j) is non-zero only when the basis functions are non-zero at least on one shared element. Due to this fact, the matrices in (26) are sparse, with very few non-zero elements.

2.6. Eigenvalue problems

Since the solution of the differential problem (5) is the stationary point of F̃(H, H^a), we derive equation (26) with respect to the unknown vector and obtain

[γ^{2} M + γ C + K - ω^{2} μ_{0} T] h = 0,

where

h = {(h_{x}^{T} h_{y}^{T} h_{z}^{T})}^{T}

. According to the known/unknown parameters, an ω-formulation or a γ-formulation can be derived [8]. In the first case, the propagation constant is provided as an input of the problem and the previous equation becomes

[γ^{2} M + γ C + K] h = ω^{2} μ_{0} T h,

where (ω, h) is the eigenvalue-eigenvector pair of the generalized eigenvalue problem (GEP) in Eq. (28). Vice versa, fixing the angular frequency ω, the Eq. (27) is a quadratic eigenvalue problem (QEP) [34, 35], where h is the eigenvector and γ its eigenvalue.

The number of rows/columns of the matrices in Eq. (28) is three times the number of the mesh nodes N̂. In the literature, different approaches have been derived in order to reduce the memory and computational efforts, e.g. the formulation in terms of the transverse magnetic field components. However, to the best of our knowledge, this formulation has been obtained only in the case of ε_xz = ε_yz = 0 [1, 8, 22, 36, 37], while its derivation is not straightforward in the more general case. Moreover, the full vectorial formulation preserves the matrix sparsity and provides a higher accuracy on the computation of the longitudinal component H_z, which is in fact directly evaluated [8]. It is worth noting that an accurate value of all field components is important especially in silicon photonics, where the concepts of pure TE and TM modes are undermined due to the high index contrast.

An easy way to solve the QEP consists in transforming it into an equivalent generalized eigenvalue problem [34]. By introducing the vector

u = γ h,

we obtain

(\begin{matrix} 0 & I \\ - K + ω^{2} μ_{0} T & - C \end{matrix}) (\begin{matrix} h \\ u \end{matrix}) = γ (\begin{matrix} I & 0 \\ 0 & M \end{matrix}) (\begin{matrix} h \\ u \end{matrix}) .

As a result, both the ω-formulation (Eq. (28)) and the γ-formulation (Eq. (30)) have been transformed into a generalized eigenvalue problem

A v = σ B v .

Note that the size of the computational problem in the case of fixed optical frequency is doubled compared to the formulation where the propagation constant is set.

The eigenvector-eigenvalue pairs, which are the numerical solutions of the electromagnetic problem, can be found using Krylov methods. Such methods allow us to compute few eigenvalues of large sparse matrices, e.g. the eigenvalues which have the largest or the smallest magnitude. If we roughly know the eigenvalues we are interested in, we can shift the spectrum close to them and then apply the method, looking for the smallest magnitude eigenvalues. Called such approximation σ₀, we subtract σ₀Bv from either term in Eq. (31)

A v - σ_{0} B v = σ B v - σ_{0} B v .

Collecting the common terms, we obtain

(A - σ_{0} B) v = (σ - σ_{0}) B v .

The generalized eigenvalue problems (31) and (33) have the same eigenvectors, but shifted eigenvalues

(v, σ) and (v, σ - σ_{0}) .

In the ω–formulation, the spectrum is shifted close to the angular frequency we are considering. In the γ–formulation, β is unknown, but its value for the propagating modes must respect the inequality

n_{\min} k_{0} \leq β \leq n_{\max} k_{0},

where n_min and n_max are the minimum and maximum refractive indices in the problem under investigation, and k₀ is the wavenumber in the vacuum (k₀ = 2π/λ). In this case we shift the spectrum close to σ₀ = n_maxk₀ when we want to compute the forward propagating modes, while we move the spectrum close to σ₀ = −n_maxk₀ when we are looking for the backward propagating modes. From Eq. (35) it is useful to define the effective index as n_eff = β/k₀ which can be easily compared with the refractive index of the waveguide.

To improve the numerical stability of the method, the QEP is properly rescaled [38, 39]. By introducing the 2-norm of the following matrices

m_{2} = {‖ M ‖}_{2}, c_{2} = {‖ C ‖}_{2}, s_{2} = {‖ K - ω^{2} μ_{0} T ‖}_{2},

we define

ξ = \sqrt{\frac{s_{2}}{m_{2}}}, δ = \frac{2}{s_{2} + c_{2} ξ} .

The eigenvalues and the matrices of the scaled eigenvalue problem are

\tilde{γ} = \frac{γ}{ξ}, \tilde{M} = ξ^{2} δ M, \tilde{C} = ξ δ C, \tilde{K} = δ K, \tilde{T} = δ T .

2.7. Lossy case: fully complex formulation

The described mode solver formulation is fully complex and the propagation loss can be directly computed. In order to evaluate the radiation loss of the waveguide, the perfectly matched layer method (PML) can be used [40], and the parameters can be easily estimated using the method presented in [14, 41].

2.8. Lossless case: real formulation

Ordinarily, the field E and H have complex components. However, the problem can be simplified considering a particular case for the permittivity tensor

ε_{r} = (\begin{matrix} ε_{x x} & ε_{x y} & - j ε_{x z} \\ ε_{x y} & ε_{y y} & j ε_{y z} \\ j ε_{x z} & - j ε_{y z} & ε_{z z} \end{matrix}),

where ε_ij are real for i, j = x, y, z. In this case, which belongs to the more general case described before, the tensor is self-adjoint (or hermitian) and the waveguide is lossless (i.e., α = 0). The problem is still quite general and the lossless magneto-optic and anisotropic materials can be described by a simplified formulation. Writing the magnetic and electric fields as

E = [E_{x} (x, y) i_{x} + E_{y} (x, y) i_{y} + j E_{z} (x, y) i_{z}] e^{j ω t - j β z},

H = [H_{x} (x, y) i_{x} + H_{y} (x, y) i_{y} + j H_{z} (x, y) i_{z}] e^{j ω t - j β z},

makes the problem (17) fully real [1, 9] and this halves the memory demand.

3. Numerical results

In this section, we are presenting some numerical results obtained by the described mode solver. The method has been fully programmed in MATLAB. The integrands in Eq.s (19) are polynomials over triangular elements and can be computed exactly [9]. The routine eigs has been used to numerically solve the eigenvalue problems. It provides the reverse communication required by the Fortran library ARPACK [42,43]. ARPACK is a collection of Fortran77 subroutines designed to solve large scale eigenvalue problems. This software is based upon an algorithmic variant of the Arnoldi process called the implicitly restarted Arnoldi method (IRAM) [44,45]. In the ω–formulation, because both of the matrices in the generalized eigenvalue problem are symmetric, the IRAM is replaced by the implicitly restarted Lanczos method (IRLM) [43].

3.1. Nonreciprocal phase shift in lossless magneto-optical waveguides

Here, we present the modal analysis results of a nonreciprocal waveguide that has been used to perform silicon-based optical isolator for the TM mode using both Mach-Zender interferometer [26, 46, 47] and micro-ring resonator configurations [48, 49]. In the former photonic circuit, the device is designed to generate constructive interference for the forward light and destructive interference for the backward light, while in the latter the different propagation constant for the clockwise (CW) and the counterclockwise (CCW) modes results in a different resonant wavelength for the two directions.

The waveguide cross-section is shown in Fig. 2. The silicon waveguide is fabricated on a silicon-on-insulator (SOI) wafer, having refractive index n_Si = 3.45 and n_SiO₂ = 1.46, respectively. In order to achieve NRPS, the top-layer of the waveguide is bonded with a cerium-substituted yttrium iron garnet (Ce:YIG) (n_Ce:YIG = 2.22, [50]) grown on a (Ca,Mg,Zr)-substituted gadolinium gallium garnet (SGGG), (n_SGGG = 1.97).

Fig. 2 Magneto-optic waveguide cross section.

Download Full Size | PDF

The Ce:YIG is a magneto-optic (MO) garnet used in silicon photonics to perform optical isolators and circulators. By applying a static magnetic field along the horizontal axis (x-direction in Fig. 2), the Ce:YIG permittivity tensor results

ε_{r} = (\begin{matrix} ε_{x x} & 0 & 0 \\ 0 & ε_{y y} & j ε_{y z} \\ 0 & - j ε_{y z} & ε_{z z} \end{matrix}),

where

ε_{x x} = ε_{y y} = ε_{z z} = n_{Ce : YIG}^{2}

, and the off-diagonal element ε_yz is responsible for the MO effect. The latter term is related to the Faraday rotation constant θ_F by

ε_{y z} = \frac{2 n_{Ce : YIG} θ_{F}}{k_{0}} .

For calculation, we assume θ_F = 4000°/cm at λ = 1550nm [48–50], and then ε_yz = 7.65 · 10⁻³. However, in high quality Ce:YIG crystal, a Faraday-rotation coefficient as large as 4500°/cm can be measured [51].

By solving the curl-curl equation for the magnetic field, as described in section 2, we have computed the three components of the magnetic field and the phase constant, for the forward and backward propagating waves, respectively. The NRPS can be directly computed as

Δ β = | β^{+} | - | β^{-} |,

where + and − refer to the forward and backward propagating directions, respectively. TE and TM forward propagating modes are shown in Fig. 3, where we assumed a 600nm × 215nm silicon waveguide cross-section, with a 380nm thick Ce:YIG layer above it. Those values maximize the NRPS for the TM mode, as it is explained in the following.

Fig. 3 Magnetic field H for the TE-mode (first row) and TM-mode (second row) at 1550nm. The fields have been normalized in term of power and the profiles are shown in H/μm.

Download Full Size | PDF

Because the value of the off-diagonal terms in the permittivity tensor are usually small (e.g., $ε_{y z} / n_{Ce : YIG}^{2} ≃ 0.015$ ), it is very important to optimize the waveguide cross-section in order to maximize Δβ and to achieve the highest isolation. In Fig. 4 we report the NRPS computed for different silicon and Ce:YIG layer thicknesses (h_MO and h_Si in Fig. 2), while we consider a 600nm wide silicon waveguide in order to guarantee the single mode regime (w_Si in Fig. 2).

Fig. 4 NRPS for the TM mode at 1550nm.

Download Full Size | PDF

Note that NRPS is maximized when the maximum of |H_x|² is located close to the boundary of the two materials, similarly to what happens using garnet double layers with opposite Faraday rotation values [52]. The same calculation for the TE mode shows a negligible split between the phase constants of the propagating and counterpropagating modes.

Considering h_MO = 380 nm, h_Si = 215 nm and w_Si = 600 nm, we have performed the modal analysis for λ within the range: [500nm, 3000nm]. Figure 5 shows the dispersion curves of the first ten modes with respect to the wavelength. As it is known, the higher the frequency (the lower the wavelength), the larger the number of modes that are confined in the structure. Moreover, the modes with an effective index higher than the refractive index of the cladding (1.97 for SGGG) are confined in the silicon-waveguide core and can propagate, while the ones which have n_eff < 1.97 are below cut-off and cannot propagate.

Fig. 5 Dispersion curves (eigenvalue spectrum) as a function of the wavelength (λ = 2πc/ω) for the waveguide under examination.

Download Full Size | PDF

3.2. Finite element method vs. perturbation method

In the literature, the most used method for investigating the nonreciprocal effect is the perturbation method [28, 47, 52–55]. Such a method can be applied when the off-diagonal elements of the permittivity tensor are rather small, which in our case means ε_yz ≪ ε_xx, ε_yy, ε_zz. When this happens, ε_yz can be assumed as a small perturbation on the permittivity tensor and the shift between forward and backward propagation constants can be estimated by applying the perturbation theory [28]. Using this approach, the electromagnetic fields are computed considering diagonal permittivity tensor and the phase constant difference is estimated like

Δ β_{p m} = \frac{2 ω ε_{0} \iint_{S} E^{*} Δ ε E d x d y}{\iint_{S} [E^{*} \times H + E \times H^{*}] \cdot i_{z} d x d y},

where Δε is a matrix which contains all the off-diagonal entries.

The previous equation can be further simplified using the semivectorial approximation (i.e., E_x = 0 and E_z = j∂_yE_y/β for the TM mode and E_y = 0 and E_z = j∂_xE_x/β for the TE mode) [47, 53, 55]. Considering our example (TM mode), Eq. (44) can be rewritten as follows

Δ β_{s v} = \frac{2 β}{ω ε_{0}} \frac{\iint_{S} ε_{y z} \frac{\partial}{\partial y} \frac{{| H_{x} |}^{2}}{n^{4}} d x d y}{\iint_{S} [E \times H^{*} + E^{*} \times H] \cdot i_{z} d x d y},

where n is the refractive index.

To conclude our analysis, we compare the NRPS computed with our model with the ones calculated using Eq. (44) and Eq. (45). In Fig. 6(a), the percentage difference between Δβ and Δβ_pm is reported for different values of layer thickness. From that figure, we can easily see that the difference between the two methods is less than 3%. Similarly, a relative difference larger than 10% is computed between Δβ and Δβ_sv, as shown in Fig. 6(b). Such a big difference can be explained considering that the E_x component is neglected for the TM mode in Eq. (45). Indeed, in a high index contrast system, like silicon-photonics, the concepts of pure TE and TM modes is undermined and all three field components become non-negligible and comparable in magnitude.

Fig. 6 Relative change of Δβ_pm 6(a) and Δβ_sv 6(b) with respect to the NRPS computed directly.

Download Full Size | PDF

Finally, in Fig. 7 we compare the three methods considering h_MO = 500nm and varying the silicon thickness only.

Fig. 7 NRPS computation assuming h_MO = 500nm and varying only the silicon thickness. In the plot we compared the three methods.

Download Full Size | PDF

3.3. Nonreciprocal loss in magneto-optical waveguides

Optical isolation based on nonreciprocal loss can be implemented by considering a magneto-optic material, which provides the nonreciprocal loss, and an optical waveguide amplifier, which compensates the forward propagation loss. Several configurations have been proposed for both indium phosphide [56–58] and silicon photonics technologies [59].

In this section, we present the modal analysis of a nonreciprocal loss waveguide for hybrid silicon photonics integrated circuit platforms. The silicon waveguide is defined on SOI wafers by etching 400nm deep trenches into a 700nm thick silicon layer [60]. The nonreciprocal loss is induced by a thin film of iron (Fe) as schematically shown in Fig. 8. The iron layer is separated from the silicon waveguide by a thin titania (TiO₂) layer to reduce and optimize the Fe-layer loss. The propagation loss for the forward light can be effectively compensated by bonding an InGaAsP multi quantum well (MQW) active layer above the silicon waveguide.

Fig. 8 NRL isolator waveguide cross section.

Download Full Size | PDF

In our simulations, we assume a 30nm-thick layer of titania and we varied the thickness of the iron between 25nm and 100nm. The refractive indices of the TiO₂ and Fe at 1550nm are assumed to be 2.1 and 3.17 + 5.27i, respectively [57]. The nonreciprocal effect can be induced on the TE mode by applying a magnetic field along the vertical direction (i.e., the y-axis in Fig. 8). For this case, the Fe permittivity tensor results

ε_{r} = (\begin{matrix} ε_{x x} & 0 & ε_{x z} \\ 0 & ε_{y y} & 0 \\ - ε_{x z} & 0 & ε_{z z} \end{matrix}),

where

ε_{x x} = ε_{y y} = ε_{z z} = n_{Fe}^{2}

, and the off-diagonal element ε_xz = 3.15 + 1.8i when the magnetization of the Fe layer is saturated [57, 61]. It is worth mentioning that, although iron is a ferromagnetic material at low frequency, at optical frequencies the magnetic susceptibility ceases to have any physical meaning and its magnetic permeability has been assumed equal to μ₀ [23, 62, 63].

The results of the modal analysis are reported in Fig. 9, where the amplitude of the magnetic field components of the TE-mode are shown.

Fig. 9 Amplitude of the magnetic field H for the TE-mode at 1550nm. The amplitude of the field components have been normalized in term of power and the profiles are shown in H/μm.

Download Full Size | PDF

In Fig. 10(a), the propagation loss for both the forward and backward modes is reported for different values of the Fe-layer thickness. The thicker the iron layer, the larger the difference between the loss in the two directions and, as a consequence, the higher the optical isolation. On the other hand, a thicker Fe-film produces larger loss that needs to be compensated by a MQW layer. With our model, we directly compute the NRL and NRPS as

Δ α = | α^{+} | - | α^{-} |,

Δ β = | β^{+} | - | β^{-} |,

where the superscript ± refer to the forward and backward waves, respectively. Their values are reported in Fig. 10(b). It is worth noting that the longer the waveguide, the larger the nonreciprocal effect. As it can be seen from the figure, a 25nm-thick iron film will produce a NRL of about −7dB/mm, while a four times thicker layer generates an isolation of almost −12dB/mm. Those values are consistent with those reported in the literature [56, 57, 59].

Fig. 10 Nonreciprocal effects on the TE-mode in the NRL-waveguide.

Download Full Size | PDF

Also in this case, a compared analysis with the perturbation method is performed. Let us note that in the lossy case, the formula (44) must be replaced by the more general one

Δ γ_{p m} = Δ α_{p m} + j Δ β_{p m} = \frac{2 ω ε_{0} \iint_{S} E^{a} Δ ε E d x d y}{\iint_{S} [E^{a} \times H - E \times H^{a}] \cdot i_{z} d x d y},

where E^a and H^a are the adjoint fields [28]. The results of this comparison are reported in table 1.

Table 1. Comparison between the two methods.

View Table

From table 1, we can observe that the results of the two methods are in good agreement. In addition, the relative difference between the two methods is smaller than the ones computed in the previous case. This can be understood because the layer of iron is much smaller than the silicon waveguide and only a small percentage of the field is confined in it. This is in agreement with the basic assumption of the perturbation method.

It is worth noting that the proposed method allows us to directly compute the NRPS and the NRL for complex structures based on new materials characterized by higher Faraday rotation constants like, for example, the bismuth iron garnet (BiIG) [64], in which the perturbation method could be rather inaccurate.

4. Conclusions

We have presented a rigorous full vectorial mode solver based on the finite element method. Such a method allows us to study a generic magneto-optic and anisotropic lossy straight waveguide. After having computed the Rayleigh-Ritz functional for the non-self-adjoint case, the finite element method has been implemented using the node based second order shape functions. Furthermore, the penalty function has been introduced to remove the spurious solutions from the eigenvalue spectral region of interest.

Considering the discretized eigenvalue problem, the γ-formulation has been derived, avoiding the time-consuming iterations necessary in previous formulations for evaluating the mode propagation constants at a given operating frequency [8]. This formulation allows speeding up simulations and improving the convergence of the eigenvalues at the same time. Moreover, all the modes and the propagation constants are computed at the same time and not one after the other as performed in the case of iterative approaches.

Concerning the magneto-optic materials the accuracy of the method has been compared with the perturbation method, the most used technique for investigating the nonreciprocal effect, providing satisfactory results. The mesh adaptivity of the finite element method allows us to efficiently compute the modes of graded index waveguides, like the anisotropic ion-implanted waveguides reported in [5], as well as the modes of high refractive index contrast slot waveguides [65,66]. The described mode solver can be effectively employed also to design microring-based optical isolators [49,67,68], where the mode of a bent waveguide with a large ring radius is approximated by a lateral shift of the straight waveguide mode [69, 70]. The presented approach provides then an important tool for integrated optical device modal analysis and design.

Acknowledgments

The author would like to thank Prof. Fabrizio Di Pasquale from Scuola Sant’Anna (Italy), Dr. Claudio Bonati from Istituto Nazionale di Fisica Nucleare in Pisa (Italy), Prof. Cornelis van der Mee, Prof. Giuseppe Rodriguez and Prof. Sebastiano Seatzu from University of Cagliari (Italy), Prof. John E. Bowers from University of California Santa Barbara (USA), and Dr. Ming-Chun (Jason) Tien from Intel for the useful discussions and their suggestions.

This work has been partially supported by the Tuscany Region through the project ARNO-T3 PAR FAS 2007–2013 and by the European Commission grant IRIS, project no. 619194 FP7-ICT.

References and links

1. J. Jin, The Finite Element Method in Electro-magnetics (Wiley, 2002), 2nd ed.

2. A. B. Fallahkhair, K. S. Li, and T. E. Murphy, “Vector finite difference modesolver for anisotropic dielectric waveguides,” J. Lightwave Technol. 26, 1423–1431 (2008). [CrossRef]

3. P. Berini and K. Wu, “Modeling lossy anisotropic dielectric waveguides with the method of lines,” IEEE Trans. Microwave Theory Tech. 44, 749–759 (1996). [CrossRef]

4. A. S. Sudbø, “Film mode matching: a versatile numerical method for vector mode field calculations in dielectric waveguides,” Pure Appl. Opt. 2, 211–233 (1993). [CrossRef]

5. S. M. Sher, P. Pintus, F. Di Pasquale, M. Bianconi, G. B. Montanari, P. De Nicola, S. Sugliani, and G. Prati, “Design of 980nm-pumped waveguide laser for continuous wave operation in ion implanted Er:LiNbO₃,” IEEE J. Quantum Electron. 47, 526–533 (2011). [CrossRef]

6. Photon Design, “Integrated optics software FIMMWAVE 4.1,” http://www.photond.com/.

7. W. N. Ye, D.-X. Xu, S. Janz, P. Cheben, M.-J. Picard, B. Lamontagne, and N. G. Tarr, “Birefringence control using stress engineering in silicon-on-insulator (SOI) waveguides,” J. Lightwave Technol. 23, 1308–1318 (2005). [CrossRef]

8. S. Selleri and M. Zoboli, “Performance comparison of finite-element approaches for electromagnetic waveguides,” J. Opt. Soc. Am. A 14, 1460–1466 (1997). [CrossRef]

9. A. Konrad, “High-order triangular finite elements for electromagnetic waves in anisotropic media,” IEEE Trans. Microwave Theory Tech. 25, 353–360 (1977). [CrossRef]

10. K. Hayata, M. Koshiba, and M. Suzuki, “Vectorial wave analysis of stress-applied polarization-mantaining optical fiber by finite element method,” J. Lightwave Technol. 4, 133–139 (1986). [CrossRef]

11. K. Hayata, K. Miura, and M. Koshiba, “Finite-element formulation for lossy waveguides,” IEEE Trans. Microwave Theory Tech. 36, 268–276 (1988). [CrossRef]

12. K. Hayata, K. Miura, and M. Koshiba, “Full vectorial finite element formalism for lossy anisotropic waveguides,” IEEE Trans. Microwave Theory Tech. 37, 875–883 (1989). [CrossRef]

13. S. Selleri and M. Zoboli, “An improved finite element method formulation for the analysis of nonlinear anisotropic dielectric waveguides,” IEEE Trans. Microwave Theory Tech. 43, 887–892 (1995). [CrossRef]

14. S. Selleri, L. Vincetti, A. Cucinotta, and M. Zoboli, “Complex FEM modal solver of optical waveguides with PML boundary conditions,” Opt.Quant. Electron. 33, 359–371 (2001).

15. J.-F. Lee, D.-K. Sun, and Z. J. Cendes, “Full-wave analysis of dielectric waveguides using tangential vector finite elements,” IEEE Trans. Microwave Theory Tech. 39, 1262–1271 (1991). [CrossRef]

16. J.-F. Lee, “Finite element analysis of lossy dielectric waveguides,” IEEE Trans. Microwave Theory Tech. 42, 1025–1031 (1994). [CrossRef]

17. J. Tan and G. Pan, “A new edge element analysis of dispersive waveguiding structures,” IEEE Trans. Microwave Theory Tech. 43, 2600–2607 (1995). [CrossRef]

18. L. Valor and J. Zapata, “An efficient finite element formulation to analyze waveguides with lossy inhomogeneous bi-anisotropic materials,” IEEE Trans. Microwave Theory Tech. 44, 291–296 (1996). [CrossRef]

19. G. Pan and J. Tan, “General edge element approach to lossy and dispersive structures in anisotropic media,” IEE Proc.-Microw. Antennas Propag. 144, 81–90 (1997). [CrossRef]

20. G. Mur, “Edge elements, their advantages and their disadvantages,” IEEE T. Magn. 30, 3552–3557 (1994). [CrossRef]

21. B. M. A. Rahman and J. B. Davies, “Penalty function improvement of waveguides solution by finite elements,” IEEE Trans. Microwave Theory Tech. 32, 922–928 (1984). [CrossRef]

22. Y. Lu and F. A. Fernandez, “An efficient finite element solution of inhomogeneous anisotropic and lossy dielectric waveguides,” IEEE Trans. Microwave Theory Tech. 41, 1215–1223 (1993). [CrossRef]

23. L. D. Landau and E. M. Lifshits, “Electrodynamics of continuous media,” in “A Course of Theoretical Physics,” vol. 8 (Pergamon, 1960).

24. H. Dötsch, N. Bahlmann, O. Zhuromskyy, H. Hammer, L. Wilkens, R. Gerhardt, P. Hertel, and A. F. Popkov, “Applications of magneto-optical waveguides in integrated optics: review,” J. Opt. Soc. Am. B 22, 240–253 (2005). [CrossRef]

25. H. Yokoi, T. Mizumoto, Y. Shoji, N. Futakuchi, and Y. Nakano, “Demonstration of an optical isolator with a semiconductor guiding layer that was obtained by use of a nonreciprocal phase shift,” Appl. Opt. 39, 6158–6164 (2000). [CrossRef]

26. Y. Shoji, T. Mizumoto, H. Yokoi, I. W. Hsieh, and R. M. Osgood, “Magneto-optical isolator with silicon waveguides fabricated by direct bonding,” Appl. Phys. Lett. 92, 071117 (2008). [CrossRef]

27. W. Zaets and K. Ando, “Optical waveguide isolator based on nonreciprocal loss/gain of amplifier covered by ferromagnetic layer,” IEEE Photonic. Tech. L. 11, 1012–1014 (1999). [CrossRef]

28. G. J. Gabriel and M. E. Brodwin, “The solution of guided waves in inhomogeneous anisotropic media by perturbation and variational methods,” IEEE Trans. Microwave Theory Tech. 13, 364–370 (1965). [CrossRef]

29. W. Chew, Waves and Fields in Inhomogenous Media (Wiley-IEEE Press, 1999). [CrossRef]

30. C. H. Chen and C.-D. Lien, “The variational principle for non-self-adjoint electromagnetic problems,” IEEE Trans. Microwave Theory Tech. 28, 878–886 (1980). [CrossRef]

31. S. R. Cvetkovic and J. B. Davis, “Self-adjoint variational formulation for lossy anisotropic dielectric waveguide,” IEEE Trans. Microwave Theory Tech. 34, 129–134 (1986). [CrossRef]

32. R. Hoffman, “Comments on “self-adjoint variational formulation for lossy anisotropic dielectric waveguide”,” IEEE Trans. Microwave Theory Tech. 34, 1227–1228 (1986). [CrossRef]

33. P. Quéffélec, M. Le Floc’h, and P. Gelin, “New method for determining the permeability tensor of magnetized ferrites in a wide frequency range,” IEEE Trans. Microwave Theory Tech. 48, 1344–1351 (2000). [CrossRef]

34. F. Tisseur and K. Meerbergen, “The quadratic eigenvalue problem,” SIAM Rev. 43, 235–286 (2001). [CrossRef]

35. A. Nicolet and C. Geuzaine, “Waveguide propagation modes and quadratic eigenvalue problems,” in “6th International Conference on Computational Electromagnetics (CEM), 2006,” (Aachen, Germany, 2006), 1–8 (2006).

36. W. C. Chew and M. Nasir, “A variational analysis of anisotropic, inhomogeneous dielectric waveguides,” IEEE Trans. Microwave Theory Tech. 37, 661–668 (1989). [CrossRef]

37. K. Hayata, M. Koshiba, M. Eguchi, and M. Suzuki, “Vectorial finite-element method without any spurious solutions for dielectric waveguiding problems using transverse magnetic-field component,” IEEE Trans. Microwave Theory Tech. 34, 1120–1124 (1986). [CrossRef]

38. H.-Y. Fan, W.-W. Lin, and P. V. Dooren, “Normwise scaling of second order polynomial matrices,” SIAM J. Matrix Anal. A. 26, 252–256 (2005). [CrossRef]

39. N. J. Higham, D. S. Mackey, F. Tisseur, and S. D. Garvey, “Scaling, sensitivity and stability in the numerical solution of quadratic eigenvalue problems,” Int. J. Numer. Meth. Eng. 73, 344–360 (2008). [CrossRef]

40. Z. S. Sacks, D. M. Kingsland, R. Lee, and J.-F. Lee, “A perfectly matched anisotropic absorber for use as an absorbing boundary condition,” IEEE Trans. Antennas Propag. 43, 1460–1463 (1995). [CrossRef]

41. A. Cucinotta, G. Pelosi, S. Selleri, L. Vincetti, and M. Zoboli, “Perfectly matched anisotropic layers for optical waveguide analysis through the finite-element beam-propagation method,” Microw. Opt. Techn. Let. 23, 67–69 (1999). [CrossRef]

42. R. B. Lehoucq, K. Maschhoff, D. C. Sorensen, and C. Yang, “ARPACK library,” http://www.caam.rice.edu/software/ARPACK/.

43. R. B. Lehoucq, D. C. Sorensen, and C. Yang, ARPACK Users’ Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods, SIAM Publications, Philadelphia (1998). [CrossRef]

44. R. B. Lehoucq and D. C. Sorensen, “Deflation techniques for an implicitly re-started Arnoldi iteration,” SIAM J. Matrix Anal. A. 17, 789–821 (1996). [CrossRef]

45. D. C. Sorensen, “Implicit application of polynomial filters in a k-step Arnoldi method,” SIAM J. Matrix Anal. A. 13, 357–385 (1992). [CrossRef]

46. T. Mizumoto, R. Takei, and Y. Shoji, “Waveguide optical isolators for integrated optics,” IEEE J. Quantum Electron. 48, 252–260 (2012). [CrossRef]

47. R. Takei and T. Mizumoto, “Design and simulation of silicon waveguide optical circulator employing nonreciprocal phase shift,” Jpn. J. Appl. Phys. 49, 052203 (2010). [CrossRef]

48. M.-C. Tien, T. Mizumoto, P. Pintus, H. Kromer, and J. E. Bowers, “Silicon ring isolators with bonded nonreciprocal magneto-optic garnets,” Opt. Express 19, 11740–11745 (2011). [CrossRef] [PubMed]

49. P. Pintus, M.-C. Tien, and J. E. Bowers, “Design of magneto-optical ring isolator on SOI based on the finite element method,” IEEE Photonic. Tech. L. 23, 1670–1672 (2011). [CrossRef]

50. H. Yokoi, “Calculation of nonreciprocal phase shift in magneto-optic waveguides with Ce:YIG layer,” Opt. Mater. 31, 189–192 (2008). [CrossRef]

51. T. Shintaku, T. Uno, and M. Kobayashi, “Magneto-optic channel waveguides in Ce-substituted yttrium iron garnet,” J. Appl. Phys. 74, 4877–4881 (1993). [CrossRef]

52. M. Wallenhorst, M. Niemöller, H. Dötsch, P. Hertel, R. Gerhardt, and B. Gather, “Enhancement of the nonreciprocal magneto-optic effect of TM modes using iron garnet double layers with opposite Faraday rotation,” J. Appl. Phys. 77, 2902–2905 (1995). [CrossRef]

53. O. Zhuromskyy, H. Dötsch, M. Lohmeyer, L. Wilkens, and P. Hertel, “Magnetooptical waveguides with polarization-independent nonreciprocal phase-shift,” J. Lightwave Technol. 19, 214–221 (2001). [CrossRef]

54. S. Yamamoto and T. Makimoto, “Circuit theory for a class of anisotropic and gyrotropic thin-film optical waveguides and design of nonreciprocal devices for integrated optics,” J. Appl. Phys. 45, 882–888 (1974). [CrossRef]

55. O. Zhuromskyy, M. Lohmeyer, N. Bahlmann, Dötsch, P. Hertel, and A. F. Popkov, “Analysis of polarization independent Mach-Zehnder-type integrated optical isolator,” J. Lightwave Technol. 17, 1200–1205 (1999). [CrossRef]

56. H. Shimizu and Y. Nakano, “First demonstration of TE mode nonreciprocal propagation in an InGaAsP/InP active waveguide for an integratable optical isolator,” Jpn. J. Appl. Phys. 43, 1561–1563 (2004). [CrossRef]

57. H. Shimizu and Y. Nakano, “Fabrication and characterization of an InGaAsP/InP active waveguide optical isolator with 14.7dB/mm TE mode nonreciprocal attenuation,” J. Lightwave Technol. 24, 38–43 (2006). [CrossRef]

58. T. Amemiya, H. Shimizu, Y. Nakano, P. Hai, M. Yokoyama, and M. Tanaka, “Semiconductor waveguide optical isolator based on nonreciprocal loss induced by ferromagnetic MnAs,” Appl. Phys. Lett. 89, 021104 (2006). [CrossRef]

59. H. Shimizu and G. Syunsuke, “InGaAsP/InP evanescent mode waveguide optical isolators and their application to InGaAsP/InP/Si hybrid evanescent optical isolators,” Opt.Quant. Electron. 41, 653–660 (2009).

60. G. Kurczveil, P. Pintus, M. Heck, J. Peters, and J. Bowers, “Characterization of insertion loss and back reflection in passive hybrid silicon tapers,” IEEE Photonics Journal 5, 6600410 (2013). [CrossRef]

61. G. Krinchik and V. Artemjev, “Magneto-optic properties of nickel, iron, and cobalt,” J. Appl. Phys. 39, 1276–1278 (1968). [CrossRef]

62. C. Kittel, “Theory of the dispersion of magnetic permeability in ferromagnetic materials at microwave frequencies,” Phys. Rev. 70, 281–290 (1946). [CrossRef]

63. P. S. Pershan, “Magneto-optical effects,” J. Appl. Phys. 38, 1482–1490 (1967). [CrossRef]

64. M. C. Sekhar, M. R. Singh, S. Basu, and S. Pinnepalli, “Giant faraday rotation in BixCe3-xFe5O12 epitaxial garnet films,” Opt. Express 20, 9624–9639 (2012). [CrossRef]

65. P. Pintus, S. Faralli, and F. Di Pasquale, “Low threshold pump power and high integration in Al₂O₃:Er³⁺ slot waveguide laser on SOI,” IEEE Photonic. Tech. L. 22, 1428–1430 (2010). [CrossRef]

66. P. Pintus, S. Faralli, and F. Di Pasquale, “Integrated 2.8μm laser source in Al₂O₃:Er³⁺ slot waveguide on SOI,” J. Lightwave Technol. 29, 1206–1212 (2011). [CrossRef]

67. P. Pintus, F. Di Pasquale, and J. E. Bowers, “Design of TE ring isolators for ultra low loss Si₃N₄ waveguides based on the finite element method,” Opt. Lett. 36, 4599–4601 (2011). [CrossRef] [PubMed]

68. P. Pintus, F. Di Pasquale, and J. E. Bowers, “Integrated TE and TM optical circulators on ultra-low-loss silicon nitride platform,” Opt. Express 21, 5041–5052 (2013). [CrossRef] [PubMed]

69. D. Jalas, A. Petrov, M. Krause, J. Hampe, and M. Eich, “Resonance splitting in gyrotropic ring resonators,” Opt. Lett. 35, 3438–3440 (2010). [CrossRef] [PubMed]

70. V. Subramaniam, G. N. De Brabander, D. H. Naghski, and J. T. Boyd, “Measurement of mode field profiles and bending and transition losses in curved optical channel waveguides,” J. Lightwave Technol. 15, 990–997 (1997). [CrossRef]

Fe layer thickness (nm)	NRL Δα (dB/mm)	NRL Δα_pm (dB/mm)	NRPS Δβ (rad/mm)	NRPS Δβ_pm (rad/mm)

25	−6.9643	−6.9550	0.4092	0.4159
50	−9.8495	−9.7799	1.1817	1.1903
75	−10.9298	−10.8404	1.4907	1.4976
100	−11.8120	−11.7101	1.6454	1.6531

Accurate vectorial finite element mode solver for magneto-optic and anisotropic waveguides

Abstract

1. Introduction

2. Mathematical model

2.1. Differential electromagnetic problem

2.2. Variational formulation: non-self-adjoint problem

2.3. Waveguide mode problem

2.4. Interpolating function

2.5. Discretized functional

2.6. Eigenvalue problems

2.7. Lossy case: fully complex formulation

2.8. Lossless case: real formulation

3. Numerical results

3.1. Nonreciprocal phase shift in lossless magneto-optical waveguides

3.2. Finite element method vs. perturbation method

3.3. Nonreciprocal loss in magneto-optical waveguides

4. Conclusions

Acknowledgments

References and links

Cited By

Figures (10)

Tables (1)

Equations (56)

Optics Express