1. Introduction

Microstructured optical fibers (MOFs) have recently attracted much interest because of their unique optical properties [2, 3, 4, 5]. Light in a MOF is confined by one of the three principal mechanisms: modified total internal reflection, photonic band gap guidance and antiresonant guidance. In all the cases, the propagation of light is strongly influenced by the geometry of the fiber cross section. To fully realize the potential of MOFs through optimal design, availability of the high performance simulation tools is essential. Although significant progress has been made in modelling of MOFs, rigorous modal analysis remains challenging, especially in the case of large number of holes or non-cylindrical hole shapes. Additional complexity comes from the MOF large refractive-index contrast, fully vectorial modal fields, and leaky nature of the modes.

A number of modelling techniques has been developed for finding the modes of MOFs - the plane wave expansion method [6], the localized function method [7] (in a scalar approximation), finite element and finite difference methods [8, 9, 10], boundary element methods [11, 12, 13, 14], and multipole expansion method [15, 16, 17, 18, 19]. The finite element and finite difference methods require volume discretization which demands considerable memory resources to achieve high accuracy. These methods, however, are well suited to treat complex shapes. The multipole method experiences convergence and accuracy problems when the holes are too closely spaced from one another or when the effective refractive index of a mode is too different from that of a cladding. Due to the nature of the expansion, multipole method is best suited for the analysis of circular inclusions, however extensions to treat complex shapes [17] and sub-inclusions [18, 19] have been demonstrated.

In this paper, we develop an accurate and efficient boundary integral method for the modal analysis of MOFs. Here we extend the boundary integral method originally developed for the modal analysis of conventional waveguides [20]. The key advantages of the proposed method are the following. For circular inclusions (similarly to the multipole method) much of the calcu- lations are done analytically resulting in a highly accurate and rapidly convergent implementation. Additionally, the method can treat inclusions of arbitrary shapes defined by a parametrical curve. The method is able to treat accurately systems that contain a large numbers of inclusions (hundreds), due to a well behaved nature of matrix elements used in the formulation. Particularly, the method is able to treat efficiently very closely spaced inclusions, which is of particular interest in the case of hollow core air-guiding fibers. Finally, symmetries of the fiber geometry can be readily taken into account by the numerical algorithm to speed up modal calculations.

 

Fig. 1. Three structures studied in the paper. (a) Hollow coreMOF with five rings of circular holes; the pitch is Λ=2.74µm, the hole diameter is d=.95Λ and the core diameter is d_c=2.5d. (b) Elliptic hollow core MOF with three layers of circular holes; the pitch is Λ=2µm, the hole diameter is d=0.9Λ and the core principal axis are a=2.3µm and b=4.6µm. (c) Solid core MOF with six silver coated elliptical holes; the outer hole principal axis are a_o=0.84µm and b_o=0.76mm, the inner hole principal axis are a_i=0.74µm and b_i=0.66µm, the pitch is Λ=1.5mm.

Download Full Size | PDF

After analyzing its performance, we believe that boundary integral method finds its place in-between the brute force finite element method and the semi-analytical multipole method. Similarly to the multipole method, boundary integral method requires small operational memory and allows almost analytical evaluation of the matrix elements. However, expansions used in the proposed boundary integral method are more stable than those used in the multipole method. This allows to study systems with a large number of arbitrary shaped inclusions, similarly to a finite element method. We validate our method on various simple structures by comparing its predictions with those of the multipole method.We then demonstrate versatility of the method by applying it to several challenging problems. First, we study confinement losses of the core guided modes of a hollow photonic crystal fiber featuring up to 5 reflector layers Fig. 1(a). Second, modal fields and birefringence of a hollow elliptical core fiber with 3 reflector layers Fig. 1(b) are characterized. Finally, birefringence in plasmonic excitation is studied in a solid core microstructured fiber with metallized elliptical holes Fig. 1(c). We show splitting in a doubly degenerate plasmonic mode when hole ellipticity is introduced. Application to pressure sensing is proposed.

2. Mathematical formulation

The schematic of a MOF geometry is shown in Fig. 2. We suppose that the fiber cross-section is located in the xy plane, while the longitudinal axis z is directed along the fiber length. Fiber cross-section consists of a finite number N_c of homogeneous inclusions of refractive index n_c embedded into a homogeneous background material of refractive index n_g. Modal field components are taken to have an exp i(β z-ωt) dependence. Here β is the propagation constant of the mode and ω is the angular frequency, which is related to the free-space wave number by ω=ck₀. Modal effective refractive index is defined as n_e=β/k₀.

 

Fig. 2. a) Schematic of aMOF cross section. b) Schematic of Re(G(s, s′)). Green’s function has a cusp when s → s′ it also exhibits oscillations $\left(\gamma =\sqrt{n_g^2-n_e^2}\right).c)$ Arbitrary shaped inclusion and a corresponding regularization circle.

Download Full Size | PDF

Given the longitudinal components E_z and H_z of the modal electromagnetic fields, all the other field components can be deduced from Maxwells equations. In each of the homogeneous dielectric regions, individual longitudinal components satisfy the Helmholtz equations:

where γ²_c,g=n²_c,g-n²_e. Moreover, on the boundaries between various dielectric regions the longitudinal and tangential components of the fields are continuous:

where $\frac{\partial }{\partial \tau }$ is the tangential derivative to the boundary contour L,$\frac{\partial }{\partial n}$ is the outer normal derivative, and ε_(c,g) is the dielectric constant of either the inclusion or the cladding. Throughout the paper H_z represents the true magnetic field scaled by the free space impedance H_z = μ₀cH^{true field}_z.

