## Abstract

We present a direct, rigorous, and fast numerical method for obtaining leaky-mode losses in optical fibers by solely solving complex propagation constants of the characteristic equation of leaky modes. Both the modified Bessel function and the Hankel function of the second kind are individually used to express the field component of leaky modes in the outermost cladding. The characteristic equation of cylindrically symmetric fiber structures, which consist of uniform and graded layers, is derived by combining the Runge–Kutta method and the exact solution of a homogeneous layer. Since complex root searching is the key technique in this method, we also present a numerical algorithm for solving the characteristic equation of optical fibers. Moreover, because for both guided and leaky modes the field distributions in the outermost cladding region have the same expression, the leaky mode can be easily obtained by choosing an improper solution, and therefore the calculation of leaky modes demonstrates the simplicity of this method. An approximation rule of branch choices for lossy material is also derived. The approach we present is consistent with the results of previously published papers.

© 2008 Optical Society of America

## 1. INTRODUCTION

The analysis of leaky modes has been applied for designing various photonic devices such as bent waveguides [1], polarizers [2], and depressed inner cladding (DIC) single-mode fibers [3, 4, 5]. Theoretically, a leaky mode is a bound mode below the cutoff frequency and should radiate energy out of the waveguide. The leaky mode has also provided an excellent approximation to radiation losses, which is thoroughly described by Snyder and Mitchell [6]. Since leaky modes are important in optimizing the performance of those specific elements, a simple and rigorous numerical method is required for analyzing leaky modes.

In an optical fiber transmission system, the most typical application of leaky-mode analysis was designing DIC single-mode fibers [7, 8, 9, 10, 11, 12, 13, 14], which can minimize dispersion effects at a desired wavelength in the $1.3\u20131.6\text{\hspace{0.17em}}\mu \mathrm{m}$ transmission window and exhibit leaky losses in the long wavelength region. In designing such DIC fibers, the perturbation theory was frequently used to analyze their leaky modes [7, 9, 10]. Kawakami and Nishida [7] used the perturbation theory to obtain an approximate formula for the leaky-mode attenuation coefficient and then applied this formula in the design of the DIC single-mode fibers. In their derivation, the inner cladding was assumed to be sufficiently thick. To find the leaky loss of the DIC fiber with a thin inner cladding layer, Maeda and Yamada [8] assumed that the waveguide structure was weakly guiding and then used the Poynting vector theorem to achieve a simplified expression of the attenuation coefficient. Cohen *et al.* [9] then proposed an approximate analytical formula, which was obtained by using a zero-order approximation, to explore the leaky losses of the fundamental mode $\left({\mathrm{LP}}_{01}\right)$ in the wide DIC fibers. Renner [10] also used the perturbation theory to derive a formula for approximating the attenuation coefficient for the ${\mathrm{LP}}_{11}$ mode of a coated DIC fiber. The above formulas can be simply applied to approach leaky loss in DIC fibers but become rather cumbersome when extended to multilayer structures with graded index profiles [11].

The leaky-mode attenuation coefficient can also be extracted using the full width at half-maximum (FWHM) of the Lorentzian peak. Shenoy *et al.* [11] calculated the leaky-mode attenuation coefficient by using the matrix approach, which is developed for analyzing the planar optical waveguide and the FWHM of the Lorentzian peak. Shenoy *et al.* transformed the cylindrically symmetric refractive index profiles into equivalent planar profiles and then used the matrix approach to determine the leakage loss of optical fibers with arbitrary refractive index profiles. Their method adopted the equivalent step index technique to approach the graded index profiles. Moreover, Pal and Priye [12] extended the technique of Shenoy *et al.* to analyze the leaky loss of a practical DIC fiber, which displayed a layered structure in the inner cladding and an axial dip in the core. The leaky loss spectrum of the ${\mathrm{LP}}_{11}$ mode of a coated DIC fiber was also estimated by using the matrix approach [13]. On the other hand, Thyagarajan *et al.* proposed a simple numerical technique for calculating the leaky loss of various types of optical fibers [14]. By using scalar multilayer approximation to deal with the graded index profile, Thyagarajan *et al.* applied the matrix method to solve the multilayer problem and then obtained the loss coefficient using the FWHM of the Lorentzian peak.

Recently, the analysis of loss or gain properties of optical fibers with arbitrary complex index profiles has attracted growing interest [15, 16]. Singh *et al.* [15] developed a rapidly converging numerical procedure for accurately evaluating the complex propagation constants in such fibers. Moreover, Qian and Boucouvalas [16] also presented a simple direct numerical procedure, which was based on the transmission line principles. This procedure was useful for designing and predicting high gain (or loss) in fibers. However, neither Singh *et al.* [15] nor Qian and Boucouvalas [16] applied their methods to solve leaky-mode problems.

