Full-vectorial finite element method in a cylindrical coordinate system for loss analysis of photonic wire bends

: This paper presents a new full-vectorial finite-element method in a local cylindrical coordinate system, to effectively analyze bending losses in photonic wires. The discretization is performed in the cross section of a three-dimensional curved waveguide, using hybrid edge/nodal elements. The solution region is truncated by anisotropic, perfectly matched layers in the cylindrical coordinate system, to deal properly with leaky modes of the waveguide. This approach is used to evaluate bending losses in silicon wire waveguides. The numerical results of the present approach are compared with results calculated with an equivalent straight waveguide approach and with reported experimental data. These comparisons together demonstrate the validity of the present approach based on the cylindrical coordinate system and also clarifies the limited validity of the equivalent straight waveguide approximation.


Introduction
The realization of high-density on-chip integration of optical components, or microphotonics, has been the subject of intense investigation in the new technological era. In particular, high index contrast waveguides with a submicron cross-section, such as silicon photonic wires [1][2][3], have attracted considerable attention. Various ultra-small functional optical components based on the photonic wires have been proposed and demonstrated [4][5][6][7][8][9].
Since most of these components contain bends, the development of the integrated circuits requires accurate recognition of the bend characteristics. Indeed, it has been demonstrated that due to the high confinement of light in the photonic wires, the bending radii can be reduced to less than several micrometers without significant losses [10,11]. In the design of the bend structures, numerical techniques play a vital role because the analytic determination of the loss properties is not always viable. To this end, progressive theoretical studies have been conducted using the finite-difference time-domain method [12][13][14][15]. However, this type of method directly computes a three-dimensional space, and designs based on this approach are often inefficient in terms of computational cost.
In contrast, modal approaches to bending waveguides, such as these based on the finite difference method [16][17][18] and the finite element method (FEM) [19][20][21][22], can be more efficient, because the solution regions decrease to two-dimensional cross sections. In addition, these modal approaches provide insights into the physical basis of bending losses in each individual eigenmode. Importantly, these approaches are divided into two categories according to coordinate systems in the formulations, i.e. approximate equivalent straight waveguide (ESW) formulations [21,22] and rigorous cylindrical coordinate system (CCS) formulations [16][17][18][19]. The ESW approximation has widely been used in the loss analysis of conventional weakly guiding structures and the validity has been confirmed in a rib waveguide geometry [17,22,23]. It is expected, however, for this approximation to become invalid in bends with smaller bending radii [23]. It has been demonstrated, in the curved rib waveguide, that the difference in the calculated effective indices between the ESW approximation and the CCS approach becomes significant when the bending radii reach several micrometers, although the difference in the calculated losses is still insignificant [17,22].
In this context, this paper proposes a full-vectorial FEM based on the CCS for efficient and accurate analysis of photonic wire bends. The analysis should be based on the CCS, because bend radii of photonic wire bends can be less than a few micrometers. The analysis should also be based on the full-vector wave equation, because modes in high index contrast waveguides are hybrid and the hybridness is further magnified due to the lack of symmetry in the bending structure. In addition, the FEM exhibits several advantages over the finite difference method in terms of geometrical adaptability and material generality, for modeling arbitrary geometries and materials of various compositions, although the application in this paper is for simple rectangular waveguides consisting of homogeneous isotropic media.
Radiation losses can be identified from the imaginary part of the complex propagation constants of the leaky modes. To deal properly with these leaky modes, the solution region must be truncated with an appropriate absorbing boundary condition. One of the most flexible and efficient absorbing boundary conditions is the perfectly matched layer (PML) [24]. Recently, the PML technique has been extended to a CCS [25] and implemented in finite difference analysis of curved waveguides [17,18]. However, this implementation procedure is based on an analytic continuation of the frequency-domain Maxwell's equations to complex space, which does not provide an easy interface for methods other than finite difference approaches [26,27]. Alternatively, so-called anisotropic PMLs [26,27], which are represented as artificial anisotropic media with properly chosen constitutive parameters, can be effectively applied in many kinds of numerical analyses regardless of the discretization methods. We therefore implement the anisotropic PML to the FEM in the local CCS, based on the formulation for the standard CCS [26].
The remainder of this paper is structured as follows. In section 2, we first describe the anisotropic PML formulation for the local CCS and then present the full-vectorial FEM formulation. The solution region is divided into hybrid edge/nodal elements [28], to solve the vector wave equation in a precise and straightforward fashion with the three components of electric or magnetic fields. In section 3, bending losses in silicon wire waveguides are evaluated using this method. The numerical results of the present approach are compared with results calculated with the ESW approach and with reported experimental data. These comparisons together demonstrate the validity of the present approach based on the CCS and also clarify the limited validity of the ESW approximation. In section 4, the important features of the FEM and the findings in the loss analysis are summarized.