Since the field components satisfy the Helmholtz equations (1), at any point r⃗ in the fiber cross section, they can be represented by the following contour integrals (also known as the single layer potentials) [20, 21]:

The contour functions e(r⃗_s) and h(r⃗_s) are called the potential densities [21] and their specification is sufficient to obtain the longitudinal field components. The function G(r⃗,r⃗_s) is the Green’s function of the Helmholtz equation in a uniform medium and it is given by:

where H⁽¹⁾0 () is the zeroth-order Hankel function of the first kind. The value of γ and the contour of integration depend on the position r⃗. If r⃗ is located inside of the kth inclusion then γ=γ_c and integration is taken along the inclusion’s boundary L=L^(k). If r⃑ is located in the fiber cladding then λ = λ_g and integration is taken along all the inclusion boundaries $L={\displaystyle \sum _{k=1}^{N_c}}{L}^{\left(k\right)}.$. Finally if r⃑ is located exactly at the inclusion boundary one can use any of the two definitions. We also note here that such a formulation satisfies the Reichardt condition [20] at infinity so it can treat both “proper” leaky modes that have fields vanishing at infinity and the “improper” leaky modes that grow at infinity. This can be possible by allowing the imaginary part of λ_g to be positive or negative.

We assume now that all the individual inclusion boundaries can be described by the parametric expressions x=x_k(s) and y=y_k(s), s being a parameter such that 0≤s≤2π, and k being the inclusion number. Inserting Eqs. (3) into the continuity conditions (2), we obtain the following equations for any point⃑r_s∈L ^(j):

where $J^{\left(k\right)}\left(s\right)={\left[{(d{x}_k⁄ds)}^2+{(d{y}_k⁄ds)}^2\right]}^{1⁄2}$ is the jacobian of the contour k. This is a system of four coupled linear integral equations for the scalar contour functions e^(k)_c, e^(k)_g, h^(k)_c, h^(k)_g, k=1…N_c. The value of n_e, for which a nontrivial solution of (5) exists, defines the effective refractive index of the fiber mode. Note, that the terms $\frac{e(s)}{2}\mathrm{and}\frac{h(s)}{2}$ appear because of the normal derivatives of the integrals in (3) (the so called double layer potentials) exhibit discontinuities when r⃗_s →L^(j) [21]. Consider now the discretization of Eqs. (5). A direct discretization would run into difficulties as Hankel functions and their tangential derivatives become singular when s′→s. Therefore, such a singularity has to be first removed from the formulation.We start with the case of circular inclusions.

2.1. Circular inclusions

For circular inclusions J^(k)(s)=a_k, where a_k is a radius of the kth inclusion. Let ψ^(k)(s′) represent any of the contour functions e^(k)_c (s′), e^(k)_g (s′), h^(k)_c (s′), or h^(k) _g (s′)and let st = $t\frac{2\pi }{2n^{(k)}}$, t=0,…,2n ^(k)-1, be an equidistant grid. Since ψ^(k)(s′) is a periodic function with a period 2π, we can choose the trigonometric interpolation to approximate ψ^(k)(s′):

As discussed in [22], the error of such an approximation decreases exponentially fast with n^(k). Expansion (6) is a key to our formulation of the boundary integral method. As it will be shown later, in comparison to the prior formulations that use Fourier expansion, using the trigonometric interpolation reduces considerably the numerical cost of resolving Eqs. (5). Consider first the integrals along the contour L^(j). These are the most problematic because of the presence of the singular point⃑r_s′=⃑r_s. When interpolations of the form (6) are inserted into the integrals along the contour L^(j) in (5), these integrals take the following form:

where Φ(s,s′) stands for G(s,s′),$\frac{\partial G(s,s\text{'})}{\partial n}$ or $\frac{\partial G(s,s\text{'})}{\partial \tau }$. In Appendix A we present details of the Green’s function derivative evaluation for the arbitrary contour shapes. The fourier transforms of Φ(s,s′) in (7) can be evaluated analytically according to the following formulas [23]:

where J_m() and H⁽¹⁾_m () are the mth order Bessel and the first kind Hankel functions, while a prime denotes the derivative with respect to the argument. The first two relations can be derived from the Graf’s addition theorem [24], while the third relation is simply the $\left(\frac{1}{a_j}\frac{d}{ds}\right)$ derivative of the first one. Thus, all the integrals in parenthesis of (7) can be evaluated analytically, and the only approximation remaining is due to the trigonometric interpolation of ψ(s′).

We now consider the integrals along the other contours L(k)k≠. Particularly, we are looking for a discrete form of these integrals as a linear function of the values ψ^(k)(st) defined on the corresponding equidistant grids $s_t=t\frac{2\pi }{2n^{\left(k\right)}}$, t=0,…,2n^(k)-1. Since ⃑r_s is not on the contour L^(k), no singularities are present in such integrals. We now distinguish two cases.

First case is when the shortest distance d^j,k_min between ⃑r_s and the contour L^(k) is relatively small d^j,k_min≲ak, which results in a sharply peaked function Φ(s,s′) with a width at a half maximum ~d^j,k_min/ak (see Fig. 2). As before, we approximate ψ^(k)(s′) with a trigonometric interpolation and obtain integrals of the form (7). In this case, however, no analytical formulas are available for the fourier transforms present in (7) and we evaluate them numerically by performing a Fast Fourier Transform (FFT). The accuracy of FFT should be higher than the accuracy of the trigonometric interpolation in order to achieve the smallest error. In fact, because of a cusp in the Green’s function the number of points used in FFT should exceed ~a_k/^dj,k_min.

