Applicability of the Effective Index Method for the Simulation of X-Cut LiNbO₃ Waveguides

Abstract

Photonic integrated circuits (PIC) find applications in the fields of microwaves, telecoms and sensing. Generally, PICs are fabricated on a base of isotropic materials such as SOI, Si₃N₄, etc. However, for some applications, anisotropic substrates such as LiNbO₃ are used. A thin film of LiNbO₃ on an insulator (LNOI) is a promising material platform for complex high-speed PICs. The design and simulation of PICs on anisotropic materials should be performed using rigorous numerical methods based on Maxwell’s equations. These methods are characterized by long calculation times for one simulation iteration. Since a large number of simulation iterations are performed during the PIC design, simulation methods based on approximations should be used. The effective index method (EIM) is an approximation-based method and is widely applied for simulations of isotropic waveguides. In this study, the applicability of EIM for simulations of anisotropic waveguides is analyzed. The results obtained by EIM are compared with the calculation results of a rigorous finite-difference frequency-domain (FDFD) method for evaluation of the EIM’s applicability limits. In addition, radiation losses in waveguides with rough sidewalls are estimated using the Payne–Lacey model and EIM. The results demonstrate the applicability of EIM for the simulation of anisotropic LNOI-based waveguides with cross-section parameters specified in this paper.

1. Introduction

The scope of photonic integrated circuits (PICs) covers optical metrology, sensors, optical signal processing, telecommunications, etc. Currently there is a wide range of materials being used for production of PICs. Each material has properties that make it suitable for a specific class of problems that can be solved using PICs. For example, Si₃N₄ is used for quantum cryptography and LIDAR applications. Low-loss telecom PICs are fabricated on various glass-based material platforms [1,2,3]. A well-known material for PIC fabrication is lithium niobate (LiNbO₃, LN). A notable property of LN is its high electro-optical coefficient, which determined its application as a material platform for the manufacture of modulators, switches and other devices with high-frequency modulation rates (up to 100 GHz) [4,5,6,7,8]. The conventional technologies for waveguides formation on LN substrate are proton exchange and titanium diffusion [9,10]. On the one hand, waveguides formed with conventional technologies are characterized with low optical losses for coupling the signal with fiber and low propagation losses. On the other hand, these optical waveguides are characterized by low refractive index contrast. The refractive index contrast (the difference between the refractive indexes of the waveguide core and those of the substrate or cladding) determines the minimum acceptable curvature radius of the waveguide. Waveguides formed with proton exchange and titanium diffusion have a typical curvature radius of about 2–3 cm, which makes it impossible to integrate dozens of basic elements of PICs on the single chip. Around the early 2000s, a technology of producing a thin film of lithium niobate on an insulator (LNOI) was developed (schematic structure of X-cut of LNOI is shown on Figure 1) [11]. It was shown that high-index-contrast waveguides with a step-index type profile could be fabricated on this material platform. These high-index-contrast waveguides based on LNOI can be used for PICs fabrication with many integrated elements with high-frequency modulation rates. Due to this, great attention has been paid to the application of LNOI in integrated photonics; this is confirmed by the large number of published papers on the subject [12,13,14,15,16,17,18,19].

The basic process of PICs development includes the simulation of single components, as well as of entire optical circuits [20]. The simulation of arbitrary waveguide structures can be performed using a group of methods based on finite differences and finite elements (FDTD, FDFD, FEMTD) [21,22,23], which are widely presented in the commercial software. Moreover, there are other high-accuracy methods for simulation of anisotropic photonic integrated elements, such as the spectral element method [24,25]. The above methods directly solve Maxwell’s equations without any approximations and can therefore provide the most accurate result. The main disadvantages of these methods are their long calculation times and their programming complexity [26,27]. During the simulation process of basic integrated optical elements, hundreds of calculation iterations are performed in order to determine the optimum parameters. Thus, the time of a single calculation iteration becomes critical, and finding ways to speed up and simplify the simulation process is an urgent task.

The purpose of this study is to evaluate the possibility of applying the effective index method (EIM) for designing an anisotropic high-contrast optical waveguide based on LNOI. EIM is a simple semi-analytical method that allows fast simulation of waveguides with different cross-sections; it could be used to simplify 3D calculations of radiation propagation in integrated optical elements (to perform so-called 2.5D simulation). In addition, EIM might be used for rapid circuit simulation with the S-matrix method [28].

