## Abstract

Modeling of microstructure fibers often involves severe computational bottlenecks, in particular when a design space with many degrees of freedom must be analyzed. Perturbative versions of numerical mode-solvers can substantially reduce the computational burden of problems involving automated optimization or irregularity analysis, where perturbations arise naturally. A basic theory is presented for perturbative multipole and boundary element methods, and the speed and accuracy of the methods are demonstrated. The specific optimization results in an elliptical-hole birefringent fiber design, with substantially higher birefringence than the intuitive unoptimized design.

©2004 Optical Society of America

## 1. Introduction

Microstructure fiber (MOF), first explored in the 1970’s [1], has recently stimulated much attention with the demonstration of supercontinuum generation [2] and proposals for a variety of fibers with useful characteristics not achievable using standard fibers. The potential for air-core data transmission is particularly interesting, since it would represent a dramatically different regime of optical communications, with drastically reduced nonlinearities [3].

Modeling of microstructure fibers has made rapid progress recently [4], including demonstrations of several accurate mode-solving algorithms [5, 6, 7, 8, 9] and investigation of a large variety of application-specific designs. While current simulation tools are general enough to provide a basic analysis of most fibers of interest, important obstacles continue to arise. One reason is that relevant fiber characteristics may require more than the basic mode-finding provided by these algorithms. For example, coupling between modes can play a critical role in analysis and design [10, 11, 12] and requires additional analysis after fiber modes have been calculated.

Another type of obstacle is the severe numerical bottleneck that arises for fibers with complex geometries, or when analysis of a large design space is needed. For this reason, most modeling work has been restricted to MOF with only two degrees of freedom. While these restricted studies are informative, much of the promise of MOF is in the design freedom of complex geometries. This freedom can only be utilized once designs spanning many degrees-of-freedom can be analyzed.

A perturbative numerical formulation would allow large speedups on important, numerically challenging problems. It is well suited to finding optimum perturbations that improve a design specification, and for analysis of structural irregularity. The latter problem can be computationally burdensome even when target design is simple, since the sensitivity analysis requires analysis over the space of all possible perturbations.

This paper presents a perturbative form of the numerical mode equations applicable to the multipole (MP) and boundary element (BE) methods. We demonstrate the formalism with examples using both simulation methods, and for both sensitivity analysis and optimization over several degrees of freedom. The perturbative expressions are carefully tested against exact solutions of the standard (non-perturbative) MP and BE equations.

## 2. Perturbative approach

A perturbative picture is very natural in describing structural irregularity of MOF. Realistic fibers have a hole geometry that deviates slightly from some ideal structure (as in Fig. 1). For useful fibers made by a reliable fabrication process, these deviations should be assumed small, but their effect on fiber performance can be large [13, 14, 15]. Perturbations can change essentially all of the mode properties: birefringence, mode area and shape, dispersion, loss, coupling lengths (of a multi-core or multi-mode fiber), etc. Analysis can involve an extreme burden for several reasons: The perturbations typically have many degrees of freedom (previous analysis of irregularity have analyzed only one or two, hopefully representative, degrees of freedom). The perturbed fiber may be non-symmetric, increasing the size of the corresponding numerical problem substantially [6]. Finally, for some reliability issues, rare events may be important, so that many random perturbations should be analyzed.

#### 2.1. Reformulating the multipole and boundary-element methods

A number of simulation methods seem sufficiently accurate and efficient for most applications. However, the multipole and boundary-element methods have the advantage of being based on very compact field representations. Because they assume a finite number of constant-dielectric regions, the entire field can be reduced to expansion coefficients on the relatively small number of “holes” in the microstructure. Because of its basic representation, the MP method is restricted to fibers with circular holes, but the BE method has no such restriction.

These methods do not have mode equations of the standard eigenvalue form * Av* =

*λ*, the most familiar starting point for perturbation theory in physics. The mode parameters

**v***λ*and

*n*

_{eff}=

*k*

_{z}/

*k*=

*λk*

_{z}/2

*π*are instead wrapped into a mode equation in a nonlinear way:

Here * v* is a vector of coefficients representing the field, and

*is a matrix imposing Maxwell’s equations for the given fiber geometry. (Other methods, such as the plane-wave or finite-element methods,*

**M***do*have the standard eigenvalue form, although it has been pointed out that care must nevertheless be taken when perturbing these equations [16]).