Second case is when the shortest distance between ⃑r_s and the contour L^(k) is relatively large d^j,k_min≳a_k, and the functions Φ(s,s′) may be considered as relatively smooth. When a large number of inclusions is present, this becomes the most common case. Here, for the integration of the function ψ^(k)(s′)Φ(s,s′), instead of (7) we apply a much simpler trapezoidal rule:

The error associated with (9) will be still exponentially small [22]. Such a discretization has a minimal cost as only a single evaluation of the kernel function is needed. The only limitation imposed on the number of points n^(k) is that it has to be large enough to resolve all the oscillations in the Green’s functions Φ(s,s′) (see Fig. 2(b)). Particularly, from the functional form of the Green’s function (4) it is straightforward to demonstrate that the number of oscillations in the interval of s-s′∈ (0,2π) is ~ k₀γ2a_k/(2π) ~ γa_k/λ (where λ is a wavelength of light). Therefore, the number of points in the boundary discretization has to be larger than the number of oscillations. In what follows, as an empirical rule, we use the trigonometric interpolation (6) when d^j,k_min<6ak, while we use trapezoidal rule (9) when d ^j,k_min≥6a_k.

When discretizations (7) and (9) are substituted back into (5), matrix equation is obtained:

The vector of unknowns X has 4N_c∑2n^(k) elements which are the values of e^(k)_c (s_t), e^(k)_g (s_t),h^(k)_c (s_t), h^(k)_g (s_t) defined on their proper discrete lattices t=0,…,2n^(k)-1, and k=1,…,N_c. The elements of matrix A(n_e) depend non-linearly on n_e. Modal effective refractive indexes are defined by the values of n_e for which the determinant of A(n_e) is zero. We note here that, as in [15], the size of the unknown vector X and the corresponding matrix A(n_e) could be cut in half by considering only the first two equations of (2) and using the other two to express e_c(s_t) and h_c(s_t) as a function of e_g(s_t) and h_g(s_t) at each inclusion. However, this results in more complex matrix elements.

2.2. Comparison of the code performance with that of a multipole method

The multipole method is not very different from the boundary integral method proposed here. In fact, one can derive the multipole method starting from the integral equations (5). If, as in [20], the coefficients of a fourier expansion of the contour functions e_g(s′) and h_g(s′) are specified and the Graf’s addition theorem is applied twice (once to derive the so-calledWijngaard expansion and once to transform the the origin of the cylindrical waves) the multipole formulation is obtained.

Particularly, the azimuthal variation of the longitudinal electric field in the vicinity of an inclusion is approximated as $E_z\left(s\text{'}\right)={\displaystyle \sum _{m=-M}^M}{c}_m{e}^{ims\text{'}},$, (and similarly for H_z) for some coefficients c_m. Here M is the multipole order and this field expansion has the same accuracy as the trigonometric interpolation (6) if M≃n is assumed. The Graf’s addition theorem involves an infinite series which for consistency must be truncated to the same order M as the field expansion. This induces some error which is independent from the error associated with the field expansion. For well separated inclusions this last error is smaller than the error in the field expansion and the method works well. As the distance between inclusions decreases, the error associated with the truncated Graf’s series increases and at some point it becomes bigger than the error associated with the fields expansion. In such case, the multipole method starts to loose accuracy unless the order of the field expansion is increased accordingly.

We now estimate the numerical cost of the multipole method for the calculation of the modes of a system of N_c identical inclusions of diameter d forming a periodic lattice of pitch λ (see Fig. 1(a), for example). In this case the total number of unknowns is 2N_c(2M+1), and the number of matrix elements Size (A(N_e))~N²_cM². One can show that to construct such a matrix one needs ~N²_cM Bessel function evaluations. In the case when inclusions are well separated from each other, the truncated Graf’s series are rapidly convergent when M ≳ k₀λ_gd. However, when inclusions become too close to each other d→Λ convergence of the truncated Graf’s series is achieved only when M≳d/(Λ-d).

We now estimate the numerical cost of the boundary integral method. Assuming 2N_p point discretization of every boundary, the number of unknowns in a system is 8N_cN_p. The number of matrix elements becomes Size(A(N_e))~N²_cN²_p. One can show that to construct such a matrix one needs ~ (const₁ ·N²_c N_p+const2 ·N_cN_FFT) Bessel function evaluations. In the case when inclusions are well separated from each other, as established in the previous subsection, we only have to resolve all the oscillations in the Green’s function resulting in (N_FFT~N_p) ≳ k₀γgd. When inclusions become too close to each other d→L the FFT order has to be high enough to resolve the cusp in the Green’s function N_FFT ≳ d/(L-d). However, due to the conception of the method, the number of boundary discretization points will still remain small N_p≳k0γ_gd resulting in a considerable performance improvement over the multipole method.

We summarize the performance of the multipole and boundary integral methods in the Table 1. From this table it is clear that both methods show comparable performances when inclusions are well separated. Boundary integral method, however, greatly outperforms multipole method in the case of closely separate inclusion both in terms of memory and simulation time.

Table 1. Performance comparison of the multipole and boundary integral methods.

View Table | View all tables in this article

