1. Introduction

Photonic crystal fibers (PCFs) have attracted enormous interest in the last decade because of their unique optical properties [1, 2, 3, 4]. PCFs are a special class of fibers characterized by a distribution of air holes running along the axis of the fiber. Light in a PCF can be confined by one of two mechanisms: total internal reflection (TIR) and photonic band gap (PBG) guidance [5, 6]. Broadly speaking, TIR is used for “solid core” PCFs, where air holes surround a high index solid core. PBG guidance is used for “hollow core” PCFs where the light is confined to a central air hole. In both cases, the propagation of light in PCFs is controlled by the precise geometry and detailed distribution of the air holes themselves. Thus, in order to fully realize the potential of PCFs through optimal design, a fast, accurate, and versatile simulation method is essential.

Although significant progress has been made, the rigorous analysis of light propagation in PCFs remains problematic because of the large index contrast, the vectorial nature of the Maxwell equations, and the complicated cross-sections of the hole geometries involved. The earliest attempt to model PCFs is the effective index method [3], which provided some useful insights. Since then, a range of methods have been developed: the plane wave expansion method [9, 10], the localized function method [7], the multipole expansion method [11, 12], finite element methods [13, 14] and boundary element methods [15, 16, 17]. Many, if not most, PCF applications involve leaky modes. Given this fact, any method that claims to accurately simulate light propagation in a PCF should provide some estimate of the confinement loss. A standard mathematical expression for this is the imaginary part of the propagation constant β in the axial direction. We should note that in principle, all the simulation methods mentioned above (with the exception of the plane wave expansion method, which implies an infinite periodic medium) are theoretically capable of yielding the confinement loss, although the computation may be expensive. Nevertheless, to the authors’ knowledge, such estimates have been published only with the multipole expansion method [11, 12] and the finite element method [19]. The multipole expansion method [11, 12] requires that the air holes be perfectly circular and that the holes be well separated from one another; while the finite element method requires a volume discretization, as well as considerable human intervention to fine tune the perfectly matched layer (PML) boundary condition and the mesh discretization. Furthermore, in order to reduce the computational burden, most simulations have also required that the holes be arranged in highly symmetric patterns, reducing the number of degrees of freedom dramatically. As indicated earlier, many questions regarding design and manufacturability can only be addressed if we can overcome these constraints while preserving high accuracy.

In this paper, we develop an integral equation (Green’s function) method for studying wave propagation in PCFs. The mathematical formulation is discussed in section 2, and is related to the method developed in [18]. We show, in section 3, that PCFs with non-circular holes and irregular hole distributions can be analyzed with high precision. The subsequent coupling of our discretization technique to fast algorithms for large-scale computation will be discussed in a companion paper [20].

2. Integral equation formulation

For the modal analysis of PCF, we follow standard convention and assume the longitudinal axis is the z-axis, and that the structure of the fiber is defined by its cross-section in the xy-plane. A typical section is shown in Fig. 1. More precisely, the fiber consists of a glass matrix denoted by V ₀ with refractive index n ₀ which is perforated by a finite number of “holes” denoted by V ₁,V ₂, …,V _M . These “holes” can be composed of any homogeneous isotropic material (including air), and their refractive indices are assumed to be given by n ₁,n ₂,…,n_M , respectively.

For a propagating mode in the fiber [21], we represent the vector electric E(x,y,z,t) and magnetic H(x,y,z,t) fields as

Here, ω is the angular frequency and β is the propagation constant of the particular mode under consideration. Modes in PCFs are typically leaky, so that β is a complex number and its imaginary part is a measure of the attenuation of the mode along the fiber.

From Maxwell’s equation, the fields can be expressed entirely in terms of the longitudinal z-components of the electric and magnetic fields. The other four components can be obtained from the relations [21]