We can still obtain systematic corrections, following the spirit of more familiar eigenvalue solutions. The geometrical change is written as proportional to a small parameter *δ*. That is, if * p* is a vector of parameters describing the geometry (hole positions, radii, etc), then

*=*

**p**

**p**_{0}+

*δ*

**p**_{1}. The matrix

*is directly dependent on the geometrical parameters, and the wavelength:*

**M**So that an expansion of * M* in orders of

*δ*can be explicitly calculated using numerical approximations to the derivatives (with small ε):

The effective index *n* and field coefficients * v* are unknowns of Eq. (1), not explicit functions of the geometry. They can also be expanded in orders of

*δ*,

The unknowns (**v**_{1}, *n*
_{1}, etc) are solved by plugging into (1),

and sorting into orders. The vectors **v**_{1}, **v**_{2}, are assumed orthogonal to **v**_{0}. In the simplest case, a first-order effective index shift can be obtained as a very simple expression,

Here **M**_{1} and * M*′

_{0}are derivatives of

*with respect to*

**M***δ*and

*n*, respectively, and

**u**_{0}and

**v**_{0}are the left and right null vectors of the unperturbed matrix

**M**_{0}(

*n*

_{0}). In the next section we discuss details of the derivation and more complicated cases.

The calculation assumes that a mode solution has already been calculated for the unperturbed problem, and requires some additional computation to generate the matrix derivatives. The advantages of this approach depend on the problem and implementation. Clearly, perturbed geometries can be analyzed using the standard (non-perturbative) version of the mode-solver. The perturbative approach solves the same problem (approximately), but in many cases it can eliminate much of the work. In Fig. 2, we see that the standard mode solver requires that a large matrix be formed and decomposed many times, until a search for singularities is completed. The perturbative approach instead requires very few large-matrix operations, once the unperturbed solution is known. In particular, for problems like sensitivity analysis, one would like to calculate the effects of many geometrical perturbations for each unperturbed geometry analyzed, and the perturbative approach offers a substantial speedup. We anticipate that for other problems speed may not be limiting, but the perturbative framework may still prove useful in understanding why perturbations to the index, loss, and fields arise.

## 3. Perturbative estimates and irregularity analysis

Detailed expressions are now presented in the context of irregularity analysis. The results included below are representative of a number of calculations comparing results of perturbative expressions with the standard, non-perturbative (multipole or boundary element) mode solvers. The following example demonstrates the simplest case, where the unperturbed solution is non-degenerate. That is, when there is only one vector **v**_{0} satisfying the unperturbed equation for *n*
_{eff} =*n*
_{0},

This could represent a mode that is physically non-degenerate. In the current example, it is instead the fundamental mode of a solid-core triangular-lattice fiber. While the mode is degenerate, it can be separated from physically degenerate modes by using symmetry [17, 18] as long as the geometric perturbations obey the relevant symmetry–in this case hole displacements preserve the sixfold rotational symmetry. The solution on the restricted space is mathematically non-degenerate.

Estimates are derived by sorting terms from Eq. (7). The zero-order equation is simply Eq. (9), the unperturbed solution. First-order terms in *δ* give an expression

in the unknown index shift *δn*
_{1} and field coefficients *δ v*

_{1}. Since

**M**_{0}(

*n*

_{0}) has exactly one right null vector, it also has a single left null vector,

**u**^{†}

_{0}

**M**_{0}(

*n*

_{0}) = 0, which we use to project out the effective index perturbation:

This readily gives Eq. (8).

Just as *δn*
_{1} is obtained by projecting out the null space of **M**_{0}, the perturbed field coefficients *δ v*

_{1}are obtained by projecting out the complimentary space. On this restricted space,

**M**_{0}is invertible, and so the solution involves a matrix

**M**_{0}

^{p}that acts as an inverse of

**M**_{0}on this subspace:

It is straightforward to obtain *δ v*

_{1},

*δ*

**u**_{1}, and

**M**_{0}

^{p}from a singular-value decomposition of

**M**_{0}(

*n*

_{0}).

First-order results will most likely prove the most useful because of their simplicity and speed. However, second-order effective index estimates are manageable, and give more accurate results. The second order equation is

