1 Introduction

Photonic crystal fibres are fibres in which the light is guided by a periodic array of air holes in a glass matrix surrounding the core. While in early work the core was characterized by the absence of an air hole, perhaps the most intriguing type of fibre is that for which the the core consists of an air hole that is larger than the others [1, 2]. In some of the modes of these fibres the energy is mostly in the air core, which immediately suggests a number of applications, including low-loss propagation, for example at wavelengths where the glass absorbs.

The guiding mechanism of the light in air-core fibres is Bragg reflection due to the periodicity of the airholes outside the core. This is a subtle effect which complicates modal calculations. To see this, recall that the existence of a mode requires the presence of a mirror, while in addition some phase condition needs to be satisfied. In conventional germanium-doped fibres the mirror is due to total internal reflection at the core-cladding interface, and the modal calculation essentially involves satisfying the phase condition. The same is true in photonic crystal fibres with a solid core. In contrast, in air-core fibres, where the mirrors are associated with Bragg reflection, a modal calculation needs to be preceded by a separate band structure calculation of the cladding. Once the regions of high-reflection are then identified, the modal calculation can proceed.

Different methods for photonic crystal fibre calculations have been discussed in the literature [2, 3, 4, 5, 6]. The most common of these is the plane wave method [2, 4], in which the field is written as a superposition of plane waves. Similarly, the calculation of the band structure of the cladding uses a plane-wave expansion. Recently, we developed a multipole method for field calculations of fibres with inclusions. In this method the field is written in terms of cylindrical harmonics centered around each hole, with consistency being enforced using a field identity [7]. The method has advantages in terms of speed and accuracy, though in its present implementation it can only deal with circular holes. In contrast to the multipole method, the plane wave method is presently implemented using periodic boundary conditions and so losses due to a finite cladding cannot be calculated in this way. The features of the multipole method are not limited to photonic crystal fibres. Indeed, it was originally developed for general diffractive structures and photonic crystals [8, 9]. Here we show its development for determining the cladding band structure, and the subsequent calculation of the air-core fibre modes, building on the multipole method with a scattering matrix approach and Bloch’s theorem [10, 11].

2 Theoretical Formulation

We commence with an outline of the multipole theory of propagation in photonic crystal fibres and proceed to relate this to the corresponding properties of the cladding. Our model [7] of the fibre, which is aligned with the z-axis, provides for a finite number (N_c ) of circular holes of refractive index n_i embedded in a finite matrix of index n_e bounded by a circular jacket, the presence of which is convenient for various applications.

Because of strong refractive index contrast of the materials, a full vector field formulation is required, necessitating field expansions for the longitudinal (z) components of both electric and magnetic fields, from which all other components may be derived. Assuming a longitudinal field dependence of exp(iβz), each of E_z and H_z satisfies a transverse (xy) Helmholtz equation (${\nabla }_{\mathit{\text{xy}}}^2$ +κ ²)V=0 with $\kappa =k_{\perp }^e=\sqrt{k^2{n}_e^2-{\beta }^2}$ in the matrix, where k is the free space wavenumber, and a corresponding form with $\kappa =k_{\perp }^i=\sqrt{k^2{n}_i^2-{\beta }^2}$ in the holes. In the vicinity of each cylinder l, we prescribe local expansions in the cylindrical harmonic basis. Using the local coordinate system centred at c_l , the exterior form of E_z is

and correspondingly for the magnetic field. Such local expansions are valid only in an annulus extending to the perimeter of the nearest scatterer, thus necessitating the introduction of a global expansion that is valid throughout the matrix

and from which a field identity determining the source coefficients ($B_m^l$ ) can be derived.

Representation (2) contains outgoing waves sourced at each cylinder (l) and a standing wave term generated at the jacket boundary, identified by superscript 0. The derivation of the field identity [7] follows from the observation that the regular, or non-singular, part of the field, characterized by coefficients A ^El =[$A_m^{\mathit{\text{El}}}$ ],¹ in the vicinity of each cylinder l arises from outgoing contributions (due to B ^Ej =[$B_m^{\mathit{\text{Ej}}}$ ]) from all other cylinders j and the standing wave originating on the jacket, A^0l=[$A_m^{E\mathit{0}}$ ]. In matrix form,

