Full-vectorial analysis of optical waveguide discontinuities using a propagation operator method based on the finite element scheme discontinuities using a propagation operator method based on the finite element scheme

: We propose an eﬃcient numerical method for the full-vectorial analysis of three-dimensional (3-D) optical waveguide discontinuities. In this method, the ﬁnite element method with higher adaptability and ﬂexibility is employed to discretize the waveguide cross section. In order to calculate the square root of the characteristic matrix, the Denman-Beavers iterative scheme is used. Applying this method to 3-D strongly guiding waveguide discontinuity problems, the modal reﬂectivities of the fundamental TE-like and TM-like modes are calculated. These results show unique vector properties and signiﬁcantly diﬀer from those of scalar analysis because various mode couplings between the ﬁeld components occur at the discontinuity facet and they cannot be ignored. that in the vector analysis due to the diﬀerence of the light conﬁnement. These results suggest that the accuracy of the scalar wave approximated analysis in the problem with relatively large lateral refractive index diﬀerence is insuﬃcient and the full vector analysis is required.


Introduction
With the recent development of communications technology, optical waveguide discontinuities have been an important subject in the design and investigation of novel optical devices and systems. These discontinuities exist in a lot of situations such as laser facets, waveguide ends, butt-joint of different types of waveguides, and cause universal problems of power loss. Therefore, an accurate and efficient method for understanding the reflection and transmission characteristics at these optical waveguide discontinuities is highly demanded.
There are a number of analytical or numerical techniques proposed in previous works for the analysis of optical waveguide discontinuity problems. In the case of weakly guiding waveguide discontinuities where mode field and refractive index mismatches between different waveguides are small, the transmission efficiency can be evaluated by using simple method based on overlap integral without considering reflected radiation wave [1]. This method is not suitable for strongly guiding waveguide problems because the reflected radiation wave cannot be necessarily ignored in such waveguides. The eigenmode expansion texhnique is often applied to three-dimensional longitudinal discontinuities and only guided modes need to be treated in the wakly guiding problem [2]. However, this technique is not easy to analyze the problems with complicated index profiles because a lot of radiation and evanescent modes must be taken into account and large computational cost is required for mode expansion. On the other hand, direct numerical approaches such as finite element method (FEM) [3,4], and finite difference time domain (FDTD) [5] method has high accuracy and flexibility in the strongly guiding waveguide problems. However, these methods require enormous computer resources, especially for discretizing three-dimensional models. Therefore, we consider that the propagation operator method (POM) [6][7][8][9] seems to be one of the best effective approaches for solving the optical waveguide discontinuity problems. POM can evaluate both reflected and transmitted field including radiation just analyzing the boundary of waveguide discontinuities. In the past studies of POM, finite difference method (FDM) is usually employed to discretize the waveguide cross section, and these recent works in [8,9], show that they are competent to analyze including evanescent modes in a full-vectorial analysis. However, its accuracy depends on the numerical discretization of cross-sectional geometries at discontinuity facets, and in order to improve the adaptability and flexibility to complicated waveguide geometries, the computational efficiency will be degraded by requiring enormous discretizing meshes. Particularly in butt-coupling between optical fibers such as multi core fiber [10][11][12] and photonic crystal fiber [13][14][15][16], which is more complicated structure than conventional single mode fibers, more accurate structural representation is considered to become important.
Therefore, we employ FEM as an alternative discretization scheme of POM to make meshes more conform for full-vectorial analysis of three-dimensional optical waveguide discontinuities. In POM analysis, the square root of a characteristic matrix which is used to determine the discontinuity condition has to be calculated. This matrix is approximated by Denman-Beavers iterative (DBI) scheme with a branch-cut rotation in the complex plane whose accuracy and stability are reported to be higher than those of Padé approximation [9]. The proposed method is applied to two types of optical waveguide facets exhibiting strong refractive index discontinuities and the calculated reflectivity is compared with that of the full-vectorial analysis method reported so far. These results of the fundamental TE-Like and TM-Like modes show the unique vector properties and significantly differ from those of the scalar analysis because various mode couplings between the field components occur at the discontinuity facet and they cannot be ignored.
The rest of this paper is organized as follows. We describe the formulation of FEM and POM in Section II, DBI scheme is mentioned in the same section. After that, numerical examples are given and discussed in Section III with concluding remarks in Section IV.