In addition, studies of leaky modes based on the planar waveguides are widely discussed [17, 18, 19, 20]. Lee *et al.* [17] employed an explicit formula for the normalized radiation mode to obtain the function of leaky modes without resorting to power integration that can simultaneously avoid normalization and orthogonality problems of leaky modes. Lin *et al.* [18] proposed a simple and intuitive approach to analyze light transmission behavior in the planar waveguide with a deformed leaky region in its midsection. Moreover, Anemogiannis *et al.* [19] presented two new nonrigorous but easy and accurate numerical methods to determine the attenuation coefficient of leaky modes in planar multilayer waveguides. Their attenuation coefficient was extracted from the FWHM of the Lorentzian peak of the reflection coefficient or the density of wave vectors of the structure. Additionally, Petracek and Singh [20] provided a simple and exact technique for finding leaky-mode solutions of the dispersion equation, which was given by using the thin film transfer matrix method.

More recently, leaky modes of microstructured optical fibers are of interest in optical fiber communication [21]. Basically, the effective geometry of microstructured optical fibers is very similar to DIC fibers. In [21], Issa and Poladian presented the adjustable boundary condition Fourier decomposition method to calculate leaky modes in arbitrary microstructured optical fibers. They devised a refinement technique to determine the radiating field outside the boundary.

Although the leaky losses of optical fibers can be obtained from the above numerical methods, their evaluation depends on various assumptions or approximations. The perturbation method gives better results only when there is a small difference between the refractive indices in each layer of the fiber. The attenuation coefficient can approach accurate results when the roots are searched in the complex plane or when the Lorentzian peaks do not overlap with each other. In other words, the FWHM of the Lorentzian peaks are less accurate in situations involving large loss or almost equal effective indices [19]. Besides, these techniques are complicated since the field components in the outermost cladding are expressed in terms of the Hankel function ${H}_{\nu}$.

In this paper, we present a direct and rigorous method for analyzing the leaky-mode loss of optical fibers without any assumptions or additional numerical works. Our method starts with solving the complex propagation constant of a leaky mode of optical fibers. Since the imaginary part of the propagation constant is the attenuation constant due to the leakage of the leaky-wave energy into the substrate and the superstrate regions, the attenuation constant of the leaky mode represents the propagation loss along the longitudinal direction, i.e., the leaky-mode loss of fibers [17, 22]. Two expressions are chosen to express the field components of leaky modes in the outermost cladding: (1) the second kind of Hankel function ${H}_{\nu}^{\left(2\right)}$ and (2) the modified Bessel function ${K}_{\nu}$. The second kind of Hankel function ${H}_{\nu}^{\left(2\right)}$ is commonly used to express an outgoing wave in the outermost region. It is applied to the perturbation method [10] and the FWHM of the Lorentzian peaks [14] to evaluate the leaky-mode loss. In this paper, the leaky mode is solved based solely on an outgoing wave expressed as the Hankel function ${H}_{\nu}^{\left(2\right)}$. Moreover, note that the Hankel function is the linear recombination of the modified Bessel function ${K}_{\nu}$. The field components of leaky modes in the outermost cladding can also be expressed as the modified Bessel function ${K}_{\nu}$. Since the Hankel function is the linear recombination of the first and second kinds of Bessel functions (${J}_{\nu}$ and ${Y}_{\nu}$), using the modified Bessel function ${K}_{\nu}$ to calculate leaky modes takes less computation time than that of the Hankel function. Because complex root searching plays a very important role in this paper, a numerical algorithm for solving the characteristic equation of optical fibers is also demonstrated.

The advantage of this method is that the proposed technique can obtain the leaky-mode loss without any assumptions or additional work. Our method is so general that even the leaky-mode loss for lossy material can be easily calculated. Note that numerous methods (including the perturbation method and the FWHM of the Lorentzian peaks) evaluate the leaky-mode loss only for the lossless material. In addition, the characteristic equation of cylindrically symmetric fiber structures, which consist of uniform and graded layers, is derived by combining the Runge–Kutta method and the exact solution of a homogeneous layer. On the other hand, our method is simple because we use the same expressions (${K}_{\nu}$ or ${H}_{\nu}^{\left(2\right)}$) in the outermost cladding layer for guided and leaky modes. By using proper branch choices, the guided and leaky modes can easily be obtained without changing the expression of the outermost cladding. A conventional rule for choosing branch cuts in planar waveguides [22, 23, 24] can be applied to calculate the leaky modes of optical fibers. Moreover, we derive an approximation rule for choosing branch cuts. This rule gives a guideline for obtaining proper and improper solutions. To check the degree of accuracy, we compare the leaky losses of DIC fibers obtained from the proposed technique with the numerical results in [9, 13]. The power loss of a six-layer dispersion flattened fiber proposed by Yokoyama *et al.* [25] is also calculated.

