Inverse Effective Index Method for Two-Dimensional Simulations of Photonic Components

The representation of three-dimensional integrated photonic structures in two-dimensional designs has been investigated through finite element analysis using COMSOL Multiphysics. The accuracy of the effective index method for various waveguide designs is studied. An alternative to the effective index method, here named the inverse effective index method, is derived in this paper. The results indicate that the commonly used effective index method has a deviation in propagation constant up to 17%, where the error is largest for low-confinement factor waveguides. The inverse effective index method vastly reduces this error to $< $0.01% for low confinement-factor geometries. This enables an accurate simulation in two-dimensions while retaining the waveguiding properties of the third dimension. As with the effective index method, the inverse effective index method does not require proprietary software.

Inverse Effective Index Method for Two-Dimensional Simulations of Photonic Components

Jens Høvik and Astrid Aksnes
Abstract-The representation of three-dimensional integrated photonic structures in two-dimensional designs has been investigated through finite element analysis using COMSOL Multiphysics.The accuracy of the effective index method for various waveguide designs is studied.An alternative to the effective index method, here named the inverse effective index method, is derived in this paper.The results indicate that the commonly used effective index method has a deviation in propagation constant up to 17%, where the error is largest for low-confinement factor waveguides.The inverse effective index method vastly reduces this error to <0.01% for low confinement-factor geometries.This enables an accurate simulation in two-dimensions while retaining the waveguiding properties of the third dimension.As with the effective index method, the inverse effective index method does not require proprietary software.
Index Terms-Effective medium theory, electromagnetic optics, photonic integrated circuits, propagation methods.

I. INTRODUCTION AND MOTIVATION
W HILE mass production of standardized microphotonic components is becoming more economical, development and testing of new designs is still time-consuming and costly.Iterative approaches with respect to fabrication of newer and better designs are therefore not ideal, and accurate characterization of components before fabrication is important to minimize both cost and time consumption.Historically a great deal of effort went into finding analytical solutions to Maxwell's equations for integrated photonics design.As geometries become more complex the derivation of closed form solutions of Maxwell's equation were not always possible.The inability to analytically solve Maxwell's equations led to the development of computational electrodynamics.
Although computational electrodynamics has become a large field over the last decades, its main methods of solution can be categorized into two distinct categories.Frequency domain methods, which most notably include the Finite Element Method (FEM) [1] and the Finite-Difference Frequency Domain (FDFD) method [2], and time domain simulations, which most notably include the finite-difference Time-Domain (FDTD) [3], and the Time-Domain Finite-Element-method (TD-FEM).Frequency-domain methods generally rely on matrix inversions, requiring large matrix solvers which put a high demand on available computational power and memory.The time domain solutions iterate solutions to Maxwell's equations in time and thus are quite computationally intensive.
In the early days of computational electrodynamics, rigorous numerical solutions to Maxwell's equations were time consuming due to the low computational power of the available hardware.This led to significant efforts towards reducing the computational cost of such numerical solutions by performing approximations by sacrificing accuracy in areas of the simulation which were not considered crucial for the design component.One such simplification is to reduce the number of computations required per unit cell, for example by assuming the electromagnetic wave only has polarization along one unit vector.An alternative is to reduce the dimensionality of the simulation space, as dimensionality scales the number of computational points with N D where D is the dimension and N is the number of grid points required along one length of the simulation space.Reducing the number of dimensions reduces both the memory and time requirements of an electromagnetic simulation.This is the main reason the effective index method (EIM) [4] is utilized so often.The effective index method was originally a simple way of estimating the propagation constants of strip-wire waveguides by solving the boundary equations of a "side view" of the component.The method was first proposed by R. M. Knox and P. P. Toulios in 1970 [4].There has since been attempts at improving the method [5] and developing new methods to rapidly find the propagation constants of waveguides such as the weighted index method [6], perturbation methods [7] and variational methods [8].As the EIM was originally intended for strip-wire waveguides, alternate methods for finding the propagation constants of other waveguide geometries such as ridge waveguides [9] were also developed, such as the spectral index method [10], [11], [12], [13].Most of these solutions were developed to circumvent the need for numerical calculations when determining the waveguiding properties of a component due to limited hardware at the time.
In this paper we consider how the effective index method is utilized today -as a convenient tool to reduce the dimensionality of a numerical simulation [14], [15].Although simple and elegant, we show how the method is not particularly accurate for a large range of waveguide geometries.Additionally, we propose our own variation of the effective index method which significantly increases the accuracy of the approximation method to the waveguide geometries where the conventional EIM fails to achieve satisfactory accuracy.The alternative method, here named the inverse effective index method (inverse EIM) is derived in this paper.