Formulation
The following description relates to a curved optical waveguide shown in Fig. 1. Instead of the standard CCS, our formulation is based on a local CCS [29], which is commonly used in the analysis of curved waveguides [17,18,23]. Modes propagate along the curved z-axis. The solution region is on the x-y plane and the x = 0 plane corresponds to the waveguide center. R represents the bend radius. The open region is truncated with the anisotropic PML in the local coordinates, formulated in the following.

Anisotropic perfectly matched layer
PMLs are equivalent to analytic continuations of Maxwell's equations to complex variable spatial domains [30]. In the local CCS, this analytic continuation is expressed in the form where s x (x), s y (y), and s z (z) are the complex stretching variables [30] and T denotes the transpose. Table 1 shows the values of s x (x), s y (y), and s z (z) in each PML region, numbered in Fig. 1(b). We assume the value of s i (i = 1, 2, 3) in Table 1 as where m is a positive integer, ρ is the distance from the border between the PML and modeled region, d PML is the thickness of the PML, and δ is the loss angle on the exterior boundary of the PML (ρ = d PML ).
From (1), we also have Understanding (6), the following identity can be verified for any vector function A(r) in the PML and in the real space: In the real space, Maxwell's equations are Inside the complex space PML, Maxwell's equations are modified as where the superscript 'c' refers to complex space PML fields [27]. Using (3) and (7), we can recast (11a) and (11b) in the real space as We here introduce a new set of field definitions: where the superscript 'a' refers to Maxwellian PML fields [27]. Substituting (13a) and (13b) into (12a) and (12b), we have with the constitutive tensors characterizing the PML

Full-vectorial finite element method
From (14a) and (14b) where k 0 is the free-space wavenumber. Dividing the waveguide cross section in Fig. 1(b) into a number of hybrid edge/nodal elements [28], we can represent the transverse field components, φ x and φ y , and the longitudinal field component, φ z , within an element as where β is the propagation constant, {φ t } e and {φ z } e are, respectively, the edge-and nodalvariable vectors for each element, {U} and {V} are the shape function vectors for the edge elements [28], {N} is the shape function vector for the nodal elements [28], and {0} is a null vector.
Substituting (18) to (16) and applying the Galerkin method, we obtain the weak form: where n is the outward unit normal vector to a boundary. We assume the Dirichlet condition and the Neumann condition on boundaries The first term on the left hand side of (24) is rewritten for an element as Similarly, the second term is rewritten as Substituting (25) and (26) to (24), we obtain the eigenvalue matrix equation for the guided modes in the curved waveguide: where [0] is a null matrix and Σ e represents the summation over all different elements. In the calculation of the integration in (31) -(34), we employ numerical integration based on the seventh order Gauss-Legendre quadrature [31]. This integration scheme samples the values of x and y at the seven sampling points in each element.
The method of inverse iteration with shifts of origin can be efficiently applied to solve this eigenvalue equation. The matrix equations to be solved in this iterative procedure are sparse complex linear systems of equations, which can be solved with appropriate techniques, such as a multifrontal method [32].

