## Abstract

Photonic crystal fibers are well-known to offer a number of unusual properties, including supercontinuum generation, large mode-areas and controllable dispersion behavior. Their manufacturability would be enhanced by a more detailed understanding of how small perturbations in the fiber’s geometric structure cause variations in the fiber’s fundamental modes. In this paper, we demonstrate that such sensitivity analysis is feasible using highly accurate boundary integral techniques.

© 2004 Optical Society of America

## 1. Introduction

Photonic crystal fibers (PCF) are optical fibers with a regular pattern of holes that run parallel to the fiber axis and radially confine light [1, 2]. Light propagates parallel to the holes in these fibers and is guided either by total internal reflection or by a photonic band gap. Although enormous progress had been made in creating a variety of PCF’s in recent years, design and “design for manufacture” remain non-trivial.

In manufacturing microstructured devices like PCF’s, it would be helpful to have a full understanding of the relationship between variations in the manufactured device and variations in the property of the device. This would illuminate the manufacturing process, aid in the development of appropriate diagnostics, aid in acceptance sampling of manufactured devices, etc. The simplest form of such an understanding would be a sensitivity analysis, which would show the linear relationship between device properties and infinitesimal changes in the PCF geometry. One can certainly perform a sensitivity analysis using experimental methods, but this would both expensive and time consuming - it would clearly be preferable to perform such analyses through numerical simulation. There has been some previous work in this direction. In [3], for example, the authors carry out a careful analysis of dispersion properties by direct computation of the dispersion curves over a range of wavelengths for a selection of material properties (such as filling fraction). Our approach to sensitivity analysis is somewhat different - we are interested here in the direct calculation of the derivative of the optical property with respect to variation in some parameter.

The particular difficulty with numerical simulation which will concern us here is the requirement for accuracy. Numerical differentiation requires dividing a potentially noisy finite difference by a small number. Unless the noise can be adequately suppressed, the derivative will have no significant digits. This problem is exacerbated when the quantity of interest is small, e.g. the imaginary part of the effective index. In this setting, the numerical method must be very accurate.

Among the most accurate numerical methods used for simulating PCF’s is the multipole expansion method [4, 5]. However, it is currently restricted to circular holes only, which obviously misses a large class of manufacturing defects. We have recently developed a method [6] for simulating PCF’s based on an integral equation formulation of Maxwell’s equations and specialized high-order quadrature rules. It is accurate and applicable to arbitrarily-shaped holes. In this note we demonstrate its application to the problem of computing the sensitivity of the fiber modes to geometric variations.

## 2. Model problem and convergence analysis

For purposes of demonstration, we will adopt a specific model of a PCF. We do not wish to imply that this model is of particular importance, only that it can be used to demonstrate the method. The model problem (Fig. 1) is a solid core PCF with confinement provided by 22 circular holes arranged in two concentric layers. The model is similar to a configuration considered for polarization maintaining PCF [7]. Two holes are omitted in the center, which provides a defect in the lattice and permits propagating modes with small imaginary parts. The index of the exterior material is 1.45 while the index of the holes is 1.0. The holes are placed in a triangular lattice with pitch Λ. The diameter *d* of the air-holes is 0.7Λ. All the simulations reported here take place at a wavelength λ = 1.1Λ.

In this paper we are concerned with the response of the two dominant modes to small changes in the geometry of the fiber. It is first necessary to understand how accurately we can compute the effective index *n _{eff}*. This can be determined through a convergence study, i.e. calculating the modes with increasing numbers of degrees of freedom. The degrees of freedom are electric and magnetic charge densities located at equispaced quadrature points on the holes[6]. Our method uses 4 charge densities at each quadrature point, i.e. with 32 quadrature points on each curve, there are 128 degrees of freedom on each curve. Because the geometric perturbations are very small, we can use the convergence study of the unperturbed hole pattern as a guide to the accuracy of the perturbed hole patterns.