$H_1^{\beta }=\frac{-i\omega n^2{\epsilon }_{\mathit{vac}}}{k^2-{\beta }^2}{\partial }_y{E}_3^{\beta }+\frac{i\beta }{k^2-{\beta }^2}{\partial }_x{H}_3^{\beta },$ $H_2^{\beta }=\frac{-i\omega n^2{\epsilon }_{\mathit{vac}}}{k^2-{\beta }^2}{\partial }_x{E}_3^{\beta }+\frac{i\beta }{k^2-{\beta }^2}{\partial }_y{H}_3^{\beta },$

$E_1^{\beta }=\frac{i\omega \mu }{k^2-{\beta }^2}{\partial }_y{H}_3^{\beta }+\frac{i\beta }{k^2-{\beta }^2}{\partial }_x{E}_3^{\beta },$ $E_2^{\beta }=\frac{-i\omega \mu }{k^2-{\beta }^2}{\partial }_x{H}_3^{\beta }+\frac{i\beta }{k^2-{\beta }^2}{\partial }_y{E}_3^{\beta }.$

 

Fig. 1. Cross-section of a model PCF, whose longitudinal axis runs in the z-direction. V ₀ denotes the glass matrix with index n ₀. In this case, there are only four “holes” in the fiber V ₁,…,V ₄ made of materials with refractive indices n ₁,…,n ₄. At each point on a material interface, ν denotes the unit normal vector and τ denotes the unit tangent vector.

Download Full Size | PDF

Here, ε_vac and µ_vac denote the permittivity and permeability of free-space, n is the given index of refraction, and the wave number k is defined by k=ωn/c, where c is the speed of light in vacuum. When the context is clear, we will denote the local wave number by k_j =ωn_j /c. In free space, we have k_vac =ω/c. We have assumed, for simplicity, that the materials are not magnetic so that the permeability is given by µ_vac throughout.

The required interface conditions are well-known, namely that the tangential components of the E and H fields be continuous. Thus,

where ν is the normal vector on the interface and [f] denotes the jump in the quantity f across an interface. For notational convenience, in the remainder of this paper, we will use the scalar unknowns E=$E_3^{\beta }$ and $H=\sqrt{\frac{{\mu }_{\mathit{vac}}}{{\epsilon }_{\mathit{vac}}}}{H}_3^{\beta }$. Using these variables, we can summarize the eigenvalue problem which determines a mode with propagation constant β as follows:

These equations can be found in standard texts and form the starting point for the analysis in [18] as well.

2.1. Representations

To derive an integral equation for the mode calculation, we begin with global representations for the E and H fields which are motivated by classical potential theory [22, 23]. The advantage of the following representations is that the fields obey Maxwell’s equations trivially everywhere except on the boundary curves Γ _j . The representations are based on the free-space Green’s function for the Helmholtz equation

$G_j(P,Q)=\frac{i}{4}{H}_0^{\left(1\right)}\left(\sqrt{k_j^2-{\beta }^2}\parallel P-Q\parallel \right),$,

where i=√-1 and $H_0^{\left(1\right)}$ is the usual Hankel function of the first kind of order zero. In the following, Γ _j denotes the boundary of the j^th hole. If P lies inside the j^th hole, we let

If P lies outside of all holes (P∈V ₀), we let

Here, M denotes the number of holes. In analogy to classical potential theory, ${\mathrm{\sigma }}_j^E$ , ${\mathrm{\sigma }}_j^H$ , are called single layer “densities” and ${\mathit{\mu }}_j^E$ ,${\mathit{\mu }}_j^H$ are called double layer “densities”, defined on the hole boundaries. Finally, the quantities S^E (1,j), S^H (1,J), S^E (2,J),S^H (2,J), d^E (1,J),d^H (1,J), and d^E (2,j)d^H (2,j) are free parameters at our disposal. We have chosen to use both single and double layer densities, as well as introducing the constants S^E (1,j) etc. in order to have sufficient freedom to obtain desirable properties for the integral equations. Their use will be described in the following.

Although Maxwell’s equations are satisfied away from the boundaries Γ _j , it remains to examine the representations Eqs. (7)–(10) upon the boundaries. When the evaluation point P ₀ is precisely on an interface, care must be taken in the calculation of appropriate limits, involving principal value integrals [24]. For example, if the limit is taken from outside the j^th hole (i.e. inside the glass), the E field takes the form

