Abstract

A boundary integral method [1] for calculating leaky and guided modes of microstructured optical fibers is presented. The method is rapidly converging and can handle a large number of inclusions (hundreds) of arbitrary geometries. Both, solid and hollow core photonic crystal fibers can be treated efficiently. We demonstrate that for large systems featuring closely spaced inclusions the computational intensity of the boundary integral method is significantly smaller than that of the multipole method. This is of particular importance in the case of hollow core band gap guiding fibers. We demonstrate versatility of the method by applying it to several challenging problems.

©2007 Optical Society of America

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 dc=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 ao=0.84µm and bo=0.76mm, the inner hole principal axis are ai=0.74µm and bi=0.66µm, the pitch is Λ=1.5mm.

Download Full Size | PPT Slide | 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 Nc of homogeneous inclusions of refractive index nc embedded into a homogeneous background material of refractive index ng. 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 ω=ck0. Modal effective refractive index is defined as ne=β/k0.

 

Fig. 2. a) Schematic of aMOF cross section. b) Schematic of Re(G(s, s′)). Green’s function has a cusp when ss′ it also exhibits oscillations (γ=ng2ne2).c) Arbitrary shaped inclusion and a corresponding regularization circle.

Download Full Size | PPT Slide | PDF

Given the longitudinal components Ez and Hz 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:

2Ez(c,g)+k02γc,g2Ez(c,g)=0
2Hz(c,g)+k02γc,g2Hz(c,g)=0,

where γ2c,g=n2c,g-n2e. Moreover, on the boundaries between various dielectric regions the longitudinal and tangential components of the fields are continuous:

Ez(c)=Ez(g)
Hz(c)=Hz(g)
Et(c)=Et(g)ik0γc2(neEz(c)τHz(c)n)=ik0γg2(neEz(g)τHz(g)n),
Ht(c)=Ht(g)ik0γc2(neHz(c)τ+εcEz(c)n)=ik0γg2(neHz(g)τ+εgEz(g)n)

where τ is the tangential derivative to the boundary contour L,n is the outer normal derivative, and ε(c,g) is the dielectric constant of either the inclusion or the cladding. Throughout the paper Hz represents the true magnetic field scaled by the free space impedance Hz = μ0cHtrue fieldz.

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]:

Ez(r)=Le(rs)G(r,rs)dls
Hz(r)=Lh(rs)G(r,rs)dls.

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:

G(r,rs)=i4H0(1)(k0γrrs),

where H(1)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=k=1NcL(k).. 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=xk(s) and y=yk(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⃑rsL (j):

1:02πec(j)(s)Gc(s,s)J(j)(s)ds=k=1Nc02πeg(k)(s)Gg(s,s)J(k)(s)ds
2:02πhc(j)(s)Gc(s,s)J(j)(s)ds=k=1Nc02πhg(k)(s)Gg(s,s)J(k)(s)ds
3:1γc2(ne02πec(j)(s)Gc(s,s)τJ(j)(s)ds02πhc(j)(s)Gc(s,s)τJ(j)(s)dshc(j)(s)2)=
1γc2(nek=1Nc02πeg(k)(s)Gg(s,s)τJ(k)(s)dsk=1Nc02πhg(k)(s)Gg(s,s)τJ(k)(s)ds+hc(j)(s)2),
4:1γc2(ne02πhc(j)(s)Gc(s,s)τJ(j)(s)ds+εc02πec(j)(s)Gc(s,s)nJ(j)(s)ds+εcec(j)(s)2)=
1γg2(nek=1Nc02πhg(k)(s)Gg(s,s)τJ(k)(s)ds+εgk=1Nc02πeg(k)(s)Gc(s,s)nJ(k)(s)dsεgec(j)(s)2)

where J(k)(s)=[(dxkds)2+(dykds)2]12 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…Nc. The value of ne, for which a nontrivial solution of (5) exists, defines the effective refractive index of the fiber mode. Note, that the terms e(s)2andh(s)2 appear because of the normal derivatives of the integrals in (3) (the so called double layer potentials) exhibit discontinuities when rsL(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)=ak, where ak 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 = t2π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′):