A. The Effective Index Method
A three-dimensional (3D) study can be decomposed into three separate simulations by removing one axis for each simulation, as indicated in Fig. 1.All three perspectives hold important information, and eliminating one axis will deviate the simulation from the physical case.Most 2D-simulations utilize the birds-eye view (Fig. 1(d)) for the component study, as it gives the best representation of the interaction between different photonic components on a single chip.
The birds-eye view leads to relatively inaccurate simulations however, as the material interactions above and below the waveguide are not accounted for.This causes the propagation constant of the guided mode to deviate significantly from the 3-dimensional scenario.The effective index method attempts to remedy this by first considering a side-view of the component (Fig. 1(b)), where the area above and below the waveguide is included.In the EIM, the guided mode is viewed from the side, and the effective index is calculated and used as the refractive index of the waveguide core when performing simulations from the birds-eye view.The material parameters of the birds-eye view simulation (Fig. 1(d)) are altered to better approximate the three-dimensional structure.When viewed from the side, the component is reduced to a planar waveguide, where the guided mode can be found by utilizing Maxwell's equations.By applying the boundary conditions the following pair of characteristic equations are obtained [16]: where h is the height of the planar waveguide, κ core is the transverse propagation constant in the core, γ sub/cl the propagation constant in the substrate and cladding, and n core/cl/sub the refractive index in the core, cladding and substrate, respectively.
) applies to out-of-plane polarization whereas (2) applies to in-plane polarization.It is important to note that the EIM only utilizes the side-view to estimate the propagation constant of the guided mode, and ( 1) is therefore conventionally defined as the TE-mode solution while (2) is defined as the TM-mode solution.The inverse effective index method utilizes the topdown view (Fig. 1(d)), thus the relative polarization direction changes.Therefore, to ensure that (1) and ( 2) are generic, the terms in-plane and out-of-plane are used rather than TE/TM.
These transcendental equations can be plotted to find the solution for the effective index of the propagating mode of the planar waveguide.This effective index is then used as the refractive index of the guiding material when simulating the waveguide from a birds-eye view (Fig. 1(d)).Altering the refractive index of the guiding material in such a way allows for the 2-dimensional simulation to somewhat compensate that one dimension is removed.The equations are solved by setting all material and geometric parameters as constant, with the effective index n ef f as the only unknown.To demonstrate the EIM, an example component is studied.A rectangular core waveguide Fig. 2. EIM is used to estimate the TE propagation constant of a slab waveguide using (1).The value of the effective index is given where the two curves intersect.
with fundamental TE mode at 1.55 μm and material values for a generic SOI component are studied, e.g.n core = 3.47, n sub = 1.44 and n cl = 1.The EIM is used to estimate the TE propagation constant and effective index.Equation 1 is used to plot the graphs for a 500 nm x 220 nm rectangular slab waveguide shown in Fig. 2. The point where the two graphs intersect is the effective index of the slab mode.This effective index is then used as the refractive index of the core in the real component study (Fig. 1(d)).In this graph the two lines intersect at n ef f = 2.825.The EIM is then essentially just an alteration of the waveguide material, from conventional silicon (n = 3.47) to a two-dimensional version of silicon appropriate for these specific waveguide geometries (n = 2.825).This greatly increases the accuracy of a two-dimensional study from the birds-eye view, but as shall be seen, is still very inaccurate, in particular for waveguide geometries with lower confinement factors or for the TM-polarization.