Numerical results and discussion
In this section, we apply the vectorial FEM in the CCS to the identification of radiation losses of fundamental modes in bent photonic wires. The results obtained with a FEM based on the ESW approximation [21,33] can demonstrate the limited validity of this approximation. The schematic representation of the curved photonic wire is shown in Fig. 1. The refractive indices of the core and cladding are set as n core = 3.5 and n clad = 1.45, respectively, assuming Si and SiO 2 . The waveguide width and height are W = 445 nm and H = 220 nm, respectively. R represents the bend radius. These parameters are set to resemble the structure fabricated in the work of Vlasov and McNab [11].
The computational domain is divided into more than 20000 linear tangential and quadratic normal (LT/QN) elements [28]. The typical computational time is ~ 50 seconds on a personal computer with a Pentium 4 processor (3.2 GHz). To evaluate the bending losses, we paid close attention in choosing the parameters for the PML: x max = 2 μ m, y min = -1 μ m, y max = 2 μ m, 6.8 < tan δ < 9.8, and d PML = 1 μ m. Each PML region is divided into 15 layers of elements in the thickness direction.
First, we analyze the bending losses spectra of the 90° bends measured in [11], using these approaches. Figure  μ m, the results obtained with the ESW approximation and the CCS approach are highly consistent, as observed in Figs. 2(b) and (c). These results also match well with the experimental data. For R = 1 μ m, however, the ESW approximation incorrectly estimates the radiation losses with lower values than those calculated with the CCS approach at longer wavelengths, as observed in Fig. 2(a). The results of the CCS approach agree substantially with the experimental data, although there are dips at longer wavelengths in the measured loss spectrum and the measured losses are slightly higher than the calculated losses at the remaining wavelengths. The origin of the dips in the experimental data is unclear, but may be due to a spatial periodicity, inadvertently introduced   in the fabrication process, as mentioned in [11]. On the other hand, the cause of the difference in the losses at the remaining wavelengths can be explained more clearly by considering the possible effect of transition losses occurring at the junction of two waveguides with different bending radii. The FEM modal solvers analyze pure bending losses due to radiation in curved waveguides with a constant curvature. In the experimental situation, however, the insertion losses of waveguide bends contain both radiation losses and transition losses, as well as losses due to fabrication errors. Recognizing that the differences between the experimental data and calculated results are around 0.1 dB/90°, we can estimate the transition loss for the bend with R = 1 μ m at approximate 0.05 dB. This value assumes that the value of the transition loss between the straight waveguide in the input side and the curved waveguide and the value of the transition loss between the curved waveguide and the straight waveguide in the transmitted side are the same. Such transition losses can be reduced by introducing, for example, lateral offsets so that the field profiles overlap better between the straight waveguide and curved waveguide. In other words, the radiation losses identify the lowest possible value of losses through a bend. In this case, it is understood that the bending losses can be reduced to negligible levels even for R = 1 μ m at wavelengths shorter than μ m, the ESW approach incorrectly estimates the radiation losses with lower values than those evaluated with the CCS approach at wavelengths where the losses increase. The difference in the losses between the two approaches increases for the smaller values of R. Compared with the case of the quasi-TE modes, this difference is more significant for each radius R. The results of the CCS approach agree substantially with the experimental data, excepting the dip at the wavelength λ ≈ 1.27 μ m in the measured loss spectrum for R = 2 μ m and the difference in possible transition losses of approximate 0.1 dB for R = 1 and 2 μ m. Next, we investigate the dependence of the radiation losses on bend radius R, to identify the limit of the validity of the ESW approximation. Figures 4(a) and (b) show the R dependence for the quasi-TE mode and for the quasi-TM mode, respectively, at λ = 1.55 μ m. For the quasi-TE mode, the calculated losses with the CCS and ESW approaches are in firm agreement for relatively large bend radii (R > 2 μ m). Similarly, for the quasi-TM mode, the results with both methods agree very well for larger bend radii (R > 20 μ m). The ESW approach has been verified for large bend radii, R > 1 mm, in a rib waveguide [21]. We can therefore validate the CCS approach from the firm agreement with the ESW approach for relatively large bend radii. In contrast, the difference in the calculated losses between the ESW approximation and CCS approach, for both modes, clearly increases as R decreases. Such significant differences have never been observed in the loss analysis of conventional weakly guiding structures, such as rib waveguides [17,22,23]. These differences suggest that, in designs of optical circuits with a number of bends with ultra-small bend radii, errors due to the misestimation of the radiation losses can damage output characteristics. For example, in the designs of add/drop filters with a wide free spectral range based on Mach-Zehnder interferometers or microring resonators [2, 6], these errors can lead to undesirable insertion losses. It is also clear from the figures that the differences for the quasi-TM mode are more significant than those for the quasi-TE mode.
This variation in the losses can be understood from the difference in the Poynting vector, S, distributions. Figures 5(a) and (b) show the Poynting vector distributions at R = 2 μ m for the quasi-TE mode and for the quasi-TM mode, respectively, where the operating wavelength is λ = 1.55 μ m. The dips in the |S z | distributions close to the boundary between air and the Si core correspond to the negative Poynting vector regions discussed in [34]. Cleary, from Fig.  5(a), the light energy in the quasi-TE mode is confined within close vicinity of the Si core. In contrast, the distribution for the quasi-TM mode spreads over a larger area, especially in the SiO 2 area, and the energy leaks radially. The curved waveguides with R > 2 μ m exhibit similar tendencies.  Finally, we examine the hybridness represented by the amplitude rates of the maximum nondominant magnetic fields to the maximum dominant magnetic fields, to demonstrate the necessity of the full-vector formulation in the analysis of photonic wire bends. The red and green lines in Fig. 6 show the hybridness for the quasi-TE mode and for the quasi-TM mode, respectively, as a function of radius R at λ = 1.55 μ m. All the values of the hybridness are higher than 0.2 in the calculated range of R. In addition, there is a tendency for the hybridness to increase as R decreases, particularly in the case of the quasi-TM mode. Such high values of non-dominant field amplitudes indicate that the scalar and semivector approximations are no longer valid in the analysis of the silicon wire waveguides.

Summary
This paper has presented a full-vectorial FEM in a local CCS for efficient and accurate analysis of photonic wire bends. The formulation is based on hybrid edge/nodal elements, to solve the vector wave equation in a precise and straightforward fashion with the three components of electric or magnetic fields. To properly truncate the solution region, the FEM was implemented with the anisotropic PML in the CCS. Using this method, we have evaluated bending losses in silicon wire waveguides. The results clarify the radiation losses of the leaky modes, which are the theoretical minimum values of the bending losses. The comparison with experimental data supports the validity of the present approach. From the comparison with the losses in ESWs, it has also been demonstrated that the ESW approximation is no longer valid for the photonic wire bends, even in the loss calculation, not only in the calculation of effective indices.