ψ(k)(s)=t=02n(k)1(12n(k)m=n(k)n(k)1eim(sst))ψ(k)(st).

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⃑rs′=⃑rs. When interpolations of the form (6) are inserted into the integrals along the contour L(j) in (5), these integrals take the following form:

I(j)=02πψ(j)(s)Φ(s,s)J(j)(s)dst=02n(j)1(12πm=n(j)n(j)1eimst02πeimsΦ(s,s)ajds)ψ(j)(st),

where Φ(s,s′) stands for G(s,s′),G(s,s')n or G(s,s')τ. 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]:

02πeimsG(s,s)ds'=iπ2Jm(k0γaj)Hm(1)(k0γaj)eims
02πeimsG(s,s')nds'=[12aj+ik0γπ2Jm(k0γaj)Hm(1)(k0γaj)eims]eims,
02πeimsG(s,s')nds'=mπ2ajJm(k0γaj)Hm(1)(k0γaj)eims

where Jm() and H(1)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 (1ajdds) 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 st=t2π2n(k), t=0,…,2n(k)-1. Since ⃑rs 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 dj,kmin between ⃑rs and the contour L(k) is relatively small dj,kminak, which results in a sharply peaked function Φ(s,s′) with a width at a half maximum ~dj,kmin/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 ~ak/dj,kmin.

Second case is when the shortest distance between ⃑rs and the contour L(k) is relatively large dj,kminak, 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:

I(k)=02πψ(k)(s)Φ(s,s)J(k)(s)dst=02n(k)1(ak2π2n(k)Φ(s,st))ψ(k)(st).

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 ~ k0γ2ak/(2π) ~ γak (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 dj,kmin<6ak, while we use trapezoidal rule (9) when d j,kmin≥6ak.

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

A(ne)·X=0.

The vector of unknowns X has 4Nc∑2n(k) elements which are the values of e(k)c (st), e(k)g (st),h(k)c (st), h(k)g (st) defined on their proper discrete lattices t=0,…,2n(k)-1, and k=1,…,Nc. The elements of matrix A(ne) depend non-linearly on ne. Modal effective refractive indexes are defined by the values of ne for which the determinant of A(ne) is zero. We note here that, as in [15], the size of the unknown vector X and the corresponding matrix A(ne) could be cut in half by considering only the first two equations of (2) and using the other two to express ec(st) and hc(st) as a function of eg(st) and hg(st) 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 eg(s′) and hg(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 Ez(s')=m=MMcmeims',, (and similarly for Hz) for some coefficients cm. Here M is the multipole order and this field expansion has the same accuracy as the trigonometric interpolation (6) if Mn 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 Nc 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 2Nc(2M+1), and the number of matrix elements Size (A(Ne))~N2cM2. One can show that to construct such a matrix one needs ~N2cM Bessel function evaluations. In the case when inclusions are well separated from each other, the truncated Graf’s series are rapidly convergent when Mk0λgd. However, when inclusions become too close to each other d→Λ convergence of the truncated Graf’s series is achieved only when Md/(Λ-d).

We now estimate the numerical cost of the boundary integral method. Assuming 2Np point discretization of every boundary, the number of unknowns in a system is 8NcNp. The number of matrix elements becomes Size(A(Ne))~N2cN2p. One can show that to construct such a matrix one needs ~ (const1 ·N2c Np+const2 ·NcNFFT) 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 (NFFT~Np) ≳ k0γ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 NFFTd/(L-d). However, due to the conception of the method, the number of boundary discretization points will still remain small Npk0γ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.

Tables Icon

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

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′):

Ψ(k)(s)=t=02n(k)1(12n(k)m=n(k)n(k)1eim(sst))Ψ(k)(st).

Assuming that →rsL (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) (st) 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:

I(j)=02πΨ(j)(s)Φ(s,s)dst=02n(j)1(12n(j)m=n(j)n(j)1eimst02πeimsΦa(s,s)ds)Ψ(j)(st),

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 aj as shown schematically in Fig. 2(c). We further distinguish two cases. In the first case, Φa(s,s′) represents the Green’s function Ga(s,s′). Then, the integral in (12) can be expressed as:

02πeimsGa(s,s)ds=02πeims[Ga(s,s')Gc(s,s)]ds+02πeimsGc(s,s)ds,

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

02πeimsΦa(s,s)ds=02πeims[Φa(s,s')ajJ(j)(s)Φc(s,s)]ds+ajJ(j)(s)02πeimsΦc(s,s)ds,

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 ΦaajJ(j)(s)Φc a FFT is performed. The value of ΦaajJ(s)Φ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)(st), are employed.

2.4. Finding the modes

Finding the propagating modes corresponds to finding values of ne for which the determinant of A(ne) 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 ne 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 ne 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 do=0.8µm, inner one di=0.7µm, pitch Λ = 1.5µm.

Download Full Size | PPT Slide | 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(ne) 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(ne).

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 ng=1.45, while the air holes have nc=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 ne; 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 (ne) 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-

Tables Icon

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)).

Tables Icon

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

clusions are coated with a silver metal having a refractive index nm=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 do=0.8µm, inner diameters di=0.7µm, and a pitch Λ=1.5µm. The glass cladding is assumed to have refractive index of ng=1.45, while the air holes have nc=1. The convergence analysis is presented in Table 4, including comparison with [19]. Again we remark that both methods have almost the same accuracy.

Tables Icon

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.

 

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 | PPT Slide | 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 dc=2.5d. The glass cladding is assumed to have refractive index of ng=1.45, while the refractive index of the air holes is nc=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 ne=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 ne 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 |Ez|, |Hz| and Sz 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: Sz fluxes for the x and y polarizations of the fundamental mode at λ=1.42µm.

Download Full Size | PPT Slide | 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 ng=1.45, while refractive index of the air holes is nc=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 Sz 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: nxe=0.93903355+6.7418E-4i and nye=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 ng=1.45, while the air holes have nc=1. The hole to hole pitch is Λ=1.5µm. Six coated elliptical inclusions are described by the outer major axis ao=0.8µm+δ, bo=0.8µm-δ and the inner major axis ai=0.7µm+δ, bi=0.7µm+δ. In our simulations we use δ=0.04µm, which defines the hole ellipticity of δ¯ =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 Sz 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: Sz fluxes for the x and y polarizations of the fundamental mode at the wavelengths of the two plasmon excitation peaks.

Download Full Size | PPT Slide | 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λ[nm]=(Δλp)δ¯, 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(xs,ys) 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 →rs′ (xs′,ys′) be an arbitrary point. Omitting the factor i4, we consider derivatives of the Hankel function H(1)0 (k0γR), where R=|→rs-→rs′|. The normal derivative at →rs is given by:

H0(1)(k0γR)n=k0γRnH1(1)(k0γR).

Considering that,

Rn=cos(rsrs',n)=(rsrs')·nR

and since n=1J(s)(ys,xs), where the prime denotes the derivative with respect to s, we obtain:

Rn=ys(xsxs)xs(ysys)J(s)R.

Thus,

H0(1)(k0γR)n=k0γxs(ysys)ys(xsxs)J(s)RH1(1)(k0γR).

In the same way, considering that →τ=τ=1J(s)(xs,ys)(xs,ys), we obtain:

H0(1)(k0γR)τ=k0γxs(xsxs)+ys(ysys)J(s)RH1(1)(k0γR).

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

H0(1)(2ak0γsinss2)n=k0γsinss2H1(1)(2ak0γsinss2),

and

H0(1)(2ak0γsinss2)τ=k0γsin(ss)2sinss2H1(1)(2ak0γsinss2).

We also note the following limits:

limss[H0(1)(k0γR)H0(1)(2ak0γsinss2)]=2iπlnJ(s)a
limss[H0(1)(k0γR)naJ(s)H0(1)(2ak0γsinss2)n]=1κ(s)J(s)iπJ(s),
limss[H0(1)(k0γR)τaJ(s)H0(1)(2ak0γsinss2)τ]=0

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

 

Fig. 7. Schematic of a coated inclusion.

Download Full Size | PPT Slide | PDF

Appendix B: Coated inclusions

Consider the jth inclusion which is coated with a material of the refractive index nm as shown schematically in Fig. 7. In this case an additional inner contour is present, and four additional contours functions must be specified: e(→rs) and h(→rs) 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 Ez, becomes:

Lo(j)eo(j)(s)Gm(s,s)J(s)ds+Li(j)ei(j)(s)Gm(s,s)J(s)ds=k=1NcL0(k)eg(k)(s)Gg(s,s)J(s)ds,

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 →rs 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:

Li(j)ec(i)(s)Gc(s,s)J(s)ds=Li(j)ei(j)(s)Gm(s,s)J(s)ds+Lo(j)eo(j)(s)Gm(s,s)J(s)ds.

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]  

References

  • View by:
  • |
  • |
  • |

  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]

2006 (1)

2005 (1)

2004 (6)

2003 (3)

2002 (5)

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]

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]

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]

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]

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]

