Analytical mode normalization and resonant state expansion for optical fibers - an efficient tool to model transverse disorder

We adapt the resonant state expansion to optical fibers such as capillary and photonic crystal fibers. As a key requirement of the resonant state expansion and any related perturbative approach, we derive the correct analytical normalization for all modes of these fiber structures, including leaky modes that radiate energy perpendicular to the direction of propagation and have fields that grow with distance from the fiber core. Based on the normalized fiber modes, an eigenvalue equation is derived that allows for calculating the influence of small and large perturbations such as structural disorder on the guiding properties. This is demonstrated for two test systems: a capillary fiber and an endlessly single mode fiber.


Introduction
Photonic crystal fibers guide light in a central defect core surrounded by a periodic cladding [1]. The guiding mechanism of the photonic crystal fiber can be a bandgap effect or modified total internal reflection in cases where the index of the core is larger than the effective cladding index. These fibers feature a high degree of light confinement, highly tunable dispersion properties [2], and single mode operation [3]. Photonic crystal fibers are extensively used in gas sensing [4], nonlinear optics such as supercontinuum generation [5], and many more applications [6][7][8][9][10][11].
In theoretical investigations, an ideal cladding is usually used to analyze such structures, while a fabricated photonic crystal fiber cladding is never truly perfect [12,13]. The fabrication process itself gives rise to shape and position disorders that influence the guiding properties. Studying that influence requires investigating many realizations [14], which is rather tedious in conventional numerical approaches. In contrast, the resonant state expansion has proven rather efficient for investigating a large set of similar three-dimensional resonator systems [15][16][17][18] and slab waveguides [19,20]. The resonant state expansion is a rigorous perturbative approach, in which the resonant states (also known as quasi-normal modes [21,22]) of a reference system are used to setup an eigenvalue equation that provides the resonant states of a perturbed system. Here, we adapt the resonant state expansion to fiber geometries, in which the core and cladding modes constitute the resonant states, and treat disorder as a perturbation of the perfect cladding system. As in any perturbation theory, the key factor in the resonant state expansion is the normalization of the resonant states. The normalization is not trivial, since the solutions of Maxwell's equations include leaky modes [23]. These modes radiate energy perpendicular to the fiber axis and have fields that grow with distance from the fiber core. This is displayed in Fig. 1 for a capillary fiber with air core and silica cladding (a) and a photonic crystal fiber with air inclusions and silica background (b). A lot of work has been devoted to the normalization of leaky modes [24][25][26][27]. The most sophisticated approach is introducing a complex coordinate transformation in the exterior that suppresses the growth [28], which is equivalent to using perfectly matched layers and extending the area of normalization to the perfectly matched layers [21]. In contrast, we derive here an analytical normalization that can be calculated without any perfectly matched layers and is valid for both guided as well as leaky modes. Our new normalization can be easily applied when using standard numerical methods for the calculation of modes.
Here, the properties of the resonant state expansion with our analytical mode normalization is demonstrated for two fiber geometries. In the first example, we use the analytical solutions for a capillary fiber as basis to model the influence of a homogeneous change of the refractive index of the fiber core on the propagation constants of the fiber modes. In the second example, we investigate the influence of diameter disorder on the modal properties of an endlessly single mode fiber [29].

