Open Access
Add to BibSonomyAbstract
In this paper we present a full-vectorial finite-difference analysis of microstructured optical fibers. A new mode solver is described which uses Yee’s 2-D mesh and an index averaging technique. The modal characteristics are calculated for both conventional optical fibers and microstructured optical fibers. Comparison with previous finite difference mode solvers and other numerical methods is made and excellent agreement is achieved.
©2002 Optical Society of America
1. Introduction
Among optical waveguides, only a few structures such as slab waveguides and step-index optical fibers can be solved analytically. For more complex waveguide structures, rigorous numerical methods have been proposed including both finite element methods and finite difference (FD) methods [1,2].
Microstructured optical fibers (MOFs) have a much more complex transverse structure compared to that of conventional optical fibers and have attracted a great deal of recent interest both in theory and experiment. Much of this interest is centered around the additional degree of freedom offered by the size and pattern of the air holes. MOFs have exhibited interesting features such as wide range of single mode operation [3,4], easily controllable dispersion [5,6] and birefringence [7], and tailorable effective modal areas [8,9]. They have shown great promise in nonlinear fiber optics such as supercontinuum generation [10] and many other novel fiber devices [11]. Since the first experimental demonstration of an MOF (also called a photonic crystal fiber or holey fiber) in 1996 [3], many modeling techniques have been applied in their characterization, including approximate effective index models [4,12], plane-wave expansion methods [13–15], localized-function methods [16,17], finite element methods [18], finite difference time domain (FDTD) methods [19,20], and multipole methods [21]. In this paper, we use finite difference frequency domain (FDFD) techniques to analyze the modal properties of MOFs.
In the FDFD mode solvers, two discretization schemes have been used. One is that first proposed by Stern [22] in which possible discontinuities lie between two adjacent mesh grids and every grid point corresponds to a unique refractive index. The wave equation in terms of transverse electric field E t (or magnetic field H t )
or
where k 0=2π/λ is the wave number in free space, ε r is the waveguide dielectric constant, and β is the propagation constant, is directly discretized by finite difference [23]. This discretization scheme is usually used in the context of the beam propagation method [23,24], which can also be used as a mode solver. Another discretization scheme is that first proposed by Bierwirth et al. [25]. In this scheme, possible discontinuities lie on the mesh grids, so that any grid point can be associated with up to four different refractive indices. The transverse magnetic components are usually used in deriving the discretization matrix [26,27].
In this paper, we present a full-vector finite difference mode solver that is based on discretization scheme first proposed by Yee [28]. Yee’s mesh is widely used in the FDTD analysis [29]. Here we use Yee’s two-dimensional mesh in our frequency domain mode solver for complex optical waveguides. In Section 2, the formulation of the method is described in detail. In order to improve the staircase approximation for curved interface, an index averaging technique is used for the cells across interfaces. Numerical results for both conventional optical fibers and MOFs are presented in Section 3. The comparison between our results and those by other methods are also made. Section 4 concludes the paper.
2. Formulation
The Yee’s 2-D mesh is illustrated in Fig. 1(a). The mesh grids for electric fields lie on possible dielectric discontinuities. Since all the transverse field components are tangential to the unit cell boundaries, the continuity conditions are automatically satisfied.
Assume the fields have dependence of exp[i(βz-ωt)]. From Maxwell’s curl equations (∇×E = -∂B/∂t, ∇×H = ∂D/∂t), after scaling E by the free space impedance , we have
Discretizing Eqs. (2) and (3) yields
where
In writing Eqs. (6a)~(6c), we have approximated the refractive indices by averaging the refractive indices of adjacent cells.
Eqs. (4) and (5) can be written in matrix form as
where I is a square identity matrix, ε rx , ε ry and ε rz are diagonal matrices determined by Eqs. (6). U x , U y , V x and V y are square matrices and depend on the boundary conditions of the rectangular computation window. For example, when the zero-value boundary condition is chosen for the computation window edges, we have
From Eqs. (7) and (8), after some algebra, we obtain an eigenvalue equation in terms of transverse electric fields:
where
Alternatively, we can obtain an eigenvalue equation in terms of transverse magnetic fields:
with
Under proper boundary conditions (zero value boundary, for example), one can get V x =- and V y =- , where the superscript T denotes transpose operation. Then it is clear from Eqs. (11) and (13) that
Other FD mode solvers [23,27] obtain eigenvalue equations similar to Eq. (10) or (12), but the coefficient matrices P and Q are different. One advantage of the present approach is that, once E t or H t is obtained, all other field components are readily obtained through Eqs. (7) and (8).
If we ignore waveguide material absorption, i.e., ε r is real, then P and Q are real sparse matrices. The sparse matrices can be stored efficiently in computers to reduce the amount of memory required for data storage. Solving the eigenvalue equation (10) or (12) using available numerical routines provides us the effective modal index neff =β/k 0 and modal fields of the guided modes.
In the finite difference analysis of waveguides with curved interfaces, the staircase approximation has to be used in the rectangular mesh. In order to improve this approximation, we use averaged refractive indices for the mesh cells across the interface. Similar techniques have been previously used in the plane expansion method [15] and FDTD analysis [30]. For example, εrz (j,l) = faεa + (1-fa )εb is used instead of Eq. (6c) for the cell that is filled with two different dielectric media εa and εb , where fa is the fraction of the cell occupied by medium εa (see Fig. 1(b)). Similarly, the averaged εrx and εry can be calculated in the cells centered on Ex and Ey , respectively. We also apply this index averaging (IA) technique in other FD mode solvers mentioned in the Introduction. The use of the averaged refractive index of interfacial cells, as will be seen in the next section, can significantly accelerate convergence and improve the modeling accuracy for waveguides with curved interfaces, such as optical fibers and MOFs.
3. Numerical Results
In this section, we first use the proposed FD mode solver to calculate the fundamental effective index of a step-index optical fiber and compare them with other FD mode solvers [23,27] and the analytical solution. Then we apply the mode solver for the characterization of MOFs, including air-hole assisted optical fibers, and, finally holey fibers with both circular air hole and elliptical air holes. Comparisons with other approaches are also given.
3.1. Step-Index Fiber
We first calculate the fundamental effective index of a step-index circular optical fiber. Since modal characteristics of a step-index circular fiber can be obtained analytically, we can illustrate the accuracy of our numerical results by comparing them with the analytical solution. We also present the results from FD mode solvers by Huang et al.[23] and by Lüsse et al.[27].
The step-index circular fiber we considered has core diameter of 6 μm, core refractive index of 1.45 at wavelength λ=1.5 μm, and air cladding with unity refractive index. The analytical solution gives the fundamental mode index neff = 1.438604. We choose such a high index difference fiber in order to illustrate the applicability of the FD mode solver in modeling waveguides with curved high-index difference interfaces. In Fig.2 we show the relative error of the calculated fundamental mode index neff from the three FD solvers (both with and without index averaging) when the number of grids along the x-axis is varied. As can be clearly seen, the index averaging technique can significantly stabilize and improve the accuracy for the three FD mode solvers. In Table 1, we list the calculated results from different approaches. Uniform meshes are used in all the numerical calculations.
We also carried out calculations on fibers with different parameters. All the comparisons reveal that the index averaging technique can significantly improve modeling accuracy. The proposed FD mode solver and the mode solver by Huang et al. [23] both give better results than that from the approach by Lüsse et al. [27] when few grids are used. Compared with that by Huang et al., the current mode solver appears to get better results even with very sparse grids.
Table 1. Calculated fundamental mode indices of a step-index fiber from different FD mode solvers. The fiber parameters are the same as in Fig. 2. The analytical solution for the fundamental mode index is 1.438604.
Fig. 2. Relative error in the calculated fundamental mode index neff of a step-index circular optical fiber. The fiber has a core diameter 6 μm, a core refractive index of 1.45 at wavelength 1.5 μm, and air cladding with unity refractive index. The calculation window is chosen to be the first quadrant of the fiber cross section with a computation window size of 6μm by 6 μm. The left boundary is magnetic wall, the bottom is electric wall; all others are zero-value boundaries.
3.2. Microstructured Optical Fiber (MOF)
With the demonstrated validity of the FD mode solvers in modeling conventional optical fiber, in this sub-section we apply them in the characterization of two types of MOFs: air-hole assisted optical fibers and holey optical fibers.
Air-hole-assisted optical fiber (AHAOF) [31] is a microstructured fiber that has air holes surrounding the raised-index guiding core (see Fig 3(a)). One interesting aspect of the AHAOF is that it has easily tailorable dispersion and polarization properties while having low transmission loss comparable to conventional optical fiber. The modal characteristics of AHAOF can be accurately obtained by using the multipole method [32]. Here we use the FD mode solvers to examine this fiber and compare the results with those by the multipole analysis. In Table 2, we list calculated fundamental mode indices. The results from both methods agree very well. In Fig. 3(b) the comparison between different mode solvers is made. In Fig. 3(c) we show the magnetic field plots of the y-polarized fundamental mode. The field amplitudes are not normalized. The relative comparison of amplitude between two cartesian components of the magnetic field indicates that Hx dominates in the y-polarized mode. In Fig. 3(d) we show the group velocity dispersion (GVD) curves calculated from the FD method and the multipole method. The agreement is excellent. In the calculation of the GVD curves, we did not consider material dispersion which is dominant at short wavelengths and assumed constant refractive index at all wavelengths investigated, i.e., n core =1.45 and n clad=1.0. While it is trivial to take into account the material dispersion by changing the refractive index at each wavelength, small differences in GVD between different approaches are better appreciated by excluding material dispersion.
Table 2. Calculated fundamental mode indices for air-hole-assisted optical fiber shown in Fig. 3(a). Core index 1.45, silica cladding index 1.42, r 0=2μm, r a=2μm, Λ=5μm, wavelength=1.5 μm, first quadrant window size 8μm by 8μm, C1: magnetic wall, C2~C4: electric wall. The multipole analysis [32] gives neff =1.4353607.
Fig. 3. (a) Schematic of an air-hole assisted optical fiber; (b) Relative errors of calculated fundamental mode index from different mode solvers when the number of grids along the x-axis is varied (fiber parameters are the same as in Table 2); (c) Magnetic field plots of the fundamental mode (y-polarization); (d) Calculated GVD curves (material dispersion not considered).
Fig. 4. (a) Calculation window of the holey fiber; (b) Comparison between different FD mode solvers (Λ=2.3 μm, r a=0.5 μm, silica refractive index of 1.45 is assumed at λ=1.5 μm, calculation window of 3Λ by 3Λ is the first quadrant bounded by C1~C4. C1: magnetic wall, C2~C4: electric wall); (c) Field plots of the y-polarized fundamental mode; (d) Calculated effective modal index β/k 0, nFSM and GVD (Silica refractive index of 1.45 is assumed at all wavelengths).
Next we analyze the holey fiber (Fig. 4(a), triangular array of circular air holes shown) with the FD mode solvers. In holey optical fibers, a mode is guided if its effective mode index satisfies nFSM < neff < nsi , where nsi is the refractive index of silica, and nFSM is the effective index of the fundamental space-filling mode (FSM) [33,34] of the infinite period cladding. nFSM is an important parameter in the effective index model of holey fibers, which can give estimation of the number of guided modes and the effective mode areas [9,12]. nFSM can be obtained through a few methods including the approximate analytical solution [33], plane wave expansion method [34], and the finite element method [18]. In this paper we use the FD methods to calculate nFSM in the elementary region (see Fig. 4(a), region bound by B1~B4) with proper boundary conditions. For example, boundary conditions for an x-polarized FSM can be: B1, B2 electric wall, B3 and B4 magnetic wall.
The effective modal indices of guided modes are calculated in the first quadrant window shown in Fig. 4(a). The results are given in Table 3 and in Fig 4(b). The field plots for the fundamental mode are shown in Fig. 4(c). As in Fig. 3(c), the field is not normalized. In Fig. 4(d) we show the calculated GVD curves from both FD mode solvers and the localized-function method. The latter seems to suffer less accuracy in the long wavelengths.
Modal birefringence between the two fundamental modes in MOFs with circular air holes has been observed both experimentally [5,10] and numerically [14,17], but remained unclear until the theoretical and numerical study by Steel et al. [35]. It is clear now that the fundamental modes are two-fold degenerate in waveguide structures with Cnv (n≥3) symmetry. The modal degeneracy of air-hole assisted fiber has been numerically confirmed by the multipole method [32], and the holey fibers by multipole method and other methods [35,36]. Here we test the FD method by numerically inspecting the degeneracy properties of MOFs. We consider both the air-hole assisted fiber (Fig.3(a)) and the holey fiber (Fig.4(a)). The parameters for the air-assisted fiber are same as in Table 2, and the holey fiber has same parameters as in Fig. 4(b), with λ=1.5 μm. Fig. 5 shows the convergence behavior of the fundamental modal birefringence B=|neff,x -neff,y |, where neff,x and neff,y are the effective index for the x-polarized and y-polarized fundamental modes, respectively. The observed birefringence is attributed to the discretization in the FD method, which introduces different symmetry from that of the ideal waveguide structure and lifts the expected degeneracy.
Fig. 5. Convergence of modal birefringence of the fundamental modes in (a) air-hole assisted optical fiber and in (b) holey optical fiber. The parameters for the two fibers are listed in Table 2 and Fig. 4(b), respectively.
Finally, we apply this method to the calculation of the fundamental mode and FSM indices, and the GVD of a holey fiber with elliptical air-holes. The structure analyzed is similar to that of Fig. 4(a), but with elliptical holes (eccentricity e=0.375) in the place of the circular holes. Steel et al. [37] analyzed these structures using the plane wave expansion method, and identified possible applications. Fig. 6 shows the numerical results for the two lowest order modes. As expected, the degeneracy of the fundamental mode is lifted, with a predicted birefringence of B=6.5×10-4 at λ=1.5 μm. The indices of the x-polarized and y-polarized FSM are also calculated and show much larger birefringence that the fundamental modes. This has important implications for the design of polarization-maintaining MOFs because it can affect both the single-mode cutoff condition and the differential loss (e.g. bending loss or leakage). As before, the material dispersion is not considered in the calculations.
Fig. 6. (a) Calculated effective index and GVD of the x- and y-polarized fundamental modes for a holey fiber with elliptical air holes. The air holes have their major axis along the y-direction. The semimajor axis is 0.8 μm; the semiminor axis is 0.5 μm. Other calculation parameters are same as in Fig. 4. The effective indices of x-polarized and y-polarized FSM are also shown. (b) Contour plots of major magnetic field components for x-polarized (left) and y-polarized (right) fundamental modes at wavelength of 1.5 μm.
4. Discussion and Conclusion
We have carried out numerical computations of modal characteristics of MOFs using different FD mode solvers together with an index averaging technique. In the calculations, we found that grid size ~λ/15 is usually enough to get fairly accurate results. However, it is wise to check the numerical accuracy and convergence by increasing the computational window size and/or decreasing the grid size. Moreover, making use of the mode symmetry can greatly reduce computational time. For a large enough computation window, the selection of closed boundary conditions (zero-value, electric wall, or magnetic wall) for window edges (except symmetry axes) seems to have negligible influence on results. However, a large computation window leads to a large matrix and thus more computational time. This is especially true when modeling higher order modes or modes near cut-off. We note that all the computations give real propagation constants for guided modes investigated. This contrasts with the multipole method which can give complex propagation constant with the imaginary part related to the leakage loss of finite MOFs [21]. It is possible, however, to apply open boundary conditions in the FD methods and model leakage loss of guided modes, though complicated formulation and iterations are needed [1].
In conclusion, we have presented a full-vectorial finite-difference mode solver based Yee’s mesh, and have applied an index averaging technique in this mode solver. Index averaging shows improvement for all of the FD mode solvers considered and is particularly effective when combined with Yee’s mesh. This results in significantly accelerated convergence and improved numerical accuracy. The validity and effectiveness of the mode solver was demonstrated in conventional step-index fibers. We investigated the modal characteristics of some MOFs with the presented method and achieved excellent agreement with previous FD mode solvers and other numerical methods.
Acknowledgement
This work is supported in part by the Laboratory for Laser Energetics of the University of Rochester and by the National Science Foundation, grant number ECS-9816251.
References and Links
1. K. S. Chiang, “Review of numerical and approximate methods for the modal analysis of general optical dielectric waveguides,” Opt. Quantum Electron. 26, s113–s134 (1994). [CrossRef]
2. C. Vassallo, “1993-1995 optical mode solvers,” Opt. Quantum Electron. 29, 95–114 (1997). [CrossRef]
3. 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, 1547–1549 (1996). [CrossRef] [PubMed]
4. T. A. Birks, J. C. Knight, and P. St. J. Russell, “Endlessly single-mode photonic crystal fiber,” Opt. Lett. 22, 961–963 (1997). [CrossRef] [PubMed]
5. J. C. Knight, J. Arriaga, T. A. Birks, A. Ortigosa-Blanch, W. J. Wadsworth, and P. St. J. Russell, “Anomalous dispersion in photonic crystal fiber,” IEEE Photon. Technol. Lett. 12, 807–809 (2000). [CrossRef]
6. A. Ferrando, E. Silvestre, J. J. Miret, and P. Andres, “Nearly zero ultraflattened dispersion in photonic crystal fibers,” Opt. Lett. 25, 790–792 (2000). [CrossRef]
7. T. P. Hansen, J. Broeng, S. E. B. Libori, E. Knudsen, A. Bjarklev, J. R. Jensen, and H. Simonsen, “Highly birefringent index-guiding photonic crystal fibers,” IEEE Photon. Technol. Lett. 13, 588–590 (2001). [CrossRef]
8. N. G. R. Broderick, T. M. Monro, P. J. Bennett, and D. J. Richardson, “Nonlinearity in holey optical fibers: measurement and future opportunities,” Opt. Lett. 24, 1395–1397 (1999). [CrossRef]
9. N. A. Mortensen, “Effective area of photonic crystal fibers,” Opt. Express 10, 341–348 (2002), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-7-341 [CrossRef] [PubMed]
10. J. K. Ranka, R. S. Windeler, and A. J. Stentz, “Visible continuum generation in air silica microstructure optical fibers with anomalous dispersion at 800nm,” Opt. Lett. 25, 25–27 (2000). [CrossRef]
11. B. J. Eggleton, C. Kerbage, P. Westbrook, R. S. Windeler, and A. Hale, “Microstructured optical fiber devices,” Opt. Express 9, 698–713 (2001), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-9-13-698 [CrossRef] [PubMed]
12. J. C. Knight, T. A. Birks, P. St. J. Russell, and J. P. de Sandro, “Properties of photonic crystal fiber and the effective index model,” J. Opt. Soc. Am. A 15, 748–752 (1998). [CrossRef]
13. A. Ferrando, E. Silvestre, J. J. Miret, P. Andres, and M. V. Andres, “Full-vector analysis of a realistic photonic crystal fiber,” Opt. Lett. 24, 276–278 (1999). [CrossRef]
14. A. Ferrando, E. Silvestre, J. J. Miret, and P. Andres, “Vector description of higher-order modes in photonic crystal fibers,” J. Opt. Soc. Am. A 17, 1333–1340 (2000). [CrossRef]
15. S. G. Johnson and J. D. Joannopoulos, “Block-iterative frequency-domain methods for Maxwell’s equations in a planewave basis,” Opt. Express 8, 173–190 (2001), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-8-3-173 [CrossRef] [PubMed]
16. D. Mogilevtsev, T. A. Birks, and P. St. J. Russell, “Localized function method for modeling defect modes in 2-D photonic crystals,” J. Lightwave Technol. 17, 2078–2081(1999). [CrossRef]
17. T.M. Monro, D. J. Richardson, N. G. R. Broderick, and P. J. Bennett, “Modeling large air fraction holey optical fibers,” J. Lightwave Technol. 18, 50–56 (2000). [CrossRef]
18. F. Brechet, J. Marcou, D. Pagnoux, and P. Roy, “Complete analysis of the characteristics of propagation into photonic crystal fibers by the finite element method,” Opt. Fiber Technol. 6, 181–191 (2000). [CrossRef]
19. G. E. Town and J. T. Lizer, “Tapered holey fibers for spot size and numerical-aperture conversion,” Opt. Lett. 26, 1042–1044 (2001). [CrossRef]
20. M. Qiu, “Analysis of guided modes in photonic crystal fibers using the finite-difference time-domain method,” Microwave Opt. Technol. Lett. 30, 327–330 (2001). [CrossRef]
21. T. P. White, R. C. McPhedran, C. M. de Sterke, L. C. Botten, and M. J. Steel, “Confinement losses in microstructured optical fibers,” Opt. Lett. 26, 1660–1662 (2001). [CrossRef]
22. M. S. Stern, “Semivectorial polarized finite difference method for optical waveguides with arbitrary index profiles,” IEE Proc. J. Optoelectron. 135, 56–63 (1988). [CrossRef]
23. W. P. Huang and C. L. Xu, “Simulation of three-dimensional optical waveguides by a full-vector beam propagation method,” IEEE J. Quantum Electron. 29, 2639–2649 (1993). [CrossRef]
24. W. P. Huang, C. L. Xu, S. T. Chu, and S. K. Chaudhuri, “The finite-difference vector beam propagation method. Analysis and Assessment,” J. Lightwave Technol. 10, 295–305 (1992). [CrossRef]
25. K. Bierwirth, N. Schulz, and F. Arndt, “Finite-difference analysis of rectangular dielectric waveguide structures,” IEEE Trans. Microwave Theory Tech. 34, 1104–1113 (1986). [CrossRef]
26. H. Dong, A. Chronopoulos, J. Zou, and A. Gopinath, “Vectorial integrated finite-difference analysis of dielectric waveguides,” J. Lightwave Technol. 11, 1559–1563 (1993). [CrossRef]
27. P. Lüsse, P. Stuwe, J. Schüle, and H. G. Unger, “Analysis of vectorial mode fields in optical waveguides by a new finite difference method,” J. Lightwave Technol. 12, 487–493 (1994). [CrossRef]
28. K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propagat. 14, 302–307 (1966). [CrossRef]
29. K. S. Kunz and R. J. Luebbers, The Finite Difference Time Domain Method for Electromagnetics, (CRC, Boca Raton, 1993).
30. S. Dey and R. Mittra, “A conformal finite-difference time-domain technique for modeling cylindrical dielectric resonators,” IEEE Trans. Microwave Theory Tech. 47, 1717–1739 (1999).
31. T. Hasegawa, E. Sasaoka, M. Onishi, M. Nishimura, Y. Tsuji, and M. Koshiba, “Hole-assisted lightguide fiber for large anomalous dispersion and low optical loss,” Opt. Express 9, 681–686 (2001), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-9-13-681 [CrossRef] [PubMed]
32. Z. Zhu and T. G. Brown, “Multipole analysis of hole-assisted optical fibers”, Opt. Commun. 206, 333–339 (2002). [CrossRef]
33. M. Midrio, M. P. Singh, and C. G. Someda, “The space filling mode of holey fibers: an analytical vectorial solution,” J. Lightwave Technol. 18, 1031–1037 (2000). [CrossRef]
34. Z. Zhu and T. G. Brown, “Analysis of the space filling modes of photonic crystal fibers,” Opt. Express 8, 547–554 (2001), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-8-10-547 [CrossRef] [PubMed]
35. M. J. Steel, T. P. White, C. M. de Sterke, R. C. McPhedran, and L. C. Botten, “Symmetry and degeneracy in microstructured optical fibers,” Opt. Lett. 26, 488–490 (2001). [CrossRef]
36. M. Koshiba and K. Saitoh, “Numerical verification of degeneracy in hexagonal photonic crystal fibers,” IEEE Photon. Technol. Lett. 13, 1313–1315 (2001). [CrossRef]
37. M. J. Steel and R. M. Osgood Jr., “Elliptical-hole photonic crystal fibers,” Opt. Lett. 26, 229–231 (2001). [CrossRef]