2000 (1)

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]

1999 (2)

1997 (1)

Abramowitz, M.

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

Alam, M. S.

Andres, M.V.

Andres, P.

Ang, W. T

Barton, G.W.

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

Bennett, P.J.

Benson, T. M.

Benson., T.M.

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]

Bjarklev, A.

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

Bjarklev, A.S.

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

Boriskina, S. V.

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]

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]

Botten, L. C.

Brechet, F.

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]

Broderick, N.G.R.

Broeng, J.

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

Burks, T.A.

Campbell, S.

Cheng, H.

Colton, D.

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

Crutchfield, W.

Cucinotta, A.

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]

Doery, M.

Ferrando, A.

Greengard, L.

Guan, N.

Habu, S.

Himeno, K.

Knight, J.C.

Koshiba, M.

Kress, R.

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

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

Kuhlmey, B. T.

Large, M.C.J.

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

Lou, J.

Lu, C.

Lu, T.

Marcou, J.

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]

Martijn de Sterke, C.

Maystre, D.

McPhedran, R. C.

Medina, F.

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]

Mesa, F.

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]

Miret, J.J.

Monro, T.M.

Nosich, A. I.

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]

Pagnoux, D.

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]

Pathmanandavel, K.

Poladian, L.

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

Renversez, G.