with c_lj =c_j -c_l , and with an analogous form for H_z . In (3), the matrices 𝓗 ^lj and 𝓙 ^l0 arise through an application of Graf’s addition theorem for Bessel functions [17] and may be regarded as change of basis transformations. For example, the matrix 𝓗 ^lj converts outgoing waves (i.e. Hankel function terms) sourced at the centre of cylinder j to standing waves (i.e. regular Bessel function series) in the coordinate system of cylinder l. Correspondingly, 𝓙 ^l0 transforms standing waves from the jacket frame of reference to standing waves in the coordinate frame of cylinder l. There exists a further identity

which together with its magnetic field analogue, expresses the total outgoing field in the vicinity of the jacket in terms of contributions from all cylinders. It involves changes of basis 𝓙^0l that convert outgoing waves from the coordinate system of cylinder l to that of the jacket.

Each of the identities (3) and (5) applies individually to single field components. However, cross-coupling of fields occurs through the application of the boundary conditions at the surface of each cylinder and the jacket. It is convenient [7] to write the boundary conditions, expressing the continuity of the tangential electric and magnetic field components, in terms of cylindrical harmonic reflection coefficients. For each cylinder l we may write

Only for in-plane incidence (β=0) do R^EH and R^HE vanish, thereby decoupling the problem into two distinct polarizations, each of which can be solved as a scalar problem. In all other cases, a full vector formulation is required and it is convenient to mirror this requirement in the notation. We thus introduce the partitioned vector notation ${\tilde{\mathbf{B}}}^l={\left[{\mathbf{B}}^{{El}^T}{\mathbf{B}}^{{Hl}^T}\right]}^T$, with the superscript T denoting the transpose. This represents the outgoing electromagnetic field for each cylinder l and we combine these in the partitioned vector 𝓑=[B ^̃l] for all cylinders. Similarly, we introduce ${\tilde{\mathbf{A}}}^l={\left[{\mathbf{A}}^{{El}^T}{\mathbf{A}}^{{Hl}^T}\right]}^T$ and 𝓐=[A ^̃l] for the regular field components and also block partitioned matrices

that enable the exterior forms of the boundary conditions at cylinder interfaces to be expressed in the compact form 𝓑=𝓡𝓐. Similarly, the boundary condition at the jacket boundary may be written as Ã=R̃⁰ B̃, involving the interior reflection coefficient R̃⁰.

In (8), 𝓙 ^B0 denotes a column vector containing matrices 𝓙 ^l0 for l=1, 2, …,N_c , while 𝓙 ^0B denotes a row vector containing matrices 𝓙 ^0l.

Manipulation of (8) and the boundary conditions then yields the field identity

In (9), the scattering operator 𝓢 comprises two essential elements: (a) $\widetilde{𝓗}$ that contains partitions 𝓗 ^lj and describes direct cylinder to cylinder interactions, and (b) 𝓙 ̃^B0R̃⁰ J $\widetilde{𝓙}$ ^0B that describes all indirect interactions between cylinders that take place via reflections (R̃₀) at the jacket interface. The matrix 𝓜=𝓜(k, β) when singular, defines a mode of the photonic crystal fibre, with the corresponding null space vector 𝓑 characterizing the associated fields. Details of the implementation of the computational methods and the search algorithm for modes are described elsewhere [7].

Locating modes of photonic crystal fibres is a numerically intensive and delicate task, guided by the requirement that they be confined to complete band gaps of the cladding. The task thus reduces to the location of fibre modes within cladding band gaps. The additional requirement that the energy be mostly confined to the central air hole can further refine the search region by requiring that the mode be located close to the light line (i.e., β≈k).

