Photonic crystal fibres: mapping Maxwell's equations onto a Schrödinger equation eigenvalue problem

We consider photonic crystal ﬁbres (PCFs) made from arbitrary base materials and introduce a short-wavelength approximation which allows for a mapping of the Maxwell’s equations onto a dimensionless eigenvalue equations which has the form of the Schrödinger equation in quantum mechanics. The mapping allows for an entire analytical solution of the dispersion problem which is in qualitative agreement with plane-wave simulations of the Maxwell’s equations for large-mode area PCFs. We offer a new angle on the foundation of the endlessly single-mode property and show that PCFs are endlessly single mode for a normalized air-hole diameter smaller than ∼ 0.42 , independently of the base material. Finally, we show how the group-velocity dispersion relates simply to the geometry of the photonic crystal cladding. [DOI: 10.2971/jeos


I n t r o d u c t i o n
Photonic crystal fibres (PCF) are a special class of optical fibres where a remarkable dielectric topology facilitates unique optical properties.Guiding of light is accomplished by a highly regular array of air holes running along the full length of the fibre, see Figure 1.In their simplest form the fibres employ a single dielectric base material (with dielectric function ε b = n 2 b ) and for the majority of fabricated fibres silica has been the most common choice [1,2].This preference has of course been highly motivated by silica's transparency in the visible to infrared regimes, but has also been strongly driven by the highly matured glass technology used in the fabrication of standard telecommunication fibres.However, a growing interest in fibre optics for other wavelengths and light sources has recently led to a renewed interest in other base materials including chalcogenide glass [3], lead silicate glass [4], telluride glass [5], bismuth glass [6], silver halide [7], teflon [8], and plastics/polymers [9].The huge theoretical and numerical effort in the PCF community has been important in understanding the dispersion and modal properties, but naturally there has been a clear emphasize on silica-based PCFs.PCFs fabricated from different base materials typically share the same overall geometry with the air holes arranged in a triangular lattice and the core defect being formed by the removal of a single air-hole.However, it is not yet fully clarified how the established understanding may be transferred to the development of PCFs made from other base materials.To put it simple; PCFs made from different base materials share the same topology, but to which extend do they also have optical characteristics in common?
When addressing this question it is important to realize that the otherwise very useful concept of scale-invariance of Maxwell's equations [10] is of little direct use in this particular context.Since PCFs made from different base materials do not relate to each other by a linear scaling of the dielectric function, ε(r) s 2 ε(r), scale invariance cannot be applied directly to generalize the results for e.g.silica to other base ma-terials.Electro-magnetic perturbation theory is of course one possible route allowing us to take small changes of the basematerial refractive index into account, but in this paper we will discuss and elaborate on an alternative approach which was put forward recently [11].
As a starting point we take the classical description based on Maxwell's equations, which for linear dielectric materials simplify to the following vectorial wave-equation for the electrical field [10] Here, ε is the dielectric function which in our case takes the material values of either air or the base material.For a fibre geometry with r along the fibre axis we have ε(r) = ε(r ⊥ ) and as usual we search for solutions of the plane-wave form e i(βr −ωt) with the goal of calculating the dispersion ω(β) relation.

S H O R T -W A V E L E N G T H A P P R O X -I M A T I O N S
Today, several approaches are available for solving Eq. ( 1) numerically, including plane-wave, multi-pole, and finiteelement methods [12]- [14].Such methods will typically be preferred in cases where highly quantitative correct results are called for.However, numerical studies do not always shine much light on the physics behind so here we will take advantage of several approximations which allow for analytical results in the short-wavelength limit λ Λ.The key observation allowing for the approximations is that typically the base material has a dielectric function exceeding that of air significantly, ε b 1.In that case it is well-known that the short-wavelength regime is characterized by having the majority of the electrical field residing in the high-index base material while the fraction of electrical field in the air holes is vanishing.
The, at first sight very crude, approximation is simply to neglect the small fraction of electrical power residing in the air-hole regions.Mathematically the approximation is implemented by imposing the boundary condition that E is zero at the interfaces ∂Ω to air.Since the displacement field D = εE is divergence free we have 0 = ε∇ • E + E • ∇ε ≈ ε∇ • E and the wave equation (1) now reduces to and since solutions are of the plane-wave form e i(βr −ωt) along the fiber axis we get where ∇ 2 ⊥ is the two-dimensional Laplacian in the transverse plane.The essentially scalar nature of this equation makes the problem look somewhat similar to the more traditional scalar approximation that has been applied successfully to PCFs in the shortwavelength regime [15].However, while that approach took the field and the dielectric function in the air holes into account we shall here solve Eq. ( 3) with the boundary condition that E is zero at the interfaces to air.We note that this of course is fully equivalent to solving Eq. ( 4) with the boundary condition that Ψ is zero at the interfaces to air.Obviously, the scalar problem posed by Eq. ( 3) is separable and formally we have that where Ω = cβ/n b is the frequency associated with the longitudinal plane-wave propagation with a linear dispersion relation (i.e.Ω ∝ β) and Ω ⊥ is the frequency associated with the transverse confinement/localization.At this point we already note how an arbitrary base-material refractive index n b enters in the frequency associated with the longitudinal plane-wave propagation.This is solely possible because we have a Dirichlet boundary condition in Eq. ( 3) which does not depend on the refractive index of the base material.

T H E S C H R Ö D I N G E R E Q U A T I O N E I G E N V A L U E P R O B L E M
In the following we rewrite Eq. ( 5) as where now γ is a dimensionless number characterizing the confinement/localization.It is of purely geometrical origin and thus only depends on the normalized air-hole diameter d/Λ.From Eq. ( 3) it follows that γ is an eigenvalue governed by a scalar two-dimensional Schr ödinger equation The same equation was recently studied in work on anti-dot lattices in a two-dimensional electron gas [16] and there will be many similarities between such quantum mechanical electron systems and the present electromagnetic problem.In the language of quantum mechanics the wavefunction ψ is subject to hard-wall boundary conditions corresponding to an infinite potential barrier in the air-hole regions.The task of calculating the optical dispersion properties has now reduced to solving the scalar two-dimensional Schr ödinger-like eigenvalue problem posed by Eq. ( 7).The strength is obviously that Eq. ( 7) is material independent and thus we may in principle solve it once and for all to get the geometry dependence of γ and thereby the optical dispersion of the PCF.In the following we solve Eq. ( 7) for different classes of modes and compare approximate results to corresponding numerically exact finiteelement solutions.

. 1 T h e f u n d a m e n t a l s p a c e -f i l l i n g m o d e
The fundamental space-filling mode is a key concept in understanding the light guiding properties of PCFs.It is the fundamental de-localized mode of Eq. ( 7) in an infinite periodic structure and the corresponding mode index corresponds to an effective material index of the artificial periodic photonic crystal lattice.In the language of band diagrams and Brillouin zones it is defined in the Γ-point where Bloch's theorem is particular simple; as illustrated in Figure 2 the wave function is subject to a simple Neumann boundary condition on the edge ∂ Ω of the Wigner-Seitz cell, i.e.
where n is a unit vector normal to the boundary ∂ Ω.We thus have to solve an eigenvalue problem on a doublyconnected domain.By a conformal mapping, leaving the Laplacian and the boundary conditions invariant, one may in principle transform the hexagonal shape with the circular hole of diameter d into a simple annular region of inner diameter d and outer diameter 2R [17] leaving us with an eigenvalue problem in the radial coordinate.Finding the exact conformal mapping may be a complicated task so here we simply use the approximate mapping with which serves to conserve the area of the Wigner-Seitz unit cell.Neglecting any small changes to the eigenvalue equation caused during the mapping we get the following eigenvalue equation in the radial coordinate, The solution is given by a linear combination of the Bessel functions J 0 and Y 0 and the eigenvalue is determined by the roots of the following expression (11) In general the equation cannot be solved exactly except for at the point {0.66, 53.26} , (12)   as is easily verified by insertion and using that α n,m is the mth zero of the nth Bessel function, i.e.J n (α n,m ) = 0.In particular we have that α 0,1 2.405 and α 1,1 3.8317.However, expanding the left and right-hand sides of Eq. ( 11) around this point to first order in both γ cl and d/Λ we get an equation from which we may isolate γ cl yielding an approximate analytical solution of the form The coefficients are given by expressions involving Bessel functions Y 0 and Y 1 evaluated at α 0,1 and α

. 2 C o r e -d e f e c t m o d e s
While the fundamental space-filling mode, discussed above, is a de-localized mode it is also possible to have strongly localized states, especially in the presence of defects such as one formed by removing an air hole from the otherwise periodic lattice of air holes, see Figures 1 and 3.In this way light can be guided along the defect thus forming the core of the fibre.The requirement for corralling the light is of course that the coredefect mode has an eigenvalue γ 2 c,1 not exceeding the corresponding eigenvalue γ 2 cl of the surrounding photonic crystal cladding.To put it slightly different the defect in the air-hole lattice effectively corresponds to a potential well of depth γ 2 cl and a radius roughly given by ρ = Λ − d/2, see Figure 3.In order to calculate the fundamental core-defect mode we first note that in the limit of d/Λ approaching 1, the problem can be approximated with that of a two-dimensional spherical infinite potential well with radius ρ = Λ − d/2, see Figure 4.For this problem the lowest eigenvalue is

∂Ω
Although this expression yields the correct scaling with d/Λ, the approximation obviously breaks down for small values of d/Λ.In that limit we follow the ideas of Glazman et al. [18], who studied quantum conductance through narrow constrictions.The effective one-dimensional energy barrier for transmission through the constriction of width Λ − d between two neighbouring holes has a maximum value of π 2 in the limit where d goes to zero.We thus approximate the d/Λ → 0 problem with that of a two-dimensional spherical potential well of height π 2 and radius Λ, where Θ(x) is the Heaviside step function which is zero for x < 0 and unity for x > 0. The lowest eigenvalue γ2 c,1 for this problem is the first root of the following equation γc,1 In the case where the potential π 2 is replaced by an infinite potential the first root is given by α 2 0,1 .Lowering the confinement potential obviously shifts down the eigenvalue and in the present case the eigenvalue is roughly γ2 c,1 ∼ π.In fact, expanding the equation to first order in γc,1 around e.g.√ π we get a straight forward, but long, expression with the numerical value γ2 c,1 3.221 which is in excellent agreement with a numerical solution of the equation.Correcting for the low-d/Λ behaviour we thus finally get For the first high-order mode we may apply a very similar analysis.This mode has a finite angular momentum of ±1 with a radial J 1 solution yielding Here, γ2 c,2 7.673 is the second eigenvalue of the twodimensional spherical potential well of height π 2 and radius Λ.Again, it can be found from an equation very similar to Eq. ( 16).

. T w o -d i m e n s i o n a l f i n i t e -e l e m e n t s o l u t i o n s
Above we have solved the geometrical eigenvalue problem analytically by means of various approximations.In this section we compare the quality of these results by direct comparison to a numerically exact solution of the twodimensional eigenvalue problem.The developments in computational physics and engineering have turned numerical solutions of partial differential equations in the direction of a standard task.Here, we employ a finite-element approach [19] to numerically solve Eq. ( 7) and calculate γ 2 versus d/Λ.We employ an adaptive mesh algorithm to provide efficient convergence.For cladding mode we implement Eq. ( 8) directly while for the defect modes we solve Eq. ( 7) on a finite domain of approximate size 10Λ × 10Λ with the defect located in the center of the domain.For d/Λ down to around 0.1 we have found this to give adequate convergence for the strongly localized defect states, thus avoiding strong proximity from the domain-edge boundary condition which for simplicity has been of the Dirichlet type.Figure 5 summarizes our numerical results for the first cladding mode γ 2 cl as well for the two first core-defect modes γ 2 c,1 and γ 2 c,2 .The dashed lines indicate the corresponding approximate expressions in Eqs. ( 13), (17), and (18) obtained by analytical means.As seen the qualitative agreement is excellent, but in order to facilitate also quantitative applications of the results we also include the results of numerical least-square error fits below.The thin solid lines show the expressions which match the numerical data within a relative error of less than 2% around the most important cut-off region d/Λ ∼ 0.42.

D E R I V E D F I B E R O P T I C A L P A -R A M E T E R S
With Eq. ( 6) at hand we have now provided a unified theory of the dispersion relation in the short-wavelength regime for PCFs with arbitrary base materials and Eq. ( 6) illustrates how geometrical confinement modifies the linear free-space dispersion relation.
In fibre optics it is common to express the dispersion properties in terms of the effective index n eff = cβ/ω versus the free-space wavelength λ = c2π/ω.From Eq. ( 6) it follows straightforwardly that which obviously is in qualitative agreement with the accepted view that n eff increases monotonously with decreasing wavelength and approaches n b in the asymptotic short-wavelength limit as reported for e.g.silica-based PCFs [2].Similarly, the group-velocity v g = ∂ω/∂β becomes while the phase velocity v p = ω/β becomes The group-velocity dispersion is most often quantified by the wave-guide dispersion parameter D wg = ∂(1/v g )/∂λ which becomes where any possible material dispersion n b (λ) has been neglected.We note that since D wg > 0 the geometry of the airhole lattice always causes a positive wave-guide dispersion parameter.Large-mode area PCFs belong to the regime with λ Λ were we predict the following general magnitude of the wave-guide dispersion of a large-mode area PCF lim Finally, the recently suggested parameter [20] can be shown to be a purely geometrically defined parameter in the large-mode area limit.From Eq. ( 6) it follows straightforwardly that lim This implies that the endlessly single-mode property [2], [20]- [22] is a wave phenomena independent of the base material refractive index of purely geometrical origin.Higher-order modes are only supported for d/Λ 0.42, as seen in Figure 5, for which V PCF π [20].

C O M P A R I S O N T O F U L L Y -V E C -T O R I A L P L A N E -W A V E S I M U L A -T I O N S
In the previous sections substantial analytical progress was made, but we still remain to address the question to which extend the basic assumption behind Eq. ( 3) holds except that it becomes exact for the fundamental modes as λ/Λ approaches zero.In Figure 6 we compare the analytical results for the effective index in Eq. ( 20) to results obtained by fully-vectorial plane-wave simulations [12] of Eq. ( 1).For the fundamental space-filling mode we have employed a basis of 2 6 × 2 6 plane waves while for the core-defect modes we have employed a super-cell configuration of size 10Λ × 10Λ in the case of d/Λ = 0.4.Panel (a) shows results for the fundamental space-filling mode for various values of the base refractive index.As clearly seen the theory agrees excellently in the short-wavelength limit while pronounced deviations occur as the wavelength increases and the field penetrates deeper into the air hole regions and also vectorial effects become important as reported in [15].Similar conclusion applies to the corresponding results in panel (c) for the fundamental core-defect mode.While d/Λ = 0.4 is a case of particular technological interest because of its endlessly single-mode property we would like to emphasize that equivalent agreement is found for other values of d/Λ (results not shown).However, since the success of the approximation in Eq. (3) really is that λ/n b 2ρ = 2Λ − d and λ/n b Λ − d for the core and cladding results, respectively, there will of course be a small d-dependence.However, this also indicates that the agreement will increase with an increasing refractive index n b of the base material as has also been observed for even higher values of n b than those studied in Figure 6 (results not shown).

D I S C U S S I O N A N D C O N C L U S I O N
It is today more than 10 years ago that the PCF was invented by Russell and co-workers [1] and there is a still growing community of researchers directing their efforts toward fabrica-tion and experimental studies of silica-based PCFs as well as quantitative modelling of the optical properties.At the same time the community is obviously facing the challenges and opportunities of new exciting fiber materials.PCFs made from different base materials share the same topology so it seems quite natural to assume that they also have at least some basic optical properties in common.Let's return to the question posed in the introduction: to which extend do PCFs made from different materials have optical characteristics in common?The present work do to some extend address this question.Most importantly we illustrate how the waveguide dispersion originates from the geometrical transverse confinement/localization of the mode and how the endlessly singlemode property arises as a sole consequence of the geometry which acts as a modal sieve.In particular we have shown how PCFs are endlessly single-mode for d/Λ 0.42 irrespectively of the base material and we have also derived an analytical expression for the wave-guide dispersion parameter applicable to large-mode area fibres.
For small-core PCFs our theory provides qualitative correct results though more quantitative insight still calls for fully vectorial simulations, see e.g.[23].However, for large-mode area PCFs our expressions not only give qualitative insight, but also quantitative correct expressions which may be used in straightforward design of large-mode area PCFs with special properties with respect to the group-velocity dispersion and the susceptibility of the fundamental mode with respect to longitudinal non-uniformities [24].

A C K N O L E D G M E N T S
C. Flindt, J. Pedersen, and A.-P. Jauho are acknowledged for stimulating discussions on the results in Section 3 and for sharing numerical results.

A p p e n d i x A T A Y L O R E X P A N S I O N O F E Q . ( 1 1 )
Taylor expanding the left and right-hand sides of Eq. ( 11) around the point

FIG. 1
FIG. 1 Optical micrograph of a silica-based PCF with a pitch Λ of order 10 microns and d/Λ ∼ 0.5.Courtesy of Crystal Fibre A/S, www.crystal-fibre.com.

FIG. 2
FIG. 2 The hexagonal Wigner-Seitz unit cell of the periodic triangular lattice with lattice constant/pitch Λ.The dashed line indicates the annular region which has the same area as the doubly-connected domain of hexagonal outer shape.Dirichlet and Neumann boundary conditions apply to the inner and outer boundaries ∂Ω and ∂ Ω, respectively.

FIG. 3
FIG. 3 The topology of a PCF formed by omitting an air hole of diameter d in a triangular lattice of air holes with pitch Λ. Dirichlet boundary conditions are applied to the boundary ∂Ω indicated by red solid lines.The red dashed line indicates twice an effective core radius ρ = Λ − d/2.

∂Ω 2ρ FIG. 4
FIG. 4 Zoom in on the core-defect in the extreme limit of d/Λ = 1.A Dirichlet boundary condition applies to the boundary ∂Ω and the dashed circle indicates the approximate core radius ρ = Λ − d/2.

2 FIG. 5
FIG. 5 Geometrical eigenvalues γ 2 versus normalized air-hole diameter d/Λ.The data points are results of finite-element simulations while the corresponding dashed lines are the approximate results in Eqs.(13),(17), and(18).The thin solid lines shows the numerical fits listed in Eqs.(19a), (19b), and (19c).Note how only a single core-defect mode is supported in the shaded region, d/Λ 0.42.

FIG. 6
FIG. 6 Effective index versus normalized wavelength λ/Λ for the fundamental spacefilling mode (a) and the fundamental core-defect mode (c) for a PCF with d/Λ = 0.4 with a base material refractive index n b varying from 1.45 to 1.6 in steps of 0.05 from below.The data points are the results of a vectorial plane-wave simulation of Eq. (1) while the dashed lines show the corresponding results based on Eq. (20) with Eqs.(19a) and (19c) for geometrical eigenvalues γ 2 cl and γ 2 c,1 .Panels (b) and (d) show the corresponding eigenfunctions ψ.