The remainder of this paper is organized as follows. Section 2 presents a derivation of the characteristic equation. The branch choices for the transverse wave vector are then discussed in Section 3. In Section 4, numerical results of the leaky losses for the typical DIC fiber are presented, and Section 4 also makes comparisons with other published results. In Section 5, we propose an approximation rule of branch choices for lossy material. A numerical algorithm for calculating the propagation constant of leaky modes is shown in Section 6. Finally, conclusions are given in Section 7.

## 2. PROBLEM FORMULATION

We begin by considering a cylinder optical waveguide structure with multiple layers (total *n* layers) consisting of graded index profiles and step index profiles. A circularly cylindrical coordinate system *r*, *ϕ*, *z* is used, where *r* is the radial, *ϕ* is the azimuthal, and *z* is the axial coordinate. All field quantities are assumed to vary sinusoidally in time with circular frequency *ω*, and the field is assumed to propagate in the *z* direction. The complex representation of field quantities is used, and the complex time factor, which is $\mathrm{exp}\left(j\omega t\right)$, is omitted throughout the study. Assume that the $i\text{th}$ layer is homogeneous with the index. The wave equation can be written as follows:

*z*components of the electric field and the magnetic field, respectively, ${n}_{i}$ represents the refractive index of the $i\text{th}$ layer, ${k}_{0}=2\pi \u2215\lambda $ is the free space wavenumber, and

*ν*is the azimuthal mode number. ${k}_{z}(=\beta -j\alpha )$ is a complex propagation constant, where the real part of ${k}_{z}$,

*β*, is the phase constant, and the imaginary part of ${k}_{z}$,

*α*, is the attenuation constant due to the leakage of leaky modes. Since the effective index is defined by ${n}_{\mathrm{eff}}=\beta \u2215{k}_{0}$, the solution of the differential equation (1) becomes

From ${E}_{z}$ and ${H}_{z}$ for the homogeneous layer, the electromagnetic tangential field distribution for a mode of azimuthal order *ν* can be written as

Assume that the index profile of the $i\text{th}$ layer is graded. According to the Maxwell equation, the relationship of the electromagnetic tangential field components can then be expressed as [26]

*i*represents the $i\text{th}$ layer. Equation (7) can be numerically solved by the Runge–Kutta method [23, 26] and can be written as

In solving the electromagnetic problem, the boundary condition should be applied to each interface of the structure. There are four types of interfaces for fiber structures. Figure 1 illustrates an example of an optical fiber consisting of these four types of interfaces, which has been studied to apply in dispersion flattened fibers [25]. The four types of interfaces of Fig. 1 are the step–step interface between layers 1 and 2, the step-graded interface between layers 2 and 3, the graded–graded interface between layers 3 and 4, and the graded-step interface between layers 4 and 5. Since the tangential field components should be continuous across the interface of layers, by matching the boundary conditions at each layer interface we can obtain the characteristic equation of an optical fiber,

where $\mathbf{W}\left({k}_{z}\right)$ is a $4\times 4$ matrix and**V**is a constant vector. The system of linear equations has a nontrivial solution when the determinant of matrix $\mathbf{W}\left({k}_{z}\right)$ equals zero. The solutions of the determinant equation,give the propagation constants ${k}_{z}$ for the modes of the structure. The propagation constants ${k}_{z}$ can be solved using Muller’s method [27]. The field distributions of cylindrically symmetric index profiles can subsequently be evaluated.

## 3. BRANCH CHOICES FOR TRANSVERSE WAVE VECTORS

Because the leaky modes radiate energy out from the waveguide, a conventional expression of an outgoing wave in the outermost cladding is expressed in terms of the Hankel function of the second kind ${H}_{\nu}^{\left(2\right)}$ [8, 10],

*r*. Therefore, the negative sign of Eq. (12) is chosen so that the transmission mode corresponds to the guided mode, while the positive sign of Eq. (12) is chosen so that the transmission mode corresponds to the leaky mode.

Since ${H}_{\nu}^{\left(2\right)}$ is a linear recombination of the modified Bessel function ${K}_{\nu}$, in this paper we also use ${K}_{\nu}$ to express the outermost cladding field distribution function. The relationship between ${K}_{\nu}$ and ${H}_{\nu}^{\left(2\right)}$ is given in Appendix A. Using this expression, the field distributions in the outermost cladding region for both guided and leaky modes can be written as

*r*. We define that the positive sign of Eq. (14) corresponds to the positive real part of ${k}_{cn}$ $(\mathrm{Re}\left({k}_{cn}\right)\u2a7e0)$ and represents a proper solution, while the negative sign of Eq. (14) refers to the negative real part of ${k}_{cn}$ $(\mathrm{Re}\left({k}_{cn}\right)\u2a7d0)$ and represents an improper solution.