However, when the limiting process is carried out from inside the j^th hole, we have

Analogous limits hold for the H field and its normal and tangential derivatives.

2.2. The integral equation

We are finally in a position to derive the integral equation for the four layer potentials. Substituting the integral representations into the boundary conditions, Eqs. (3)–(6), we arrive at

where

$b_t^E(1,j)=\frac{\beta }{k_0^2-{\beta }^2},b_n^E(1,j)=-\frac{k_{\mathit{vac}}}{k_0^2-{\beta }^2},b_t^H(1,j)=\frac{\beta }{k_0^2-{\beta }^2},b_n^H(1,j)=\frac{k_0^2}{k_{\mathit{vac}}}{k}_j^2-{\beta }^2,$ $b_t^E(2,j)=\frac{\beta }{k_j^2-{\beta }^2},b_n^E(2,j)=-\frac{k_{\mathit{vac}}}{k_j^2-{\beta }^2},b_t^H(2,j)=\frac{\beta }{k_j^2-{\beta }^2},b_n^H(2,j)=\frac{k_j^2⁄k_{\mathit{vac}}}{k_{\mathit{vac}}-{\beta }^2},$,

are a set of constants associated with the j^th hole. Note now that the last terms in Eqs. (19) and (20) consist of the differences of complicated expressions involving the normal derivatives of double layer potentials. These are generally referred to as hypersingular and are problematic both from the point of view of quadrature and because of adverse effects on the behavior of iterative methods. Fortunately, by choosing the coefficients d^E (1,j)d^E (2,j),d^H (1,j) and d^H (2,j) so that

we are able to cancel the hypersingular terms (at least to leading order) and remove the associated difficulties. We would also like to preserve as much diagonal dominance as possible in the linear system. For this, note that the fifth terms in Eqs. (19) and (20) consist of some constant multiple of the densities ${\sigma }_j^H$ and ${\sigma }_j^E$ , respectively. Since these are diagonal operators, we would like to make sure their influence doesn’t vanish. Thus, we choose the coefficients S^E (1,j),S^E (2,j),S^H (1,j) and S^H (2,j) so that

These choices are intended to make the integral equation as close to what is known as a “Fredholm equation of the second kind” as possible [25]. Ignoring the mathematical details, such equations tend to lead to more robust discretizations and are better suited for iterative solution techniques. Unfortunately, our system of equations is not quite in this class because of the terms involving tangential derivatives

$\frac{1}{2}\left(b_t^E(1,j)d^E(1,j)+b_t^E(2,j)d^E(2,j)\right)\frac{\partial {\mu }_j^E\left(P_0\right)}{\partial \tau }$

and

$\frac{1}{2}\left(b_t^H(1,j)d^H(1,j)+b_t^H(2,j)d^H(2,j)\right)\frac{\partial {\mu }_j^H\left(P_0\right)}{\partial \tau }$

For algebraic reasons, these terms cannot be eliminated while simultaneously eliminating the hypersingular terms. It is, perhaps, worth noting here how our integral equation differs from that proposed in [18]. In [18], instead of representing the electric and magnetic fields via single and double layer densities as in Eqs. (17)–(20), the authors chose to represent the field via two distinct single layer densities, one for the field inside the hole and the other for the glass matrix. They arrived at a mixture of first and second kind integral equations, which generally requires a preconditioner for satisfactory convergence behavior.

It is important to note that in our formulation, one is no longer seeking eigenvalues of a differential equation. Instead, as is common in optical waveguide analysis, one is seeking nontrivial solutions to a homogeneous linear system in which β plays a nonlinear role. The advantage gained is that we avoid having to discretize the cross-sectional structure, only the material interfaces themselves. To find such solutions numerically, we first have to discretize all the integrals in Eqs. (17)–(20). Using techniques like those in [27], we have designed specialized quadrature rules which integrate, up to a given precision, all functions which appear. Moreover, they do so with “spectral” or “superalgebraic” accuracy [26]. This means that the error decreases faster than any finite power of h, where h is the mesh spacing on an interface. The impact on simulations will be seen in the numerical examples. The details of these quadratures will be given in [20].