Where $\mathcal{M}$
_{2} includes all of the second-order matrix terms:

Multiplying by *δ u*

_{1}

^{†}eliminates

*δ*

^{2}

*δ*

**v**_{2}. The remaining scalar equation can easily be solved in the remaining unknown,

*δ*

^{2}

*n*

_{2}.

Results of an example calculation are shown in Figure 3. Mode effective index and loss were calculated for fifty pseudorandom displacements of the holes. Perturbative estimates are compared with results of the standard mode-soler. Errors in the effective-index shift are quite small in all cases, and errors roughly follow the expected accuracy trends (dashed lines indicating Error ~ *풪*{Δ*n*
^{2}} and Error ~ 풪{Δ*n*
^{3}} for first-order and second-order estimates, respectively), down to an apparent “noise floor” in the 10^{-6} range. Loss estimates are also fairly accurate, and impressive when we consider that the loss is represented by an imaginary part of the effective index, with unperturbed values around *𝓘*{*n*
_{eff}} ≈ 6×10^{-7}. Naturally, the imaginary part is much more difficult to estimate than the real part, and can easily be buried in noise arising from a careless implementation (poor choice of *ε* in the numerical derivatives, etc.). The current results are improved from a previous implementation [19]: they allow accurate first-order estimates of loss for geometrical deviations up to *σ*_{x}
= *σ*_{y}
= .01Λ. In the Fig. 3, we see that second-order estimates can cover even larger deviations, where confinement loss is a factor of two above or below its unperturbed value.

Excess confinement loss due to irregularity is calculated here for the first time, to my knowledge. While it is intuitively obvious that some perturbations of a fiber can cause changes in confinement loss (for example, particular perturbations amount to a change in hole pitch), detailed studies are needed to indicate when this excess loss is of practical relevance. Since the standard multipole method has provided a number of accepted loss calculations, the perturbative version (along with test calculations) is an ideal tool for further studies. Of course, real fibers have other sources of loss, and so excess confinement loss should be interpreted as a lower bound.

#### 3.1. Degenerate modes and birefringence

Irregularity-induced birefringence is a key application motivating a perturbative approach, and is a good example of a problem where degenerate modes are unavoidable. The formulas for first and second-order perturbations have fairly straightforward generalizations to the degenerate case, and are outlined below. The expressions have previously been tested against standard multipole birefringence values [19], again showing excellent agreement. The results also qualitatively confirms other published results [13, 14, 20, 15], that for fibers with large *d*/Λ even percent-level geometrical imperfections can cause prohibitively large birefringence. A systematic index difference of 10^{-4} means that picoseconds of skew due to polarization mode dispersion can accumulate in just a few meters, if mode-mixing is small. Similarly, the birefringent beat-length is very short (~1 cm), and the practicality of spinning such fibers effectively is questionable.

For degenerate solutions, there is still some unperturbed solution *δ v*

_{1}6 [Eq. (5)] near the exact solution

**v**, but now even the unperturbed solution is not completely determined

*a priori*. It could be any linear combination of the

*K*degenerate modes, and the correct combination will depend on the perturbation. This has a simple interpretation in terms of birefringence: for a symmetric unperturbed structure, we cannot even guess the axis of birefringence, let alone its magnitude, until we know something about the perturbation. Similarly, even if the unperturbed problem has degenerate modes that can be separated by symmetry, a randomly perturbed structure typically will not have the same symmetry, and a separation into symmetry classes will not be possible.

We can in general represent the unperturbed solution using a small *K*-dimensional vector of superposition coefficients * x*:

Here, * B* is a tall, thin

*N*×

*K*matrix satisfying

**M**_{0}(

*n*

_{0})

*= [0]*

**B**_{N×K}. (Typically

*N*is in the hundreds or thousands and

*K*is two; accidental degeneracies can larger, but symmetry-induced degeneracies are always of order two [21]). Similarly, the projection we need is a short, wide

*K*×

*N*matrix

**A**^{†}, satisfying

**A**^{†}

**M**_{0}(

*n*

_{0}) = [0]

_{K×N}. The projection of the first-order Eq. (10) gives an expression for the effective index shift,

with a generalized-eigenvalue form (* Px* =

*n*

_{1}

*). Solving this for*

