Efficient optimization of hollow-core photonic crystal fiber design using the finite-element method

We employ a ﬁnite-element (FE) solver with adaptive grid reﬁnement to model hollow-core photonic crystal ﬁbers (HC-PCFs) whose core is formed from 19 omitted cladding unit cells. We optimize the complete ﬁber geometry for minimal ﬁeld intensity at material/air interfaces, which indicates low loss and high damage threshold, using multidimensional optimization. The optimal design shows a 99.8 % power fraction within the air and an overlap with a Gaussian mode of 96.9 %. [DOI: 10.2971/jeos.2006.06011]


I n t r o d u c t i o n
Photonic crystal fibers (PCFs) can exhibit a photonic bandgap that enables light to be guided predominantly in air. Recently, hollow-core PCFs (HC-PCFs) have attracted attention due to their low loss and nonlinearity [1] and their ability to convey laser pulses at very high peak power [2]. In order to predict the loss and optimize fiber production, accurate numerical modeling of the fiber is crucial. The oldest and still frequently applied numerical method for calculating mode solutions in these fibers is the plane-wave expansion (PWE) method [1,3,4]. Subsequently, a large variety of simulation methods have been devised; see [5] for a review.
One of the most promising solution methods for mode field profiles in HC-PCFs is the finite-element method (FEM) [6]- [11]. In this work, we use the commercial FEM program JCMmode [12]- [14] which employs higher-order vectorial elements and automatic triangulation with adaptive iterative grid refinement, resulting in the application of an optimized basis to the fiber geometry. Adaptive triangulation, being controlled by a residuum-based error estimator, strongly improves the accuracy of the computed field and results in the glass/air interfaces becoming particularly closely sampled by triangle edges. This is because fast field variations occur near the interfaces, and the normal component of the electric field suffers a discontinuity across the interfaces. The increased accuracy is particularly important when computing loss due to mode coupling induced by interface roughness [15,16], as we will demonstrate in this work.
To solve the algebraic equations appearing in the FEM discretization of the Maxwell eigenproblem, JCMmode uses an inverse shifted Arnoldi subspace iteration method with a multi-level preconditioner. For matrix inversion the sparse LU solver PARDISO is employed [17]. The exact description of the solver is beyond the scope of this paper and the reader is referred to [12] for details. Computing a single fiber mode solution with two refinement iterations typically requires only about 40 seconds on a fast workstation, making the method fast enough for automatic multidimensional parameter optimization. Our purpose is to optimize the design of a HC-PCF whose hollow core corresponds to 19 omitted cladding cells with a target minimum attenuation of less than 30 dB/km at 589 nm, for use in high-power pulsed laser sodium guide star beam relays [18]. We chose 19-cell core HC-PCFs since we expect 7-cell (or smaller) core fibers to guide a larger fraction of light in glass with a larger field intensity existing at the material/air interfaces. These smaller-core fibers will therefore show higher loss [19] and a lower power threshold for the onset of nonlinear processes such as stimulated Brillouin scattering (SBS). On the other hand, fibers with larger cores, such as 37-cell, do not feature quasi-single mode behavior and suffer from higher macrobending losses [16].
The contribution of the current work is to obtain a realistic design of a 19-cell core HC-PCF that has been optimized to reduce the field intensity at the glass/air interfaces, and as such is expected to show lower-loss than has so far been demonstrated. The obtained design also yields enhanced field exclusion from the silica and a high overlap with a Gaussian incident mode, indicating an improved power handling capability. The improved design was found using a robust automatic multidimensional optimization procedure for the complete HC-PCF geometry. To our knowledge this is the first demonstration of automatic optimization of a realistic HC-PCF geometry.
We describe the simulation details in Section 2, the simulation validation and optimization of a HC-PCF in Section 3, and we conclude in Section 4.