This paper considers the simulation of the optical anisotropic waveguide on an X-cut LNOI substrate, directed along the Y-crystallographic axis. As far as the authors know, EIM has not been used before for anisotropic waveguides analysis on LNOI substrate. As a result of the simulation, the range of etching depths and slope angles of the sidewalls were determined, for which the difference between the results of calculation by the finite-difference method and by EIM does not exceed 1%. This makes it possible to use EIM to speed up and simplify the simulation process of optical waveguides.

Since EIM is applicable for the LNOI waveguide analysis, optical radiation losses could be evaluated by the Payne–Lacey model. The model uses the root-mean-square (RMS) roughness and the correlation length for the radiation losses calculation. In the paper, RMS is determined from the line edge roughness (LER) obtained from experiments performed earlier in [29].

2. Materials and Methods

2.1. Calculations of Optical Waveguides with EIM

Optical waveguides based on the X-cut of LNOI are anisotropic; consequently it is necessary to solve the vector wave equation for analysis of radiation propagation along these waveguides. Additional difficulties appear during the modes simulation if the waveguide axis does not coincide with one of the crystallographic axes of LN, since the permittivity tensor included in the wave equation will not be diagonal [12,13,14,15,16,17,18,19].

The problem of solving the vector wave equation for anisotropic waveguides might be simplified by considering the optical waveguide directed along one of the crystallographic axes, when the permittivity tensor becomes diagonal. In this case, if we assume that the optical modes in the waveguide are polarized strictly along the transverse axes, then the electric displacement vector will be parallel to the electric field vector. Then the vector wave equation should be [30]:

$${\nabla }^{2}E+k^2{n}^2·E=\nabla \left(\nabla ·E\right)$$

(1)

where E is the electric field vector, k is the wave vector in vacuum and n is the diagonal refractive index tensor. Equation (1) might be represented as a simplified vector wave equation:

$$\left({\nabla }_t^2+k^2{n}_t^2-{\beta }^2\right)e_t=-{\nabla }_t\left(e_t·{\nabla }_t\ln n_t^2\right)$$

(2)

where ${\nabla }_t^2$ is the transverse nabla operator, n_t = n_x is xx component of the refractive index tensor, e_t = e_x is x component of the electric field for modes polarized along the x axis, n_t = n_y is yy component of the refractive index tensor, e_t = e_y is y component of the electric field for modes polarized along the y axis and β is the propagation constant. Equation (2) describes the propagation of quasi-TE and quasi-TM modes and could be represented as two scalar wave equations with polarization corrections characterizing the structural anisotropy of the waveguide (Equation (3)):

$$\begin{gathered} \frac{{\partial }^2{e}_x}{\partial x^2}+\frac{{\partial }^2{e}_x}{\partial y^2}+\left(k^2{n}_x^2-{\beta }^2\right)e_x=-\frac{\partial }{\partial x}\left(\frac{1}{n_x^2}\frac{\partial n_x^2}{\partial x}{e}_x\right), \\ \frac{{\partial }^2{e}_y}{\partial x^2}+\frac{{\partial }^2{e}_y}{\partial y^2}+\left(k^2{n}_y^2-{\beta }^2\right)e_y=-\frac{\partial }{\partial y}\left(\frac{1}{n_y^2}\frac{\partial n_y^2}{\partial y}{e}_y\right). \end{gathered}$$

(3)

In addition, in Equation (2), the term responsible for the interaction between the transverse components of the electric field is omitted, which follows from the assumption of a strong polarization of the TE and TM modes along the transverse axes [30].

Since it is assumed that the electric fields of the TE and TM modes are strongly polarized along the x and y axes, the solution of Equation (3) can be obtained using methods of separation of variables. According to [31], the wave Equation (3), which takes into account polarization effects, could be solved using EIM. EIM is a well-developed semi-analytical method for solving wave equations based on separation of variables. Widespread application of EIM is explained by its simplicity even among semi-analytical methods. The idea of EIM is that the cross-section of a real 3D waveguide is sequentially substituted by equivalent slab 2D waveguides. It is known that EIM works well when the propagating fundamental TE- and TM-modes are far from the cutoff condition [31,32,33]. To apply EIM, the waveguide must be slab-like. A waveguide is slab-like when it is strongly structurally anisotropic, for example, when the waveguide width is greater than its height. In addition, if a shallow-etched waveguide is considered, then a smaller etching depth will correspond to greater slab-likeness, since the refractive index gradient along the horizontal axis of the waveguide cross-section will be less expressed. However, the slab-likeness will vary for each configuration of waveguide materials and cross-sections.