2.3. Arbitrary shaped inclusions

For the arbitrary shaped inclusions, the discretization procedure is similar to that used for the circular inclusions. Particularly, the trigonometric interpolation is now used for the functions ψ^(k)(s′)=y^(k)(s′)J^(k)(s′):

Assuming that →r_s ∈ L ^(j), discretization of the integrals along the contours L^(k)_k≠j is done exactly as in the previous section, obtaining linear equations in terms of the ψ^(k) (s_t) by employing

either a FTT transform or a trapezoidal integration. Now consider evaluation of the integrals along the contour L^(j). When interpolations (11) are substituted into (5), the following expressions are obtained:

where index a shows that Φ(s,s′) is evaluated along the boundary of an arbitrary shaped inclusion. In this case no analytical formulas are available for the Fourier transforms of Φ^a(s,s′). Moreover, we can not use FFT efficiently because of the singularity at s=s′. To circumvent this problem we introduce a regularization circle following [20]. Particularly, for every jth inclusion, we consider a circle of a comparable diameter a_j as shown schematically in Fig. 2(c). We further distinguish two cases. In the first case, Φ^a(s,s′) represents the Green’s function G^a(s,s′). Then, the integral in (12) can be expressed as:

where index c denotes the regularization circle. Second integral in the righthand side can be evaluated analytically. Function G^a-G^c is not singular any more, and its fourier transform can be evaluated by using FFT. The value of G^a-G^c when s′→s is given in the Appendix A. In the second case, Φ^a(s,s′) represents the normal $\frac{\partial G(s,s\text{'})}{\partial n}$ or the tangential $\frac{\partial G(s,s\text{'})}{\partial \tau }$ derivative of the Green’s function. Then, the integral in (12) can be expressed as:

where Φ^c denotes the corresponding normal or tangential derivative of the Green’s function on the circle, and J^(j)(s) is the jacobian of the inclusion contour at s. Again, the second integral on the righthand side is evaluated analytically while for the non singular function ${\Phi }^a-\frac{a_j}{J^{\left(j\right)}\left(s\right)}{\Phi }^c$ a FFT is performed. The value of ${\Phi }^a-\frac{a_j}{J\left(s\right)}{\Phi }^c$ when s′ → s is given in the Appendix A.

We note that our discretization scheme, although closely related, is significantly different from the one used in [20]. In that work, instead of a trigonometric polynomial, Fourier expansion of the potential densities ψ^(k)(s′) is used. In our case, the matrix elements are given by the single integrals, while in [20] they are given by the double Fourier integrals which are numerically more intensive to evaluate. Finally, the numerically efficient trapezoidal rule (9) used to calculate the majority of the matrix elements in our method, can not be used in [20] as in their case Fourier coefficients, rather than the equidistant values ψ^(k)(s_t), are employed.

2.4. Finding the modes

Finding the propagating modes corresponds to finding values of n_e for which the determinant of A(n_e) in Eq.(10) is zero. This can be done in several ways. One approach is the integral root searching technique [25]. Another approach [15] is to compute a map of the modulus of the determinant over the region of interest and then use the local minima of this map as initial values for further root refinement. We perform a search for the minima of the determinant along a single line in the complex plane and then use these local minima as initial values for further root refinement. If the imaginary component of n_e is expected to be much smaller

than its real component a search line can be typically chosen along the real axis. Finally, as a refinement algorithm we use Newton method. Determinant derivative with respect to n_e is calculated numerically using first order finite difference scheme.

 

Fig. 3. Convergence analysis and comparison with the multipole method for the three simple test structures: (a) six circular holes; diameter d=5µm, pitch L=6.75µm (b) six elliptic holes; axis a=2.5µm b=1.5µm, pitch L=6.75µm (c) six metal coated cylinders; outer diameter d_o=0.8µm, inner one d_i=0.7µm, pitch Λ = 1.5µm.

Download Full Size | PDF

In fact, we consider the determinant only for simple structures with less than 10 inclusions. For the structures with a higher number of inclusions, instead of the determinant, we consider the matrix A(n_e) smallest eigenvalue. When the determinant goes to zero, so does the smallest eigenvalue. We find that the smallest eigen value possess a wider convergence zone compared to that of the determinant. Furthermore, by working with the smallest eigenvalue, one avoids very large values typical for the determinants. Thus, for complex structures, we find the effective refractive index of a mode by performing root searching for the smallest eigenvalue of A(n_e).

An important part of the method is inclusion of symmetries. According to the symmetry, the modes of a fiber are separated into distinct classes. Calculations are then performed separately for every class. Considerable reduction of the overall computational cost is achieved as only a small part of a structure has to be used in a simulation.

3. Study of the code accuracy for the simple test structures

In order to validate the method, we perform convergency analysis and accuracy comparison with a multipole method for the three simple test structures shown in Fig. 3. First, we consider the structure presented in [15] and shown in Fig. 3(a). Is consists of a single ring of six equally spaced circular holes with diameters d=5µm and a pitch Λ = 6.75µm. The glass cladding is assumed to have refractive index of n_g=1.45, while the air holes have n_c=1. The wavelength is λ=1.45µm. In the Table 3 we present convergence study of the effective refractive index of one of the low order modes as calculated by the boundary integral method, as well as comparison with a multipole method [15]. We remark that for the same number of unknowns both methods have similar accuracy. As mentioned earlier, this is to be expected in this case of well separated inclusions. We also remark that in our implementation choosing 24 discretization points per hole results in a ten digit accuracy for n_e; this is to be compared with the results of [14] where 50 points per hole were needed to achieve the same accuracy.

