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

In this work, a dielectric waveguide mode solver is presented considering a general nonreciprocal permittivity tensor. The proposed method allows us to investigate important cases of practical interest in the field of integrated optics, such as magneto-optical isolators and anisotropic waveguides. Unlike the earlier developed mode solver, our approach allows for the precise computation of both forward and backward propagating modes in the nonreciprocal case, ensuring high accuracy and computational efficiency. As a result, the nonreciprocal loss/phase shift can be directly computed, avoiding the use of the perturbation method. To compute the electromagnetic modes, the Rayleigh-Ritz functional is derived for the non-self adjoint case, it is discretized using the node-based finite element method and the penalty function is added to remove the spurious solutions. The resulting quadratic eigenvalue problem is linearized and solved in terms of the propagation constant for a given frequency (i.e., γ−formulation). The main benefits of this formulation are that it avoids the time-consuming iterations and preserves the matrix sparsity. Finally, the method is used to study two examples of integrated optical isolators based on nonreciprocal phase shift and nonreciprocal loss effect, respectively. The developed method is then compared with the perturbation approach and its simplified formulation based on semivectorial approximation. © 2014 Optical Society of America OCIS codes: (130.0130) Integrated optics; (000.3860) Mathematical methods in physics; (000.4430) Numerical approximation and analysis; (160.3820) Magneto-optical materials; (230.7370) Waveguides; (230.3240) Isolators. 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). 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). 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). 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:LiNbO3,” IEEE J. Quantum Electron. 47, 526–533 (2011). 6. Photon Design, “Integrated optics software FIMMWAVE 4.1,” http://www.photond.com/. #205244 $15.00 USD Received 22 Jan 2014; revised 31 Mar 2014; accepted 28 Apr 2014; published 20 Jun 2014 (C) 2014 OSA 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015737 | OPTICS EXPRESS 15737 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). 8. S. Selleri and M. Zoboli, “Performance comparison of finite-element approaches for electromagnetic waveguides,” J. Opt. Soc. Am. A 14, 1460–1466 (1997). 9. A. Konrad, “High-order triangular finite elements for electromagnetic waves in anisotropic media,” IEEE Trans. Microwave Theory Tech. 25, 353–360 (1977). 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). 11. K. Hayata, K. Miura, and M. Koshiba, “Finite-element formulation for lossy waveguides,” IEEE Trans. Microwave Theory Tech. 36, 268–276 (1988). 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). 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). 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). 16. J.-F. Lee, “Finite element analysis of lossy dielectric waveguides,” IEEE Trans. Microwave Theory Tech. 42, 1025–1031 (1994). 17. J. Tan and G. Pan, “A new edge element analysis of dispersive waveguiding structures,” IEEE Trans. Microwave Theory Tech. 43, 2600–2607 (1995). 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). 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). 20. G. Mur, “Edge elements, their advantages and their disadvantages,” IEEE T. Magn. 30, 3552–3557 (1994). 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). 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). 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). 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). 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–071117–3 (2008). 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). 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). 29. W. Chew, Waves and Fields in Inhomogenous Media (Wiley-IEEE Press, 1999). 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). 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). 32. R. Hoffman, “Comments on “self-adjoint variational formulation for lossy anisotropic dielectric waveguide”,” IEEE Trans. Microwave Theory Tech. 34, 1227–1228 (1986). 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). 34. F. Tisseur and K. Meerbergen, “The quadratic eigenvalue problem,” SIAM Rev. 43, 235–286 (2001). 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). 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). #205244 $15.00 USD Received 22 Jan 2014; revised 31 Mar 2014; accepted 28 Apr 2014; published 20 Jun 2014 (C) 2014 OSA 30 June 2014 | Vol. 22, No. 13 | DOI:10.1364/OE.22.015737 | OPTICS EXPRESS 15738 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). 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). 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). 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). 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). 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). 45. D. C. Sorensen, “Implicit application of polynomial filters in a k-step Arnoldi method,” SIAM J. Matrix Anal. A. 13, 357–385 (1992). 46. T. Mizumoto, R. Takei, and Y. Shoji, “Waveguide optical isolators for integrated optics,” IEEE J. Quantum Electron. 48, 252–260 (2012). 47. R. Takei and T. Mizumoto, “Design a


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][9][10][11][12][13][14] and edge elements (also called Nédélec elements) [1,8,[15][16][17][18][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 crosssection is in the xy-plane and the light propagates only along the z-direction.
Using Cartesian coordinate system, we have considered the following relative permittivity tensor 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 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 timeconsuming 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].

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 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 where J and ρ are the electric current and charge density, respectively, the constants ε 0 and μ 0 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 B(E, H).
By computing E(r) from (4d) and replacing it into (4c), the result is a differential problem for H(r) where V is the domain and S is the boundary.In Eq. ( 5), we used the relation c = 1/ √ ε 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.

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 = ε T r ) nor hermitian (i.e., ε r = ε † r ).As a result, the differential operator [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 where H a is the adjoint magnetic fields, J a is the adjoint source distribution, and B a (H a ) = 0 are the adjoint boundary conditions.Note that J a and B a (H a ) = 0 are chosen according to the direct problem, while the relative tensors ε a r and μ a r are equal to ε T r and μ T r when we consider the real inner product, or ε † r and μ † r in case we use the complex inner product [30].
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 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 ε T r and μ T r in the previous formula [30].A similar formulation can be derived for E and the electric adjoint field E a .

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 where γ = jβ + α is the propagation constant, β the phase constant and α the attenuation.
Similarly, the adjoint fields are 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][31][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 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 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 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).

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 where m is the number of nodes in the single element, N e i (x, y) is the basis function for the node i = 1, . . ., m, and (h e xi , h e yi , h e zi ) are the values of H e n at the interpolation node i.The approximate field H n in the whole cross-section is then 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 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 ) 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.