Richardson, D.J.

Rodriguez-Berral, R.

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]

Roy, P.

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]

Russell, P.

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

Russell, P.S.J.

Saitoh, K.

Selleri, S.

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]

Sewell, P.

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]

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]

Silvestre, E.

Skorobogatiy, M.

Stegun, I. A.

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

Takenaga, K.

van Eijkelenborg, M.A.

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

Vincent, L.

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]

Wada, A.

Wang, X.

White, T. P.

Yevick, D.

Zhao, C. L.

Zoboli, M.

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]

IEEE J. Quantum Electron. (1)

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]

IEEE J. Sel. Top. Quantum Electron. (1)

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]

IEEE Photon. Technol. Lett. (1)

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]

Int. J. RF Microw. Comp. Eng. (1)

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]

J. Lightwave Technol. (3)

J. Opt. Soc. Am. A (1)

J. Opt. Soc. Am. B (4)

Opt. Express (3)

Opt. Fiber Technol. (1)

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]

Opt. Lett. (3)

Science (1)

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

Other (6)

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

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

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

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

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

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

Cited By

OSA participates in Crossref's Cited-By Linking service. Citing articles from OSA journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (7)

Fig. 1.
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 dc =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 ao =0.84µm and bo =0.76mm, the inner hole principal axis are ai =0.74µm and bi =0.66µm, the pitch is Λ=1.5mm.
Fig. 2.
Fig. 2. a) Schematic of aMOF cross section. b) Schematic of Re(G(s, s′)). Green’s function has a cusp when ss′ it also exhibits oscillations ( γ = n g 2 n e 2 ) . c ) Arbitrary shaped inclusion and a corresponding regularization circle.
Fig. 3.
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 do =0.8µm, inner one di =0.7µm, pitch Λ = 1.5µm.
Fig. 4.
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.
Fig. 5.
Fig. 5. Birefringence of the fundamental mode of a PCF with elliptic hollow core. (b) Outset: Sz fluxes for the x and y polarizations of the fundamental mode at λ=1.42µm.
Fig. 6.
Fig. 6. Loss dispersion curves for the two polarizations of the fundamental mode of a MOF with one ring of metallized elliptic holes. Outset: Sz fluxes for the x and y polarizations of the fundamental mode at the wavelengths of the two plasmon excitation peaks.
Fig. 7.
Fig. 7. Schematic of a coated inclusion.