Now, we consider a solid core fiber with all the structural parameters as in the previous case, except with elliptical holes instead of the circular ones (see Fig. 3(b)). The major axis of the ellipses is taken as a=5µm, while the minor one is taken as b=3µm. Convergence data for one of the low order modes of this fiber is presented in [17]. Our results are presented in Table 3, including comparison with [17]. The methods agree well, with a consistency of up to five significant digits in the real part of (n_e) and two significant digits in the imaginary part.

Our next comparison is with a coated MOF studied in [19]. The details of how the boundary integral method equations are modified in such a case are given in Appendix B. The in-

Table 2. Effective refractive index of a mode (of a symmetry class p=1 as defined in [15]) of a solid core MOF featuring one ring of six holes (see Fig. 3(a)).

View Table | View all tables in this article

Table 3. Effective refractive index of a mode of a solid core MOF featuring one ring of six elliptic inclusions (see Fig. 3(b)). The results are for the fundamental mode where the nodal line of the E_z field is horizontal. For the other polarization the value 1.446429072+ 2.9898E-6i is obtained by us and 1.446427235+2.9601E-6i by [17].

View Table | View all tables in this article

clusions are coated with a silver metal having a refractive index n_m=0.433480043837 + i8.70529497278 at λ=1.45µm. Geometry of the structure is shown in Fig. 3(c); it consists of six coated cylinders with outer diameters d_o=0.8µm, inner diameters d_i=0.7µm, and a pitch Λ=1.5µm. The glass cladding is assumed to have refractive index of n_g=1.45, while the air holes have n_c=1. The convergence analysis is presented in Table 4, including comparison with [19]. Again we remark that both methods have almost the same accuracy.

Table 4. Effective refractive index of a mode of a solid core MOF with one ring of six coated holes (see Fig. 3(c)). Results are for the fundamental core guided mode.

View Table | View all tables in this article

 

Fig. 4. Hollow core MOF with 5 rings of holes in the reflector. (a) Dispersion curve of the fundamental mode. (b) Loss as a function of the number of reflector layers.

Download Full Size | PDF

4. Demonstration of the code potential for the study of complex structures

4.1. Loss of the hollow core PCF featuring a large number of reflector layers

As a first example we consider finding the fundamental core guided mode of a hollow-core PCF presented in Fig. 1(a). The PCF consists of five rings of circular holes arranged on a hexagonal lattice and surrounding a hollow core formed by the two missing rings in a fiber center. Overall, there are 120 holes in the cladding. The hole to hole pitch is Λ=2.74µm, the hole diameter is d=0.95Λ, and the core diameter is d_c=2.5d. The glass cladding is assumed to have refractive index of n_g=1.45, while the refractive index of the air holes is n_c=1. Dispersion curve for the fundamental core guided mode of this fiber is shown in Fig. 4(a) (dashed line). The minimum imaginary part, corresponding to the center of the bandgap, is obtained at λ=1.51mm where we find n_e=0.98451599741954+3.434721E-8i. For these calculations, the number of discretization points (2n) per hole is chosen according to the following distribution: for the central hole n=32, for the five rings of holes starting from the inner one we have taken respectively n={16,16,14,12,10}. When n is increased by 2, that is n={34,18,18,16,14,12} the change in the value of n_e equals to 2.6E-9+5.4E-10i, signifying convergence of the small imaginary part. Convergence analysis is also performed for the order of FFT, which must be at least 256 to guarantee the low value of the overall error. This is a relatively high number and it increases the computational cost of the matrix elements, however note that we are dealing with a difficult case where the spacing between the inclusions is small. We also present in Fig. 4(b) radiation loss of the hollow core PCF as a function of the number of rings in the reflector. All the loss calculations are performed at a single wavelength λ=1.51µm. At this wavelength for the case of five rings we also show the |E_z|, |H_z| and S_z of the fundamental mode in the outset of Fig. 4.

 

Fig. 5. Birefringence of the fundamental mode of a PCF with elliptic hollow core. (b) Outset: S_z fluxes for the x and y polarizations of the fundamental mode at λ=1.42µm.

Download Full Size | PDF

4.2. Large birefringence of a hollow elliptical core PCF

Next, we consider the structure presented in Fig. 1(b). It consists of three layers of circular holes arranged on a hexagonal lattice. There are four missing holes at the center of the fiber replaced by a central elliptic core. The hole pitch is λ=2µm, the hole diameter is d=0.9Λ, and the elliptical core has axis a=2.3µm and b=4.6µm. The glass cladding is assumed to have refractive index of n_g=1.45, while refractive index of the air holes is n_c=1. This structure is similar to the structure of a highly birefringent fiber proposed in [26], except for the shape of a central hole. In Fig. 5 we present birefringence of the fundamental mode of this fiber. Interestingly, the birefringence changes its sign around λ=1.41µm. At the outset of this figure we show the S_z fluxes for the x and y polarizations of the fiber fundamental mode at λ=1.42µm. The values for the effective refractive indices at this wavelength are as follows: n^x_e=0.93903355+6.7418E-4i and n^y_e=0.93816250+2.2133E-3i.

4.3. Loss birefringence of a MOF containing metal coated elliptical inclusions