**Qx***n*

_{1}is a negligible computation since

*and*

**P***are small (usually 2×2).*

**Q**## 4. Application: Optimization

In the above sections, we saw that a perturbative calculation can speed up the analysis of geometrical irregularity in microstructure fibers. A number of other applications are worth looking into: group velocity and dispersion calculations require derivatives with wavelength; corrections to analytical approximate solutions might provide more useful agreement with exact results; a similar approach could perhaps be useful in untangling nearly-degenerate modes, a significant issue for important microstructure fibers such as aircore and multi-core couplers. The diversity of successful perturbative approaches in the physics literature offers encouragement to look for motivation beyond simple accounting of FLOPS.

In this section, we look at automated design optimization, another obvious application where fast calculation of perturbations addresses an important bottleneck. Gradients are obtained by intentionally perturbing the geometry independently along each degree of freedom. These gradients can then be incorporated into a smart high-level optimizing algorithm.

The following example demonstrates this optimization approach, using the design of a high-birefringence fiber with elliptical air holes, similar to those discussed in [22, 23, 8]. It has two nice features:

- It highlights the flexibility of the boundary element method to handle arbitrary non-circular hole shapes (in contrast to the multipole method), and
- It uses a design problem simple enough that the optimum solution might seem intuitively obvious.

High birefringence is a feature of polarization-maintaining fiber, a product with several important uses that could potentially be improved by the unique properties of microstructure fibers. We fix the size and shape of the holes (hole area = .3Λ^{2}; *η* = 2 = major axis/minor axis) but allow some of the holes to rotate so that they are not necessarily aligned along a common axis. For problems such as this (and aircore fibers [24]), the ability to directly simulate non-circular holes is an obvious advantage. We can intuitively guess that birefringence will tend to be maximized by aligning the holes to a single axis.

Figure 4 shows results of an optimization where two holes adjacent to the core were fixed, and the remaining four inner holes were allowed to independently rotate. The initial and final geometries are shown along with the birefringence as a function of iteration step in the optimization. As expected, optimization aligns the four holes that are initially randomly oriented. At each optimization step, a gradient was calculated using first-order perturbative estimates of the index shifts for each hole orientation. The vector of geometric parameters in this case is just a vector of orientation angles for the holes * p* = [

*ϕ*

_{1}

*ϕ*

_{2}…]

^{T}, and the

*j*th element of this gradient corresponds to hole

*j*:

Here the birefringent index difference *n*
_{Bi} has been defined as the difference in index between the “slow” and “fast” nearly-degenerate modes (*n*
_{slow} > *n*
_{fast}). Once a gradient has been calculated for step *m*, the geometry for the next step **p**^{m+1} can be determined:

For simplicity, a crude gradient-crawling update with fixed perturbation magnitude *δ* was used,

The perturbative approach simply provides fast gradients, and can be incorporated equally well into any optimization algorithm that uses gradients (see, for example [25, section 10.5]). Selecting suitable smart optimizers is an area for further work. The boundary element mode-solver used follows the description in [26], and is similar to that described in [8].

The modal dispersion curve was essentially the same as the similar fiber in [23] for wavelengths below λ ~ 2Λ. The fields are essentially linearly polarized, dominated by either x-or y-directed components of the magnetic field, as expected. Gradient estimates were tested against results of the non-perturbative boundary-element method. In Fig. 5, changes in birefringence for each step in the optimization are shown, with first-order estimates plotted against non-perturbative results. The agreement assures us that the implementation is correct and that the change at each step is not too large (better agreement is obtained for smaller changes; in the example, *δ* = .06 limits the rotations per step to Σ_{j}Δ${\varphi}_{j}^{2}$ = (3.4degrees)^{2}).

Interestingly, if all of the holes are allowed to freely rotate, the optimum structure is not the intuitive “perfectly aligned” formation. The results of an unconstrained optimization are shown in Figs. 6 and 7, and the final birefringence is seen to be appreciably larger than the initial perfectly-aligned geometries. The *x*- and *y*-aligned initial geometries result in equivalent final optimized structures, although the *y*-aligned is rotated 120 degrees and has not completely converged after 50 iterations. As one might expect, the optimization arrives at a structure with *C*
_{2v} symmetry (reflection on two orthogonal axes).