B. Inverse Effective Index Method
State-of-the-art simulation methods are still being improved and methods like the beam propagation method, beam envelope method [17], 2.5D FDTD simulation methods [18] and the eigenmode expansion method [19] offer ways of achieving higher accuracies or faster simulations of 3D structures.These methods have in common that they require proprietary software.The EIM is still widely used for simulating photonic structures [20], [21], [22] because it can be solved using simple equations and graph plotting software.In addition, the EIM is compatible with most simulation approaches such as FDTD or FEM, meaning that it can be utilized regardless of the simulation software of the individual user, unlike many other approaches that improve the simulation accuracy.
A fully-vectorial 3-dimensional simulation of a photonic component can take hours to days depending on the computational power available as well as the complexity and size of the geometry under study.However, the supported modal fields and propagation constants of a 3-dimensional waveguide are inherently two-dimensional problems.These can be solved by considering the cross-section of the waveguide under study Fig. 3. Inverse EIM is utilized to find the refractive index of the core material using (2).The refractive index is given where the two curves intersect.
(Fig. 1(c)) and utilizing e.g. the finite-difference method.In the finite-difference method, the wave equation is discretized which yields an eigenvalue matrix across the simulation grid.This matrix can be solved using linear algebra algorithms which in turn gives the supported mode profiles and propagation constants of the cross-section of the waveguide [23], [24].On modern hardware such a computation only takes a few seconds.
The inverse effective index method presented here utilizes that modal fields and propagation constants are inherently 2D problems.As explained above, the eigenvalue matrix is solved to calculate β, the propagation constant of the 3D-waveguide.β is therefore a known in (1) and (2) (as opposed to an unknown for the EIM) and the required index of the core (n core ) is the unknown.In summary, (1) and ( 2) are used to solve for the refractive index which must be given to the core to ensure that the 2D-simulation best approximates the 3D study.It is imperative that the correct equations are utilized, as unlike the EIM, which solves for the side-view of the waveguide (Fig. 1(b)), the inverse EIM must solve for the birds-eye view (Fig. 1(d)).With a change in perspective, the relative polarization direction also changes.Equation (2) must therefore be utilized for TE-modes and (1) for TM-modes.Assuming a 500 nm × 220 nm waveguide propagating the fundamental TE-mode, the required refractive index of the core material can be found by plotting (2) with n core as the unknown and utilizing β as found by the finite-difference method.The resulting graph can be seen in Fig. 3.Although the left-hand side of (1) and ( 2) are the same for both equations, since κ core = n 2 core k 2 0 − n 2 ef f k 2 0 , plotting as a function of n core (first term) as opposed to n ef f (second term) inverts the graph and alters its shape relative to the value on the x-axis.In Fig. 3 the two graphs intersect at n core = 2.763, as opposed to n ef f = 2.825 for the EIM.Although this comparison gives useful insight at how these two methods differ, these are simply material parameters to be used for subsequent simulations.The difference in material parameters must not be confused with the alteration of propagation constant in the subsequent simulation when utilizing these material parameters.The result of interest is how these material parameters affect the propagation constant Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
of the guided mode in the two-dimensional simulation relative to the three-dimensional study.This is examined in more detail in Section III.

III. METHOD VERIFICATION
The accuracy of 2D simulations utilizing either no approximation (material parameters identical to the 3-dimensional case), the EIM or the inverse EIM is analyzed through two sets of studies.The first study simulates all geometries that belong to the single-mode domain of silicon-on-insulator strip-wire waveguides.For each simulation the propagation constant of the guided mode in two dimensions is compared to the propagation constant in the 3D case.This study is intended to probe the simulation accuracy for a wide array of waveguide geometries.Each simulation in this study utilizes the standard telecommunications wavelength of 1.55 μm with material parameters n core = 3.47, n sub = 1.44 and n cl = 1.The second study simulates ring resonators (5 μm radius) to investigate the accuracy of the approximation methods when representing a common photonic component of certain complexity.All simulations are performed using the FEM software COMSOL Multiphysics and are conducted on a 3.6 GHz quad-core computer with 32 GB of RAM.For each simulation the waveguide geometry is fed into a script which either calculates the refractive index of the core using the EIM or the inverse EIM, and feeds this data to COMSOL Multiphysics which performs the finite-element analysis.Except for the calculations done to approximate the refractive index of the core, the simulation run in COMSOL Multiphysics is identical for each approximation method.

A. Propagation Constant Accuracy
The entire single-mode domain of silicon-on-insulator (SOI) strip-wire waveguides was simulated iterating the width and height of the core by 10 nm increments for a total of 3496 individual simulations.The propagation constants of the fundamental TE-and TM-modes of the 3D waveguides were then compared to the propagation constants found in 2D-simulations.A matrix plot of the results are shown in Fig. 4 for TE-modes and Fig. 5 for TM modes.