Finally, we consider loss birefringence of the fundamental mode in a MOF containing six elliptical air holes coated with a thin silver layer. The geometry is given in Fig. 1(c). The inclusions are coated with silver. The refractive index of silver is calculated from the interpolation of measured data like in [19]. The glass cladding is assumed to have a refractive index n_g=1.45, while the air holes have n_c=1. The hole to hole pitch is Λ=1.5µm. Six coated elliptical inclusions are described by the outer major axis a_o=0.8µm+δ, b_o=0.8µm-δ and the inner major axis a_i=0.7µm+δ, b_i=0.7µm+δ. In our simulations we use δ=0.04µm, which defines the hole ellipticity of $\overline{\delta }$ =2|a-b|/(a+b)=10%. We now characterize losses of the two fundamental mode polarizations as a function of the wavelength.

When ellipticity parameter is taken to zero δ=0 (circular inclusions), both polarizations are degenerate. Here, loss curve of the fundamental mode is presented as dashed in Fig. 6. The wavelength of maximal loss ~1.41 corresponds to the point of phase matching of a core guided mode with a plasmon propagating on the interface between silver and glass. When ellipticity is introduced, wavelengths of phase matching of a fundamental core guided mode with a plasmon become somewhat different for the two polarizations. For example, for δ=0.04µm, dispersion

curves for the losses of the two polarizations are presented in Fig. 6 (solid curves). Plasmonic resonances for both polarizations are clearly identifiable as maxima in the modal losses. For the x-polarization the maximum of losses is at λ_m=1.419µm, while for the y-polarization it is at λ_m=1.407µm. The corresponding S_z fluxes are shown in the outset of Fig. 6. From the flux distributions it is clear that at the wavelengths of phase matching with a plasmon, core guided mode is well mixed with a plasmonic wave propagating on the silver-glass interface.

 

Fig. 6. Loss dispersion curves for the two polarizations of the fundamental mode of a MOF with one ring of metallized elliptic holes. Outset: S_z fluxes for the x and y polarizations of the fundamental mode at the wavelengths of the two plasmon excitation peaks.

Download Full Size | PDF

In principle, by measuring spectral splitting in the plasmon excitation peaks Δλ_p for the two polarizations of a fundamental core guided mode, one can envision detection of the hole ellipticity δ̄. This principle can be used in pressure sensors. Thus, by starting with a fiber containing circular metallized inclusions and by compressing the fiber uniaxially one will induce ellipticity in the hole structure. Such an ellipticity can then detected by measuring splitting in the plasmon excitation wavelengths. To characterize sensitivity of a pressure sensor we define sensitivity as Sλ[nm]=$S_{\lambda }\left[nm\right]=\frac{\partial \left(\Delta {\lambda }_p\right)}{\partial \overline{\delta }}$, which in our case gives S=120nm. Assuming that 0.1nm shift between two plasmonic peaks can be resolved, ellipticity detection limit is estimated at 8·10^-4.

5. Conclusion

A high performance boundary integral method for the modal analysis of MOFs is presented. The method can treat a large number of arbitrary shaped inclusions with boundaries defined by the individual parametric curves. For circular inclusions, in particular, majority of the calculations are done analytically, ensuring high accuracy and rapid convergence. Both solid core and hollow core fibers can be treated; multilayer (coated) inclusions can be easily accommodated. The method was tested on several simple problems, and it shows excellent agreement with simulations performed by other groups. Moreover, we demonstrated that the numerical cost of the boundary integral method is considerably smaller than that of the multipole method for the case of closely spaced inclusions, which is of particular importance for the case of hollow core photonic crystal fibers. We have established that unlike the multipole method, when spacing between inclusions decreases no convergence problems arise; only the computational cost of some of the matrix elements (the order of FFT) increases. To demonstrate the above mentioned advantages of the boundary integral method we applied it to several challenging problems. First, we studied confinement loss of a core guided mode of a hollow photonic crystal fiber featuring large number of reflector layers (5 layers with 120 holes). Second, modal birefringence of a hollow elliptical core fiber with 3 reflector layers was characterized. Finally, birefringence of the plasmon-coupled core modes of a solid core microstructured fiber with metallized elliptical holes was studied, and application to pressure sensing was proposed.

Appendix A: Normal and tangential derivatives of the Green’s functions

Consider a point →rs(x_s,y_s) on a contour L described by the parametric expressions: x=x(s) and y=y(s), s being a parameter such that 0≤s≤2p. Let →r_s′ (x_s′,y_s′) be an arbitrary point. Omitting the factor $\frac{i}{4}$, we consider derivatives of the Hankel function H⁽¹⁾₀ (k₀γR), where R=|→r_s-→r_s′|. The normal derivative at →r_s is given by:

$\frac{\partial H_0^{\left(1\right)}\left(k_0\gamma R\right)}{\partial n}=-k_0\gamma \frac{\partial R}{\partial n}{H}_1^{\left(1\right)}\left(k_0\gamma R\right).$

Considering that,

$\frac{\partial R}{\partial n}=\cos \left({\overrightarrow{r}}_s-{\overrightarrow{r}}_{s\text{'}},\overrightarrow{n}\right)=\frac{\left({\overrightarrow{r}}_s-{\overrightarrow{r}}_{s\text{'}}\right)·\overrightarrow{n}}{R}$

and since $\overrightarrow{n}=\frac{1}{J\left(s\right)}\left(y_s^{\prime },-x_s^{\prime }\right)$, where the prime denotes the derivative with respect to s, we obtain:

$\frac{\partial R}{\partial n}=\frac{y_s^{\prime }\left(x_s-x_{s\prime }\right)-x_s^{\prime }\left(y_s-y_{s\prime }\right)}{J\left(s\right)R}.$

Thus,

(15)$$\frac{\partial H_0^{\left(1\right)}\left(k_0\gamma R\right)}{\partial n}=k_0\gamma \frac{x_s^{\prime }\left(y_s-y_{s\prime }\right)-y_s^{\prime }\left(x_s-x_{s\prime }\right)}{J\left(s\right)R}{H}_1^{\left(1\right)}\left(k_0\gamma R\right).$$

In the same way, considering that →τ=$\overrightarrow{\tau }=\frac{1}{J\left(s\right)}(x_s^{\prime },y_s^{\prime })$(x′_s,y′_s), we obtain:

(16)$$\frac{\partial H_0^{\left(1\right)}\left(k_0\gamma R\right)}{\partial \tau }=-k_0\gamma \frac{x_s^{\prime }\left(x_s-x_{s\prime }\right)+y_s^{\prime }\left(y_s-y_{s\prime }\right)}{J\left(s\right)R}{H}_1^{\left(1\right)}(k_0\gamma R).$$

When the contour L is a circle with radius a and →r_s′ is a point on that contour, the equations (15) and (16) take the following forms:

(17)$$\frac{\partial H_0^{\left(1\right)}\left(2a{k}_0\gamma \mid \sin \frac{s-s\prime }{2}\mid \right)}{\partial n}=-k_0\gamma \mid \sin \frac{s-s\prime }{2}\mid H_1^{\left(1\right)}\left(2a{k}_0\gamma \mid \sin \frac{s-s\prime }{2}\mid \right),$$

and

(18)$$\frac{\partial H_0^{\left(1\right)}\left(2a{k}_0\gamma \mid \sin \frac{s-s\prime }{2}\mid \right)}{\partial \tau }=-k_0\gamma \frac{\sin \left(s-s\prime \right)}{2\mid \sin \frac{s-s\prime }{2}\mid }{H}_1^{\left(1\right)}\left(2a{k}_0\gamma \mid \sin \frac{s-s\prime }{2}\mid \right).$$

We also note the following limits:

(19)$$\underset{s\prime \to s}{lim}\left[H_0^{\left(1\right)}\left(k_0\gamma R\right)-H_0^{\left(1\right)}\left(2a{k}_0\gamma \mid \sin \frac{s-s\prime }{2}\mid \right)\right]=\frac{2i}{\pi }\ln \frac{J\left(s\right)}{a}$$
$$\underset{s\prime \to s}{lim}\left[\frac{{\partial H}_0^{\left(1\right)}\left(k_0\gamma R\right)}{\partial n}-\frac{a}{J\left(s\right)}\frac{\partial H_0^{\left(1\right)}\left(2a{k}_0\gamma \mid \sin \frac{s-s\prime }{2}\mid \right)}{\partial n}\right]=\frac{1-\kappa \left(s\right)J\left(s\right)}{i\pi J\left(s\right)},$$
$$\underset{s\prime \to s}{lim}\left[\frac{{\partial H}_0^{\left(1\right)}\left(k_0\gamma R\right)}{\partial \tau }-\frac{a}{J\left(s\right)}\frac{\partial H_0^{\left(1\right)}\left(2a{k}_0\gamma \mid \sin \frac{s-s\prime }{2}\mid \right)}{\partial \tau }\right]=0$$

with κ(s) being the curvature of L at s.

 

Fig. 7. Schematic of a coated inclusion.

Download Full Size | PDF

Appendix B: Coated inclusions

Consider the jth inclusion which is coated with a material of the refractive index n_m as shown schematically in Fig. 7. In this case an additional inner contour is present, and four additional contours functions must be specified: e(→r_s) and h(→r_s) at both contours on the coating side. When the boundary conditions are considered at the outer contour of this inclusion, Eqs. (5) should be modified. The first one, corresponding to the continuity of E_z, becomes:

(20)$$\underset{L_o^{(j)}}{\int }{e}_o^{(j)}(s\prime )G_m(s,s\prime )J(s\prime )ds\prime +\underset{L_i^{(j)}}{\int }{e}_i^{(j)}(s\prime )G_m(s,s\prime )J(s\prime )ds\prime =\sum _{k=1}^{N_c}\underset{L_0^{(k)}}{\int }{e}_g^{(k)}(s\prime )G_g(s,s\prime )J(s\prime )ds\prime ,$$

where L^(j)_o denotes the outer contour, L^(j)_i denotes the inner one, e^(j)_o denotes the potential density at the outer contour and e^(j)i denotes the potential density at the inner one (both densities are defined on the coating side). Suppose that the contours are circular. Since →r_s is on the outer boundary, the first integral in (20) is discretized by using the analytical formulas in (8) while the second integral is discretized by performing a FFT.

When the same boundary condition is considered at the inner contour we obtain:

(21)$$\underset{L_i^{(j)}}{\int }{e}_c^{(i)}(s\prime )G_c(s,s\prime )J(s\prime )ds\prime =\underset{L_i^{(j)}}{\int }{e}_i^{(j)}(s\prime )G_m(s,s\prime )J(s\prime )ds\prime +\underset{L_o^{(j)}}{\int }{e}_o^{(j)}(s\prime )G_m(s,s\prime )J(s\prime )ds\prime .$$