Note that choosing the correct signs in Eqs. (12, 14) is critical in obtaining the correct complex roots of the determinant (10). The sign choice determines whether the field distribution of the mode will be proper or improper. Since guided modes in the outermost cladding are expressed in terms of the modified Bessel function ${K}_{\nu}$, we focus on introducing a conventional rule in choosing branch cuts for ${K}_{\nu}$ in this paper. The same rule can be applied to selecting branch choices for ${H}_{\nu}^{\left(2\right)}$.

Consider that the field component in the outermost cladding is the modified Bessel function ${K}_{\nu}$. In lossless material, proper or improper solutions correspond to the location of the effective index. If the effective index of a mode is greater than ${n}_{n}$, the calculated propagation constant represents a proper mode $({k}_{cn}=+\sqrt{{k}_{z}^{2}-{k}_{0}^{2}{n}_{n}^{2}})$, whereas the solution refers to an improper mode $({k}_{cn}=-\sqrt{{k}_{z}^{2}-{k}_{0}^{2}{n}_{n}^{2}})$ if ${n}_{\mathrm{eff}}$ is less than ${n}_{n}$. For guided modes, the real part of ${k}_{cn}$ is positive:

In lossy material, the branch choice is more complicated than that of lossless material. Note that the indices of each layer are expressed as complex numbers. The conventional rule of branch choices requires that the loci of ${k}_{z}$ in the complex plane must be traced continuously to produce smooth curves as a function of a defined physical parameter. Therefore, the rules of the branch cut for lossy material depend on the loci of ${k}_{z}$. Figure 2 shows the loci of proper and improper solutions. We notice that both improper and proper attenuation constants have identical values at ${\omega}_{0}$. That means $\omega ={\omega}_{0}$ corresponds to a double root of the propagation constant. Because the value of the propagation constant should vary continuously with a chosen variable *ω*, the propagation constant ${k}_{z}$ should change from a guided mode to a leaky mode at $\omega ={\omega}_{0}$ as *ω* varies from $\omega >{\omega}_{0}$ to $\omega <{\omega}_{0}$.

When the complex propagation constants ${k}_{z}$ of leaky modes are obtained, for the transmission distance *L* the power losses of modes can be easily calculated from

## 4. RESULTS

Consider the DIC fiber structure analyzed by Cohen *et al.* [9]. The inset of Fig. 3 shows the refractive index profile. The refractive index ${n}_{1}$ of the core is slightly larger than the index ${n}_{3}$ of the outer cladding, while the refractive index ${n}_{2}$ of the inner cladding is smaller than that of the outer cladding. The radii of the core and inner cladding are indicated by *a* and *b*, respectively. The fiber parameters used in the calculations are $a=3.75\text{\hspace{0.17em}}\mu \mathrm{m}$, ${n}_{1}={n}_{3}(1+\Delta )$, ${n}_{2}={n}_{3}(1+{\Delta}^{\prime})$, and $\Delta -{\Delta}^{\prime}=0.5\%$. *Δ* and ${\Delta}^{\prime}$ are the core-outer cladding index difference and inner–outer cladding index difference, respectively, which are expressed by $\Delta =({n}_{1}-{n}_{3})\u2215{n}_{3}$ and ${\Delta}^{\prime}=({n}_{2}-{n}_{3})\u2215{n}_{3}$. Figure 3 illustrates the leaky losses of the fundamental mode as a function of the wavelength for different values of *Δ* and $b\u2215a$. Figure 3a shows that the curves are applied to different values of $\Delta =0\%$, 0.2%, 0.23%, 0.25%, and 0.27% for $b\u2215a=6$ and 7. Figure 3b demonstrates that the curves are applied to different values of $b\u2215a$
$(=5.5\u20138)$ for $\Delta =0$. The solid curves correspond to the results obtained by using the approximate analytical formula of [9], while the dashed curves correspond to the present method. Clearly, the two results are in good agreement with the present method and the formula of [9]. The maximum relative error is less than 3%.

Consider a lossy optical fiber with a material loss of $0.176\text{\hspace{0.17em}}\mathrm{dB}\u2215\mathrm{km}$ in each layer of the fiber. Figure 4 shows the attenuation as a function of wavelength for a DIC optical fiber. The radii of the core and inner cladding are $a=3.75\text{\hspace{0.17em}}\mu \mathrm{m}$ and $b=22.5\text{\hspace{0.17em}}\mu \mathrm{m}$, respectively. As shown in Fig. 4, a turning point for the proper and improper solutions occurs at a wavelength of $1.4071\text{\hspace{0.17em}}\mu \mathrm{m}$. The transverse wave vector ${k}_{cn}$ must be chosen properly at this turning point so that the locus of propagation constant ${k}_{z}$ in the complex plane can be traced continuously. The power loss is $\sim 0.18\text{\hspace{0.17em}}\mathrm{dB}\u2215\mathrm{km}$ when the wavelength is less than $1.4071\text{\hspace{0.17em}}\mu \mathrm{m}$. Note that the loss increases rapidly when the wavelength is larger than $1.4071\text{\hspace{0.17em}}\mu \mathrm{m}$, consistent with the property of the leaky mode.