B. Ring Resonator Simulations
A silicon ring resonator of radius 5 μm was modelled and simulated in a 3D finite element study.The wavelength was swept from 1.56 μm to 1.635 μm and the transmission measured.The sweeping range was chosen to encompass several resonance orders while being within the range of a typical tunable laser source around the communications wavelength at 1.55 μm.An identical resonator geometry was then simulated in a 2D domain using no approximation, the EIM and the inverse EIM.The propagating mode was for all simulations the fundamental TE-mode.The results are presented in Fig. 6.For ease of comparison, the resonance orders of the various resonances are included in the plot (labelled Mxx where xx is the order number).V. DISCUSSION Simulations comparing the propagation constants of twoand three-dimensional waveguides confirm that for low confinement-factor modes the error can be as high as 100% when not utilizing any approximation method.This error is reduced to <17% when utilizing the effective index method.Using the inverse effective index method reduces this margin of error to <0.01%.The error present in this case is most likely numerical noise and uncertainty in finding the exact interpolation of the two graphs as shown in Fig. 3.The same statistical noise is likely to be present during the effective index approach, but not visible when the error margin is already large.It is apparent from Figs. 4 and 5 that the errors are largest for waveguides with a small cross section where the confinement factor is low.This is evident as a low confinement factor waveguide is less accurately depicted when one dimension is disregarded.As the waveguide geometry nears the multimode regime and most of the field is situated inside the core, the error decreases as the confinement factor decreases.For the inverse effective index method, the confinement factor and specific waveguide geometry have no impact on the simulation accuracy as all dimensions are taken into consideration when performing the inverse effective index method.
Simulating ring resonators using no approximation yielded an error of almost 350 nm (or 19 resonance orders).Due to the high refractive index of the silicon it also yields spurious higher-order modes in the waveguide, even though such modes should not exist for the specified geometries.For the effective index method the resonances are erroneous by three resonance modes or a total of 50 nm, while for the inverse effective index approach the resonances of the three-and two-dimensional studies overlap to within 50 pm.On a 32 GB intel quad-core computer running at 3.6 GHz, the 2D simulations took approximately 2.5 hours whereas the 3D-simulations took 30 hours.
Ring resonators are in essence one-dimensional problems where the response can be accurately predicted by the effective index of the guided mode alone.Therefore, it is not surprising that a 2D study where the guided modes match with the 3D case is accurate to within 50 pm.Nonetheless, the behaviour of more complicated structures that also mainly rely on the propagation constants such as Vernier enhanced Mach-Zehnder interferometers, long period grating couplers or ring resonator arrays for frequency filtering can be more challenging to predict analytically.For such components, the accuracy of the simulation is paramount.The ring resonator was chosen as an example due to its compact geometry making it easy to simulate in 3D allowing for the comparison of the two-and three-dimensional simulations.

VI. CONCLUSION
The inverse effective index method has been developed to perform 2D simulations that accurately mimic 3D simulations when considering the effective index of the guided mode.Compared to the effective index method our method is a large improvement for low confinement-factor waveguides as the error decreases to less than 0.01 percent from as high as 17 percent for the EIM and over 100 percent when utilizing no approximation method.The inverse EIM was also applied to the simulation of a ring resonator, and was found to closely mimic the 3D result.For simulations, the two-and three-dimensional resonance modes overlap to within 50 pm.One main advantage of both EIM and inverse-EIM is their compatibility with any simulation software, making them highly flexible as approximation methods.Utilizing this method requires a mode solver, which only takes seconds and is freely available online [25].

Fig. 1 .
Fig. 1.(a) The component under study illustrated in three dimensions.The bottom three figures are 2D representations of the same component, with (b) side-view, x-axis removed.(c) Front view (waveguide cross-section), y-axis removed and (d) top-view, z-axis removed.

Fig. 4 .
Fig. 4. Deviation in the propagation constants of the fundamental TE modes of simulated 2D waveguides with varying width and height compared to the 3D case for (a) no approximation method used, (b) effective index method and (c) the inverse effective index method.

Fig. 5 .
Fig. 5. Deviation in the propagation constants of the fundamental TM modes of simulated 2D waveguides with varying width and height compared to the 3D case for (a) no approximation method used, (b) effective index method and (c) the inverse effective index method.

Fig. 6 .
Fig. 6.Simulations of a 5 µm radius ring resonator propagating the fundamental TE-mode using (a) no approximation, (b) EIM and (c) inverse EIM.The resonance order is specified below each resonance.