The simulation accuracy depends on the direction of polarization of the electric field. In order to determine the polarization fraction of radiation along the transverse axes, Equation (4) was used [19]:

$${\gamma }_{TE}=\frac{\iint {\left|e_x\right|}^{2}dxdy}{\iint \left({\left|e_x\right|}^2+{\left|e_y\right|}^2\right)dxdy}.$$

(4)

A finite-difference frequency-domain (FDFD) method was chosen as the reference method for the effective indexes and the mode polarization fraction calculation since it solves the full-vector wave equation without any approximation, except for the finite difference partial derivative approximation. The results obtained with EIM were compared with the results obtained with FDFD and the relative errors of the effective indexes of modes were calculated. The FDFD calculations were carried out using Difra Lab software, which allows the simulation of diagonally anisotropic waveguides.

The size of the calculation window was 4.1 μm in width and 3.6 μm in height; a uniform grid with 300 cells for each axis was used for the simulation. The refractive indexes for LN and silicon oxide (SiO₂) were calculated using the dispersion equations:

$$\begin{gathered} n_{{LiNbO}_3}^o=\sqrt{1+\frac{{2.6734\lambda }^2}{{\lambda }^2-0.018}+\frac{{1.229\lambda }^2}{{\lambda }^2-0.059}+\frac{12.614{\lambda }^2}{{\lambda }^2-474.6}}, \\ {n}_{{LiNbO}_3}^e=\sqrt{1+\frac{{2.9804\lambda }^2}{{\lambda }^2-0.0205}+\frac{{0.5981\lambda }^2}{{\lambda }^2-0.066}+\frac{8.9543{\lambda }^2}{{\lambda }^2-416.08}}, \\ {n}_{Si{O}_2}=\sqrt{1+\frac{{0.6961663\lambda }^2}{{\lambda }^2-{0.0684}^2}+\frac{{0.4079426\lambda }^2}{{\lambda }^2-{0.1162}^2}+\frac{8.9543{\lambda }^2}{{\lambda }^2-{9.8962}^2}}, \end{gathered}$$

(5)

where $n_{{LiNbO}_3}^o$ and $n_{{LiNbO}_3}^e$ are the ordinary and extraordinary refractive indexes of LN, respectively, and $n_{Si{O}_2}$ is the refractive index of silicon oxide [34]. The simulation was carried out for a wavelength range from 1.5 μm to 1.6 μm.

2.2. The Selection of Waveguide Cross-Section Parameters for Calculations with EIM

In paper [29], the etching technology of waveguides based on bulk LN was developed in order to further its transfer to LNOI substrate. To determine the parameters of the waveguide cross-section, calculations of the mode composition of optical radiation were carried out. However, slope angle of sidewalls and anisotropy were not taken into account in [29]. The waveguide cross-section might be adjusted, taking in consideration the slope angle of sidewalls and anisotropy. Therefore, the target structure is a shallow-etched waveguide with a width of 1.1 μm, an etching depth of 0.3 μm, and a slope angle of sidewall of 75° (Figure 2).

In this paper, the applicability of EIM for calculation of the waveguide, shown in Figure 2, as well as waveguides with varying parameters of the cross-section was evaluated. The limits of applicability of the EIM to the calculation of anisotropic waveguides based on the X-cut of LNOI for a set of etching depths from 0.1 μm to 0.6 μm were studied by comparison of the calculation results of the optical mode characteristics obtained with FDFD and EIM methods. The waveguide width for each etching depth was chosen as slightly less than the width, at which the first-order modes appear. Additionally, the variation of the slope angle of the waveguide depending on the technological parameters should be considered. Therefore, the numerical study of the mode characteristics of waveguides with slope angles from 60° to 90° in steps of 10° was also implemented. The parameters of the cross-section of the studied waveguides are given in

Table 1. The parameters of the cross-section of the studied waveguides based on X-cut LNOI.