Basic equation
We consider a three-dimensional (3-D) optical waveguide connected to a discontinuity cross section which is shown in Fig. 1, where x and y are the transverse directions, and z is the propagation direction. From Maxwell's equations, the vectorial wave equation is expressed as follows: where k 0 is the free space wavenumber, Φ denotes either the electric field E or the magnetic field H, p and q are given by where n is the refractive index.

Finite element discretization
Dividing the waveguide cross section into hybrid edge/nodal elements [17], and considering light wave propagating along the z-direction with the propagation constant β, then the field Φ in each element is expanded as where {φ t } and {φ z } are, respectively, the edge and nodal variables in each element, {U} and {V } are the shape function vectors for edge elements and {N} is the shape function vector for nodal elements. By using third-order nodal elements and QT/CuN edge elements [17], we can reduce unknown variables of φ and expect to calculate efficiently compared to the FEM with second-order elements.
Substituting (4) into (1) and applying FEM based on the Galerkin method to the transverse xy plane, we obtain the following matrix eigenvalue equation: where {0} is a null vector, and the transverse field components {φ}, the finite element matrices , are given by where e extends over all the different elements.

POM analysis
Now, we assume that there is no longitudinal refractive index variation except for discontinuity facets as shown in Fig. 1. When the electromagnetic field Φ is regarded as an arbitrary field which is not certain eigenmode, z-dependence of light propagation cannot be expressed as e −jβz . Thus, {φ} in (5) is replaced by {Φ(z)} considering z direction dependence, and β is replaced by jd/dz. After that, we obtain the following matrix form second-order ordinary differential equation: As a solution to the differential equation (9), the electromagnetic field can be formally expressed as represent column vectors standing for eigenmode amplitude of forward and backward waves, respectively. Considering there is no backward wave in the output side waveguide, the boundary condition of a certain electromagnetic field at the discontinuity boundary can be expressed as follows: where the subscripts 1 and 2 indicate the input waveguide side and the output waveguide side, respectively. Next, the longitudinal boundary condition is rigorously established by the continuity of the transmission power at the discontinuity boundary that is expressed using the transverse field components. Defining Ψ as a pair of the electromagnetic field for Φ in (1), it is given by where c is light velocity. Ψ is also regarded as H or E when Φ indicates E or H. Therefore, the power transmitted by the i-th eigenmode can be calculated as follows: In addition, utilizing (5)-(8), {φ z,i } can be written by Therefore, (13) is transformed to And now, we have to consider transmission power as an arbitrary propagation wave to account for all eigenmodes. Thus, β i in (15) Since the power flow in z-direction has to be continuous at the discontinuity boundary, the continuity relation at the waveguide discontinuity facet is expressed as Therefore, by utilizing (11) and (17), we can derive the reflection amplitude φ (−) t,1 as follows: where [F i ](i = 1, 2) can be regarded as an impedance operator or an admittance operator depending on whether φ is an electric field or magnetic field. The transmission amplitude can also be obtained by using the incident amplitude and the reflection amplitude, from the boundary condition (11).