Theory
Maxwell's Equations can be summarized in real space and frequency domain with time dependence exp(−iωt) by the compact operator form [30] with electric and magnetic fields E and H, respectively, permittivity and permeability tensors ε and µ, respectively, and k 0 = ω/c. The right-hand side contains the electric source term J E = −4πij/c with current density j, and the magnetic source term J H that has been introduced for the sake of symmetry. For optical fibers, the permittivity and permeability tensors are translationally symmetric along the direction of propagation, which we choose as the z direction of our coordinate system. Defining the Fourier transform in this direction aŝ with r | | being the projection of r to the xy plane and the hat denoting Fourier transformed quantities, the Fourier transform of Eq. (1) yields The Green's dyadic [31] of Eq. (3) satisfies the relation and provides the solutions¿ of Eq. (3) for a given sourceˆ aŝ The Green's dyadic can be expanded in terms of the resonant states [15][16][17]30,[32][33][34], which are solutions of Eq. (3) in the absence of sources for outgoing boundary conditions with eigenvectorŝ ¿ n and eigenvalues β n :ˆ 0 (r | | ; β n )¿ n = 0.
with ⊗ denoting the outer vector product, and N n being the normalization constant in order to assign the appropriate weight to the resonant states, since Eq. (6) provides the resonant field distribtuions only up to a constant factor. The factor −1/2 has been introduced for later convenience. The superscript R denotes the reciprocal conjugate resonant state, which is a solution of Eq. (6) at −β n . Note that Eq. (7) is only valid within the regions of spatial inhomogeneities of the fiber [36], where the leaky modes do not exhibit any growth. Furthermore, ∆ˆ cuts denotes cut contributions due to branch cuts in the involved analytical functions. In the following, we will focus on the contribution of the resonant states, keeping in mind that we can treat the cut contributions in a similar manner in numerical calculations [17,32]. The derivation of the normalization constant is described in detail in Appendix A. The resulting normalization can be split into two terms comprising of a surface and a line integral that are evaluated on a circle with radius R outside the region of inhomogeneities, which yields with the surface term which is proportional to the integral over the z component of the real-valued Poynting vector, and the line term where the subscript R indicates that the integrand is evaluated at radius R, and For truly guided modes, the resonant states decay outside the regions of spatial inhomogeneities, so that the line term vanishes in the limit of R → ∞. This results in the rather well-known normalization of resonant states by the integral over the z component of the Poynting vector [28]. For leaky modes, both the line and the surface term diverge. However, their sum countervails this divergence, resulting in a normalization constant independent of the radius of normalization, see Fig. 1. Hence, it is possible to calculate the normalization constant for a small area surrounding the regions of spatial inhomogeneities, without the need of including perfectly matched layers [21] or, equivalently, complex coordinates [28]. Furthermore, it should be noted that this approach also simplifies the normalization of truly guided modes in numerical calculations, since it allows to restrict the normalization integrals, and, thus, the computational domain, to a small area.
Using the normalization to gauge the correct weight of the resonances, it is possible to determine the resonant states of a perturbed system (denoted by subscript ν) with perturbation ∆ε and ∆µ that exhibits the same translational symmetry as ε and µ and vanish outside the regions of spatial ingomogeneities. The Maxwell operatorˆ of the perturbed system can be separated into the operatorˆ 0 of the unperturbed system and the deviation ∆ˆ asˆ =ˆ 0 + ∆ˆ , with Thus, we can recast Eq. (6) in the formˆ Using Eq. (5), we therefore obtain Next, we construct the resonant states of the perturbed system as a linear combination of the normalized resonant states of the unperturbed system: Using this ansatz in Eq. (14) and equating it for each¿ n independently, we obtain where The above equations describe a linear eigenvalue problem with β ν as the eigenvalue. Note that the sum in Eq. (15) is carried out over all resonant states of the unperturbed system, but in real calculations, a truncated basis is used to expand¿ ν . The choice of the basis size has to be taken large enough to accurately account for the perturbations in the system.