For a well confined photonic crystal fibre with a large number of holes, the jacket is negligible and (9) is well approximated by (I-𝓡𝓗)𝓑=0. For a periodic lattice, the Bloch condition, namely that $V\left(\mathbf{r}+{\mathbf{c}}_j\right)=e^{i{\mathbf{k}}_0{\mathbf{c}}_j}V\left(\mathbf{r}\right)$ where k₀=(k0_x, k0_y ) denotes the Bloch vector, relates the field coefficients of every cylinder through the quasiperiodicity constraint ${\tilde{\mathbf{B}}}^j=\tilde{\mathbf{B}}{e}^{i{\mathbf{k}}_0{\mathbf{c}}_j}$. Here, B̃ denotes multipole field coefficients associated with a central cell, centred at the coordinate origin. The field identity (9) thus reduces to

In (10), R̃ is the exterior reflection matrix given in (6) for any of the identical cylinders and S ^A =[$S_{\mathit{\text{nm}}}^A$ ]=[$S_{n\mathit{-}m}^A$ ] is a Toeplitz matrix whose bands contains array lattice sums [13]. These follow from the quasiperiodicity condition and the structure of the field identity, and are given by

where the summation excludes the central lattice point. Lattice sums underpin the application of multipole methods to periodic structures and a description of them, and techniques for their computation, may be found in the review of McPhedran et al. [13].

Returning to (10), the matrix M̃=M̃(k, k ₀ , β) and thus the location of complete band gaps requires the scanning of a subset of k-β parameter space, locating all modes by searching for (effectively) zero determinants of M̃ for all possible orientations k ₀. In practice, this task is undertaken by traversing the perimeter of the irreducible part of the first Brillouin zone of the particular lattice. This method has been successfully used in our earlier computation of in-plane (β=0) band diagrams [14] for 2D photonic crystals. However, in the present three-dimensional case, the additional requirement of scanning over a β interval is a computationally demanding task.

Instead, it is advantageous to regard the lattice as a stack of one-dimensional diffraction gratings, a reformulation that is equivalent [11, 15] to the array multipole treatment. The method relies on the use of plane wave scattering matrices that characterize diffraction by an elemental cylinder grating, and the subsequent formulation of an eigenvalue problem from which we may compute the modes of the photonic crystal cladding.

The derivation of the scattering matrices outlined below is the conical diffraction generalization [16] of the in-plane formulation [8, 9]. We consider a single grating of period d that gives rise to a set of upward and downward propagating plane waves exp i(α_sx±χ_sy) where α_s =2πs/d+k_0x and $\chi s=\sqrt{k_{\perp }^{e^2}-{\alpha }_s^2}$. In what follows, plane wave fields will be characterized by vectors f of plane wave coefficients that comprise both TE($f_s^E$ ) and TM($f_s^H$ ) entries. The plane wave expansions

represent the transverse electric and magnetic plane wave fields above the grating (Fig. 1), with δ- denoting an arbitrary incoming field and f ⁺ denoting an upward, outgoing field. An analogous expansion for the field below the grating may also be written, this time in terms of an incoming field δ ⁺ and an outgoing field f ^-. In (12) and (13), ξ_s =χ_s/k, while the TM and TE plane wave modes are given respectively by ${\mathbf{R}}_s^H$ =q _s /q_s exp(iq _s ·r), ${\mathbf{R}}_s^E$ =ŷ×q _s /q_s exp(iq _s ·r), with q _s =(α_s, β) and r=(x, z).

Around each cylinder of the grating, fields are again expressed in a local, multipole series (1). Proceeding in an analogous manner to that described above, we are led to an identity for the field coefficients Ã=S̃ ^G B̃+Ã^ext, which is the analogue of (8) for gratings. Here, S̃ ^G is a matrix of lattice sums for the grating, the definition of which follows from that of the array sum (11) by restricting the summation to the central line of cylinders [8, 9]. The exterior driving field Ã^ext is due to the incoming plane waves δ ^± and it is straightforward to show that Ã^ext=J̃+ξ ^δ̃⁺+J̃-$\tilde{\xi }$ $\tilde{\delta }$ -, where $\tilde{\xi }$ is a transformation that converts TE/TM plane wave coefficients to the E_z-H_z system, in which the multipole scattering problem is most naturally expressed. The matrices J ^̃± perform changes of basis converting upward and downward going plane waves into series of cylindrical harmonic terms, with J ^̃+=diag [J ⁺, J ⁺], J ⁺=[$J_{\mathit{\text{ns}}}^{\mathit{+}}$ ]=[exp(inθ_s )], and J̃-=diag [J-, J-], J-=[${\mathbf{J}}_{\mathit{\text{ns}}}^{\mathit{-}}$ ]=[(-1) ⁿ exp(-inθ_s )] and exp(iθ_s )=(χ_s +iα_s )/k _⊥. Here, the diag function assembles its arguments into diagonal, or block diagonal, form.

 

