Open Access
Add to BibSonomy
Abstract
We propose an efficient numerical method for the full-vectorial analysis of three-dimensional (3-D) optical waveguide discontinuities. In this method, the finite element method with higher adaptability and flexibility 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 reflectivities of the fundamental TE-like and TM-like modes are calculated. These results show unique vector properties and significantly differ from those of scalar analysis because various mode couplings between the field components occur at the discontinuity facet and they cannot be ignored.
© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. 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–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–12] and photonic crystal fiber [13–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$\rm \acute {e}$ 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.
2. POM based on FEM
2.1 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, ${\boldsymbol{\Phi}}$ denotes either the electric field $\boldsymbol {E}$ or the magnetic field $\boldsymbol {H}$, $p$ and $q$ are given by where $n$ is the refractive index.2.2 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 $\beta$, then the field ${\boldsymbol{\Phi} }$ in each element is expanded as
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 $\{\phi \}$, the finite element matrices $[K]$ and $[M]$, are given by with2.3 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 ${\boldsymbol{\Phi}}$ is regarded as an arbitrary field which is not certain eigenmode, $z$-dependence of light propagation cannot be expressed as $e^{-j\beta z}$. Thus, $\{\phi \}$ in (5) is replaced by $\{\Phi (z)\}$ considering $z$ direction dependence, and $\beta$ is replaced by $jd/dz$. After that, we obtain the following matrix form second-order ordinary differential equation:
where $[Q]$ is the propagation operator matrix which is defined by $\sqrt {[M]^{-1}[K]}$. As a solution to the differential equation (9), the electromagnetic field can be formally expressed as where $\left \{\phi _t^{(+)}\right \}$ and $\left \{\phi _t^{(-)}\right \}$ 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: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 $\boldsymbol{\Psi }$ as a pair of the electromagnetic field for $\boldsymbol{\Phi }$ in (1), it is given by
where $c$ is light velocity. $\boldsymbol{\Psi }$ is also regarded as $\boldsymbol {H}$ or $\boldsymbol {E}$ when $\boldsymbol{\Phi }$ indicates $\boldsymbol {E}$ or $\boldsymbol {H}$. Therefore, the power transmitted by the $i$-th eigenmode can be calculated as follows: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
2.4 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 $[A]_0$ and $[B]_0$ represent initial values of the iterations, $[I]$ is the unit matrix and $k$ is a number of iterations. $[A]_k$ and $[B]_k$ converge $[Q]$ and $[Q]^{-1}$ respectively for sufficiently large $k$. 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
Fig. 2. Calculation error of propagation operator by using DBI method and Pad$\rm \acute {e}$ approximation as a function of the iterative number of DBI method or Pad$\rm \acute {e}$ order.
3. 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.
3.1 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~\mu \textrm {m}$. The fundamental $E^x$ or $E^y$ mode with operating wavelength $\lambda = 1.3~\mu \textrm {m}$ is launched. In our analysis, considering the structural symmetry of the waveguide geometry, the quarter region ($3\times 2~\mu \textrm {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.
Fig. 4. Modal reflectivity of the fundamental modes as a function of the core width corresponding to the waveguide structure shown in Fig. 3(a).
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~\mu \textrm {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~\mu \textrm {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~\mu \textrm {m}$ is about two times larger than that in the case of $w=2.0~\mu \textrm {m}$.
Fig. 5. Field distribution of the incident and reflected field components of the fundamental $E^x$ modes for waveguide structure which core width $w = 2.0~\mu \textrm {m}$ shown in Fig. 3(a). Incident field: (a)$E_x$ and (b)$E_y$; Reflected field: (c)$E_x$ and (d)$E_y$.
Fig. 6. Field distribution of the incident and reflected field components of the fundamental $E^x$ modes for waveguide structure which core width $w = 0.4~\mu \textrm {m}$ shown in Fig. 3(a). Incident field: (a)$E_x$ and (b)$E_y$; Reflected field: (c)$E_x$ and (d)$E_y$.
3.2 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~\mu \textrm {m}$, $t_1=2~\mu \textrm {m}$, and $t_2=1.3~\mu \textrm {m}$. Considering the structural symmetry, the half region ($3\times 4~\mu \textrm {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.
Fig. 8. Modal reflectivity of the fundamental modes as a function of the core width corresponding to the waveguide structure shown in Fig. 7(a).
Figure 9 shows the convergence behavior of the calculated modal reflectivity of the fundamental $E^x$ mode corresponding to $w = 2~\mu \textrm {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.
Fig. 9. Modal reflectivity of the fundamental modes and calculation time, as a function of the number of unknown variables in entire analysis domain.
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~\mu \textrm {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.
Fig. 10. Modal reflectivity of the fundamental modes and matrix error calculated by Eq. (23), as a function of the number of iteration of DBI method.
The electromagnetic field distributions of the incident, reflected and transmitted waves at the discontinuity facet between the rib waveguide with $w = 2~\mu \textrm {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~\mu \textrm {m}$ shown in Fig. 7(a). Incident field: (a)$E_x$ and (b)$H_x$; Reflected field: (c)$E_x$ and (d)$H_x$; Transmitted field: (e)$E_x$ and (f)$H_x$.
Fig. 12. Field distribution of the incident, reflected and transmitted field components of the fundamental $E^x$ modes and $E^y$ modes for scalar wave approximation analysis with the waveguide structure which core width $w = 2.0~\mu \textrm {m}$ shown in Fig. 7(a). Incident field: (a)$E_x$ and (b)$H_x$; Reflected field: (c)$E_x$ and (d)$H_x$; Transmitted field: (e)$E_x$ and (f)$H_x$.
4. 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).
References
1. M. Eguchi and Y. Tsuji, “Influence of reflected radiation waves caused by large mode field and large refractive index mismatches on splice loss evalution between elliptical-hole lattice core holey fibers and conventional fibers,” J. Opt. Soc. Am. B 30(2), 410–420 (2013). [CrossRef]
2. H.-P. Nolting and G. Sztefka, “Eigenmode matching and propagation theory of square meandertype couplers,” IEEE Photon. Technol. Lett. 4(12), 1386–1389 (1992). [CrossRef]
3. M. Koshiba, Optical Waveguide Theory by the Finite Element Method (KTK Scientific/Kluwer Academic, 1992).
4. Y. Tsuji and M. Koshiba, “Finite element method using port truncation by perfectlymatched layer boundary conditions for optical waveguide discontinuity problems,” J. Lightwave Technol. 20(3), 463–468 (2002). [CrossRef]
5. K. S. Ye, “Numerical solution of initial boundary value problems involving Maxwell’s equations,” IEEE Trans. Anlennas Propag. AF-14, 302–307 (1966).
6. Y.-P. Chiou and H.-C. Chang, “Analysis of optical waveguide discontinuities using the Padé approximants,” IEEE Photon. Technol. Lett. 9(7), 964–966 (1997). [CrossRef]
7. H. El-Refaei, I. Betty, and D. Yevick, “The application of complex Padé approcimants to reflection at optical waveguide facets,” IEEE Photon. Technol. Lett. 12(2), 158–160 (2000). [CrossRef]
8. H. A. Jamid and Md. Zahed M. Khan, “3-D Full Vectorial Analysis of Strong Optical Waveguide Discontinuities Using Pade Approximants,” IEEE J. Quantum Electron. 43(4), 343–349 (2007). [CrossRef]
9. S. Wu, J. Xiao, and X. Sun, “Full-vectorial analysis of optica waveguide discontinuities using denman-beavers iterative scheme,” J. Lightwave Technol. 33(2), 511–518 (2015). [CrossRef]
10. S. Inao, T. Sato, H. Hondo, M. Ogai, S. Sentsui, A. Otake, K. Yoshizaki, K. Ishihara, and N. Uchida, “High density multicore-fiber cable,” Proc. Int. Wire & Cable Symp., pp. 370–384, 1979.
11. G. Le Noane, D. Boscher, P. Grosso, J. C. Bizeul, and C. Botton, “Ultra high density cables using a new concept of bunched multicore monomode fibers: A key for the future FTTH networks,” Proc. Int. Wire & Cable Symp., pp. 203–210, 1994.
12. T. Hayashi, T. Taru, O. Shimakawa, T. Sasaki, and E. Sasaoka, “Design and fabrication of ultra-low crosstalk and low-loss multi-core fiber,” Opt. Express 19(17), 16576–16592 (2011). [CrossRef]
13. J. C. Knight, T. A. Birks, P. St. J. Russell, and D. M. Atkin, “All-silica single mode optical fiber with photonic crystal cladding,” Opt. Lett. 21(19), 1547–1549 (1996). [CrossRef]
14. T. A. Birks, J. C. Knight, and P. St. J. Russell, “Endlessly single-mode photonic crystal fiber,” Opt. Lett. 22(13), 961–963 (1997). [CrossRef]
15. J. Broeng, D. Mogilevstev, E. Barkou, and A. Bjarklev, “Photonic crystal fibers: a new class of optical waveguides,” Opt. Fiber Technol. 5(3), 305–330 (1999). [CrossRef]
16. A. R. Bhagwat and A. L. Gaeta, “Nonlinear optics in hollow-core photonic bandgap fibers,” Opt. Express 16(7), 5035–5047 (2008). [CrossRef]
17. S. Kawai, A. Iguchi, and Y. Tsuji, “Study on high precision and stable finite element beam propagation method based on incomplete third order Hybrid Edge/Nodal Element,” J. Lightwave Technol. 36(11), 2278–2285 (2018). [CrossRef]
18. J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. 114(2), 185–200 (1994). [CrossRef]
19. M. Koshiba and Y. Tsuji, “Curvilinear hybrid edge/nodal elements with triangular shape for guidedwave problems,” J. Lightwave Technol. 18(5), 737–743 (2000). [CrossRef]
20. T. Q. Tran and S. Kim, “Stability condition of finite-element beam propagation methods in lossy waveguide,” IEEE J. Quantum Electron. 50(10), 808–814 (2014). [CrossRef]
21. K. Morimoto and Y. Tsuji, “Analysis of butt coupling of optical waveguides using propagation operator method based on finite element method,” IEICE TRANSACTIONS on Electronics (Japanese Edition) J101-C(5), 210–216 (2018).