In addition, this work calculates the leakage losses for the ${\mathrm{HE}}_{21}$ mode in a coated DIC fiber, as discussed in [13], where the parameters of each layer are core radius $a=4.15\text{\hspace{0.17em}}\mu \mathrm{m}$, inner cladding radius $b=24.9\text{\hspace{0.17em}}\mu \mathrm{m}$, and outer cladding radius $c=62.5\text{\hspace{0.17em}}\mu \mathrm{m}$. Assume that ${n}_{1}$, ${n}_{2}$, ${n}_{3}$, and ${n}_{4}$ are the refractive indices of the core, the inner cladding, the outer cladding, and the coating, respectively. The index of each layer is ${n}_{3}=1.445$, ${n}_{1}-{n}_{2}=5.344\times {10}^{-3}$, ${n}_{3}-{n}_{2}=1.732\times {10}^{-3}$, and ${n}_{4}=1.50$. By using the proposed method, the power loss (17) of the ${\mathrm{HE}}_{21}$ mode as a function of the wavelength is plotted as a solid curve in Fig. 5 , and the dots correspond to the leaky loss of the ${\mathrm{LP}}_{11}$ mode by using the matrix method [13]. Table 1 shows the results obtained from [13] and our method. The effective index ${n}_{\mathrm{eff}}$ and the power loss ${P}_{\text{loss}}$ of the ${\mathrm{HE}}_{21}$ mode obtained by using our method are tabulated in columns 3 and 5 of Table 1, respectively. The corresponding effective index ${n}_{e}^{m}$ and leaky loss ${\alpha}_{L}^{m}$ obtained by the matrix method of the ${\mathrm{LP}}_{11}$ are listed in columns 2 and 4, respectively (the superscript *m* indicates the matrix method). Compared with the result listed in column 6 of Table 1, the relative error $(=\mid {\alpha}_{L}^{m}-{P}_{\text{loss}}\mid \u2215{P}_{\text{loss}})$ is less than 10%.

Next, we consider a multiple-cladding structure discussed in [25], in which the refractive index profile is shown as Fig. 1. The radii of each layer are indicated by ${r}_{1}$, ${r}_{2}$, ${r}_{3}$, ${r}_{4}$, and ${r}_{5}$. Their values used in the calculations are 1.3754, 5.98, 8.74, 11.5, and $23\text{\hspace{0.17em}}\mu \mathrm{m}$, respectively. The refractive indices ${n}_{1}$, ${n}_{2}$, ${n}_{3}$, and ${n}_{4}$ are 1.4689, 1.4593, 1.4651, and 1.4593, respectively. The relative refractive index between the center core and outermost cladding is defined by $\Delta =({n}_{1}-{n}_{5})\u2215{n}_{5}$. Figure 6 shows the leaky losses of the fundamental mode as a function of the wavelength for different values of $\Delta =0\%$, 0.52%, 0.53%, 0.535%, and 0.54%. As shown in Fig. 6, for these different values of *Δ*, the cutoff wavelengths of the fundamental mode correspond to 1.337, 1.463, 1.534, and $1.608\text{\hspace{0.17em}}\mu \mathrm{m}$, respectively. Consider the material loss of $0.176\text{\hspace{0.17em}}\mathrm{dB}\u2215\mathrm{km}$ in each layer of the fiber. Figure 7 shows the attenuation as a function of wavelength for this lossy fiber with $\Delta =0.53\%$. As shown in Fig. 7, a turning point between the guided mode and the leaky mode is found near the wavelength of $1.4607\text{\hspace{0.17em}}\mu \mathrm{m}$. This means that the turning point of the guided mode and the leaky mode is at the wavelength of $1.4607\text{\hspace{0.17em}}\mu \mathrm{m}$. Similar to Fig. 6, the attenuation of the guided modes is $\sim 0.18\text{\hspace{0.17em}}\mathrm{dB}\u2215\mathrm{km}$, where it is close to the material loss. The attenuation is steeply increased when the wavelength is larger than the turning point $(\lambda =1.4607\text{\hspace{0.17em}}\mu \mathrm{m})$. At the wavelength of $1.55\text{\hspace{0.17em}}\mu \mathrm{m}$, the power loss is $\sim \mathrm{10,000}$ times greater than that in Fig. 4.