DBI method
In order to obtain the propagation operator matrix, the square root matrix has to be calculated. This matrix can be obtained by using the eigenvalue decomposition method. However, it suffers from a large computational cost, therefore, an efficient and highly accurate approximation method is required. In our research, DBI method [9] which is an iterative scheme based on a recurrence formula is employed and is given as follows: [ where Here, when solving the square root of the real matrix by the DBI method, the solution often does not converge because the square root matrix obtained by this iterative scheme encounter only real solutions if the initial matrix is real. Therefore, the branch cut rotation technique [9] which transforms the real matrix into the complex matrix by multiplying a complex coefficient is often used to permit a complex matrix as a solution matrix. We show the approximation error of the propagation operator as a function of the iteration number of DBI method in Fig. 2. As a test matrix, we use the matrix constructed in a later numerical example of a rib waveguide. The approximation error is defined by where N is the matrix size, a ij is matrix components of [Q] 2 andã ij is matrix components obtained by squaring [Q], and [Q] is the square root matrix calculated by DBI method. Figure 2 also indicates the calculation results by Padé approximation which widely used in the past study [8], however, it has numerical instability and does not necessarily converge in the case of the propagation operator matrix constructed by finite element matrices. On the other hand, in the DBI method, the error reduces to about 10 −15 in 12 iterations and converges. Moreover, since (17) and (18) are independent of each other, it is possible to improve the calculation efficiency by parallel calculation.

Numerical simulation results
The proposed method is utilized to calculate the modal reflectivities of the three-dimensional waveguide discontinuities where the waveguide is assumed to be terminated by the air facet.

Buried waveguide
First, we consider the buried waveguide as shown in Fig. 3(a) where the refractive index of core and cladding are n 1 = 3.54, n 2 = 3.17, respectively, and the core height is h = 0.4 µm. The fundamental E x or E y mode with operating wavelength λ = 1.3 µm is launched. In our analysis, considering the structural symmetry of the waveguide geometry, the quarter region (3 × 2 µm) of the waveguide cross-section is discretized by 66 third-order triangular elements, as shown in Fig. 3(b). In this case, the number of unknown variables is 696, and it is confirmed that sufficient analysis accuracy can be obtained by using QT/CuN element. Figure 4 shows the modal reflectivities of the fundamental E x and E y modes as a function of the core width w. These results of the proposed method agree well with the results obtained by three-dimensional FEM analysis and POM analysis based on FDM which is demonstrated in [8]. For the sake of comparison, the modal reflectivities obtained by scalar wave approximation analysis of the present POM are also indicated in the same figure. It is seen that the results  by scalar wave analysis get closer to those by the full-vectorial analysis when the core width becomes relatively large. However, a large deviation of the modal reflectivities occurs when the core width becomes small. This is because scalar wave analysis ignores the minor field components, although the influence of the minor field components becomes large when the core width decreases.
The field distribution of the incident and reflected waves in the fundamental E x mode incidence case where the core width is w = 2.0 µm are illustrated in Fig. 5. We note that these amplitudes are normalized by the maximum amplitude of each component. In terms of the reflected modes, a spurious reflection occurrs from the end of the analysis region can be seen because absorption boundary condition such as perfectly matched layer (PML) [18] is not placed at the end of the analysis region. Figure 6 shows the field distribution of the incident and reflected waves in the fundamental E x mode and E y mode incidence cases when the core width is w = 0.4 µm. Comparing the results in Figs. 5 and 6, the influence of the minor field component becomes large that the maximum value of the minor reflected field amplitude with w = 0.4 µm is about two times larger than that in the case of w = 2.0 µm.

Rib waveguide
Next, the proposed method is utilized to analyze the typical rib waveguide shown in Fig. 7(a). The refractive indices of the core and the substrate are set as n 1 = 3.44 and n 2 = 3.34, and the cover region is assumed to be air whose refractive index is 1. The structural parameters are set as h = 1.1 µm, t 1 = 2 µm, and t 2 = 1.3 µm. Considering the structural symmetry, the half region  (3 × 4 µm) of the waveguide cross-section is discretized by 56 third-order triangular elements, as shownin Fig. 7(b). The number of unknown variables is 573 with good convergence property. Figure 8 shows the modal reflectivity of the fundamental E x and E y modes as a function of the core width. These results of the proposed method and the three-dimensional FEM show a similar behavior while the results of the scalar wave approximation analysis greatly differ from those of the vectorial analysis when the core width is small. Figure 9 shows the convergence behavior of the calculated modal reflectivity of the fundamental E x mode corresponding to w = 2 µm as a function of the number of unknown variables. In this case, the iteration number of DBI method is fixed to 10 times. It is known that convergence depends on the order of the element used in FEM. Figure 9 also shows the results obtained by   using the LT/QN elements and CT/LN elements [19] in addition to the QT/CuN elements. In terms of the highest order QT/CuN element, the reflected power almost converges when the number of unknown variables is about 600. However, more unknown variables are required for ensuring the convergence in the case of using LT/QN element, moreover, it is understood that the lowest order CT/LN element further degrades convergence. It can be seen that CPU time rapidly increases with respect to the total number of unknown variables because the propagation operator matrix is dense. Therefore, QT/CuN element which can obtain convergence value in sufficiently small unknown variables can achieve the best efficient calculation. Figure 10 shows the modal reflectivity of the fundamental E x mode and the matrix error of the propagation operator as a function of the number of iterations in DBI method. Here, w = 2 µm is assumed and the matrix error is defined by (23). In this example, QT/CuN elements are used and the number of unknown variables is 573. From Fig. 10, it can be seen that the modal reflectivity almost converges on about 6 iterations before the matrix error sufficiently converges. The CPU time of the whole analyzing process is 34 seconds at 6 times iteration of DBI method, and the CPU time for calculating the propagation operator matrices on the incident side and the output side is 31 seconds. Hence, it is found that the total calculation time greatly depends on the calculation time of the propagation operator matrix. We also note that the calculation time of the propagation operator matrix increases linearly with the number of iterations of the DBI method. Thus, reducing the number of iteration of DBI method directly leads to efficient computation. The electromagnetic field distributions of the incident, reflected and transmitted waves at the discontinuity facet between the rib waveguide with w = 2 µm and air are shown in Figs. 11 and 12, when the fundamental E x mode or E y mode are launched. Figures 11 and 12 show the amplitudes of the electric and magnetic field of the vector wave analysis and the scalar wave analysis, respectively. In these figures, the amplitudes are normalized by the maximum value of each field component for us to be easy to observe the smaller radiated field. We can see that the fields obtained by the scalar the vector analysis are different from each other especially around the material boundary between the rib and air. As a result, the radiated field in the scalar analysis seems to be larger than that in the vector analysis due to the difference of the light confinement. These results suggest that the accuracy of the scalar wave approximated analysis in the problem with relatively large lateral refractive index difference is insufficient and the full vector analysis is required. Fig. 11. Field distribution of the incident, reflected and transmitted field components of the fundamental E x modes and E y modes for full-vectorial wave analysis with the waveguide structure which core width w = 2.0 µm shown in Fig. 7

Conclusion
In order to analyze the three-dimensional waveguide discontinuity problems, an efficient numerical approach of the full-vectorial propagation operator method based on the finite element scheme is proposed. In order to perform vector wave analysis based on FEM, the POM is reformulated, and the square root matrix required to obtain the characteristic matrix is efficiently calculated by using DBI method.
The numerical results of the reflected fundamental mode power in the problem of two types of waveguide facets clearly differ from those of the scalar wave analysis and the necessity of the full-vectorial analysis is demonstrated. Furthermore, the results based on vector wave analysis are in reasonable agreement with those of the full-vectorial three-dimensional FEM which is a more rigorous analysis method and in which PML boundary condition is applied in every direction. In these numerical examples by POM, PML to suppress spurious reflection from the artificial computational boundary is not imposed because PML may sometimes degrade the computational stability [20]. Although in these problems, the obtained fundamental mode power does not seem to be affected to a great extent whether PML is applied or not, the accuracy may be somewhat degraded especially when the relatively large diffraction occurs and the lateral light propagation at the discontinuity cross-section cannot be ignored [21]. Therefore, the investigation of the stable absorption boundary conditions will make this method greatly useful approach for various applications and we would like to study them in the future work.

Funding
Japan Society for the Promotion of Science (JSPS) (18K04276).