Etching Depth, μm	Waveguide Width, μm	Sidewall Angles, °
0.1	2.1	60–90
0.2	1.5
0.3	1.1
0.4	0.9
0.5	0.8
0.6	0.7

. In each simulation iteration, when the slope angle of sidewalls varied, the rest of the parameters were fixed. This was done in order to take into account the influence of the slope angle of sidewalls on mode hybridization and the additional error in the accuracy of calculating the effective refractive indexes using the EIM.

2.3. The Estimation of Optical Radiation Losses with EIM and the Payne–Lacey Model

As it shown in [29], the cross-sections of real (manufactured) integrated optical waveguides after the reactive ion etching stage differ from ideal ones (simulated). For the complex analysis of optical characteristics of real waveguides, the sidewall roughness was taken into account and optical radiation losses were evaluated. An RMS roughness and a correlation length $L_c$ typically characterize a random waveguide sidewall perturbation and are used for estimation of the radiation losses along the waveguide [35]. The LER parameter is commonly used to characterize the technological processes of waveguide fabrication. RMS might be calculated as LER/3.

Commonly, the sidewall roughness of the waveguide is inherited during reactive ion etching from the edge roughness of the mask [36]. For this reason, in the paper, the LER of the waveguide was determined as the LER of the mask. As shown in [29], the chromium hard mask was used for the fabrication of waveguides on an LNOI substrate. For LER analysis, a top view image of the chromium mask was obtained using a scanning electron microscope (SEM) Tescan Mira3 (Figure 3). In addition, the same analysis of the sidewall roughness of the fabricated waveguide on the LNOI substrate is not possible due to its charging during SEM investigation.

SEM images were obtained at 100 kx magnification at three different points for one chromium strip. Then, pixel-by-pixel image processing was carried out with the OpenCV2-Python module. The chromium mask edges were defined as the intensity maxima in each pixel line of the grayscale image obtained from the secondary electron detector. Such image processing eliminates the human factor in determining the edge of the chromium mask. The illustrative example of the SEM image and the plot of mask edges of the chromium strip with high LER value is shown on Figure 4.

For the LER calculation, a middle line was built for the left and right edges of the chromium strip. The LER was calculated as [37]:

$${LER}_k=3·\sqrt{\frac{\sum _{i=1}^n{\left(x_{ik}-\overline{x_k}\right)}^2}{n}},$$

(6)

where $i$—point at which deviation is calculated, $k$ = 1, 2 are the left and right edge lines of the chromium mask. It is assumed that the LER value is within the range of three standard deviations from the profile midline (calculated separately for the left and right edges).

For the entire mask edges of the chromium strip, LER is calculated using Equation (7):

$$LER=\sqrt{\frac{{LER}_1^2+{LER}_2^2}{2}}$$

(7)

For conversion of pixels to nanometers, the final LER value was multiplied by a factor of 2.5, which corresponded to ratio of 2.5 nm/pixel for images at 100 kx magnification and 1024 × 1024 resolution. The LER values of the hard chromium mask, formed with the process described in [29], ranged from 20 to 40 nm.

The method of photolithography used determines the mask edge roughness. In [29] the contact lithography was applied for waveguide pattern transfer to the LNOI substrate. According to [38], the typical correlation length $L_c$ for the contact lithography patterning is 160 nm.

After the LER and RMS determination, optical radiation losses could be evaluated with the Payne–Lacey model, which was initially developed for symmetric 2D slab waveguide [39,40]. The Payne–Lacey model may be used for studying optical losses in 3D waveguides while applying EIM to obtain a 2D slab equivalent to the original waveguide and approximating the 3D radiative modes with 2D planar radiative modes [35,38].

The accuracy of estimation of radiation losses by the Payne–Lacey model directly depends on the calculation accuracy of the effective refractive index by the EIM; as a result, it can also be used only for waveguides with those cross-sectional parameters for which the EIM has an acceptable calculation error.

3. Results

The FDFD simulation showed that fundamental TE- and TM-modes in the waveguide depicted on Figure 2 are strongly polarized along x and y axes, which is expressed in γ_TE values exceeding 98% (γ_TE is around 0 for TM-mode). Thus, the assumption that an anisotropic optical waveguide could be approximated by an isotropic waveguide with the refractive index corresponding to a certain polarization is correct for the considered cross-section (shown on Figure 2) with a known accuracy. Furthermore, for the specified cross section on Figure 2, simulations were carried out using the EIM and FDFD methods in the wavelength range from 1.5 to 1.6 μm.