2.3. Finding modes

After the integrals are discretized by high-order accurate quadratures, we obtain a rather involved linear system which, for a fixed value of β, we write in the generic form

The unknown vector x represents point values of the single and double layer densities on all the hole boundaries. Thus, if N points are used to discretize all interfaces, there are 4N unknowns. Inspection of the full equations in the preceding section shows that the entries of the matrix A are analytic, nonlinear functions of β. Finding the propagating modes corresponds to finding values of β for which the system Eq. (21) has nontrivial solutions.

Our mode searching strategy is based on defining a new function of β:

$f\left(\beta \right)=\frac{1}{(\overrightarrow{\varphi },{\mathbf{A}}^{-1}\left(\beta \right)\overrightarrow{\psi })},$,

where ϕ⃗ and ψ⃗ are two fixed random vectors. It is straightforward to verify that the function f(β) is an analytic function of its argument. Moreover, since ‖A ^-1(β)_‖→∞ when β corresponds to a mode, we have that f(β)=0. In short, the singular matrix problem has been turned into a complex root finding process for the function f(β). Many tools are available for finding complex roots of scalar equations; we have chosen to work with Muller’s method [26] because it is both efficient and robust.

3. Numerical results

3.1. Comparison with previous results

For our first example, we use a problem described in [11]. The cross-sectional geometry of the example is given in Fig. 2. The geometry consists of a single ring of six equally spaced circular holes with hole diameter d=5µm and hole pitch Λ=6.75µm. The glass matrix is assumed to have a refractive index of n ₀=1.45 while that of the air holes is n_i =1.0. The wavelength for the example is λ=1450 nm.

In numerical simulations of complex physical phenomena, a priori estimates of accuracy are useful guidelines, but one does not know the accuracy with confidence until a convergence study is done. Figure 3 presents a convergence analysis of the effective index of the second highest mode (listed below in Table 1) as a function of the number of discretization (quadrature) points per hole. The errors are approximated by comparison with the effective indices obtained from 70 quadrature points per hole. The figure demonstrates the spectral convergence of our scheme.

 

Fig. 2. Six circular holes.

Download Full Size | PDF

 

Fig. 3. Convergence study for the second mode

Download Full Size | PDF

Table 1 lists our results for the propagating modes. The results are obtained by discretizing each hole using 50 quadrature points (leading to a linear system of size 1200). As we showed earlier in our convergence study, 50 quadrature points per hole yields 11 significant digits of accuracy, of which we show 10 digits. For comparison, we also list, in columns 5 and 6, the corresponding results of [11] using the multipole expansion method. The results of [11] are seen to be very accurate: their effective index of the fundamental mode agrees with our result to seven digits of accuracy. This is consistent with the results in our convergence study, Figure 3: the data in [11] are obtained by using multipole expansions of order M=5, which is equivalent to using 21 discretization points in our scheme. Following [11], the modal confinement loss (in dB/m) listed in column 4 are obtained from the imaginary part of the effective index via the formula

with the vacuum wavelength λ.

Table 1. Effective index of the PCF depicted in Fig. 2

View Table | View all tables in this article

3.2. Perturbations on the boundary of holes

In our second example, we demonstrate the method on non-circular boundaries. The cross-sectional geometry is shown in Fig. 4. This geometry is a perturbation of that in the previous example: the change is that the circles are made into “cookies” by replacing the parameterization

for circles with

$\parallel \mathbf{r}\left(\theta \right)-{\mathbf{c}}_i\parallel =r\left[1+\frac{h}{100}\sin \left(7\theta \right)\right],$,

where r(θ) is the position vector of points on the ith hole boundary, 0≤θ<2π, r=d/2 is the average radius and c_i is the ith hole center location. The cookie shape is not meant to be a realistic model for manufacturing defects in PCFs.