These results suggest some preliminary directions for fabricating real polarization-sensitive fibers. For one thing they indicate that more flexible fabrication methods, such as sol-gel casting [27, 28], drilling [29], and extrusion [30] may have important advantages in achieving design specifications, even when the fine-tuning of the geometry is fairly subtle. This issue should be pursued further along with studies of how hole shape can be controlled in the preform and during draw.

## 5. Conclusions

We have presented a perturbative numerical approach and demonstrated its application for the multipole and boundary-element mode-solvers. We have shown that it is applicable to two important types of problems, sensitivity analysis and design optimization, where solutions must span a large design space, involving an acute computational bottleneck. A large speedup was found compared to a non-perturbative implementation.

Specifically, we have looked at the problem of geometrical irregularity using a perturbative multipole method. Previous results were extended by looking at perturbations with many degrees of freedom, and looking at irregularity-induced excess confinement loss. Optimization with many degrees of freedom was also demonstrated: elliptical-hole highly birefringent fibers were designed using a perturbative boundary-element method. In the process, it was discovered that an optimized elliptical-hole fiber can achieve appreciably higher birefringence than the intuitive design with all holes oriented on the same axis.

Perturbative expressions have been tested and confirmed by comparison with the non-perturbative solutions. While this approach offers improved efficiency that may be important for some problems, the framework can hopefully be extended in the future to provide qualitatively new modeling abilities.

## References and links

**1. **P. Kaiser and H. W. Astle. “Low-loss single-matrial fibres made from pure fused silica,” Bell Syst. Tech. J. **53**, 1021–39 (1974).

**2. **S. A. Diddams and D. J. Jones, et al. “Direct link between microwave and optical frequencies with a 300THz femtosecond laser comb,” Phys. Rev. Lett. , **84**, 5102–5 (2000). [CrossRef] [PubMed]

**3. **S. G. Johnson and M. Ibanescu, et al. “Low-loss asymptotically single-mode propagation in large-core OmniGuide fibers,” Opt. Express **9**, 748–79 (2001). [CrossRef] [PubMed]

**4. **John Fini and Ryan Bise. “Progress in fabrication and modeling of microstructured optical fiber,” Jap. J. App. Phys. **43**, 5717–5730 (2004). [CrossRef]

**5. **A. Ferrando and E. Silvestre, et al. “Full vector analysis of a realistic photonic crystal fiber,” Opt. Lett. **24**, 276–8 (1999). [CrossRef]

**6. **T. P. White, R. C. McPhedran, L. C. Botten, G. H. Smith, and C. M. deSterke. “Calculations of air-guiding modes in photonic crystal fibers using the multipole method,” Opt. Express, **9**, 721–32 (2001). [CrossRef]

**7. **F. Brechet and J. Marcou, et al. “Complete analysis of the characteristics of propagation into photonic crystal fibers, by the finite element method,” Opt. Fiber Tech. **6**, 181–191 (2000). [CrossRef]

**8. **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–92 (2003). [CrossRef]

**9. **S. G. Johnson and J. D. Joannopoulos. “Block-iterative frequency-domain methods for Maxwell’s equations in a planewave basis,” Opt. Express **8**, 173–190 (2001). software available at http://ab-initio.mit.edu/mpb. [CrossRef] [PubMed]

**10. **T. A. Birks, J. C. Knight, and P. S. J. Russell. “Endlessly single-mode photonic crystal fiber,” Opt. Lett. **22**, 961–3 (1997). [CrossRef] [PubMed]

**11. **T. Hasegawa, T. Saitoh, D. Nishioka, E. Sasaoka, and T. Hosoya. “Bend-insensitive single-mode holey fibre with SMF compatibility for optical wiring applications,” In European Conference on Optical Communications, paper We2.7.3, (2003).

**12. **James A. West, Charlene M. Smith, Nicholas F. Borrelli, Douglas C. Allan, and Karl W. Koch. “Surface modes in air-core photonic band-gap fibers,” Opt. Express, **12**, 1485–96 (2004). [CrossRef]

**13. **A. Cucinotta, S. Selleri, L. Vincetti, and M. Zoboli. “Perturbation analysis of dispersion properties in photonic crystal fibers through the finite element method,” J. Lightwave Technol. **20**, 1433–42 (2002). [CrossRef]