To compare the computation time of using ${K}_{\nu}$ with that of using ${H}_{\nu}^{\left(2\right)}$, a three-layer DIC fiber with $\Delta =0.23\%$ for $b\u2215a=6$ and a six-layer dispersion flattened fiber with $\Delta =0.53\%$ are considered. The 100 discrete wavelengths between 1.5 and $1.65\text{\hspace{0.17em}}\mu \mathrm{m}$ are calculated. A personal computer with an $800\text{\hspace{0.17em}}\mathrm{MHz}$ Pentium III processor is used to calculate the leaky losses of these cases. In a lossy three-layer DIC fiber, the computing time is $0.7\text{\hspace{0.17em}}\mathrm{s}$ when ${K}_{\nu}$ is used. By using ${H}_{\nu}^{\left(2\right)}$, on the other hand, the computing time is $0.95\text{\hspace{0.17em}}\mathrm{s}$. The results show that ${K}_{\nu}$ saves $\sim 26\%$ of the computation time. Similarly, for a lossy six-layer dispersion flattened fiber, the computing time of ${K}_{\nu}$ is $1.52\text{\hspace{0.17em}}\mathrm{s}$, while the computing time of ${H}_{\nu}^{\left(2\right)}$ is $1.63\text{\hspace{0.17em}}\mathrm{s}$. The calculation of the transfer matrices of a six-layer dispersion flattened fiber takes more time than that of a three-layer DIC fiber. The computing time can save only 7% in this case.

## 5. APPROXIMATION RULE OF BRANCH CHOICES FOR LOSSY MATERIAL

Consider that the modified Bessel function ${K}_{\nu}$ is used to express the field components in the outermost cladding. Assume that the complex index in the outermost cladding is ${n}_{n}={n}_{nr}-j{n}_{ni}$. In general, due to $\beta \u2aa2\alpha $ and ${n}_{nr}\u2aa2{n}_{ni}$, ${k}_{cn}^{2}$ can be expressed as

The sign of the transverse wavenumber ${k}_{cn}$ of ${H}_{\nu}^{\left(2\right)}$ is opposite from that of ${K}_{\nu}$. The approximation rule of branch choice for ${H}_{\nu}^{\left(2\right)}$ is similar to that of ${K}_{\nu}$, but with a different sign. Consider the field components in the outermost cladding in terms of the Hankel function ${H}_{\nu}^{\left(2\right)}$. The transverse wavenumber square ${k}_{cn}^{2}$ can be defined as

## 6. ALGORITHM FOR SEARCHING SOLUTIONS OF THE CHARACTERISTIC EQUATION

In this paper, we use Muller’s method to solve the determinant equation (10) of an optical fiber. Three initial guesses are required by using Muller’s method, since the three initial guesses are highly sensitive for complex root searching. Here we present an algorithm to determine the initial guesses.

First, we choose an initial interval $({n}_{l},{n}_{h})$ as a region where roots may exist; ${n}_{l}$ and ${n}_{h}$ are the low and high bound intervals of the effective indices, respectively. Then we evaluate the characteristic function associated with this region. When an effective index approaches the root $(\mathrm{det}\left[W\left({k}_{z}\right)\right]=\mathrm{det}\left[W({n}_{\mathrm{eff}}\cdot {k}_{0})\right]\approx 0)$, the absolute value of the function drops. We can use those effective indices around the dips as initial guesses. Since Muller’s method is a stable and fast numerical method, we can obtain the solution within 15 iteration times with a relative error less than ${10}^{-8}$.

Consider the case of Fig. 3. Assume that ${n}_{1}$, ${n}_{2}$, and ${n}_{3}$ are 1.46, 1.4527, and 1.46, respectively. Figure 8 shows the absolute value of the determinant equation (10) as a function of the effective index. The radii of core and inner cladding are 3.75 and $22.5\text{\hspace{0.17em}}\mu \mathrm{m}$, respectively. The modified Bessel function ${K}_{\nu}$ is applied to express the function of leaky modes. The effective index ranges from 1.4527 to 1.46, which is divided by 1000 points. As shown in Fig. 8, a dip occurs around the effective index of 1.45636022. We use 1.45635292, 1.45636022, and 1.45636752 as three initial guesses. After three iteration times, the solution of the effective index $({n}_{\mathrm{eff}}=1.45635951-j2.784186\times {10}^{-10})$ is obtained, while the relative error is less than ${10}^{-8}$.