Tables (4)

Tables Icon

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

Tables Icon

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)).

Tables Icon

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

Tables Icon

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.

Equations (35)

Equations on this page are rendered with MathJax. Learn more.

2 E z ( c , g ) + k 0 2 γ c , g 2 E z ( c , g ) = 0
2 H z ( c , g ) + k 0 2 γ c , g 2 H z ( c , g ) = 0 ,
E z ( c ) = E z ( g )
H z ( c ) = H z ( g )
E t ( c ) = E t ( g ) i k 0 γ c 2 ( n e E z ( c ) τ H z ( c ) n ) = i k 0 γ g 2 ( n e E z ( g ) τ H z ( g ) n ) ,
H t ( c ) = H t ( g ) i k 0 γ c 2 ( n e H z ( c ) τ + ε c E z ( c ) n ) = i k 0 γ g 2 ( n e H z ( g ) τ + ε g E z ( g ) n )
E z ( r ) = L e ( r s ) G ( r , r s ) d l s
H z ( r ) = L h ( r s ) G ( r , r s ) d l s .
G ( r , r s ) = i 4 H 0 ( 1 ) ( k 0 γ r r s ) ,
1 : 0 2 π e c ( j ) ( s ) G c ( s , s ) J ( j ) ( s ) d s = k = 1 N c 0 2 π e g ( k ) ( s ) G g ( s , s ) J ( k ) ( s ) d s
2 : 0 2 π h c ( j ) ( s ) G c ( s , s ) J ( j ) ( s ) d s = k = 1 N c 0 2 π h g ( k ) ( s ) G g ( s , s ) J ( k ) ( s ) d s
3 : 1 γ c 2 ( n e 0 2 π e c ( j ) ( s ) G c ( s , s ) τ J ( j ) ( s ) d s 0 2 π h c ( j ) ( s ) G c ( s , s ) τ J ( j ) ( s ) ds h c ( j ) ( s ) 2 ) =
1 γ c 2 ( n e k = 1 N c 0 2 π e g ( k ) ( s ) G g ( s , s ) τ J ( k ) ( s ) d s k = 1 N c 0 2 π h g ( k ) ( s ) G g ( s , s ) τ J ( k ) ( s ) d s + h c ( j ) ( s ) 2 ) ,
4 : 1 γ c 2 ( n e 0 2 π h c ( j ) ( s ) G c ( s , s ) τ J ( j ) ( s ) d s + ε c 0 2 π e c ( j ) ( s ) G c ( s , s ) n J ( j ) ( s ) ds + ε c e c ( j ) ( s ) 2 ) =
1 γ g 2 ( n e k = 1 N c 0 2 π h g ( k ) ( s ) G g ( s , s ) τ J ( k ) ( s ) d s + ε g k = 1 N c 0 2 π e g ( k ) ( s ) G c ( s , s ) n J ( k ) ( s ) d s ε g e c ( j ) ( s ) 2 )
ψ ( k ) ( s ) = t = 0 2 n ( k ) 1 ( 1 2 n ( k ) m = n ( k ) n ( k ) 1 e i m ( s s t ) ) ψ ( k ) ( s t ) .
I ( j ) = 0 2 π ψ ( j ) ( s ) Φ ( s , s ) J ( j ) ( s ) d s t = 0 2 n ( j ) 1 ( 1 2 π m = n ( j ) n ( j ) 1 e i m s t 0 2 π e i m s Φ ( s , s ) a j d s ) ψ ( j ) ( s t ) ,
0 2 π e i m s G ( s , s ) d s ' = i π 2 J m ( k 0 γ a j ) H m ( 1 ) ( k 0 γ a j ) e i m s
0 2 π e i m s G ( s , s ' ) n d s ' = [ 1 2 a j + i k 0 γ π 2 J m ( k 0 γ a j ) H m ( 1 ) ( k 0 γ a j ) e i m s ] e i m s ,
0 2 π e i m s G ( s , s ' ) n d s ' = m π 2 a j J m ( k 0 γ a j ) H m ( 1 ) ( k 0 γ a j ) e i m s
I ( k ) = 0 2 π ψ ( k ) ( s ) Φ ( s , s ) J ( k ) ( s ) d s t = 0 2 n ( k ) 1 ( a k 2 π 2 n ( k ) Φ ( s , s t ) ) ψ ( k ) ( s t ) .
A ( n e ) · X = 0 .
Ψ ( k ) ( s ) = t = 0 2 n ( k ) 1 ( 1 2 n ( k ) m = n ( k ) n ( k ) 1 e i m ( s s t ) ) Ψ ( k ) ( s t ) .
I ( j ) = 0 2 π Ψ ( j ) ( s ) Φ ( s , s ) d s t = 0 2 n ( j ) 1 ( 1 2 n ( j ) m = n ( j ) n ( j ) 1 e i m s t 0 2 π e i m s Φ a ( s , s ) d s ) Ψ ( j ) ( s t ) ,
0 2 π e i m s G a ( s , s ) d s = 0 2 π e i m s [ G a ( s , s ' ) G c ( s , s ) ] d s + 0 2 π e i m s G c ( s , s ) d s ,
0 2 π e i m s Φ a ( s , s ) d s = 0 2 π e i m s [ Φ a ( s , s ' ) a j J ( j ) ( s ) Φ c ( s , s ) ] d s + a j J ( j ) ( s ) 0 2 π e i m s Φ c ( s , s ) d s ,
H 0 ( 1 ) ( k 0 γ R ) n = k 0 γ x s ( y s y s ) y s ( x s x s ) J ( s ) R H 1 ( 1 ) ( k 0 γ R ) .
H 0 ( 1 ) ( k 0 γ R ) τ = k 0 γ x s ( x s x s ) + y s ( y s y s ) J ( s ) R H 1 ( 1 ) ( k 0 γ R ) .
H 0 ( 1 ) ( 2 a k 0 γ sin s s 2 ) n = k 0 γ sin s s 2 H 1 ( 1 ) ( 2 a k 0 γ sin s s 2 ) ,
H 0 ( 1 ) ( 2 a k 0 γ sin s s 2 ) τ = k 0 γ sin ( s s ) 2 sin s s 2 H 1 ( 1 ) ( 2 a k 0 γ sin s s 2 ) .
lim s s [ H 0 ( 1 ) ( k 0 γ R ) H 0 ( 1 ) ( 2 a k 0 γ sin s s 2 ) ] = 2 i π ln J ( s ) a
lim s s [ H 0 ( 1 ) ( k 0 γ R ) n a J ( s ) H 0 ( 1 ) ( 2 a k 0 γ sin s s 2 ) n ] = 1 κ ( s ) J ( s ) i π J ( s ) ,
lim s s [ H 0 ( 1 ) ( k 0 γ R ) τ a J ( s ) H 0 ( 1 ) ( 2 a k 0 γ sin s s 2 ) τ ] = 0
L o ( j ) e o ( j ) ( s ) G m ( s , s ) J ( s ) d s + L i ( j ) e i ( j ) ( s ) G m ( s , s ) J ( s ) d s = k = 1 N c L 0 ( k ) e g ( k ) ( s ) G g ( s , s ) J ( s ) d s ,
L i ( j ) e c ( i ) ( s ) G c ( s , s ) J ( s ) d s = L i ( j ) e i ( j ) ( s ) G m ( s , s ) J ( s ) d s + L o ( j ) e o ( j ) ( s ) G m ( s , s ) J ( s ) d s .

Metrics