**14. **A. Peyrilloux, T. Chartier, L. Berthelot, A. Hideur, G. Mélin, S. Lempereur, D. Pagnoux, and P. Roy. “Thoeretical and experimental study of the birefringence of a photonic crystal fiber,” J. Lightwave Technol. **21**, 536–9(2003). [CrossRef]

**15. **I. K. Hwang, Y. J. Lee, and Y. H. Lee. “Birefringence induced by irregular structure in photonic crystal fiber,” Opt. Express **11**, 2799–2806 (2003). [CrossRef] [PubMed]

**16. **Steven G. Johnson, M. Ibanescu, M. A. Skorobogatiy, O. Weisberg, J. D. Joannopoulos, and Y. Fink. “Perturbation theory for Maxwell’s equations with shifting material boundaries,” Phys. Rev. E **65**, 066611 (2002). [CrossRef]

**17. **M. J. Steel, T. P. White, C.Martijn de Sterke, R. C. McPhedran, and L. C. Botten. “Symmetry and degeneracy in microstructured optical fibers,” Opt. Lett. **26**, 488–91 (2001). [CrossRef]

**18. **J. M. Fini. “Improved symmetry analysis of many-moded microstructure optical fibers,” J. Opt. Soc. Am. B, **21**, 1431–6 (2004). [CrossRef]

**19. **J. M. Fini. “Perturbative modeling of irregularities in microstructure optical fibers,” In *Conference on Lasers and Electro-Optics (CLEO)*, TOPS vol. 96, paper CThX6, (Optical Society of America, Washington, D.C.,2004).

**20. **M. Koshiba and K. Saitoh. “Polarization-dependent confinement losses in actual holey fibers,” Photon. Technol. Lett. **15**, 691–3 (2003). [CrossRef]

**21. **Paul R. McIsaac. “Symmetry-induced modal characteristics of uniform waveguides-I: Summary of results,” Microwave Theory and Techniques **23**, 421–9 (1975). [CrossRef]

**22. **M. J. Steel and R. M. Osgood. “Elliptical-hole photonic crystal fibers,” Opt. Lett. **26**, 229–31 (2001). [CrossRef]

**23. **M. J. Steel and R. M. Osgood. “Polarization and dispersive properties of elliptical-hole photonic crystal fibers,” J. Lightwave Technol. **19**, 495–503 (2001). [CrossRef]

**24. **Charlene M. Smith, Natesan Venkataraman, Michael T. Gallagher, Dirk Müller, James A. West, Nicholas F. Borrelli, Douglas C. Allan, and Karl W. Koch. “Low loss hollow-core silica/air photonic bandgap fibre,” Nature, **424**657–9, (2003). [CrossRef]

**25. **William H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. *Numerical recipes in C, the art of scientific computing*. (Cambridge University Press, New York,1992).

**26. **C. C. Su. “A surface integral equations method for homogeneous optical fibers and coupled image lines of arbitrary cross sections,” Microwave. Theory and Technol. , **33**, 1114–9, (1985). [CrossRef]

**27. **R. Bise and R. S. Windeler, et al. “Tunable photonic band gap fiber,” In *Optic al Fiber Communications Conference (OFC)*, TOPS vol. 70, paper ThK3, (Optical Society of America, Washington, D.C., 2002).

**28. **R. Bise and D. Trevor. “Sol-gel-derived microstructured fibers: fabrication and characterization,” To appear in *Optical Fiber Communications Conference (OFC)*, (Optical Society of America, Washington, D.C.,2005).

**29. **T. Hasegawa, E. Sasaoka, M. Onishi, M. Nishimura, Y. Tsuji, and M. Koshiba. “Hole-assisted lightguide fiber for large anomalous dispersion and low optical loss,” Opt. Express **9**, 681–6, (2001). [CrossRef] [PubMed]

**30. **Alexander Argyros, Ian M. Bassett, Martijn A. van Eijkelenborg, M.C.J. Large, Joseph Zagari, N.A.P. Nicorovici, Ross C. McPhedran, and C.Martijn de Sterke. “Ring structures in microstructured polymer optical fibres,” Opt. Express **9**, 813–20, (2001). [CrossRef] [PubMed]