. 1 G e o m e t r y
The accurate description of the fiber geometry is essential to obtain accurate simulation results. We assume that the fiber cross section consists of undoped silica and air, separated by smooth interfaces. The HC-PCFs under study contain a dodecagonal core surrounded by a glass core wall, in turn surrounded by seven rings of hexagonal cladding cells. 1 (a) The simulation domain, which equals one quarter of the fiber cross section.
The boundary conditions are shown along the four edges, where E and H denote the electric and magnetic tangential field components, respectively. In addition, the intensity of the HE 11 -like mode is overlaid in color. (b) Enlarged fiber detail, showing the four fiber geometry parameters: pitch Λ, strut thickness w, cladding meniscus radius r, and core surround thickness t.
The core replaces the central 19 cladding cells, and the core surround vertices are located on a circle. Figure 1a shows a plot of a fiber cross section where silica is black. Figure 1b details the following four fiber geometry parameters: pitch Λ (cladding lattice constant, measured between the centers of opposite struts in a cell), strut thickness w, cladding meniscus radius r [20], and core surround thickness t [21].
Fibers of the kind shown in Figure1a possess six rotational reflection planes in addition to the point symmetry about the origin and hence belong to the point group C 6ν [22]- [25]. Nondegenerate eigenmodes of the fiber can be found by computing the field in just one 30 o symmetry sector and continuing the field by repeated reflection and rotation. The most relevant (fundamental) mode is HE 11 -like and so is two-fold degenerate. It can be conveniently computed in a 90 o sector, around which we choose the following boundary conditions: short circuit (zero tangential electric field) on one inner edge and the outer edges and open circuit (zero tangential magnetic field) on the other inner edge, as shown in Figure 1a. We note that this choice of boundary conditions excludes numerous modes such as TE 01 , HE 21 , TM 01 and many higher-order modes, which are not the focus of our present study since they are not so well concentrated in the central core and consequently experience much higher attenuation and nonlinearity than HE 11 .
Low-loss HC-PCFs contain numerous air holes separated by delicate glass struts and typically have air filling fractions of well over 90%. The exact thickness and shape of the strut intersections have a large influence on the mode properties [20,21]. JCMmode automatically discretizes the fiber simulation domain by a triangular grid. Expressing the geometry in a vector format rather than in a bitmap ensures the glass/air interfaces are arbitrarily precisely reproduced in the triangulation. To automatically optimize fiber parameters, it is desirable to express the interface description as a polygon list that can be directly read into the mode solver. We have implemented the interface renderer as a Matlab script that is controlled by the geometry parameters of Figure 1b.

. L o s s E s t i m a t i o n
The overall attenuation of light power in the HE 11 mode of many current HC-PCFs is dominated by scattering loss due to roughness at glass/air interfaces [15,16]. The origins of the surface roughness are mainly due to thermal surface capillary waves (SCWs) that become frozen in when the silica solidifies during the fiber draw. The scattering will be decreased as the field intensity at the interfaces is reduced. We use a relative measure of the scattering strength provided by the quantity F defined as where ε 0 and µ 0 are the vacuum permittivity and permeability, respectively,ẑ is the unit vector along the fiber axis, E and H are the electric and magnetic fields within the fiber crosssection, and × denotes the vector cross product. In the numerator, the E-field is evaluated along a path just inside each hole interface, i.e. within the air. The integral in the denominator runs over the fiber simulation domain of area A c . Since reduced field strength at the interfaces is closely associated with a lowering of the power residing within the glass [21], designs which minimize F will also show an enhanced damage thresh-old in high-power pulsed laser applications [2]. Consequently, we employ F as the figure of merit in the simulation.
The outer boundary conditions in Figure 1a imply a perfect conductor. In reality, the cladding is surrounded by bulk fused silica with a thickness of about 50 microns. This situation can be better modeled applying so-called transparent boundary conditions to the outer computational boundaries, emulating an infinitely wide silica coating. Transparent boundary conditions are realized with the perfectly matched layer (PML) method in JCMmode [6,12]. The field decaying in the PML along rays pointing outwards is discretized with typically around ten finite elements of second polynomial order on each ray [12]. The PML leads to complex eigenvalues and leaky eigenmodes which radiate to the exterior. We have run extensive comparisons and found that in the parameter regimes considered throughout this paper, both the mode field solution of the HE 11 -like mode and the real part of its eigenvalue do not change to numerical accuracy when PMLs are introduced within JCMmode [26]. Moreover, SCW-induced loss in dB scales like the inverse cube of the fiber design wavelength [16], while confinement loss in dB scales only like the inverse wavelength, hence the relative importance of confinement loss decreases quickly with reduced wavelength. Attenuation of the HE 11 -like mode in fiber structures close to our considered 19-cell cores design is already dominated by SCW-induced loss at 1550 nm [19], and hence we will neglect confinement loss in this study.
Given that HC-PCFs with a similar (but not fully optimized) design to the ones we are considering in the present work have been produced with an attenuation of 1.2 dB/km at around 1550 nm [15], we expect that the loss at 589 nm can in principle be lowered to below 22 dB/km. Achieving this goal requires extreme production precision due to the reduced feature sizes within the fiber cross-section that are required at the shorter wavelength.