Figure 5 shows the convergence analysis for this example. We can observe from the figure that spectral accuracy is again achieved after the Maxwell equations are fully resolved. The number of points per hole required for full accuracy is larger than in the previous example, but this is expected because of the need to fully resolve the seven bumps on the boundary of the holes.

 

Fig. 4. Six cookie-shaped holes.

Download Full Size | PDF

Table 2 shows how the effective indices of the fundamental mode polarizations change as we gradually raise the height of the cookie bumps. From the data, we find that the imperfection of the boundary does not have a large effect on the propagation constants: the slight change shown in the table can be attributed as much to the change of the air hole area as to the change in the shape itself. The gradual splitting of the two polarizations is due to that fact that the geometry in Fig. 4 no longer has perfect six-fold symmetry.

 

Fig. 5. Convergence study of the PCF depicted in Fig.4 (h=6%)

Download Full Size | PDF

Table 2. Effective index of the PCF depicted in Fig. 4

View Table | View all tables in this article

3.3. Irregular hole shapes and locations

Our last example deals with PCFs whose cross-sectional geometries are shown in Figures 6 and 7. Fig. 6 shows a PCF with three layers of circular holes arranged in a triangular lattice pattern. The circular holes have the same diameter d=1.6µm and the hole-to-hole pitch is Λ=3µm. The refractive indices and the wavelength λ=1450nm are the same as those in Example 1. The PCF shown in Fig. 7 is obtained from Fig. 6 by changing the circular holes into ellipses and by displacing the centers of the holes in random directions for short distances (the actual distance by which each hole is moved is also random, with a maximum distance of 0.2µm). The ellipses have a long axis equal to the radius of the circles and a short axis 15% shorter than the long one. The orientations of the ellipses are also random.

 

Fig. 6. A PCF with circular holes.

Download Full Size | PDF

 

Fig. 7. A PCF with irregular shaped holes.

Download Full Size | PDF

The results of our computation for this example are presented in Figures 8, 9, 10 and Table 3. Figure 8 is the convergence study for the PCF shown in Fig. 7. Once again, our method achieves spectral accuracy. Fewer points are required per hole to fully resolve the Maxwell equations in this example because the hole sizes are small. Figures 9 and 10 show the dispersion curves of the fundamental modes.

The data in Table 3 show how the effective index of the fundamental mode changes when we increase the number of hole layers from one to three. For comparison, the last row shows the effective index of the PCF with three regular layers of circular holes as in Fig. 6. The data show the confinement loss of the fundamental modes drops significantly from the one layer case to the case with three layers. The last line of the table also shows that the changes in the effective indices are not large when we change hole shape from regular to irregular, despite the fact that the irregular holes have an area 15% smaller than the circular holes.

 

Fig. 8. Convergence study of the PCF depicted in Fig.7

Download Full Size | PDF

Table 3. Effective index of the fundamental mode of PCFs with one to three layers of irregular holes and irregular locations (Fig.7). The results in last row are for the PCF with three layers of regular holes (Fig.6).

View Table | View all tables in this article

4. Conclusions

We have presented a new integral equation method for the calculation of the modes in a PCF, which is capable of fully-resolved solutions of Maxwell’s equations in arbitrary PCF geometries. The method discretizes the material interfaces only, is robust even when the index contrasts are large, and requires a small number of degrees of freedom because of the underlying high order (spectral) accuracy.

 

Fig. 9. Dispersion curve (real part).

Download Full Size | PDF

 

Fig. 10. Dispersion curve (imaginary part).

Download Full Size | PDF

With this method (and the acceleration techniques described in [20]), one can accurately calculate physical quantities over wide dynamic ranges, such as the real and imaginary parts of a slightly leaky mode. High accuracy also makes it possible to compute finite difference approximations of derivatives of physical quantities with respect to changes in the geometry. This is an essential element in design optimization.

The present implementation does not consider a finite outer cladding nor does it allow for holes with complex internal structures. These restrictions are straightforward to remove.