Figure 2 shows the results of the convergence study for the unperturbed pattern of holes. The vertical axis is the absolute value of the difference with the result computed from 40 quadrature points per hole. In this particular model problem the error in the imaginary part of the effective index is always smaller than the error in the real part, so the error shown in Fig. 2 is dominated by the real part. The authors know of no theorem that this will true for all models. In the following analysis, we will be conservative and not assume that the error is dominated by its real part.

The convergence in Fig. 2 is approximately exponential in the number of quadrature points. This is expected for the method of [6] where the order of the quadrature increases as quadrature points are added. In the remainder of this paper we will use 32 quadrature points per hole and it will be understood that *n _{eff}* will have an error less than 10

^{-10}. Table 1 shows the effective indices of the two dominant modes. In the following we will refer to the mode with the larger effective index as mode 1 and the other as mode 2. We will refer to the effective index of mode 1 as

*n*

_{1}and the other as

*n*

_{2}. The modal confinement loss (in dB/m) is related to the imaginary part of the effective index via the formula

with the vacuum wavelength *λ* given in nanometers.

## 3. Sensitivity analysis

In the following, we will analyze the sensitivity of the effective indices *n*
_{1} and *n*
_{2} to geometric variations. We will individually vary the diameter, position, and eccentricity of the holes and calculate the response of the effective indices to that variation. It should be emphasized that although the unperturbed configuration shown in Fig. 1 has two planes of reflection symmetry, the perturbations considered here break the reflection symmetries. If the perturbations had been symmetric, it would have been possible to reduce the number of degrees of freedom by a factor of four.

We will compute derivatives of the effective indices by the application of a standard two point centered finite difference formula. It is important to understand the sources of error in the finite difference calculation so that we may estimate the accuracy of the sensitivity analysis. There are two sources of error in the finite difference approximation: numerical noise in the effective index and higher order derivative terms. More precisely, the centered difference expression yields the true derivative plus two error terms:

where *ε* is the numerical error in the computation of *n*, *h* is the size of the perturbation, and *x* is any geometric parameter of interest. In applying this formula to the current problem, we know that *ε* has a magnitude less than 10^{-10}, but its phase may be random.

Let us now choose *h* = 10^{-3}. Then the last term in Eq. (2) is bounded by 10^{-7} and it remains only to estimate the second to last term in Eq. (2) to be confident that our computation of $\frac{\partial n}{\partial x}$ is accurate. We have computed the third derivative of the effective indices with respect to the diameter of holes 1 and 11, based on a four point finite difference formula. Table 2 shows the results. If we choose *h* to be 10^{-3}, then for hole 1 the absolute value of the estimated error of *∂n*/*∂x* is approximately 10^{-7}. The corresponding uncertainties for hole 11 are also 10^{-7} as the last term in Eq. (2) dominates when the third derivative is small. We did not repeat this third derivative analysis for other holes or parameters - the errors appear to be safely smaller than the derivative estimates which we report below.

Table 3 shows the sensitivity of the modal indices to variations in the diameters of the airholes. Because of the two symmetry planes of the unperturbed configuration, we restrict our attention to 7 of the original 22 holes.

The sensitivity vector displayed in Table 3 confirms certain intuitions. For example, all of the imaginary parts are negative, i.e. increasing hole diameters will reduce mode leakage. Holes in certain positions have more effect on the real parts than holes in other positions. For example, holes 1 and 12 are located on the long side of the confining defects. They have more effect on the real parts of the second mode than the first. Hole 3, which is located on the short side, affects the real part of mode 1 more than mode 2. Some observations are less easy to explain: hole 12 has the strongest effect upon the imaginary part of both modes.

An important question is how small perturbations must be in order for a linear sensitivity analysis to be applicable. The effective index may be expanded in a Taylor series of the model parameters about the unperturbed model. The sensitivity vector in Table 3 gives the first order deviation. As the perturbation grows larger, second and higher order terms in the Taylor series become appreciable. The point at which the linear terms are dominated by higher order terms is determined by the size of the derivatives. The derivatives can, of course, be evaluated by application of finite differences.

To answer this question fully is beyond the scope of the current paper because it involves all possible cross-interactions between parameters. However, we may compute a crude estimate by considering only a single parameter, the diameter of hole 1, and restricting consideration to only linear and quadratic terms in the Taylor expansion. In the same computation that computed the centered difference approximation to *∂n _{i}*/