The obtained results of the EIM simulations are in good agreement with the results obtained with the accurate FDFD method for the considered waveguide cross-section, and the relative error of the effective refractive index calculation does not exceed 1% for both polarizations (Figure 5).

In simulations of single mode waveguides, the dependencies of the effective refractive index and polarization fraction along the x and y axes from the etching depth were studied at the wavelength 1.55 μm (see the waveguides cross-section parameters in Table 1). The simulation showed that with increase of etching depth, the effective refractive indexes of the TE- and TM-modes decrease to values close to the substrate refractive index (SiO₂), and the polarization fraction of the TE-mode along the transverse axis decreases while the value for the TM-mode increases. However, the smaller the slope angle of sidewalls, the faster the effective refractive index and the fraction of polarization of the TE-mode decrease. In the case of the TM-mode, the slope angle of sidewalls affects the polarization fraction, limiting its maximum value (Figure 6).

Moreover, at an etching depth of about 0.4 μm, the values of the effective refractive indexes of the TE and TM modes become almost the same, which means that inter-mode hybridization occurs.

According to results obtained with FDFD method, when the etching depth increases further, the effective refractive index of the TM-mode becomes larger than the effective refractive index of the TE-mode (Figure 7a). The same is observed in the results obtained by EIM, however, in case of EIM, the inter-mode hybridization region is shifted to an etching depth of 0.5 μm (Figure 7b).

To evaluate the accuracy of EIM, the relative error of calculating the effective refractive indexes of the TE- and TM-modes by EIM and FDFD method was analyzed for the waveguides described in Table 1. As can be seen from Figure 8, the relative calculation error increases with an increase of the waveguide etching depth and with decrease of the sidewall slope angle. The plot shows that the calculation error of the effective refractive index by the EIM does not exceed 1% at the etching depths less than 0.3 μm and at the sidewall slope angle at least 70°.

After the EIM accuracy estimation, optical radiation loss could be calculated with the Payne–Lacey method. The calculation results are shown on Figure 9. As seen from the plot, the increase of LER from 20 to 40 nm leads to growth of radiation optical losses: from 0.5 to 1.8 dB/cm for the TE-mode, and from 1.8 to 7.1 dB/cm for the TM-mode at 1.55 μm. The longer the wavelength, the less sensitive it is to the sidewall roughness, resulting in lower optical radiation losses.

All simulations were run on an Intel core i7-6700K processor. The time and memory consumption were 15 s and about 10 MB for the FDFD simulation iteration and 1 s and about 15 KB for the EIM iteration, respectively.

4. Discussion

Using the EIM and FDFD methods, the effective refractive indexes of the fundamental TE- and TM-modes in the waveguide with a cross section shown in Figure 2 were calculated for the wavelengths 1.5–1.6 μm. The results demonstrate that as the wavelength increases, the mode localization in the waveguide core decreases. It is known that in this case the mode field is more sensitive to the waveguide–environment interface, as a result the hybridity of the mode becomes more expressed (the amplitude of longitudinal field components increases). Wave equations (Equation (3)) take into account the refractive index only along the axis parallel to direction of field polarization and do not consider the hybrid mode nature, which is expressed in the possible interaction of transverse electric field components. This leads to an increase in the calculating error in the effective refractive indexes with increase of wavelength (Figure 5b) by EIM. Nevertheless, the value of the relative error of calculations of the effective refractive index with EIM does not exceed 1.0% in wavelength range 1.5–1.6 μm for the considered waveguide cross section.

The executed calculations in this paper allow the estimation of the applicability limits of EIM for a simulation of anisotropic waveguides based on LNOI (Figure 8). We assume that if the calculating error of the effective refractive index does not exceed 1%, the EIM could be used for the rapid estimation of the parameters of anisotropic waveguides based on LNOI.

It is known that the accuracy of EIM depends on the proximity to the cutoff condition and slab-likeness of the waveguide for the propagating mode. Slab-likeness could be defined as the polarization fraction of radiation along one of the transverse axes of the waveguide, as well as the absence of a sharp refractive index gradient. The higher the polarization fraction of radiation along the one of the axes, the smaller the real part of the amplitudes of electromagnetic field $E_z$ and $H_z$ components, and the closer the propagating modes to TEM-modes.