Results and discussion
We first consider as our unperturbed system a capillary fiber with core index 1 and cladding index 1.44 having a core radius of 8 µm. The values of the propagation constant, and hence, the effective index of the fundamental HE 11 mode as well as those of higher-order modes have been determined analytically by solving their characteristic equation [24, 37, 38] at a wavelength of 1 µm. The fields of the fiber are proportional to Bessel functions inside the core and outgoing Hankel functions in the cladding region. A homogeneous perturbation of ∆n is introduced inside the core of the fiber changing the core index to n core +∆n. As our perturbation is azimuthally symmetric, we only require modes of the same symmetry as the fundamental core mode to set up our eigenvalue problem of Eq. (16). The comparison of the resonant state expansion with the exact analytical solution for the fundamental and higher order modes of azimuthal order m = 1 is shown in Fig. 2 for (a) ∆n = 0.07 and (b) ∆n = 0.17. We can see that there is a good agreement not only for the fundamental mode (indicated by the arrow) but also for the higher order modes of the system. The number of modes used is 154, with pairs of modes with propagation constants β n and −β n . The relative error given by |1 − n RSE eff /n exact eff | is on the order of 10 −6 for the fundamental mode for a core perturbation of ∆n = 0.17. The average relative error of the higher-order modes is on the order of 10 −3 for the same perturbation. As a second example, we consider a silica-air photonic crystal fiber of air holes with radius r 0 = 0.25 µm in four cladding rings with pitch 2.3 µm around a single defect core. We numerically solve for modes of the fiber [39-42] with a perfect cladding structure and use them as basis for a perturbed system, in which we introduce diameter disorder in each and every inclusion in the cladding region. The range of the diameter disorder is determined by the disorder parameter ∆ as r 0 ± ∆. Within that radius range of width 2∆, a uniform distribution of disorder is used. The probability density for a uniform distribution is given as, We set up our eigenvalue problem with 190 modes with effective indices relatively close to the fundamental mode. The comparison of the real and imaginary part of the effective index obtained from the resonant state expansion (red crosses) and full numerical calculations (black circles) can be seen in Fig. 3 (c) and (d), respectively, for 20 realizations and ∆ = 0.1 µm at a wavelength of 1.55 µm. Evidently, there is a good agreement between the two methods for the shown realizations.
In Fig. 4 (a) and (b), we display the real and imaginary parts of the effective index averaged over 200 realizations for disorder parameters ranging from ∆ = 0 to 0.11 µm. More specifically, we generate 200 sets of random numbers between 0 and 1 for each air hole and multiply them with different values of ∆ in order to generate the disordered fibers. The standard deviation of the effective index is plotted as error bars that grow with increasing ∆. Interestingly, the average Re(n eff ) has a linear dependence with ∆ while the Im(n eff ) exhibits a more quadratic behavior.

Conclusion
We have derived an analytical normalization for modes in fiber geometries that is valid not only for guided but also for leaky modes. We have shown that the normalization constant is independent of the radius of integration even for leaky modes with fields that grow with distance from the fiber core. Thus, it is possible to set up an eigenvalue equation that allows us to calculate the effective refractive indices of modes in a perturbed system. The accuracy of this so-called resonant state expansion has been demonstrated for capillary-type and photonic crystal fibers. For the latter, we have studied diameter disorder in the cladding of a silica-air photonic crystal fiber for different disorder parameters averaged over many realizations. Here, the resonant state expansion is clearly superior compared to full numerical simulations, since it does not require to repeatedly solve Maxwell's equations, while the numerical effort for solving the eigenvalue equation is rather low. Thus, it is possible to derive the influence of disorder on the guiding properties such as propagation constant and loss efficiently.
Here, σ n (r ) is chosen to vanish outside the region of spatial inhomogeneities. Taking the source term and convoluting with the Green's dyadic in the limit β → β n , we get ¿ n (r ) = lim This can be only fulfilled for ∫ dr ′ ¿ R n (r ′ )σ n (r ′ ) = −2N n .
To derive the normalization equation, we multiply Eq. (19) with¿ R n (r ) and subtract a zero in the form of 0 =¿(r ; β) ·ˆ 0 (r ; −β n )¿ R n (r ), to obtain, ¿ R n (r ) ·ˆ 0 (r ; β)¿(r ; β) −¿(r ; β) ·ˆ 0 (r ; −β n )¿ R n (r ) = (β − β n )¿ R n (r ) · σ n (r ). (23) Dividing by β − β n and integrating over the spatial inhomogeneities in the limit β → β n , we get − 2N n = lim The subscript z indicates the integration of the z component in the second term which results in the surface integral of Eq. (9) when using that, due to symmetry, the in-plane components of the electric field and the z component of the magnetic field of resonant states with eigenvalues β n and −β n are identical, while we have to multiply all other components with −1 in order to convert ¿ R n into¿ n . The first term can be converted to a line integral by using the divergence theorem. The curve of integration is taken as a circle of radius R outside the region of inhomogeneities. For evaluating the limit β → β n , we carry out a Taylor expansion around β n aŝ