Again, for circular contours, the integrals along the inner contour are evaluated with the help of analytical formulas in (8) while the one along the outer contour is discretized by performing a FFT. The rest of the equations in (5), are modified in a similar way.

References and links

1. Matlab implementation of the code is available at http://www.photonics.phys.polymtl.ca/codes.html

2. P. Russell“Photonic crystal fibers,” Science 299, 358–362 (2003). [CrossRef]   [PubMed]  

3. A. Bjarklev, J. Broeng, and A.S. Bjarklev “Photonic crystal fibers,” Kluwer Academic Publishers, Boston, (2003).

4. T.A. Burks, J.C. Knight, and P.S.J. Russell “Endlessly single-mode photonic crystal fibers,” Opt. Lett. 22, 961–963 (1997). [CrossRef]  

5. M.C.J. Large, L. Poladian, G.W. Barton, and M.A. van Eijkelenborg, “Microstructured Polymer Optical Fibres,” Springer, Sydney, (2007)

6. 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]  

7. T.M. Monro, D.J. Richardson, N.G.R. Broderick, and P.J. Bennett “Holey optical fibers: an efficient modal model,” J. Lightwave Technol. 17, 1093–1102 (1999). [CrossRef]  

8. 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]  

9. K. Saitoh and M. Koshiba, “Full-Vectorial Imaginary-Distance Beam PropagationMethod Based on a Finite Element Scheme: Application to Photonic Crystal Fibers,” IEEE J. Quantum Electron. 38, 297 (2002). [CrossRef]  

10. A. Cucinotta, S. Selleri, L. Vincent, and M. Zoboli, “Holey fiber analysis through the finite element method,” IEEE Photon. Technol. Lett. 14, 1530–1532 (2002). [CrossRef]  

11. X. Wang, J. Lou, C. Lu, C. L. Zhao, and W. T Ang, “Modeling of PCF with multiple reciprocity boundary element method,” Opt. Express 12, 961–966 (2004), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-5-961 [CrossRef]   [PubMed]  

12. N. Guan, S. Habu, K. Takenaga, K. Himeno, and A. Wada “Boundary element method for analysis of holey optical fibers,” J. Lightwave Technol. 21, 1787–1792 (2003). [CrossRef]  

13. T. Lu and D. Yevick, “A vectorial boundary element method analysis of integrated optical waveguides,” J. Lightwave Technol. 21, 1793–1807 (2003). [CrossRef]  

14. H. Cheng, W. Crutchfield, M. Doery, and L. Greengard, “Fast, accurate integral equation methods for the analysis of photonic crystal fibers I: Theory,” Opt. Express 12, 3791–3805 (2004),http://www.opticsinfobase.org/abstract.cfm?URI=oe-12-16-3791 [CrossRef]   [PubMed]  

15. T. P. White, B. T. Kuhlmey, R. C. McPhedran, D. Maystre, G. Renversez, C. Martijn de Sterke, and L. C. Botten “Multipole method for microstructured optical fibers. I. Formulation,” J. Opt. Soc. Am. B 19, 2322–2330 (2002). [CrossRef]  

16. B. T. Kuhlmey, T. P. White, G. Renversez, D. Maystre, L. C. Botten, C. Martijn de Sterke, and R. C. McPhedran “Multipole method for microstructured optical fibers. II. Implementation and results,” J. Opt. Soc. Am. B 19, 2331–2340 (2002). [CrossRef]  

17. S. Campbell, R. C. McPhedran, and C. Martijn de Sterke “Differential multipole method for microstructured optical fibers,” J. Opt. Soc. Am. B 21, 1919–1928 (2004). [CrossRef]  

18. M. Skorobogatiy, K. Saitoh, and M. Koshiba, “Coupling between two collinear air-core Bragg fibers,” J. Opt. Soc. Am. B 21, 2095–2101 (2004). [CrossRef]  

19. B. T. Kuhlmey, K. Pathmanandavel, and R. C. McPhedran, “Multipole analysis of photonic crystal fibers with coated inclusions,” Opt. Express 14, 10851–10864 (2006). [CrossRef]   [PubMed]  

20. S. V. Boriskina, T.M. Benson., P. Sewell, and A. I. Nosich “Highly efficient full-vectorial integral equation solution for the bound, leaky and complex modes of dielectric waveguides,” IEEE J. Sel. Top. Quantum Electron. 8, 1225–1231 (2002). [CrossRef]  

21. D. Colton and R. Kress “Integral equation methods in scattering theory,” John Wiley & Sons, New York, (1983).

22. R. Kress “Linear integral equations,” Springer-Verlag, New York, (1989).

23. S. V. Boriskina, P. Sewell, and T. M. Benson “Accurate simulation of two-dimensional optical microcavities with uniquely solvable boundary integral equations and trigonometric Galerkin discretization,” J. Opt. Soc. Am. A 21, 393–402 (2004). [CrossRef]  

24. M. Abramowitz and I. A. Stegun “Handbook of mathematical functions,” Dover, New York, (1965).

25. R. Rodriguez-Berral, F. Mesa, and F. Medina, “Systematic and efficient root finder for computing the modal spectrum of planar layered waveguides,” Int. J. RF Microw. Comp. Eng. 14, 73–83 (2004). [CrossRef]  

26. M. S. Alam, K. Saitoh, and M. Koshiba “High group birefringence in air-core photonic bandgap fibers,” Opt. Lett. 30, 824–826 (2005). [CrossRef]   [PubMed]