R E S U L T S A N D D I S C U S S I O N
In this section, we present the convergence and scaling behavior of the solutions and finally report the results of optimization runs.
3 . 1 V a l i d a t i o n , c o n v e r g e n c e a n d s c a l i n g In order to validate the solutions of our FE-solver, we have investigated its convergence behavior with the number of included triangular elements and grid refinement steps. We have further compared its mode field solutions with those of the PWE solver used in [16,21,27] which operates on a bitmap approximation to the considered geometry.
The computational procedure starts by producing an initial geometry list. JCMmode then creates the initial triangular grid. Next, this grid is recursively refined twice, using adaptive grid refinement based on a local error criterion. The final solution is the mode field whose eigenvalue n eff is closest to the user-provided estimate n eff,est , in our case the effective re-fractive index of the HE 11 -like mode. In a postprocessing step, JCMmode computes the two integrals appearing in Eq. (1), yielding the value of F. Finally, the mode field diameter d MFD of the field and the fraction of light intensity in glass η can be determined. The refractive index of silica at a wavelength of 589 nm equals 1.4585. Figure 2 shows the intensity |E| 2 of the optimized HE 11 -like mode along the x-axis cutting though the center of the fiber along the lower edge in Figure 1a, on a logarithmic scale (the optimized structure will be defined in the Section 3.3). The gray vertical bars represent the glass along the slice. In order to reduce pixilation effects when using the PWE method, we linearly interpolate the refractive index at Fourier samples that are intersected by the glass/air interfaces, displayed as different shades of gray. In the FEM simulation, such interpolation is not necessary. In order to obtain a symmetric mode profile from the FEM and PWE solutions, we used the symmetrized electric mode field E = E vert + i E horiz , where E vert equals the vertically polarized mode field found with the boundary conditions shown in Figure 1a, E horiz is the corresponding horizontally polarized mode where the boundary conditions on the inner edges are exchanged, and i is the imaginary unit.
The modes E vert and E horiz form a 2-fold degenerate pair of HE 11 -like modes which belong to the same eigenvalue, and |E| 2 exhibits the full symmetry of the considered fiber structure.
The intensity profile is close to a Gaussian within the fiber core. As expected, the FEM curve is closer to the PWE 1024 curve than PWE 512 along most of the slice. This behavior reverses near the right end of the plot, where the PWE 1024 curve deviates slightly from the FEM curve (however note the logarithmic scale). We attribute this deviation to the different boundary conditions -the PWE simulation was computed on a rectangular supercell with an x range of ±5Λ and y range of ±3 √ 3Λ and naturally has periodic boundary conditions, while the FEM simulation domain is not periodic and hence allows the field to decay faster. 1 The relative difference in |E| 2 between the FEM and the PWE 1024 and PWE 512 curves, evaluated on the outer glass/air interface of the core surround near x = 2.34Λ, where the contribution to F in this slice is maximal, is 8% and 24%, respectively. In the limit of large N, number of Fourier basis-waves and computational domain size, the curves should converge. The FE-method yields FΛ 0 = 0.0646 and light intensity in glass η = 0.00171, and the PWE 1024 method gives FΛ 0 = 0.0650 and η = 0.00169, corresponding to 0.7 % and −1.3 % deviation from FEM, respectively. In order to obtain the PWE result of FΛ 0 , one has to rescale the E-component of the PWE solution normal to the interfaces according to the intersection fraction of the pixel with the interface. Without this rescaling procedure, we obtain FΛ 0 = 0.049 (24% deviation from FEM). We have also compared the FEsolver to our PWE solver for a variety of other parameter combinations and found similar results. Our PWE algorithm uses a targeted eigensolver, but is otherwise not well optimized for speed. The PWE 1024 solution requires about 8 hours of CPU time, and we expect more optimized solvers such as MIT Photonic-Bands [3] to be only about a factor of 2 faster, while attaining the same accuracy. This large discrepancy is only in part due to the employed PWE method not fully taking advantage of the fiber symmetry, but is mostly due to finite elements being much more efficient at sampling complicated, discontinuous structures. Figure 3a shows the dependence of the relative error of the effective index, ∆n eff = |n eff − n eff,qe |/n eff,qe on the number of unknowns N of the finite-element expansion of the solution on the triangular grid. The effective refractive index of the quasi-exact solution, n eff,qe is obtained from a well-converged solution computed on a fine grid with about N = 6 × 10 6 . The four curves depict the convergence behavior with linear and quadratic field interpolating functions (elements) in each triangle and for adaptive grid refinement as well as for regular, non-adaptive grid refinement. Linear finite elements are parameterized by 6 unknowns, quadratic elements by 14. The solver computed the HE 11 -like mode that is displayed in Figure 1a in color. With adaptive grid refinement (♦, •), an automatic a-posteriori error estimator decides which triangles of the finite-element grid are further refined. In contrast, with regular grid refinement (×, ) every triangle is refined in each refinement step. Adaptive grid refinement results in dramatically reduced computation times for a specified accuracy, even though the asymptotic convergence rate for large N depends only on the polynomial order of the elements. On the other hand, it is easy to see that quadratic elements ( , •) have a much better convergence rate than linear elements (×, ♦). Consequently, the quadratic adaptive method with quadratic finite elements (•) was used in the optimization study presented in the next subsection. Additionally, the computation times for this simulation are displayed as numbers next to the graph (solution of the algebraic problem, without triangulation and without post-processing, computed on a 64-bit workstation with 4 dual-core AMD Opteron processors, 2.2 GHz, and 32 GB RAM).  The quadratic-element adaptive mesh curve is very well approximated by a straight line with the function 9.469×10 −5 × N 1.0185 , i.e. the method behaves approximately as an "order N" method. The JCMmode eigensolver exploits the sparseness of the wave equation in the FE basis to enable the fast evaluation. The number of unknowns N rivals or betters that attainable using an FFT-based plane-wave scheme, but the FE basis is inherently superior in dealing with complicated geometries and the electromagnetic field boundary conditions which exist at interfaces. The peak RAM requirement scales linearly with N and equals 4,527 Bytes per unknown or approximately N/237,000 in GB. All simulations in this work except those shown in Figure 3 had about N = 340,000 unknowns and consumed 1.5 GB RAM (third data point on the red curves with circles in Figure 3a; initial triangulation followed by two adaptive refinement steps). Figure 4 shows the HE 11 -like mode in a small region of the fiber rendered as a 3D-plot. 1 The E || = 0 boundary conditions that we choose on the outer boundary imply reflections. However, our FEM simulation domain contains more cladding layers (7) than the PWE supercell, so that between x = 0 and about x = 4Λ the decay of |E| 2 is found to remain the same if PMLs are used instead. Height is proportional to the optical intensity |E| 2 . The intensity of the light in air is also color-coded and shaded, while the intensity in glass in shown in gray. Superimposed is the triangular grid of the final refinement step. The triangles tend to be smallest near the glass/air interfaces, in particular near the glass nodes where three struts are joined. The field discontinuities lead to gaps in the 3D intensity plots and are most pronounced where the electric field vector is oriented almost normal to the interfaces. Since the major contributions of the line integral in the numerator of Eq. (1) stem from these regions, the FE-method is particularly well suited to the computation of F and hence of surface roughness induced attenuation.