This paper shows an evaluation of the accuracy of EIM applied to waveguides based on LNOI through an estimation of the calculating accuracy of the effective refractive indexes. That is why for the simulations with variation of the etching depth, the waveguide width was chosen to be slightly less than the width corresponding to the occurrence of the first-order propagating modes, i.e., studied modes were far from the cutoff condition. If the waveguide width is reduced, then the calculation accuracy of the effective refractive indexes will decrease as the propagating mode approaches the cutoff condition.

The results obtained with the FDFD method show that for the TE-mode the etching depth of 0.1 μm corresponds to the maximum polarization fraction along the horizontal axis (${\gamma }_{TE}$ = 99.8%), while the TM-mode is highly hybridized (60% ≤ ${\gamma }_{TM}$ ≤ 69%). Nevertheless, the calculating accuracy of the effective refractive index of the TM-mode with EIM remains comparable to the accuracy of the FDFD method for all studied sidewall slope angles of the waveguide. As the waveguide etching depth increases, the TE-mode becomes more hybrid, since the refractive index gradient becomes sharper along the horizontal axis. However, since in each calculation iteration the mode remains far from the cutoff condition, the polarization fraction drops only to values of about 96%. In turn, the polarization fraction of the TM-mode increases with increasing the waveguide etching depth to values close to 98%.

In addition, an increase of the etching depth leads to a decrease in the effective refractive indexes of both modes, and in the range of etching depths from 0.4 μm to 0.5 μm, a hybridization region of the TE- and TM-modes appears when their effective refractive indexes become equal. With a further increase of the etching depth, the effective refractive index of the TM-mode begins to exceed the effective refractive index of the TE-mode. These effects are possibly associated with the influence of structural anisotropy on the amplitudes of the electromagnetic field components.

Considering all of the above, the design of the waveguide must be carried out using precise methods (such as FDFD), since this determines the required cross section of the waveguide. The estimation of effective refractive indexes over a wide range of wavelengths could be performed using the EIM. This approach is valid because the typical difference of characteristics of simulated and fabricated waveguides often exceeds 1%.

5. Conclusions

The simulation results obtained by FDFD and EIM showed that EIM is applicable to the analysis of anisotropic waveguides on X-cut of LNOI directed along the Y crystallographic axis when the considered modes are strongly polarized along the transverse axes. The relative error is less than 1% for the set of cross sections with sidewall slope angles from 60° to 90° and etching depths from 0.1 to 0.3 μm for 1.55 μm wavelength. Typically, the waveguide is designed for the specific central wavelength. It can be assumed that if the mode is far from cutoff condition at the certain wavelength, then the EIM will provide accurate effective refractive index values over the wide wavelength range, as can be seen in Figure 5. However, if the waveguides are considered with sidewall angles less than 70° and etching depths of more than 0.3 μm, the relative error increases. The increase of the relative error is due to the fact that in the work, the sidewalls slope angle was not taken into account in the EIM simulations. However, the general EIM approach demonstrated in [29] could provide better accuracy when the sidewalls slope angle is approximated by the set of planar waveguides of different heights. In addition, the EIM can be applied for simulations of higher order modes. In these cases, the calculating accuracy of the effective indexes will also be determined by the proximity of the mode to the cutoff condition, the sidewall slope angle, and the etching depth as well as for the fundamental modes’ simulations. Our research group is developing the fabrication process of ridge waveguides based on LNOI and our preliminary results show that it is possible to obtain waveguides with a sidewall slope angle of about 80° (we intend to publish these results in the near future). Thus, we intend to apply the EIM for the fast estimation of mode characteristics in fabricated waveguides.

Moreover, EIM is a method based on the assumption that x and y electromagnetic field components are separable. It means that a three-dimensional problem of field propagation in a photonic integrated element can be reduced to a two-dimensional problem with the real refractive indexes replaced by the effective refractive indexes of the media, if the considered photonic element is parallel to one of the crystallographic axes. In addition, as shown in [28], EIM could be integrated with an S-matrix simulator. Thus, it is possible to reduce calculation time for fast circuit simulation.