Fig. 1. Geometry of the unit cell (defined by the fundamental translation vectors e₁ and e₂. Phase origins P ₁ and P ₂ and representations of the incoming (δ ^±) and outgoing (f ^±) plane wave trains is depicted.

Download Full Size | PDF

Proceeding as in other work [8, 9, 16], the vector forms of the outgoing, plane wave diffracted fields may be shown to be f ^±=δ ^±+2w$\tilde{\xi }$^T K ^̃±B̃ where w=k/(${\mathit{\text{dk}}}_{\perp }^e$ ²). In this expression, the terms δ ^± denote the non-diffracted part of the incident field, while the second term represents the diffracted field and involves changes of basis K ^̃± from cylindrical harmonic terms to upward and downward propagating plane waves. Here, K ^̃+=diag [K ⁺,K ⁺], K ⁺=[${\mathbf{K}}_{\mathit{\text{sn}}}^{+}$ ]=[exp(-inθ_s )] and K ^̃-=diag [K ^-,K ^-], K ^-=[$K_{\mathit{\text{sn}}}^{-}$ ]=[(-1) ⁿ exp(inθ_s )] Applying the field identities with the cylinder boundary condition B̃=R̃Ã, we form

where F=[f ^T -, ${\mathrm{f}}_{+}^{\mathit{T}}$] ^T , Δ=[δ^T -, ${\delta }_{+}^T$] ^T , K̃=[K ^T -,${\mathrm{K}}_{+}^T$] ^T , and J̃=[J _-, J ₊].

The calculation of this scattering matrix W ⁽⁰⁾ is aided by the up-down symmetry of the grating. Such symmetry enables the diffraction problem to be decomposed into E-symmetric and H-symmetric problems that lead to the factorization of the dense multipole scattering matrix K̃(Ĩ-R̃S̃ ^G )^-1 R̃J̃ in (14) into a block diagonal form, each block of which is half the dimension of the original.

The structure and significance of W ⁽⁰⁾, a matrix that arises in the subsequent eigenvalue problem, may be revealed by setting δ _- and δ ₊ in turn to zero. We thus deduce that relative to the natural phase origin (denoted by a superscripted (0)) at the centre of the central cylinder

where (${\mathbf{R}}_j^{\left(0\right)}$ ,${\mathbf{T}}_j^{\left(0\right)}$ ) denote reflection and transmission scattering matrices for incidence respectively from above (j=1) and below (j=2) the grating. Naturally, ${\mathbf{R}}_1^{\left(0\right)}$=${\mathbf{R}}_2^{\left(0\right)}$≡R ⁽⁰⁾, together with an analogous relationship for the transmission matrices, for a vertically symmetric structure such as a cylinder grating.

The formation of the band structure of the infinite cladding now requires that the grating be embedded in a hexagonal lattice that is characterized by basis vectors e₁=d(1, 0) and e₂=(s_x, d_y )=d(1,√3)/2 which define the unit parallelogram cell. In order to apply the Bloch condition that defines the eigenvalue problem for the lattice, it is necessary to relocate the phase origins of the scattering matrices (R _j ,T _j ) to points P_j =±(s_x, d_y )/2 located at the centres of the upper and lower edges of the parallelogram cell and separated by a lattice constant e₂. As in [11], this is achieved by a transformation W=𝓠𝓟W ⁽⁰⁾ 𝓟𝓠, where 𝓠=diag(Q̃, Q̃^-1, and 𝓟=diag (P̃,P̃) respectively phase the scattering matrices to perform the the lateral and vertical shifts of origin. Both Q̃ and P̃ involves operations on both TE and TM partitions and so each must be defined accordingly with Q̃=diag (Q,Q), with Q=exp(iα_ps_x /2), and P̃=diag (P,P), with P=exp(iχ_pd_y /2).

When the grating is embedded in the infinite lattice of Fig. 1, the Bloch condition requires that f-=µδ- and f⁺=µ ^-1 δ ⁺, where the eigenvalue, or Bloch factor, is given by µ=exp(-ik₀·e₂). For the particular symmetry of centred rectangular and hexagonal lattices, in which each grating layer is translated by half a period relative to the surrounding layers, there exists a significant simplification that enables the formulation of the eigenvalue problem in a stable and compact form-one that is of half the dimension and which is not afflicted by the stability problems of 𝓣-matrix methods. Here, Q=exp(iα₀d/4)Q₀ with Q₀=diag[exp(ipπ/2)], enabling the the eigenvalue problem to be recast in the form

with T′=Q${̃}_0^{-1}$ P̃T ⁽⁰⁾ P̃${\mathbf{Q}}_0^{~-1}$, R′=${\mathrm{Q}}_0^{-1}$ P̃R ⁽⁰⁾ P̃${\mathrm{Q}}_0^{~-1}$, µ′=µexp(iα₀d/2). Then, with

which has the same eigenvalues. In (17), c=(µ+µ ^-1)/2 and s=-i(µ-µ ^-1)/2. Some straightforward manipulation then leads to two algebraic eigenvalue equations which are equivalent to one another

Since the eigenvalues of (18) occur in the form 1/(2c), it follows that µ′ and 1/µ′ occur as a matched pair, corresponding to upward and downward travelling states that carry energy provided |µ′|=1.

We conclude by observing that plane wave methods compute the frequency as the eigenparameter given the Bloch vector. They therefore require an iterative process to handle material constants that are a function of frequency. In contrast, the process described above inverts the parameter order and is thus well suited to the computation of band diagrams for dispersive media. In the scattering matrix technique, the method takes a component of the Bloch vector and a given frequency as input data, embeds the material constants for that frequency during the scattering matrix calculation, and then generates the remaining component of the Bloch vector as the eigenparameter.

3 Air-guided modes

The location of complete band gaps requires traversing k-β parameter space, assessing everywhere whether or not propagating states (i.e., modes with |µ|=1) exist. We do this for each β by traversing the perimeter of the Brillouin zone for all required frequencies. Here this is achieved efficiently since the scans of ΓM and ΓXM in the Brillouin zone may be generated from the solution of (18) for the normal incidence $\left(k_{0x}=0\right)$ operation [11] of two diffraction problems, respectively corresponding to gratings aligned with e1 and e2, and with periods of d and √3d.

 

Fig. 2. Dispersion diagram for a hexagonal microstructured optical fibre overlaid with the light line (magenta curve) and the modal dispersion curve (red curve). The fibre data are: hole diameter d=4.026 µm, hole spacing Λ=5.7816 µm, central hole diameter d_c =13.1 µm.

Download Full Size | PDF

The background of Fig. 2 is the result of such scans over a 101×101mesh in k-β space, requiring approximately 50 minutes of computation on a 600 MHz Pentium III system. The method is rapidly convergent and a well stabilized solution was obtained from the solution of eigenvalue problems (18) with matrices of dimension 14, corresponding to only 7 plane wave orders (-3,-2, …, 0, …, 3). The color density is proportional to the number of propagating states, and thus complete band gaps are represented by the white fingers in which the mirror condition, referred to in Sec. 1 is satisfied. The magenta curve is the light line where β=k, and the red curve gives the dispersion relation for the fundamental mode. It is seen to extend into the region where cladding modes do exist; this is so because the mode calculation was performed for a structure with finite number of holes, for which the cladding reflectivity is a smooth function of wavelength. This should be contrasted to the band diagram calculations, in which the use of suitable lattice sums implies an infinite system.

The dispersion curve for the mode given in Fig. 2 shows it to be located on the high frequency side of the light line, so that n _eff ∑β/k<1. This means that the guiding mechanism cannot possibly be total internal reflection. An effective indexless than unity of course is not unphysical since this relates to a phase index rather than a group index.

The absorption loss of a mode may be calculated in two ways which we have found to be numerically equivalent. The first involves the use of a complex value of the matrix refractive index n_e and the direct determination of the complex propagation constant β, by making the determinant of 𝓜 in (9) zero. The second involves taking n_e to be real, calculating β and the associated modal fields, and then calculating the Ohmic loss of these fields in the matrix material from Poynting’s theorem, using the complex n_e . In doing so, note that the current induced by the electric field J=σE, where σ=ω∊₀ℑ(n ²).

The results of such calculations for the mode in Fig. 2 is shown in Fig. 3. It shows the variation with wavelength of the loss reduction factor αrf, defined to be the absorption in the bulk matrix material divided by the absorption loss of the mode (calculated from the imaginary part of the propagation constant β). The reduction in loss, varies between 3 and 5.2, reaching its maximum in the middle of the gap, where the confinement of radiation in the photonic crystal is strongest, and decaying near the edges of the gap, as the evanescent waves in the crystal become propagating.

In our actual calculations, we used 4 rings of holes, which constitutes a considerable computational effort. In this structure, the confinement losses [12] are of the same order of magnitude as those due to absorption. By comparing the results of calculations in which n_e is real to those in which it is complex, it is straightforward to determine the contributions of the absorption and confinement losses to the total attenuation. The results in Fig. 3 are obtained from the absorption losses only; this is justified since ℑ(n_e )≪ℜ(n_e ), and so the confinement losses are unaffected by the absorption of the matrix material. Moreover as long as ℑ(n_e )≪ℜ(n_e ) the curve in Fig. 3 does not on the actual absortion coefficient of of the material. It should be noted that the confinement losses can be reduced arbitrarily by increasing the number of holes. In contrast, the absorption losses do not decrease in this way, and absorption is therefore the effect ultimately determining the attenuation.

 

Fig. 3. Wavelength variation of the absorption reduction factor for the fibre of Fig. 2.

Download Full Size | PDF

The spatial variations of the magnitudes of the longitudinal components of the electric and magnetic fields are shown in Fig. 4, together with the z-component of the Poynting vector. Since the latter of these is quadratic in the fields, it is better confined than either E_z or H_z , though there is nonetheless significant fluxlea kage into the regions between the central hole and the first ring of smaller holes, and the region between the first and second ring of holes. Further work is needed to design structures with a larger fraction of the energy propagating in air, and therefore higher values of α_rf.

 

Fig. 4. Longitudinal components of the electric field, the magnetic field and the Poynting vector at the wavelength, λ=3.428µm, where the MOF of Fig. 2 has maximum reduction in material absorption. At this wavelength, the matrix index is n_e =1.39+i0.0003.

Download Full Size | PDF

4 Conclusions

The multipole method [7] outlined here has been used in studies of microstructured optical fibers MOFs and has proven to be a useful tool in modelling propagation and in estimating confinement losses [12]. As noted in the introduction, it handles finite collections of scatterers exactly and is not affected by numerical convergence problems and the use of periodic boundary conditions that afflict plane wave and other methods, and which, in turn, lead to fallacious outcomes such as the numerical birefringence of modes. By virtue of solving for the propagation constant β as a function of k, the method’s handling of material dispersion is demonstrably superior to techniques (such as plane wave eigenvalue methods) which solve for k as a function of β. Furthermore, the method is capable of yielding the geometric loss of finite solid core MOFs with estimates of ℑ(n _eff) down to 10^-14 being accessible.

The quality of the numerical results depends on the truncation of the multipole expansions- which, for computational purposes contain 2M+1 harmonics (-M,-M+1, …, 0, 1, …M) - and is tested by assessing the stability of the null vector(s) 𝓑 of the matrix Min (9 with increasing 𝓜, and also the convergence of the global and local field expansion at the boundaries of the inclusions. For systems which exhibit high symmetry (e.g. C_6v symmetry), the classification of the mode structure may be determined from group theoretic considerations, leading to substantial reductions in the computational requirements. In such cases [7], the computational demands are quite modest and, for example, a complete set of modes for a singled ring of six hexagonally spaced holes each of diameter 5 µm and pitch of 6.75 µm at a wavelength λ=1.45 µm in glass of refractive index n_e =1.45 may be computed for 1.4<ℜ(N _eff)<1.45 on a 733 MHz Pentium III system in under 3 minutes using less than 2 MByte of memory. The corresponding calculation for a three ring hexagonal structure of 55 holes requires under an hour of computation on a 500 MHz Compaq Alphastation and approximately 15 Mbyte of memory. Indeed, on such a workstation, we have undertaken studies of solid core MOFs with up to 180 holes. In all such calculations, the ratio of the diameter to the pitch (δ/Λ) is a significant factor in the choice of series truncation limits. While M=5 is suitable in many circumstances, it is necessary to use higher values with increasing proximity (i.e. increasing δ/Λ) of the cylinders.

The only significant drawback to the present technqique is its restriction to circular holes. However, it is possible to lift this restriction through the numerical computation of the coefficients that comprise the multipole reflection matrices of 𝓡 of (7). In this case, the partitions that comprise 𝓡 are now dense rather than diagonal.

The results reported here require the solution of multipole coefficients for each cylinder, with the accurate characterization of the field around the large central hole requiring additional terms M ₀ in the multipole series. The resulting matrix 𝓜 would be expected to have dimension 2[N_c (2M+1)+(2M ₀+1)] for an N_c +1 hole system. Here M ₀ is the truncationorder of the central hole and M is that for all others. For the four-ring structure considered here, with M=5 and M ₀=19, the matrix 𝓜 would have dimension 2012. The dimension of the matrix can be reduced, however, by using symmetry considerations, by a factor that depends on the type of mode that being considered. For nondegenerate modes the matrix dimension reduces by a factor 6, while for degenerate modes it reduces by a factor of roughly 3.5 [7]. Nonetheless, even with these reductions in the matrix size, the computational demands are near the limits of standard workstations, and the quality of the band gap confinement is still somewhat less than desirable.

Indeed, the optimization of the absorption reduction factor is likely to push us towards systems where the number of rings of air holes is large, and the holes come close together. At this point, the present implementation of the multipole treatment would become numerically cumbersome, with many coefficients required for each of a large number of holes. This argues for the extension of the multipole treatment in two theoretical directions and the consideration of a more advanced computational implementation. The first is an asympotic analysis where the outer rings of holes in a large system are treated perturbatively, rather than by direct solution. The second is an asymptotic analysis valid when air holes come close together, in which their nearest neighbour interaction is treated analytically, rather than by a multipole expansion that converges slowly in this case. At the same time, the structure of the algorithm is highly amenable to parallelization using the OpenMP protocol (available on symmetric multi-processor, shared memory systems) to implement highly efficient parallelisation in the loop structure and in the numerical linear algebra routines.

A general feature of the multipole method is that the matrix 𝓜 needs to be carefully conditioned to avoid under- and overflows in calculating the determinant. In addition, the minima of the determinant need to be carefully examined. Since the modes of structures with sixfold symmetry are either nondegenerate and doubly degenerate [7], the determinantal minima must be due to one or two eigenvalues of 𝓜 being small. These genuine minima can be distinguished from false ones, in which many eigenvalues are small simultaneously, by a singular value decomposition [7].

We conclude be recalling that the Bloch coeffient |µ| describes the spatial decay of the field over a lattice period. Hence, we suggest that structures with smaller values of |µ| require fewer rings of holes to achieve a given level of confinement losses. This may lead to a powerful tool for the design of photonic crystal fibres.

This work benefited from support from the Defence Science Technology Organization under the PA-10 project; helpful discussions with Dr. David Lancaster are acknowledged.