. 2 M u l t i d i m e n s i o n a l o p t i m i z a t i o n
To efficiently find local or even global minima of F, we use the multidimensional optimization implemented in the Matlab routine fminsearch, which is based on the Nelder-Mead Simplex Method [28]. This algorithm approaches a minimum in n dimensions by enclosing it in an (n+1)-dimensional simplex (an interval in 1D, a triangle in 2D, a tetrahedron in 3D, etc.) and hence does not rely on local derivatives, making it more robust to fast-varying objective functions. As a scan or an optimization proceeds, the eigenvalue of the HE 11 -like mode changes. Hence, even when tracking the eigenvalue by setting n eff,est to the eigenvalue of the previously computed mode, it occasionally happens that the solver locks onto a mode different from HE 11 -like. This situation can be fixed by re-running the solver for multiple modes if it is found that the change in F exceeds a certain threshold and proceeding with the mode that has the lowest value of F.
Although fminsearch implements an unconstrained optimization method, the accessible parameter space can be effectively bounded by having the objective function return a large value (e.g. 10 10 times a typical value of F) instead of F if one of the trial optimization parameters is out of bounds [27].

. O p t i m i z a t i o n r u n s
We have run various optimizations with different starting parameters. All of these optimizations finally terminated at the same fiber geometry parameter set. A typical optimization was started at the parameter set w/Λ = 0.035, r/Λ = 0.2, t/Λ = 0.1 and normalized wavevector kΛ = 14.0 and ended after 55 simplex iterations involving 102 evaluations of F at the parameter set w/Λ = 0.028629, r/Λ = 0.19084, t/Λ = 0.11444 and kΛ = 16.0364 ( Figure 1a shows the corresponding fiber cross section). We imposed the constraints w/Λ ≥ 0.028 and an air filling fraction ≤ 93%, since fibers outside these limits are extremely hard to fabricate. At the wavelength of λ 0 = 2π/k = 589 nm the final optimized geometrical lengths are Λ 0 = 1503 nm, w = 43.2 nm, r = 285 nm, and t = 172 nm. The air filling fraction equals 93%, n eff = 0.997749, η = 0.00171.
If the fiber is irradiated with a perfect Gaussian beam [29] whose waist lies at the fiber end, a maximum coupling efficiency [30] of 96.9% can be achieved with a mode field diameter of 2.822Λ 0 = 4.24 µm (1/e 2 diameter). We note that the fraction of light traveling in glass, η ≈ 0.2%, is extremely low and we expect it to lead to a high damage threshold of the optimized fiber (for comparison, other studies have considered HC-PCFs with η ≈ 2% [31]). In principle, the damage threshold at a given wavelength can be assumed to scale inversely proportional to η, but the exact value depends critically on the fiber geometry, launch conditions, pulse format, and repetition rate of the input laser [2].  Two regions of low FΛ 0 and hence low predicted loss are visible, the wider and deeper one having a minimum of FΛ 0 = 0.0646 at kΛ 0 = 16.0364 and a width which corresponds to 13.2 nm at Λ 0 = 1503 nm (value of FΛ 0 10 % above minimum at the most). The peaks within the bandgap which separate the low-loss regions correspond to anti-crossing events [32]. We note that the proportionality factor arising in the approximate relation between F and the fiber attenuation depends on k [16], and hence the graph in Figure 5 cannot be directly regarded as a wavenumber dependent loss profile.
It is well known that minimization algorithms are prone to becoming stuck in local minima. This problem is particularly serious if the objective function fluctuates strongly as a function of some parameters, as F obviously does as a function of kΛ 0 . In addition, we also found that F fluctuates strongly as a function of t. It is therefore beneficial to scan F as a function these two variables independently, and then to start the Nelder-Mead optimization at the parameters producing the lowest F. This procedure can also be repeated iteratively. We have employed it to find the optimization starting parameters mentioned in the beginning of this subsection. Subsequent extensive parameter scans have produced no lower value of F than FΛ 0 = 0.0646. Hence we consider it likely that the parameter set given above points to a global minimum in F, for the rather broad class of practical 19-cell core HC-PCF geometries that is spanned by the geometrical parameterization we have employed, with w/Λ ≥ 0.028, an air filling fraction ≤ 93%, and a refractive index of glass near 1.46. We will report on detailed studies of the parameter space in a forthcoming publication.
Although the optimization is performed at a single target wavelength, the obtained optimized structure will necessarily show a reasonably broad bandwidth for low-loss transmission (in the example above, the bandwidth is 13.2 nm centered near 589 nm). This is because some coupling occurs between a core mode and a core interface mode of compatible symmetry, which would lead to a non-optimum interface field intensity if the associated mode-crossing is nearby in wavelength, as detailed in [27]. Indeed, the bandwidth for the optimized design is found to be close to a local maximum; small changes in the fiber geometry from this design lead to a decrease in useable bandwidth. Note, however, that the separation between mode-crossing events for substantially thinner core walls (smaller t) is typically larger than for the design optimized for minimum interface field intensity, as discussed in [21].
The cost function to be minimized in the optimization can be changed to impose different desired aspects to the fiber performance. For example, by introducing an integration over a target wavelength range, a broader bandwidth for low-loss propagation can be enforced. Of course the minimum loss will be compromised compared to a single wavelength optimization, but the benefits of an increased bandwidth may in some circumstances outweigh the cost in loss. Other aspects can be naturally be built into the cost function, such as a required dispersion characteristic. Automatic optimization of a variety of HC-PCF attributes is currently under investigation.

C o n c l u s i o n s
We compute the electric and magnetic fields of modes in a 19-cell core HC-PCF using the commercial FEM solver JCMmode [12]- [14]. Taking advantage of the high computational speed of the solver, we drive it by a multidimensional optimization algorithm to obtain a realistic fiber geometry in which the field intensity at the material interfaces and hence the expected loss is minimized. The glass strut thickness, meniscus radius, core surround thickness and cladding lattice constant are all allowed to vary during the optimization procedure. As a result, we compute optimal parameters for currently producible 19-cell core HC-PCFs, which to our knowledge is the first complete geometry optimization in HC-PCFs. The optimal design shows a 99.8 % power fraction within the air and an overlap with a Gaussian mode of 96.9 %. The FEmethod is mature, reliable, and very efficient in modeling HC-PCFs, to the degree that it allows us to study the fiber parameter space using multidimensional optimization.