Consider a lossy optical fiber in Fig. 7. The radii of ${r}_{1}$, ${r}_{2}$, ${r}_{3}$, ${r}_{4}$, and ${r}_{5}$ are 1.3754, 5.98, 8.74, 11.5, and $23\text{\hspace{0.17em}}\mu \mathrm{m}$, respectively. The refractive index of each layer is ${n}_{1}=1.4689-j5.0\times {10}^{-12}$, ${n}_{2}=1.4593-j5.0\times {10}^{-12}$, ${n}_{3}=1.4651-j5.0\times {10}^{-12}$, ${n}_{4}=1.4593-j5.0\times {10}^{-12}$, and ${n}_{5}=1.4689-j5.0\times {10}^{-12}$. Figure 9 shows the absolute value of the determinant equation (10) as a function of the effective index. The effective index ranges from 1.4593 to 1.4689 and is divided by 1000 points. As shown in Fig. 9, two dips occur around the effective indices of 1.461120116 and 1.459896304. Three initial guesses of 1.461110516, 1.461120116, and 1.461129716 can obtain the solution of $1.461122504-j1.913361764\times {10}^{-7}$ in seven iteration times with a relative error less than ${10}^{-8}$. This solution refers to the effective index of the ${\mathrm{HE}}_{11}$ mode. By using the same procedure, the solution of the second dip is $1.45989599888-j5.870866346\times {10}^{-6}$ in four iteration times, which refers to the effective index of the ${\mathrm{EH}}_{11}$ mode.

## 7. CONCLUSION

In this paper, we present a direct and rigorous numerical method for calculating leaky-mode losses in optical fibers. The modified Bessel function ${K}_{\nu}$ and the Hankel function of the second kind ${H}_{\nu}^{\left(2\right)}$ are used to express the outgoing expression in the outermost cladding. By matching the boundary condition, the complex propagation constants are evaluated, and the leaky-mode loss can be obtained from the imaginary part of the propagation constant (*i.e.*, the attenuation constant). A numerical algorithm for solving the characteristic equation of optical fibers is also presented. By properly choosing branch cuts, we can easily obtain the guided and the leaky modes of optical fiber. Additionally, an approximation rule for choosing branch cuts is developed for solving improper solutions, which correspond to leaky modes.

This method simplifies the procedure used for solving leaky modes and reduces the computation time, since the process used for branch selection is similar to that used in planar dielectric waveguides. As a result, the analysis developed in this paper agrees with the results of the previously published DIC optical fiber analyzed by perturbation theory and the matrix approach. Moreover, since ${K}_{\nu}$ is the linear recombination of ${H}_{\nu}^{\left(2\right)}$, the modified Bessel function ${K}_{\nu}$ saves the computation time when compared with the Hankel function ${H}_{\nu}^{\left(2\right)}$. This method can be applied to fibers with cylindrically symmetric structures including lossy material.

## APPENDIX A

The relationship between the modified Bessel function ${K}_{\nu}\left({k}_{cn}r\right)$ and the Hankel function of the second kind ${H}_{\nu}^{\left(2\right)}$ at an angle of ${k}_{cn}r$ between $-\pi \u22152$ and *π* is expressed as [28]

*A*is a coefficient,

*j*is $\sqrt{-1}$, and ${k}_{cn}$ is defined by Eq. (14). The field distribution of leaky modes in the outermost cladding is expressed in terms of the Hankel function of the second kind ${H}_{\nu}^{\left(2\right)}$, according to Eq. (11),where the angle of $-j{k}_{cn}$ is between 0 and $\pi \u22154$. Therefore, when the modified Bessel function ${K}_{\nu}\left({k}_{cn}r\right)$ is used to express the outgoing wave, the angle of ${k}_{cn}$ in Eq. (19) is set to $\pi \u22152\sim 3\pi \u22154$.

## ACKNOWLEDGMENT

The authors thank I-Shou University, Taiwan, for financially supporting this research under contract ISU96-01-02.

**1. **K. Thyagarajan, M. R. Shenoy, and A. K. Ghatak, “Accurate numerical method for the calculation of bending loss in optical waveguides using a matrix approach,” Opt. Lett. **12**, 296–298 (1987). [CrossRef]

**2. **K. Thyagarajan, S. Diggavi, and A. K. Ghatak, “Analytical investigation of leaky and absorbing planar structures,” Opt. Quantum Electron. **19**, 131–137 (1987). [CrossRef]

**3. **S. Kawakami and S. Nishida, “Characteristics of a doubly clad optical fiber with a low-index inner cladding,” IEEE J. Quantum Electron. **QE-10**, 879–887 (1974). [CrossRef]

**4. **Y. Akasaka, R. Sugizaki, S. Arai, Y. Suzuki, and T. Kamiya, “Dispersion flat compensation fiber for dispersion shifted fiber,” in 22nd European Conference on Optical Communication (IEEE, 1996), pp. 221–224.

**5. **T. Tsuda, Y. Akasaka, S. Sentsui, K. Aiso, Y. Suzuki, and T. Kamiya, “Broad band dispersion slope compensation of dispersion shifted fiber using negative slope fiber,” in *24th European Conference on Optical Communication* (IEEE, 1998), pp. 233–234.

**6. **A. W. Snyder and D. J. Mitchell, “Leaky mode analysis of circular optical waveguides,” Opto-electronics (London) **6**, 287–296 (1974). [CrossRef]

**7. **S. Kawakami and S. Nishida, “Perturbation theory of a doubly clad optical fiber with a low-index inner cladding,” IEEE J. Quantum Electron. **QE-11**, 131–138 (1975).