*∂d*

_{1}, we have computed the second derivitive for

*n*

_{1}to be -8.7 × 10

^{-2}+ 7.4 × 10

^{-5}

*i*and for

*n*

_{2}to be -8.9 × 10

^{-2}+ 1.2 × 10

^{-4}

*i*. Clearly the crossover between linear and quadratic domination will come when

Now let us substitute the values for the derivatives of *n*
_{1} which we have computed into Eq. (3). It is interesting to consider the real and imaginary parts separately. We find that the real part of *n*
_{1} will be nearly linear as long as the perturbation in *d*
_{1} is small compared to Λ. But the imaginary part puts a stronger constraint on the perturbation: the perturbation should be small compared to 0.045Λ. The constraints are roughly similar for *n*
_{2}, and for the crossover between linear and cubic terms in the Taylor expansion. If Λ is one micron, then the linear sensitivity analysis appears to be valid for diameter excursions up to several tens of nanometers.

Tables 4 and 5 shows the sensitivity of the effective indices to shifting the centers of individual holes in the positive x and y direction respectively, where x and y are the horizontal and vertical directions in Fig. 1. Some of the reported derivatives are essentially zero in agreement with the symmetry of the unperturbed model. For these modes only the order of magnitude are shown. On the average, the holes in the second layer have less influence on the modes.

Table 6 shows the sensitivity to changing the shape from a perfectly circular shape to an ellipsoidal shape. In this case we wish to separate the shape dependence of the curve from the area dependence of the curve. We shall therefore chose a parametrization of the ellipses which preserves the enclosed area. The ellipse is given as the solution to

In other words, the sizes of the semi-major and semi-minor axes are inversely proportional. Positive ε corresponds to an elongation in the x-direction and negative ε creates an elongation in the y-direction. On the whole, shape deformation creates less change in the effective indices than changing the area of the holes. One interesting regularity is that the real parts of the derivatives always have the same sign for modes 1 and 2. The imaginary parts frequently have opposite signs.

## 4. Conclusions

We have demonstrated the feasibility of sensitivity analysis in a simple but non-trivial photonic crystal fiber geometry. Obviously the same approach could compute the sensitivity with respect to the index of the material in the holes, although we have not done so here. The crucial requirement for the calculation is the accuracy of the method and the ability to handle all relevant geometric perturbations. We will consider more complex design problems in hollow-core (photonic band gap) fiber at a later date.

## References and links

**1. **P.S.J. Russel, “Photonic crystal fibers,” Science **299**, 358–362 (2003). [CrossRef]

**2. **B.J. Eggleton, C. Kerbage, PS. Westbrook, R.S. Windeler, and A. Hale, “Microstructured optical fiber devices,” Opt. Express **9**, 698–713 (2001), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-9-13-698 [CrossRef] [PubMed]

**3. **A. Ferrando, E. Sylestre, and P. Andres, “Designing the properties of dipsersion-flattened photonic crystal fibers,” Opt. Express **9**, 687–697 (2001), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-9-13-687 [CrossRef] [PubMed]

**4. **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]

**5. **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]

**6. **H. Cheng, W.Y. 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.opticsexpress.org/abstract.cfm?URI=OPEX-12-16-3791 [CrossRef] [PubMed]

**7. **T.P. Hansen, J. Broeng, S.E.B. Libori, E. Knudsen, A. Bjarklev, J.R. Jensen, and H. Simonsen, “Highly birefrin-gent index-guiding photonic crystal fibers,” IEEE Photonics Technology Letters **13**, 588–90 (2001). [CrossRef]