Discretized functional
Using the approximation of the previous paragraph and introducing the tensor p = ε −1 0 ε −1 r , 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 . . .
and the matrices R e 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 e i (x, y) and N e j (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 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 γ 2 , the sum of the matrices multiplied by γ, and the sum of the matrices multiplied by ω 2 , 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 For this purpose, let us extend the previous notation defining three vectors with unknown coefficients where N is the number of the mesh nodes.Therefore, the discrete functional computed in the cross-section is where M, C, K and T are 3 N × 3 N 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.

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 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 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 we obtain 0 As a result, both the ω-formulation (Eq.( 28)) and the γ-formulation (Eq.( 30)) have been transformed into a generalized eigenvalue problem 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 σ 0 , we subtract σ 0 Bv from either term in Eq. ( 31) Collecting the common terms, we obtain The generalized eigenvalue problems (31) and ( 33) have the same eigenvectors, but shifted eigenvalues 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 where n min and n max are the minimum and maximum refractive indices in the problem under investigation, and k 0 is the wavenumber in the vacuum (k 0 = 2π/λ ).In this case we shift the spectrum close to σ 0 = n max k 0 when we want to compute the forward propagating modes, while we move the spectrum close to σ 0 = −n max k 0 when we are looking for the backward propagating modes.From Eq. ( 35) it is useful to define the effective index as n e f f = β /k 0 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 The eigenvalues and the matrices of the scaled eigenvalue problem are

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

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 where ε i j 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 makes the problem (17) fully real [1,9] and this halves the memory demand.

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

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 2 = 1.46, respectively.In order to achieve NRPS, the top-layer of the waveguide is bonded with a ceriumsubstituted yttrium iron garnet (Ce:YIG) (n Ce:Y IG = 2.22, [50]) grown on a (Ca,Mg,Zr)substituted gadolinium gallium garnet (SGGG), (n SGGG = 1.97).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 where ε xx = ε yy = ε zz = n 2 Ce:Y IG , and the off-diagonal element ε yz is responsible for the MO For calculation, we assume θ F = 4000 o /cm at λ = 1550 nm [48][49][50], and then ε yz = 7.65 • 10 −3 .However, in high quality Ce:YIG crystal, a Faraday-rotation coefficient as large as 4500 o /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 600 nm × 215 nm silicon waveguide cross-section, with a 380 nm thick Ce:YIG layer above it.Those values maximize the NRPS for the TM mode, as it is explained in the following.Because the value of the off-diagonal terms in the permittivity tensor are usually small (e.g., ε yz /n 2 Ce:Y IG 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 600 nm wide silicon waveguide in order to guarantee the single mode regime (w Si in Fig. 2).
Note that NRPS is maximized when the maximum of |H x | 2 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: [500 nm, 3000 nm]. 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 e f f < 1.97 are below cut-off and cannot propagate.

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][53][54][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 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∂ y E y /β for the TM mode and E y = 0 and E z = j∂ x E x /β for the TE mode) [47,53,55].Considering our example (TM mode), Eq. ( 44) can be rewritten as follows 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.Finally, in Fig. 7 we compare the three methods considering h MO = 500nm and varying the silicon thickness only.

Nonreciprocal loss in magneto-optical waveguides
Optical isolation based on nonreciprocal loss can be implemented by considering a magnetooptic 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][57][58] and silicon photonics technologies [59].In the plot we compared the three methods.
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 400 nm deep trenches into a 700 nm 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 2 ) 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.In our simulations, we assume a 30 nm-thick layer of titania and we varied the thickness of the iron between 25 nm and 100 nm.The refractive indices of the TiO 2 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 where ε xx = ε yy = ε zz = n 2 Fe , and the off-diagonal element ε xz = 3.15 + 1.8i when the magnetization of the Fe layer is saturated [57,61] 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 μ 0 [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.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   that in the lossy case, the formula (44) must be replaced by the more general one where E a and H a are the adjoint fields [28].The results of this comparison are reported in table 1.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.

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

Fig. 1 .
Fig. 1.Cross section of a uniform waveguide loaded with inhomogeneity.
TM-mode at 1550 nm.

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

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

Fig. 9 .
Fig. 9. Amplitude of the magnetic field H for the TE-mode at 1550 nm.The amplitude of the field components have been normalized in term of power and the profiles are shown in H/μm.
(b).It is worth noting that the longer the waveguide, the larger the nonreciprocal effect.As it can be seen from the figure, a 25 nm-thick iron film will produce a NRL of about −7 dB/mm, while a four times thicker layer generates an isolation of almost −12 dB/mm.Those values are consistent with those reported in the literature[56,57,59].
(a) Propagation loss of forward and backward waves.(b) NRL and NRPS of the TE mode.
− p zx J eT p yx J e − p yy N e p xz J e − p zy N eT p zx N eT − p xz N e p xy N e − p xx J e p xy J eT − p yy N eT p yx N eT − p xx J eT 0 . It is worth mentioning that, although iron is #205244 -$15.00USD Received 22 Jan 2014; revised 31 Mar 2014; accepted 28 Apr 2014; published 20 Jun 2014

Table 1 .
Comparison between the two methods.