**8. **M. Maeda and S. Yamada, “Leaky modes on W-fibers: mode structure and attenuation,” Appl. Opt. **16**, 2198–2203 (1977). [CrossRef]

**9. **L. G. Cohen, D. Marcuse, and W. L. Mammel, “Radiating leaky-mode losses in single-mode lightguides with depressed-index claddings,” IEEE J. Quantum Electron. **QE-18**, 1467–1472 (1982). [CrossRef]

**10. **H. Renner, “Leaky-mode loss in coated depressed-cladding fibers,” IEEE Photon. Technol. Lett. **3**, 31–32 (1991). [CrossRef]

**11. **M. R. Shenoy, K. Thyagarajan, and A. K. Ghatak, “Numerical analysis of optical fibers using matrix approach,” J. Lightwave Technol. **6**, 1285–1291 (1988). [CrossRef]

**12. **B. P. Pal and V. Priye, “The effect of an axial dip and ripples in the inner cladding on the leakage loss of ${\mathrm{LP}}_{01}$ mode in depressed index clad fibre,” IEE Proc.: Optoelectron. **137**, 311–314 (1990). [CrossRef]

**13. **B. P. Pal, R. L. Gallawa, and I. C. Goyal, “${\mathrm{LP}}_{11}$-mode leakage loss in coated depressed clad fibers,” IEEE Photon. Technol. Lett. **4**, 376–378 (1992). [CrossRef]

**14. **K. Thyagarajan, S. Diggavi, A. Taneja, and A. K. Ghatak, “Simple numerical technique for the analysis of cylindrically symmetric refractive-index profile optical fibers,” Appl. Opt. **30**, 3877–3879 (1991). [CrossRef]

**15. **R. Singh and E. K. Sharma, “Propagation characteristics of single-mode optical fibers with arbitrary complex index profiles: a direct numerical approach,” IEEE J. Quantum Electron. **37**, 635–640 (2001). [CrossRef]

**16. **X. Qian and A. C. Boucouvalas, “Propagation characteristics of single-mode optical fibers with arbitrary complex index profiles,” IEEE J. Quantum Electron. **40**, 771–777 (2004). [CrossRef]

**17. **S. L. Lee, Y. Chung, L. A. Coldren, and N. Dagli, “On leaky mode approximations for modal expansion in multilayer open waveguides,” IEEE J. Quantum Electron. **31**, 1790–1802 (1995). [CrossRef]

**18. **Y. Z. Lin, J. H. Zhan, and S. M. Tseng, “A new method of analyzing the light transmission in leaky and absorbing planar waveguides,” IEEE Photon. Technol. Lett. **9**, 1241–1243 (1997). [CrossRef]

**19. **E. Anemogiannis, E. N. Glytsis, and T. K. Gaylord, “Determination of guided and leaky modes in lossless and lossy planar multilayer optical waveguides: reflection pole method and wavevector density method,” J. Lightwave Technol. **17**, 929–941 (1999). [CrossRef]

**20. **J. Petracek and K. Singh, “Determination of leaky modes in planar multilayer waveguides,” IEEE Photon. Technol. Lett. **14**, 810–812 (2002). [CrossRef]

**21. **N. A. Issa and L. Poladian, “Vector wave expansion method for leaky modes of microstructured optical fibers,” J. Lightwave Technol. **21**, 1005–1012 (2003). [CrossRef]

**22. **N. H. Sun, C. C. Chou, H. W. Chang, J. K. Butler, and G. A. Evans, “Radiation loss of grating-assisted directional couplers using the Floquet-Bloch theory,” J. Lightwave Technol. **24**, 2409–2415 (2006). [CrossRef]

**23. **R. E. Collin and F. J. Zucker, eds., *Antenna Theory* (McGraw-Hill, 1969), p. 203.

**24. **N. H. Sun, J. K. Butler, G. A. Evans, L. Pang, and P. Congdon, “Analysis of grating-assisted directional couplers using the Floquet-Bloch theory,” J. Lightwave Technol. **15**, 2301–2315 (1997). [CrossRef]

**25. **Y. Yokoyama, T. Kato, M. Hirano, M. Onishi, E. Sasaoka, Y. Makio, and M. Nishimura, “Practically feasible dispersion flattened fibers produced by VCD technique,” in *24th European Conference on Optical Communication* (IEEE, 1998), pp. 131–132.

**26. **M. O. Vassell, “Calculation of propagating modes in a graded-index optical fibre,” Opto-electronics (London) **6**, 271–286 (1974). [CrossRef]

**27. **R. L. Burden and J. D. Faires, *Numerical Analysis* (PWS-KENT, 1984).

**28. **M. Abramowitz and I. A. Stegun, *Handbook of Mathematical Functions* (Dover, 1972).